#  Elastic meson form factors in a unified scheme

Roger J. Hernández-Pinto 111Speaker<EMAIL_ADDRESS>Facultad de Ciencias
Físico-Matemáticas, Universidad Autónoma de Sinaloa, Ciudad Universitaria, CP
80000 Culiacán, México  M. Bedolla-Hernández  Facultad de Ciencias en Física y
Matemáticas, Universidad Autónoma de Chiapas,Tuxtla Gutiérrez, Chiapas, México
L. X. Gutiérrez-Guerrero  CONACyT-Mesoamerican Centre for Theoretical Physics,
Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México A. Bashir
Instituto de Física y Matemáticas, Universidad Michoacana de San Nicolás de
Hidalgo, Morelia, Michoacán, México

###### Abstract

The extraction of elastic form factors for mesons in the context of the
contact interaction model is revisited in this manuscript. The dressed masses
of quarks and mesons are determined through the gap and Bethe-Salpeter
equations. The generic elastic scattering process $M\gamma\to M$ is studied
for the meson $M$ formed of two differently flavored quarks. The charge radii
of scalar, pseudoscalar, vector and axial-vector mesons are also extracted by
virtue of explicit calculation of the meson elastic form factors.

###### pacs:

12.38.-t, 12.38.Aw, 14.40.-n, 13.40.Gp, 24.85.+p

day month yearday month year Quantum Chromodynamics, Mesons, Elastic Form
Factors

##  1 Introduction

The correct and precise description of hadron properties is a prohibitively
difficult task. Perturbative quantum chromodynamics (QCD) helps us understand
with great precision the Standard Model of particle physics in the regime of
asymptotic freedom. However, hadrons are manifestations of other emergent
properties of QCD, namely, the dynamical chiral symmetry breaking (DCSB) and
confinement. A standard method to investigate the internal structure of
hadrons is by the extraction of their form factors (FF). The prediction of
perturbative QCD for meson FF, [1], shows that there exists
$Q_{0}>\Lambda_{\rm QCD}$, such that

$Q^{2}F_{0^{-}}(Q^{2})\stackrel{{\scriptstyle
Q^{2}>Q_{0}^{2}}}{{\approx}}16\pi\alpha_{s}(Q^{2})f_{0^{-}}^{2}w_{0^{-}}^{2}(Q^{2}),$
(1)

where: $f_{0^{-}}$ is the meson’s leptonic decay constant; $\alpha_{s}(Q^{2})$
is the leading-order strong running-coupling and

$w_{0^{-}}(Q^{2})=\frac{1}{3}\int_{0}^{1}dx\,\frac{1}{x}\,\varphi_{0^{-}}(x;Q^{2})\,,$
(2)

where $\varphi_{0^{-}}(x;Q^{2})$ is the meson’s dressed-valence-quark parton
distribution amplitude. In this work, the elastic scattering process
$M\gamma\to M$ between a meson $M$ and a photon $\gamma$ is analyzed, with the
purpose of extracting information on the internal structure of mesons.

We start by providing the key ingredients for the calculation of the FF of
mesons in the contact interaction (CI) model. Firstly, the dressed mass of the
constituent quarks is calculated in this model. Then, the solution of the
Bethe-Salpeter equation (BSE) to compute the meson masses is detailed.
Finally, we extract the FF of mesons and determine their charge radii in this
unified scheme.

##  2 Contact Interaction and the gap equation

To study the internal structure of mesons in the $M\gamma\to M$ process, we
first need to determine the dressed mass of the quarks which constitute the
probed meson. In the CI model, it is calculated through solving the gap
equation for the quark propagator, given by

$\displaystyle S(p,M_{q})^{-1}=\imath\gamma\cdot p+m_{q}$
$\displaystyle+\int\frac{d^{4}q}{(2\pi)^{4}}g^{2}D_{\mu\nu}(p-q)\frac{\lambda^{a}}{2}\gamma^{\mu}S(q,m)\frac{\lambda^{a}}{2}\Gamma_{\nu}(q,p)\,,$
(3)

where $S(p,m)=(p^{2}+m^{2})^{-1}$ is the quark propagator, $M_{q}$ is the
dressed quark mass, $m_{q}$ is its current value, $\gamma^{\mu}$ and
$\lambda^{a}$ are the Dirac and Gell-Mann matrices respectively,
$\Gamma_{\nu}$ is the quark-photon vertex and $D_{\mu\nu}$ is the gluon
propagator. In perhaps the simplest truncation scheme, we consider
$\Gamma_{\nu}\equiv\gamma_{\nu}$. To solve Eq. (3) analytically, we use the
CI. Hence the gluon propagator takes the form [2],

$\displaystyle g^{2}D_{\mu\nu}(p-q)=\delta_{\mu\nu}\frac{1}{m_{G}^{2}}\,,$ (4)

with $m_{G}$ being the gluon mass parameter. The solution of Eq. (3) gives
dressed masses of the quarks. They are found by solving the transcendental gap
equation,

$\displaystyle
M_{q}=m_{q}+\frac{M_{q}}{3\pi^{2}m_{G}^{2}}\mathcal{C}^{iu}(M_{q}^{2})\,,$ (5)

where $\mathcal{C}^{iu}(x)/x\equiv\Gamma(-1,x\,\tau_{\rm
UV})-\Gamma(-1,x\,\tau_{\rm IR})$, with $\Gamma(a,b)$ the upper incomplete
gamma function, and $\tau_{\rm IR}$ and $\tau_{\rm UV}$ are the parameters
that regularize the integral [3, 4].

##  3 Bethe-Salpeter equation

Meson masses are computed through the BSE. This equation, in the CI model,
reads as [5],

$\displaystyle\Gamma(k;P)=-\frac{4}{3}\frac{1}{m_{G}^{2}}\int\frac{d^{4}}{(2\pi)^{4}}\gamma_{\mu}\chi(q;P)\gamma^{\mu}$
(6)

where $\chi(q;P)=S(q+P)\Gamma(q;P)S(q)$ and $\Gamma(q;P)$ is the Bethe-
Salpeter amplitude (BSA), $k$ is the external relative momentum and $P$ the
four momentum of the meson. Eq. (3) can be represented diagrammatically in
Figure 3.  Diagrammatic representation of the BSE. Blue (solid) circles
represent dressed quark propagators $S$, red (solid) circle is the meson BSA
and the blue (solid) rectangle is the dressed-quark-antiquark scattering
kernel $\mathcal{K}$. We remark that the BSA are necessary for the calculation
of the $M\gamma\to M$ process in our approach. They encode the interaction
properties between the meson and its constituent quarks. In the context of the
CI, the BSA can be decomposed as,

$\displaystyle\Gamma^{H}(P)=A^{H}(P)E^{H}(P)+B^{H}(P)F^{H}(P)\,,$ (7)

with $H={\rm S,PS,V,AV}$, representing the four kinds of mesons: Scalar (S),
Pseudoscalar (PS), Vector (V) and Axial-Vector (AV).

BSA | $A(P)$ | $B(P)$
---|---|---
$\Gamma^{\rm S}$ | $\imath\,I_{D}$ | –
$\Gamma^{\rm PS}$ | $\imath\,\gamma_{5}$ | $\frac{1}{2M_{R}}\gamma_{5}(\gamma\cdot P)$
$\Gamma^{\rm V}_{\mu}$ | $\gamma_{\mu}^{\rm T}$ | –
$\Gamma^{\rm AV}_{\mu}$ | $\gamma_{5}\,\gamma_{\mu}^{\rm T}$ | –

Table 1: BSA for scalar, pseudoscalar, vector and axial-vector mesons. In
here, $\gamma_{\mu}^{\rm T}$ satisfies $\gamma_{\mu}^{\rm T}P^{\mu}=0$.

The coefficients $A_{H}$ and $B_{H}$ are functions of Dirac matrices, $P$ and
the reduced mass of the constituent quarks,
$M_{R}=M_{q}M_{\bar{q}\prime}/(M_{q}+M_{\bar{q}\prime})$. The tensor structure
used in this work is shown in Table 1.

In the case of pseudoscalar mesons the BSE can be written in terms of $E^{\rm
PS}$ and $F^{\rm PS}$ as

$\displaystyle\left[\begin{array}[]{c}E^{\rm PS}(P)\\\ F^{\rm
PS}(P)\end{array}\right]=\frac{4\hat{\alpha}_{\rm
IR}}{3\pi}\left[\begin{array}[]{cc}\mathcal{K}^{\rm PS}_{EE}&\mathcal{K}^{\rm
PS}_{EF}\\\ \mathcal{K}^{\rm PS}_{FE}&\mathcal{K}^{\rm
PS}_{FF}\end{array}\right]\left[\begin{array}[]{c}E^{\rm PS}(P)\\\ F^{\rm
PS}(P)\end{array}\right]$ (14)

where $\hat{\alpha}_{\rm IR}=\alpha_{\rm IR}/m_{g}^{2}=1/(4\pi m_{G}^{2})$,
with $m_{g}=500$ MeV the mass of the gluon generated dynamically in QCD [6]
and,

$\displaystyle\mathcal{K}^{\rm PS}_{EE}$
$\displaystyle=\int_{0}^{1}d\alpha\left\\{\mathcal{C}^{iu}(\omega_{1})\right.$
$\displaystyle\left.+(M_{q}M_{\bar{q}\prime}-\alpha(1-\alpha)P^{2}-\omega_{1})\,\overline{\mathcal{C}}_{1}^{iu}(\omega_{1})\right\\},$
(15) $\displaystyle\mathcal{K}^{\rm PS}_{EF}$
$\displaystyle=\frac{P^{2}}{2M_{R}}\int_{0}^{1}d\alpha((1-\alpha)M_{\bar{q}\prime}+\alpha\,M_{q})\overline{\mathcal{C}}_{1}^{iu}(\omega_{1}),$
(16) $\displaystyle\mathcal{K}^{\rm PS}_{FE}$
$\displaystyle=\frac{2M_{R}^{2}}{P^{2}}\mathcal{K}^{\rm PS}_{EF},$ (17)
$\displaystyle\mathcal{K}^{\rm PS}_{FF}$
$\displaystyle=\frac{1}{2}\int_{0}^{1}d\alpha((\alpha-1)M_{\bar{q}\prime}^{2}-M_{q}M_{\bar{q}\prime}-\alpha
M_{q}^{2})$ $\displaystyle\times\overline{\mathcal{C}}_{1}^{iu}(\omega_{1}),$
(18)

where $\omega_{1}=M^{2}_{\bar{q}\prime}+\alpha
M_{q}^{2}+\alpha(1-\alpha)P^{2}$ and
$\overline{\mathcal{C}}_{1}^{iu}(z)=\mathcal{C}_{1}^{iu}(z)/z=\Gamma(0,z\,\tau_{\rm
UV})-\Gamma(0,z\,\tau_{\rm IR})$. Finally, the eigenvalue equation, Eq. (14),
has a solution at $P^{2}=-m_{H}^{2}$. In the case of vector, axial-vector and
scalar mesons respectively, the BSE read as

$\displaystyle\mathcal{K}_{\rm V}(P^{2})$
$\displaystyle=\frac{2\hat{\alpha}_{\rm
IR}}{3\pi}\int_{0}^{1}d\alpha\mathcal{L}_{\rm
V}(P^{2})\overline{\mathcal{C}}_{1}^{iu}(\omega_{1}),$ (19)
$\displaystyle\mathcal{K}_{\rm AV}(P^{2})$
$\displaystyle=\frac{2\hat{\alpha}_{\rm
IR}}{3\pi}\int_{0}^{1}d\alpha(\mathcal{C}_{1}^{iu}(\omega_{1})+\mathcal{L}_{\rm
G}(P^{2})\overline{\mathcal{C}}_{1}^{iu}(\omega_{1})),$ (20)
$\displaystyle\mathcal{K}_{\rm S}(P^{2})$
$\displaystyle=\frac{4\hat{\alpha}_{\rm
IR}}{3\pi}\int_{0}^{1}d\alpha\left[\mathcal{L}_{\rm
G}(P^{2})\overline{\mathcal{C}}_{1}^{iu}(\omega_{1})\right.$
$\displaystyle-\left.\left(\mathcal{C}^{iu}(\omega_{1})-\overline{\mathcal{C}}_{1}^{iu}(\omega_{1})\right)\right],$
(21)

where

$\displaystyle\mathcal{L}_{\rm V}(P^{2})$
$\displaystyle=M_{\bar{q}\prime}M_{q}-(1-\alpha)M^{2}_{\bar{q}\prime}$
$\displaystyle-\alpha M^{2}_{q}-2\alpha(1-\alpha)P^{2},$ (22)
$\displaystyle\mathcal{L}_{\rm G}(P^{2})$
$\displaystyle=M_{q}M_{\bar{q}\prime}+\alpha(1-\alpha)P^{2}.$ (23)

The normalized solutions of Eqs. (19)-(21) provide masses of the vector
($m_{\rm V}$), axial-vector ($m_{\rm AV}$) and scalar ($m_{\rm S}$) mesons
when,

$\displaystyle 1+(-1)^{i}\mathcal{K}_{j}(-m_{j}^{2})=0$ (24)

with $i=\\{0,1,0\\}$ for $j=\\{\rm S,V,AV\\}$ respectively.

The procedure described above provides good phenomenological results in the
case of vector mesons and flavour non-singlet pseudoscalar meson ground-states
[7, 8, 10, 9, 11]. For the parity partners, spin-orbit repulsion has to be
taken into account [12, 13, 14, 15, 16]. In the CI, a phenomenological
coupling $g_{SO}<1$ is introduced as a multiplicative common factor to the
kernels. The numerical values for mesons with $J=0^{+},1^{+}$ are [17],

$\displaystyle g_{SO}^{0^{+}}=0.32\,,\quad g_{SO}^{1^{+}}=0.25\,.$ (25)

To conclude this section, it is important to remark that the solutions of Eq.
(14) and Eqs. (19)-(21) constrain the values of $E^{H}$ and $F^{H}$. These
quantities in turn are crucial in subsequently computing the elastic FF of
mesons.

##  4 Elastic Form Factors of mesons

Elastic meson FF are calculated from the triangle diagram for $M\gamma\to M$,
with a quark and an anti-quark circulating in the loop, see Figure 4. We point
out that in this article, the meson is considered to be made of two
differently flavoured quarks.  Feynman diagram for the $M\gamma M$ vertex
which permits the extraction of elastic meson FF of all mesons. Labelling of
internal momenta is included in the diagram.

Considering that the photon interacts with the quark while the anti-quark is a
spectator, Feynman rules permit us to write this process as,

$\displaystyle\Lambda_{(\nu)\mu(\alpha)}$ $\displaystyle=2\,N_{c}\,{\rm
Tr}\left[\imath\Gamma^{M}_{(\nu)}(-p_{2})S_{F}(t+p_{2},M_{q})\Gamma_{\mu}\right.$
$\displaystyle\times\left.S_{F}(t+p_{1},M_{q})\imath\Gamma^{M}_{(\alpha)}(p_{1})S_{F}(t,M_{\bar{q}\prime})\right]$
(26)

where $\Gamma_{\mu}\equiv iP_{T}(Q^{2})\gamma^{\mu}$ is the photon-quark
dressed vertex, $P_{T}(Q^{2})$ is the dressing function [18], and the
subscript between parentheses indicates that the meson may or may not have
Lorentz indices, i.e. $\Lambda_{\mu}$ shall be understood for scalar and
pseudoscalar mesons and $\Lambda_{\nu\mu\alpha}$ shall be used for vector and
axial-vector mesons. The dressing function is given by,
$P_{T}(Q^{2})=(1+K_{\gamma}(Q^{2}))^{-1}$, where

$\displaystyle
K_{\gamma}(Q^{2})=\frac{1}{3\pi^{2}m_{G}^{2}}\int_{0}^{1}d\alpha\,\alpha(1-\alpha)Q^{2}\overline{\mathcal{C}}_{1}(\omega)$
(27)

with $\omega=M_{q}^{2}+\alpha(1-\alpha)Q^{2}$.

In the case of scalar and pseudoscalar mesons, there is only one
electromagnetic FF, $F_{\rm S,PS}^{\rm em}$. On the other hand, vector and
axial-vector mesons have three form factors due to their tensor structure,
$F_{\rm V,AV}^{\rm 1,2,3}$ [19]. In terms of these FFs, the electric
($G_{E}$), magnetic ($G_{M}$) and quadrupole ($G_{Q}$) FFs are

$\displaystyle G_{E}$ $\displaystyle=F^{1}+\frac{2}{3}\eta G_{Q}\,,\quad
G_{M}=-F^{2}\,,$ (28) $\displaystyle G_{Q}$
$\displaystyle=F^{1}+F^{2}+(1+\eta)F^{3}\,.$ (29)

where $\eta\equiv Q^{2}/(4m_{H}^{2})$ and $m_{H}$ is the meson mass. After
detailed analytical calculations, all FF can be written as,

$\displaystyle F^{j}_{i}(Q^{2})$
$\displaystyle=\int_{0}^{1}d\alpha\,d\beta\,\alpha\,$
$\displaystyle\times\left(\mathcal{A}^{j}_{i}\,\overline{\mathcal{C}}_{1}(\omega_{2})+(\mathcal{B}^{j}_{i}-\mathcal{A}^{j}_{i}\,\omega_{2})\,\overline{\mathcal{C}}_{2}(\omega_{2})\right)$
(30)

with $i=\\{\rm S,PS,V,AV\\}$ and $j=\\{\rm em\\}$ for the scalar and
pseudoscalar mesons scenario and $j=\\{1,2,3\\}$ for the vector and axial-
vector meson case. In addition,
$\overline{\mathcal{C}}_{2}(z)=\left(\exp(-z\,\tau^{2}_{{\rm
UV}})-\exp(-z\,\tau_{\rm IR})\right)/(2z)$ and
$\omega_{2}=\alpha\,M_{q}^{2}+(1-\alpha)M_{q\prime}^{2}-\alpha(1-\alpha)\,M^{2}_{H}+\alpha^{2}\,\beta\,(1-\beta)\,Q^{2}$.
Naturally, the calculation of the FF also yields information of the magnetic
and quadrupole moments for vector and axial-vector mesons, defined as

$\displaystyle G^{i}_{M}(Q^{2}=0)\equiv\mu^{i}\,,\qquad
G^{i}_{Q}(Q^{2}=0)\equiv\mathcal{Q}^{i}\,,$ (31)

with $i=\\{\rm V,AV\\}$, and these results can be also compared with
experimental observation. Note that in the case of mesons composed of
differently flavored quarks, we define the meson FF, $\mathcal{F}_{i}^{j}$, as
the sum of two contributions since the photon can interact either with the
quark, $F_{i}^{j}$ or with the anti-quark, $\bar{F}_{i}^{j}$. Then [20],

$\displaystyle\mathcal{F}_{i}^{j}(Q^{2})=e_{q}F_{i}^{j}(Q^{2})+e_{\bar{q}}\bar{F}_{i}^{j}(Q^{2})$
(32)

with $e_{q}$ and $e_{\bar{q}}$ the quark and anti-quark charges respectively.
The procedure described up can readily be adapted for neutral meson with same
flavored quarks.

Finally, once the analytic expression of the coefficients $\mathcal{A}$ and
$\mathcal{B}$ are obtained, the comparison between the CI model and the
experimental results can also be carried out through the charge radii of
different mesons, computed through the relation

$\displaystyle
r^{2}_{i}\equiv\Big{|}6\frac{d}{dQ^{2}}F_{i}(Q^{2})\Big{|}_{Q^{2}=0}\,.\quad
i=\\{\rm S,PS,V,AV\\}$ (33)

This observable can constrain several models for mesons. An important feature
of the FF is the dependence on $Q^{2}$. In Ref [2], it is pointed out that the
behaviour of the pion electromegnetic FF at large $Q^{2}$ tends to a constant
value in the CI. For the rho meson, it turns negative for $Q^{2}\sim 6$ GeV.
In this spirit, it would be interesting to delve into the behaviour of the FF
of all kind of mesons in order to understand the predictions of the internal
structure through the CI model.

##  5 Conclusions

In this document, the extraction of the elastic FF of scalar, pseudoscalar,
vector and axial-vector mesons in the framework of the CI model has been
presented in detail. The determination of the dressed quark and meson masses
according to the gap and BSE helps constrain the parameters of the model. The
study of the $M\gamma\to M$ process is presented in order to extract the
elastic FF of mesons formed with two differently flavored quarks. Once the
general expression is found, the charge radii of scalar, pseudoscalar, vector
and axial-vector mesons can also be computed and compared with experimental
results. Finally, interesting information such as the magnetic and quadrupole
moments can also also be obtained in order to establish the extent of validity
of the CI model more realistically.

## Acknowledgments

This research was supported by CONACyT through the Project No. A1- S-33202
(Ciencia Basica), Ciencia de Frontera 2021-2042 and Sistema Nacional de
Investigadores. The work is also supported by PROFAPI 2022 Grant No.
PRO_A1_024 (Universidad Autónoma de Sinaloa), the COST Action CA16201
(PARTICLEFACE) and MCIN/AEI/10.13039/501100011033, Grant No.
PID2020-114473GB-I00 and the CIC project 4.10 of UMSNH.

## References

  * $1$. G. P. Lepage and S. J. Brodsky, Phys. Lett. B 87, 359-365 (1979) doi:10.1016/0370-2693(79)90554-9
  * $2$. L. X. Gutierrez-Guerrero, A. Bashir, I. C. Cloet and C. D. Roberts, Phys. Rev. C 81, 065202 (2010) doi:10.1103/PhysRevC.81.065202 [arXiv:1002.1968 [nucl-th]].
  * $3$. G. Krein, C. D. Roberts and A. G. Williams, Int. J. Mod. Phys. A 7, 5607-5624 (1992) doi:10.1142/S0217751X92002544
  * $4$. C. D. Roberts, Prog. Part. Nucl. Phys. 61, 50-65 (2008) doi:10.1016/j.ppnp.2007.12.034 [arXiv:0712.0633 [nucl-th]].
  * $5$. E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232-1242 (1951) doi:10.1103/PhysRev.84.1232
  * $6$. F. Gao, S. X. Qin, C. D. Roberts and J. Rodriguez-Quintero, Phys. Rev. D 97, no.3, 034010 (2018) doi:10.1103/PhysRevD.97.034010 [arXiv:1706.04681 [hep-ph]].
  * $7$. C. Chen, L. Chang, C. D. Roberts, S. Wan and D. J. Wilson, Few Body Syst. 53, 293-326 (2012) doi:10.1007/s00601-012-0466-3 [arXiv:1204.2553 [nucl-th]].
  * $8$. S. x. Qin, L. Chang, Y. x. Liu, C. D. Roberts and D. J. Wilson, Phys. Rev. C 85, 035202 (2012) doi:10.1103/PhysRevC.85.035202 [arXiv:1109.3459 [nucl-th]].
  * $9$. S. x. Qin, L. Chang, Y. x. Liu, C. D. Roberts and D. J. Wilson, Phys. Rev. C 84, 042202 (2011) doi:10.1103/PhysRevC.84.042202 [arXiv:1108.0603 [nucl-th]].
  * $10$. P. Maris, AIP Conf. Proc. 892, no.1, 65-71 (2007) doi:10.1063/1.2714348 [arXiv:nucl-th/0611057 [nucl-th]].
  * $11$. I. C. Cloet, A. Krassnigg and C. D. Roberts, eConf C070910, 125 (2007) [arXiv:0710.5746 [nucl-th]].
  * $12$. R. Bermudez, L. Albino, L. X. Gutiérrez-Guerrero, M. E. Tejeda-Yeomans and A. Bashir, Phys. Rev. D 95, no.3, 034041 (2017) doi:10.1103/PhysRevD.95.034041 [arXiv:1702.04437 [hep-ph]].
  * $13$. A. Bashir, R. Bermudez, L. Chang and C. D. Roberts, Phys. Rev. C 85, 045205 (2012) doi:10.1103/PhysRevC.85.045205 [arXiv:1112.4847 [nucl-th]].
  * $14$. L. Chang, Y. X. Liu and C. D. Roberts, Phys. Rev. Lett. 106, 072001 (2011) doi:10.1103/PhysRevLett.106.072001 [arXiv:1009.3458 [nucl-th]].
  * $15$. L. Chang and C. D. Roberts, Phys. Rev. C 85, 052201 (2012) doi:10.1103/PhysRevC.85.052201 [arXiv:1104.4821 [nucl-th]].
  * $16$. L. Chang and C. D. Roberts, AIP Conf. Proc. 1361, no.1, 91-114 (2011) doi:10.1063/1.3622689 [arXiv:1003.5006 [nucl-th]].
  * $17$. L. X. Gutiérrez-Guerrero, G. Paredes-Torres and A. Bashir, Phys. Rev. D 104, no.9, 094013 (2021) doi:10.1103/PhysRevD.104.094013 [arXiv:2109.09058 [hep-ph]].
  * $18$. P. Maris and P. C. Tandy, Phys. Rev. C 61, 045202 (2000) doi:10.1103/PhysRevC.61.045202 [arXiv:nucl-th/9910033 [nucl-th]].
  * $19$. M. S. Bhagwat and P. Maris, Phys. Rev. C 77, 025203 (2008) doi:10.1103/PhysRevC.77.025203 [arXiv:nucl-th/0612069 [nucl-th]].
  * $20$. P. T. P. Hutauruk, I. C. Cloet and A. W. Thomas, Phys. Rev. C 94, no.3, 035201 (2016) doi:10.1103/PhysRevC.94.035201 [arXiv:1604.02853 [nucl-th]].


# Dynamic heterogeneity in active glass-forming liquids is qualitatively
different compared to its equilibrium behaviour

Kallol Paul1 Saroj Kumar Nandi1 Smarajit Karmakar1<EMAIL_ADDRESS>1
TIFR Center for Interdisciplinary Science, Tata Institute of Fundamental
Research, 36/P Gopanpally Village, Serilingampally Mandal, RR District,
Hyderabad, 500075, Telangana, India

###### Abstract

Activity driven glassy dynamics is ubiquitous in collective cell migration,
intracellular transport, dynamics in bacterial and ants colonies as well as
artificially driven synthetic systems such as vibrated granular materials,
etc. Active glasses are hitherto assumed to be qualitatively similar to their
equilibrium counterparts at a suitably defined effective temperature,
$T_{eff}$. Combining large-scale simulations with analytical mode-coupling
theory for such systems, we show that, in fact, an active glass is
qualitatively different from an equilibrium glassy system. Although the
relaxation dynamics can be similar to an equilibrium system at a $T_{eff}$,
effects of activity on the dynamic heterogeneity (DH), which has emerged as a
cornerstone of glassy dynamics, is quite nontrivial and complex. In
particular, active glasses show dramatic growth of DH, and systems with
similar relaxation time and $T_{eff}$ can have widely varying DH. Comparison
of our non-equilibrium extended mode-coupling theory for such systems with
simulation results show that the theory captures the basic characteristics of
such systems. Our study raises fundamental questions on the supposedly central
role of DH in controlling the relaxation dynamics in a glassy system and can
have important implications even for the equilibrium glassy dynamics.

Dynamic Heterogeneity $|$ Activity $|$ Mode Coupling Theory $|$ Four-point
Susceptibility

Active glass, characterizing the slow dynamics in a dense-active system of
self-propelled particles (SPP), with a self-propulsion force, $f_{0}$, and a
persistence time, $\tau_{p}$, of their motion, remains an intriguing problem
of statistical physics. A number of recent experiments show that glass-like
dynamics, that is the extreme dynamical slow-down of a system without any
discernible structural change or phase transition Berthier and Biroli (2011);
Angelini  _et al._ (2011), is important for the dynamics in several biological
and biology-inspired systems, such as a cellular monolayer Angelini  _et al._
(2011); Garcia  _et al._ (2015), cell-cytoplasm Parry  _et al._ (2014); Sadati
_et al._ (2014); Zhou  _et al._ (2009); Nishizawa  _et al._ (2017), colonies
of bacteria Takatori and Mandadapu (2020) and ants Gravish  _et al._ (2013,
2015), vertically vibrated rods Dijksman  _et al._ (2011); Dauchot  _et al._
(2005); J.Deseigne  _et al._ (2010), and several other artificial active
systems Dreyfus  _et al._ (2005); Palacci  _et al._ (2010); Klongvessa  _et
al._ (2019). These systems can be conveniently modeled via a dense active
matter of SPPs Ramaswamy (2010); Marchetti  _et al._ (2013). Motivated by the
experimental findings, such systems have been studied via simulations Ni  _et
al._ (2013); Berthier (2014); Mandal  _et al._ (2016); Flenner  _et al._
(2016); Mandal  _et al._ (2017, 2020) as well as extended analytical theories
of equilibrium glasses, such as mode-coupling theory (MCT) Berthier and
Kurchan (2013); Nandi and Gov (2017); Liluashvili  _et al._ (2017); Feng and
Hou (2017); Nandi (2018); Szamel (2016) and random first order transition
theory Nandi (2018). Studies on single-particle dynamics Chaki and Chakrabarty
(2020) and jamming transition Ni  _et al._ (2013); Merrigan  _et al._ (2020);
Bi  _et al._ (2016) have also been extended for active systems. The current
consensus in the field is that activity pushes the glass and jamming
transitions to lower temperature or higher density, though the nature of an
active glass seems to be qualitatively similar to an equilibrium glass and the
difference lies in the quantitative description, such as in the form of long-
ranged velocity correlation Szamel (2016); Caprini  _et al._ (2020) or the
evolving effective temperature, $T_{eff}$, much like a sheared system Nandi
and Gov (2017); Berthier and Kurchan (2013). However, this apparent similarity
of glassy dynamics in an active and equilibrium system is somewhat puzzling as
activity is known to exhibit nontrivial behavior, such as flocking Vicsek  _et
al._ (1995) or giant number fluctuation Narayan  _et al._ (2007); Ramaswamy
_et al._ (2003); Marchetti  _et al._ (2013); Ramaswamy (2010) in a dilute
system, and calls for rigorous theoretical exploration of the problem.

Apart from developing a coherent framework for the slow dynamics in these
biological systems, such theoretical studies can also provide deeper insights
into the physics of equilibrium glasses as they extend the scope and extent of
the problem in the form of control parameters as well as emergent behaviors.
In this work, through a combination of detailed large scale molecular dynamics
simulations and an analytical study via extended active MCT Nandi and Gov
(2017), we show that, in fact, an active glass is qualitatively different from
its equilibrium counterparts: whereas the relaxation dynamics is similar to
that of equilibrium glasses, dramatic effects of activity is observed in the
dynamic heterogeneity (DH). DH refers to the coexistence of heterogeneous
transient local dynamic behavior and has emerged as a salient feature of
glassy dynamics Ediger (2000); Berthier (2011); Berthier and Biroli (2011).
Confirmed both in simulations Yamamoto and Onuki (1998); Franz and Parisi
(2000); Karmakar  _et al._ (2009) and experiments Weeks  _et al._ (2000);
Dauchot  _et al._ (2005), existence of DH in glassy systems provide the much
sought after growing length scale to interpret the problem as a critical
phenomena with an underlying transition where the DH length scale controls the
relaxation dynamics in a glassy system. Compared to equilibrium systems,
active glasses show dramatic growth in DH with increasing activity and systems
with similar relaxation times can have varying DH suggesting a possible
decoupling from relaxation dynamics.

Although most works till date have focussed on the relaxation dynamics in an
active glass, there do exist a handfull of studies Avila  _et al._ (2014);
Malinverno  _et al._ (2017); Klongvessa  _et al._ (2019); Cerbino  _et al._
(2021) investigating different aspects of DH revealing nontrivial effects of
activity: Avila et al found Avila  _et al._ (2014) strong DH and large dynamic
length scale in a driven granular fluid compared to those in equilibrium
systems, experiments on active cellular monolayers Malinverno  _et al._
(2017); Cerbino  _et al._ (2021) show higher DH with increasing activity
though the relaxation time becomes smaller, a trend opposite to what one
expects for equilibrium glasses. Our theoretical study rationalizes these
puzzling results. We first show that the relaxation dynamics of an active
glass can be described similar to that of an equilibrium system at an
effective temperature, $T_{eff}$ (Fig. 1), however, the behavior of DH is
drastically different (Fig. 2). We then present reliable computations of
length-scales associated with DH in an active system via three different
methods (Fig. 3) and finally, illustrate the nontrivial aspects of activity
(Fig. 4). In particular, we show that systems with the same relaxation time
can have different DH; thus the dynamic length scale may not always control
the relaxation dynamics in a glassy system, at least when energy is pumped at
a microscopic length scale.

In this work, we study two different, but well-known, three-dimensional models
of glassy systems: the Kob-Andersen binary mixture (3dKA) with $80:20$ ratio
of two types of particles and another binary mixture (3dHP) with $50:50$ ratio
of two types of particles with diameter ratio of $1.4$. Results of 3dKA model
are presented in this article and some of the results of 3dHP model are shown
in SM, Fig. S10. Activity in our systems are introduced by designating $10\%$
of the particles, chosen at random, as active. These active particles are then
given constant forces of magnitude $f_{0}$ in random directions with the
constraints that vector sum of these applied active forces adds up to zero;
this ensures momentum conservation of the centre of mass in our simulations.
The applied active forces are reshuffled every $\tau_{p}$ time. In this work,
we keep density and $\tau_{p}=1.0$ constant and study the glassy properties as
functions of temperature, $T$, and $f_{0}$. We also extend active MCT,
presented in Ref. Nandi and Gov (2017), to study DH following the IMCT
formalism presented in Biroli  _et al._ (2006). More details of the models and
the simulations can be found in the materials and methods section and in the
SM. The details of the theoretical calculations will be presented elsewhere.

## Results

Figure 1: Characterization of relaxation dynamics. (a) Decay of the overlap
function, $Q(t)$, and growth of MSD, $\langle\Delta r(t)^{2}\rangle$, (inset)
becomes faster with incresing activity. (b) Symbols are $\tau_{\alpha}$ as a
function of $T$ for a passive system and the line is a fit with the MCT
prediction $\tau_{\alpha}\sim(T-T_{MCT})^{-\gamma}$ that gives $T_{MCT}=0.423$
and $\gamma=2.36$. (c) $\tau_{\alpha}$ as a function of $f_{0}$ at different
$T$, symbols are simulation data and the lines are active MCT prediction (not
individual fits). (d) Equation (1) predicts a master curve when
$\tau_{\alpha}(T-T_{MCT})^{\gamma}$ is plotted as a function of
$f_{0}^{2}/(T-T_{MCT})$ at different $T$ and $f_{0}$: the simulation data
(symbols) agrees well with the theoretical prediction (line).

### The relaxation dynamics

We first show that the relaxation dynamics of an active glass is similar to
that of an equililibrium system at a suitably chosen $T_{eff}$. The overlap
function, $Q(t)$, and the mean-squared displacement (MSD), $\langle\Delta
r(t)^{2}\rangle$, are shown in Fig. 1(a) and in the inset respectively, as a
function of time $t$ (see Materials and Methods for definitions). At a
particular $T=0.45$, both $Q(t)$ and MSD for the passive system ($f_{0}=0$)
show the plateau, characteristic of glassy dynamics. As activity increases
from $0$, decay of $Q(t)$ and growth of MSD become faster, as has also been
found in earlier simulations Berthier (2014); Flenner  _et al._ (2016); Mandal
_et al._ (2016) and experiments Klongvessa  _et al._ (2019). The relaxation
time, $\tau_{\alpha}$, is defined as $Q(\tau_{\alpha})=1/e$. MCT for the
active systems Nandi and Gov (2017) predicts

$\tau_{\alpha}\sim|T-T_{MCT}+Kf_{0}^{2}|^{-\gamma}$ (1)

where $\gamma$ is the same exponent as in a passive system, $T_{MCT}$ is the
mode-coupling transition temperature in the absence of activity, and $K$, a
constant. For the particular model of activity that we are using in this work,
defined as SNTC in Ref. Nandi and Gov (2017), we have
$K=H\tau_{p}/(1+G\tau_{p})$ with $H$ and $G$ being two constants. We show the
data of $\tau_{\alpha}(T)$ for the equilibrium 3dKA model in Fig. 1 (b) (all
simulation data, presented here are for this model, results for the 3dHP model
are shown in SM Fig. S10). A fit of the MCT prediction,
$\tau_{\alpha}\sim|T-T_{MCT}|^{-\gamma}$, with the simulation data gives
$\gamma\simeq 2.36$ and $T_{MCT}\simeq 0.423$. These numbers are similar to
previously reported results for the passive system Kob and Andersen (1995); S
and Sastry (2003). Using these values, we fit one set of data for the active
system and obtain $K\simeq 0.015$. We then plot Eq. 1 along with the
simulation data in Fig. 1(c) and show that the theory agrees remarkably well
with the simulation data. Note that once $K$ is obtained from one set of data,
there are no other free parameters in the theory. Moreover, Eq. 1 shows plot
of $\mathcal{Y}=\tau_{\alpha}(T-T_{MCT})^{\gamma}$ as a function of
$\mathcal{X}=f_{0}^{2}/(T-T_{MCT})$ at different values of activity should
follow a master curve: $\mathcal{Y}\sim(1+K\mathcal{X})^{-\gamma}$. Fig. 1(d)
shows excellent agreement between simulation data and the MCT prediction.

From Eq. (1), we also find that both the MCT transition temperature, $T_{C}$,
for the active system as well as the calorimetric glass transition
temperature, $T_{g}$, defined via $\tau_{\alpha}(T=T_{g})=10^{6}$, behave as
$T_{C,g}=a-bf_{0}^{2}$, with $a$ and $b$ being two constants. As shown in SM
Fig. S7, this prediction also agrees well with simulation data. The excellent
agreement between the theory, Eq. 1, and simulation data with $\gamma$ being
the same exponent as for the passive system can also be interpreted as the
relaxation dynamics of the active system being similar to that of the
equilibrium glass at $T_{eff}=T+Kf_{0}^{2}$. However, as we show below, an
active glass behaves quite differently from an equilibrium glass. The
nontrivial effects of activity manifest via the DH and the dynamic length
scale, characterized via four-point correlation function, $\chi_{4}(t)$. We
now present the extension of active MCT to study DH via $\chi_{4}(t)$ along
with the simulation results for DH and dynamic length scale.

### Active In-homogeneous Mode Coupling Theory (Active-IMCT)

To the best of our knowledge, MCT for an active system has not been extended
to study the DH, the detailed derivation will be presented elsewhere, here we
briefly outline the main result. We follow the IMCT formalism, developed in
Ref. Biroli  _et al._ (2006), and obtain the four-point correlation function
in terms of a corresponding susceptibility via linear response theory. Since
an active system is inherently out of equilibrium, one must study both the
correlation and (integrated) response functions, $Q(t)$ and $F(t)$
respectively, as well as their corresponding susceptibilities. We impose a
weak external field, $\bar{\varepsilon}$, derive the equations of motion for
$Q(t)$ and $F(t)$ in the presence of $\bar{\varepsilon}$, and obtain the
susceptibilities, $\chi_{Q}(t)$ and $\chi_{F}(t)$, defined as

$\chi_{Q}(t)=\frac{\partial
Q(t)}{\partial\bar{\varepsilon}}\bigg{|}_{\bar{\varepsilon}\to
0}\,\,\,\text{and}\,\,\,\chi_{F}(t)=\frac{\partial
F(t)}{\partial\bar{\varepsilon}}\bigg{|}_{\bar{\varepsilon}\to 0}.$ (2)

The equations of motions for $\chi_{Q}(t)$ and $\chi_{F}(t)$ are

$\displaystyle\frac{\partial\chi_{Q}(t)}{\partial t}$
$\displaystyle=\mu(t)+(1+\zeta)Q(t)-(T-p)\chi_{Q}(t)$ $\displaystyle-$
$\displaystyle\int_{0}^{t}m(t-s)\frac{\partial\chi_{Q}(t)}{\partial
s}\mathrm{d}s-\int_{0}^{t}\Sigma(t-s)\frac{\partial Q(s)}{\partial
s}\mathrm{d}s$ (3a) $\displaystyle\frac{\partial\chi_{F}(t)}{\partial t}$
$\displaystyle=(1+\zeta)F(t)-(T-p)\chi_{F}(t)$ $\displaystyle-$
$\displaystyle\int_{0}^{t}m(t-s)\frac{\partial\chi_{F}(t)}{\partial
s}\mathrm{d}s-\int_{0}^{t}\Sigma(t-s)\frac{\partial F(s)}{\partial
s}\mathrm{d}s,$ (3b)

$\displaystyle\text{where, }\mu(t)$
$\displaystyle=-\int_{t}^{\infty}\Delta(s)\frac{\partial\chi_{F}(s-t)}{\partial
s}\mathrm{d}s,$ (4) $\displaystyle\zeta=\int_{0}^{\infty}$
$\displaystyle\Delta(s)\frac{\partial\chi_{F}(s)}{\partial
s}\mathrm{d}s,\,\,\,p=\int_{0}^{\infty}\Delta(s)\frac{\partial F(s)}{\partial
s}\mathrm{d}s$ (5) $\displaystyle m(t)$
$\displaystyle=2\lambda\frac{Q(t)^{2}}{T_{eff}(t)},$ (6)
$\displaystyle\Sigma(t)$
$\displaystyle=2\lambda\frac{Q(t)\chi_{Q}(t)}{T_{eff}(t)}+2\lambda\frac{Q(t)^{2}}{T_{eff}(t)^{2}}\kappa
f_{0}^{2},$ (7) $\displaystyle\text{and}\,\,\,T_{eff}(t)$
$\displaystyle=\frac{\partial Q(t)}{\partial t}\bigg{/}\frac{\partial
F(t)}{\partial t},$ (8)

where $\Delta(s)$ is the active force statistics:
$\Delta(s)=f_{0}^{2}\exp[-s/\tau_{p}]$. $\kappa$ is a system-dependent
constant, set to $\kappa=2.0$ in this work and a more detailed analysis will
be presented elsewhere. As shown in Ref. Biroli  _et al._ (2006),
$\chi_{Q}(t)\equiv\chi_{4}(t)$ is the desired four-point correlation function.
These equations must be simoultaneously solved along with the equations for
$Q(t)$ and $F(t)$, as presented in Ref. Nandi and Gov (2017).

Figure 2: Characterization of $\chi_{4}(t)$ as a function of $T$ and activity.
(a) Simulation data for $\chi_{4}(t)$ at different $f_{0}$ and $T=0.45$.
$\chi_{4}^{P}$, the peak value of $\chi_{4}(t)$, decreases with increasing
activity. (b) Active-IMCT results for $\chi_{4}(t)$ at different $f_{0}$ and
$\lambda=1.95$; theory shows similar results as in simulation. (c) Symbols are
simulation data for $\chi_{4}^{P}$ as a function of $T$ at different $f_{0}$
and lines are fits with a power-law where the exponent depends on $f_{0}$ (see
text). (d) $\chi_{4}^{P}$ as a function of $f_{0}$ at different $T$; fitting
the simulation data (symbols) with a power-law form gives the exponent $\mu$
that depends on $T$ (inset).

### Dynamic heterogeneity in active glasses

We now present the simulation results showing the effect of activity on DH,
characterized via $\chi_{4}(t)$, defined in the Materials and Methods. Fig.
2(a) shows $\chi_{4}(t)$ at $T=0.45$ and different values of $f_{0}$;
$\chi_{4}^{P}$, the peak value of $\chi_{4}(t)$, decreases with increasing
$f_{0}$. Similar behavior is also found within MCT as shown in Fig. 2(b) where
we have used $\lambda=1.95$ (the MCT transition is found at $\lambda=2.0$ in
our notation). Fig. 1(c) shows that $f_{0}$ decreases $\tau_{\alpha}$ and
drives the system away from the glassy regime, therefore, decrease of
$\chi_{4}^{P}$ with increasing $f_{0}$ is expected; but, the trend is
different from equilibrium-like behavior. Fig. 2(c) shows the simulation data
of $\chi_{4}^{P}$ as a function of $T$ at two different $f_{0}$ as well as the
data for $f_{0}=0$. The lines are fit with a function
$\chi_{4}^{P}\sim(T-T_{C})^{-\nu}$. If the glassy dynamics of the active
system resembles the equilibrium dynamics at a suitably defined $T_{C}$, one
expects $\nu$ to remain same at different activity and identical to that of
the passive system. However, we find the exponent $\nu$ depends on activity:
$0.92$ for the passive system, whereas $\nu=1.03$ and $1.13$ for $f_{0}=2.0$
and $2.5$ respectively. Irrespective of the detailed values of $\nu$, the
important point of this analysis is that effect of activity on DH is quite
different from the equilibrium-like dynamics. To probe this point from another
perspective, we now study $\chi_{4}^{P}$ as a function of $f_{0}$ at different
$T$ (Fig. 2d). A fit of the function $\chi_{4}^{P}\sim(1+bf_{0}^{2})^{-\mu}$,
where $b$ is a constant, with the simulation data shows $\mu$ depends on $T$
(inset, Fig. 2d). Thus, varying activity not-only changes $T_{C}$, it also
affects the exponents. The time $\tau_{peak}$, when $\chi_{4}(t)$ attains its
maximum, provides a measure of a relaxation time. As shown in the SM Fig. S8,
$\tau_{peak}$ is proportional to $\tau_{\alpha}$. Thus, the relaxation
dynamics, characterized via either of $Q(t)$ or $\chi_{4}(t)$, can be
described in terms of a $T_{eff}$ (Figs. 1c and 1d). These results show quite
distinctive effects of activity on relaxation dynamics and DH; whereas effect
of activity can simply be treated in terms of a $T_{eff}$ for the former, it
cannot be done so for the latter.

Within MCT, the nontrivial effect of activity on the dynamic heterogeneity is
contained in the last term in Eq. (7). This term is a direct consequence of
the nonequilibrium nature of the system and absent in equilibrium. An
important characteristic of activity is the giant number fluctuation (GNF)
Ramaswamy  _et al._ (2003); Narayan  _et al._ (2007), where number
fluctuation, $\delta N$, is related to the mean, $\langle N\rangle$, as
$\langle\delta N^{2}\rangle\sim\langle N\rangle^{\alpha}$ with $\alpha>1$. GNF
has been experimentally measured in Ref. Giavazzi  _et al._ (2017) for the
system of Ref. Malinverno  _et al._ (2017) and it has been argued that the
nontrivial behavior of $\chi_{4}(t)$ with activity might be related to GNF
Giavazzi  _et al._ (2017); Cerbino  _et al._ (2021). The detailed mechanism of
how GNF affects DH and why its effect on the relaxation dynamics remains
equlibrium-like remains an intriguing question. The distinctive natures of the
effects of activity on the relaxation dynamics and DH length scales imply they
decouple from each other. Such a decoupling, even in equilibrium, has been
suggested in Karmakar  _et al._ (2009), where it has been convincingly shown
that the finite size effects of $\chi_{4}^{P}$ and $\tau_{\alpha}$ are not
controlled by the same length scale. In fact, this is one of the first direct
signatures of the existence of yet another length scale, the static length
scale of amorphous order Karmakar  _et al._ (2015), related to structural
relaxation time. It remains to be seen whether the strong decoupling of DH and
structural relaxation in active glasses is an enhanced manifestation of the
same equilibrium phenomena or distinct from it.

Figure 3: Estimation of dynamic heterogeneity length scale. (a) Schematic of
block analysis method showing how the whole system is subdivided into smaller
blocks to do the finite size scaling analysis. (b) The block size dependence
of peak height of $\chi_{4}(t)$, $\chi_{4}^{P}$, for the passive system. (c)
The finite size scaling collapse of the data in panel (b) to obtain the
dynamic heterogeneity length scale, $\xi_{D}$ (Eq. 11). (d) Similar scaling
collapse, as in (c), of $\chi_{4}^{P}$ for an active system with $f_{0}=2.5$.
(e) Block size dependence of the self part of the van Hove correlation
function, $G_{s}(x,\tau_{\alpha})$, computed at structural relaxation time. At
small $L_{B}$, $G_{s}(x,\tau_{\alpha})$ shows non-Gaussian behaviour and
becomes progressively Gaussian with increasing $L_{B}$.
$G_{s}(x,\tau_{\alpha})$ becomes Gaussian when $L_{B}$ exceeds a crossover
length scale that gives $\xi_{D}$. (f) Binder cummulant of the block averaged
$G_{s}(x,\tau_{\alpha})$ at different $T$. (g) Scaling collapse of the Binder
cummulant computed in panel (f) to obtain $\xi_{D}$. (h) $\xi_{D}$ as a
function of scaled temperature, $T/T_{K}$, where $T_{K}$ is the Kauzmann
temperature at various activities. Note that the same length scale collapses
both the $\chi_{4}^{P}$ and Binder cummulant of van Hove function data. Inset:
Comparison of different estimates of $\xi_{D}$, obtained via block analysis of
$\chi_{4}^{P}$ and using spatial displacements correlation function,
$g_{uu}(r,\tau_{\alpha})$ (see text for details). The agreement is very good
suggesting the robustness of the methods as well as the estimate of $\xi_{D}$.

### Dynamic Heterogeneity Length Scale

Having discussed the non-trivial effects of activity on the higher order
correlation functions like $\chi_{4}(t)$, we now turn our focus to compute the
dynamic length scale, $\xi_{D}$, associated with DH itself. Reliable estimate
of $\xi_{D}$ in model active glassy systems is rare Ghoshal and Joy (2020); in
this work, we compute and compare $\xi_{D}$ via three different methods - (i)
finite size block scaling of $\chi_{4}$, (ii) block analysis of van Hove
function and (iii) spatial correlation of displacements at structural
relaxation time (see SM for further details): the agreement among them shows
reliability of our results. We have used the method of block analysis that is
developed to study the growing dynamical correlations in equilibrium
supercooled liquids in Chakrabarty  _et al._ (2017); Bhowmik  _et al._ (2018).
The simplicity and the novelty of the method makes it immediately useful for
our current study as we need to compute $\xi_{D}$ not only as a function of
$T$ but also as a function of changing $f_{0}$. Block analysis made it much
easier for us to do these measurements for different values of $f_{0}$ with
much less computational difficulty. Briefly, in block analysis method, a
somewhat bigger system size with linear size $L$ is simulated at the desired
state points and then the whole system is divided into smaller blocks of
length $L_{B}=L/n$, where $n$ is the number of blocks chosen, as schematically
shown in Fig. 3(a). We typically choose largest blocks to have the linear size
$L/3$, to ensure that other boundary effects do not influence the results. It
is important to note that in this analysis, the blocks can be thought of in a
Grand-Canonical ensemble, thus ensuring all possible fluctuations e.g. number
of particles, density, and, $T$ contribute in $\chi_{4}(t)$ properly
Chakrabarty  _et al._ (2017); Flenner and Szamel (2010). Now, the dynamic
susceptibility is computed for each such blocks as

$\chi_{4}(L_{B},t)=\frac{NL_{B}^{3}}{L_{0}^{3}}\left<[Q(L_{B},t)-\left<Q(L_{B},t)\right>]^{2}\right>,$
(9)

where

$Q(L_{B},t)=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\left<w(\mid{\textbf{r}_{j}(t)-\textbf{r}_{j}(0)}\mid)\right>.$
(10)

$N_{B}$ is the number of blocks with size $L_{B}$, $n_{i}$ is the number of
particles in the $i^{th}$ block at time $t=0$, and the window function
$w(x)=\Theta(a-x)$ where $\Theta$ is the Heaviside step function and the value
of the parameter $a$ is chosen to remove the decorrelation arising from
vibrations of particles inside the cages formed by their neighbours. The
dynamic length scale is obtained from a detailed finite-size scaling (FSS)
analysis of the dependence of the four-point dynamic susceptibility on the
block size by assuming the following scaling ansatz,

$\chi_{4}^{P}(L_{B},T)=\chi_{0}(T)\mathcal{F}\left(\frac{L_{B}}{\xi_{D}}\right)$
(11)

where, $\chi_{0}(T)=\lim_{L_{B}\to\infty}\chi_{4}^{P}(L_{B},T)$. The data for
all temperatures can be collapsed to a master curve using the two parameters,
$\chi_{4}^{P}(\infty,T)$ and $\xi_{D}$.

Fig. 3(b) shows $\chi_{4}^{P}(L_{B},T)$ as a function of block length $L_{B}$
for various $T$ for the passive system ($f_{0}=0$) and Fig. 3(c) shows the
data collapse, Eq. (11), to obtain $\xi_{D}$ as shown in Fig. 3(h). In Fig.
3(d), we show similar data collapse for an active system with $f_{0}=2.5$.
Similar analyses for all other values of $f_{0}$ are reported in SM along with
other relevant discussion. The excellent quality of the data collapse observed
in the analyses suggests the reliability of the extracted $\xi_{D}$. However,
in the absence of previous data to compare with, we reconfirm the reliability
of the estimates of $\xi_{D}$ via two other methods as mentioned before.

It was demonstrated in Ref. Bhowmik  _et al._ (2018) that $\xi_{D}$ can be
obtained by analyzing the non-Gaussian behaviour of the self part of the van
Hove function, $G_{s}(x,\tau_{\alpha})$ (defined in Method section and
elaborated in SM). In this method, one looks at the behaviour of the van Hove
function as a function of spatial coarse-graining scale, $L_{B}$, similar to
the block analysis with the main difference, that here one computes the van
Hove function of the displacement fields of the centre of mass (CoM) of the
individual blocks. Intuitively, as the coarse-graining volume gradually
increases, the system will eventually become spatially homogeneous once the
scale of coarse-graining exceeds the typical dynamic heterogeneity length
scale. Fig. 3(e) shows $G_{s}(x,\tau_{\alpha})$ with increasing $L_{B}$
illustrating the change of nature of $G_{s}(x,\tau_{\alpha})$, from non-
Gaussian to Gaussian form, as $L_{B}$ increases. $G_{s}(x,\tau_{\alpha})$
becomes Gaussian at a crossover scale $L_{B}$ which is different at each $T$
and the $T$-dependence of this crossover scale has the same $T$-dependence of
$\xi_{D}$ Bhowmik  _et al._ (2018). To estimate this crossover length scale in
an unbiased way, we computed the Binder cumulant of the distribution as shown
in Fig. 3(f). It is clear that with decreasing $T$, the Binder cummulant
increases strongly as $L_{B}$ decreases. As $L_{B}\to\infty$, Binder cumulant
should reach 0 for at all $T$. We can do a one parameter scaling collapse
(Fig. 3g) using the scaling ansatz
$B(L_{B},T)=\mathcal{G}\left(L_{B}/\xi_{D}\right)$ where $\xi_{D}$ is the
crossover length scale, same as the one used in collapsing the
$\chi_{4}^{P}(L_{B},T)$ data in the top panels. The quality of the data
collapse shows that these two methods yield the same $\xi_{D}$.

Finally, we have extracted the same length scale via yet another way used
earlier in (Tah and Karmakar, 2020), by studying the spatial correlation of
the displacement fields of the particles over structural relaxation time scale
characterized by the correlation function $g_{uu}(r,\tau_{\alpha})$, where $r$
is the radial distance (see the definition in the Method section and further
details in SM). In the inset of Fig. 3(h), the comparison of the length scales
are shown. The excellent agreement again highlights the reliability of the
obtained length scale. Fig. 3(h) shows the growth of $\xi_{D}$ for various
activities as a function of temperature. Note the dramatic growth of
heterogeneity length scale with increasing activity, which is one of the main
results of this work and we are going to discuss the salient features of this
observation in the next section.

Figure 4: Nontrivial effects of activity on glassy dynamics. (a-b) Qualitative
measures of CRR showing the fraction of fast relaxing particles, CRR provides
a typical estimate of DH, for a passive (a) and an active system (b) at same
$\tau_{\alpha}=10^{3}$ and shows larger DH in an active system. (c) $\xi_{D}$
as a function of $\tau_{\alpha}$ at different $f_{0}$. The data shows that for
a given $\tau_{\alpha}$, a system with higher $f_{0}$ has larger $\xi_{D}$.
(d) Fit of simulation data (symbols) for $\xi_{D}$ with the power-law
prediction of MCT (lines); the exponent $\gamma_{D}$ is expected to be
constant for an effective equilibrium-like scenario, however, we find that
$\gamma_{D}$ almost linearly increases with $f_{0}$ (inset). (e) In
simulation, we choose a set of $T$ and $f_{0}$ such that $T_{eff}$ remains
similar. (f) $\chi_{4}(t)$ for the parameters as in (e); system with larger
$f_{0}$ has higher $\chi_{4}^{P}$. Inset: All these systems have similar
$\tau_{\alpha}$ as seen from the plot of $Q(t)$. (g) $\chi_{4}(t)$, obtained
from active-IMCT, also has similar behavior as in simulation, i.e., higher
$\chi_{4}^{P}$ for larger $f_{0}$ when $Q(t)$, and hence $\tau_{\alpha}$, are
similar (inset). (h) Variation of van Hove function for different activity;
$T$ is choosen such that the structural relaxation time is same. The
increasing non-Gaussianity signifies the drastic increase in dynamic
heterogeneity with increasing activity. (i) Increasing spatial correlation, as
measured by the excess displacement-displacement correlation function,
$\Gamma(r,\tau_{\alpha})$, with increasing activity. These correlation
functions are again computed at the same relaxation time to highlight the
growth of spatial heterogeneity with increasing activity (see text for
details).

### Dramatic Effect of Activity on Dynamical Heterogeneity

We now discuss the dramatic effect of the non-equilibrium active force on the
DH in an active glass. Within the DH picture of equilibrium glassy dynamics,
it is beleived that $\xi_{D}$ controls the relaxation dynamics Berthier and
Biroli (2011); Ediger (2000). We have already shown that activity affects DH
in a way that is different from equilibrium-like behavior (Fig. 2c, 2d).
Active-IMCT also suggests that, the memory kernel within MCT has a nontrivial
effect of activity via the presence of the second term in Eq. (7). It is
easier to demonstrate the nontrivial effects of activity by considering
different systems whose parameters are chosen such that the typical relaxation
times of the systems remain same. In an equilibrium system, DH, as it controls
relaxation dynamics, should remain unchanged when $\tau_{\alpha}$ is same.
However, the scenario for an active glass, surprisingly, turns out to be
drastically different.

DH essentially refers to the coexistence of dynamic fast and slow moving
regions in the same system; a set of adjuscent particles with similar
$\tau_{\alpha}$ relaxes collectively and is known as cooperatively rearranging
region (CRR) that gives a measure of DH Royall and Kob (2017). Figs. 4(a) and
4(b) show typical estimates of CRRs of the fast moving particles in an
equilibrium system and an active system respectively, both having the same
$\tau_{\alpha}=10^{3}$ (see SM, Sec. 5 for details). Whereas the CRR for the
former is disconnected at this particular $\tau_{\alpha}$, that for the latter
becomes system-spanning. This qualitative comparison of typical CRRs in the
two systems show greater DH in an active system at the same $\tau_{\alpha}$.
Now, for a more quantitative analysis, we show $\xi_{D}$ as a function of
$\tau_{\alpha}$ in Fig. 4(c). As is evident, at a given $\tau_{\alpha}$, a
system with higher activity (that is larger $f_{0}$) has larger $\xi_{D}$.
Note the drastic increase of $\xi_{D}$ for the largest activity case,
$f_{0}=2.5$, $\xi_{D}$ increases by a factor of $\sim 30$ compared to its high
temperature value, whereas in passive system, the growth is hardly by a factor
of $4-5$ in the comparable window of $\tau_{\alpha}$. Our results are
consistent with the large dynamic length scales found in simulations of a
driven granular fluid Avila  _et al._ (2014).

MCT predicts power-law divergence of $\xi_{D}$ at $T_{C}$, where $T_{C}$ is
the mode-coupling transition temperature for an active system, obtained from
fitting the data of $\tau_{\alpha}$. Fig. 4(d) shows simulation data of
$\xi_{D}$ as a function of $(T-T_{C})/T_{C}$ and the lines are fits with the
function $\xi_{D}\sim[(T-T_{C})/T_{C}]^{-\gamma_{D}}$; the data agree well
with the power-law form in the range of $T$ and $f_{0}$ where the simulations
are carried out. An effective equilibrium-like description implies same
$\gamma_{D}$ at different activity, as for the case of relaxation dynamics,
Eq. (1). Consistent with the $T$-dependence of $\chi_{4}^{P}$ (Fig. 2d), we
find $\gamma_{D}$ depends on activity and almost linearly increases with
$f_{0}$, as shown in the inset of Fig. 4(d); this result, again, highlights
the nontrivial activity-dependence of $\xi_{D}$.

To further illustrate the decoupling of DH and $\tau_{\alpha}$ in the presence
of activity, we now show simulation data for $Q(t)$ and $\chi_{4}(t)$ at
different $T$ and $f_{0}$ in Fig. 4(f). The parameters are chosen in such a
way that all the systems have very similar $\tau_{\alpha}$, as is avident from
the plots of $Q(t)$ shown in the inset of Fig. 4(f). Analysis of
$\tau_{\alpha}$ in terms of $T_{eff}$ implies all the systems have similar
$T_{eff}$, this is indeed true as shown in Fig. 4(e). If DH controls the
relaxation dynamics, as in an equilibrium system, one expects $\chi_{4}(t)$
for all the systems to overlap. However, as shown in Fig. 4(f), this is not
the case; $\chi_{4}(t)$ increases with increasing $f_{0}$. Similar behavior is
also found within our active-IMCT as shown in Fig. 4(g) with the inset showing
$Q(t)$ overlapping with each other for these chosen parameters. Thus, the
effects of activity on the DH and relaxation dynamics are different and the
unique relation between DH and $\tau_{\alpha}$ breaks down in an active glass.
Similar conclusions can be drawn by studying the effect of activity on other
dynamical quantities as well. Fig. 4(h) shows the increasing non-Gaussian
behaviour in $G_{s}(x,\tau_{\alpha})$ with increasing activity while the
structural relaxation time is kept constant, indicating growing DH with
increasing activity even when the structural relaxation time is same.
Similarly, Fig. 4(i) once more demonstrates that spatial correlation, computed
by the excess part of the displacement-displacement correlation function,
$\Gamma(r,\tau_{\alpha})$, increases markedly with increasing activity,
confirming the strong decoupling of DH and structural relaxation dynamics in
active glasses.

## Conclusions

To conclude, we have studied the dynamic heterogeneity in a glass-forming
system in the presence of non-equilibrium active forcing. We find that the
effect of activity on relaxation dynamics of an active glass can be described
in terms of a $T_{eff}$, much like an equilibrium system. However, the notion
that the effect of activity simply fluidizes the system with no other
significant changes in its dynamical properties is not accurate and an active
glass is fundamentally different from an equilibrium glassy system, the
nontrivial effect of activity manifests via the DH. Through the intensive
research of the last couple of decades, DH has emerged as one of the salient
features of glassy systems Ediger (2000); Berthier and Biroli (2011); Berthier
(2011) and is believed to control the relaxation dynamics. In the presence of
activity, systems with the same $\tau_{\alpha}$ can have different DH and
$\xi_{D}$, thus, masking the unique relation between DH and the relaxation
dynamics in an active glass and casting doubt on the presumed central role of
DH in glassy dynamics. It is possible that the giant number fluctuations,
which is a feature of active systems Narayan  _et al._ (2007); Ramaswamy  _et
al._ (2003); Ramaswamy (2010), affects DH and $\tau_{\alpha}$ differently. In
any case, our results expose the incompleteness of our understanding on how DH
affects relaxation dynamics and the fundamental properties of DH in glassy
systems in general.

We have shown that four-point susceptibility that quantifies DH in a glassy
system has a nontrivial behavior in an active glass. Considering glass
transition as a critical phenomenon, activity-dependence of the exponents of
the power-law behavior of both $\chi_{4}^{P}$ and $\xi_{D}$ posits a serious
challenge in developing a microscopic theory for such a system. Yet, the
active-IMCT does predict the nontrivial effects of activity on DH, which, the
theory suggests, comes from an activity-specific term within the memory
kernel, Eq. (7). To the best of our knowledge, exploration of the effect of
GNF on the glassy dynamics is scant Giavazzi  _et al._ (2017). Studying GNF in
a dense glassy system may not be straightforward, however, the vestige of the
effect should be there in the system. Due to its relevance in several
biological systems, active glass is an exciting system to investigate this
effect as an interesting future direction.

An important result of the current work is the dramatic growth of DH in an
active system, consistent with simulation results Avila  _et al._ (2014) on
driven granular fluids, which is a prototype active system Narayan  _et al._
(2007); Dijksman  _et al._ (2011); Ramaswamy (2010). In particular, we show
that dynamic heterogeneity grows dramatically with increasing activity even
when $\tau_{\alpha}$ remains constant. Recent experiments on confluent
cellular monolayer have found increased $\chi_{4}^{P}$ via the application of
GTPase RAB5A, even when the system fluidizes, i.e., $\tau_{\alpha}$ decreases
Malinverno  _et al._ (2017); Cerbino  _et al._ (2021). Application of RAB5A
may affect a confluent monolayer in several ways, such as affecting the
junctional proteins, and thus, changing inter-cellular interaction Sadhukan
and Nandi (2021) as well as affecting the motor proteins Palamidessi  _et al._
(2019). Via a combination of theory and experiment, it was convincingly shown
in Ref. Giavazzi  _et al._ (2017) that at least one important effect of RAB5A
on the monolayer is increasing cellular motility. Thus, these exciting
results, though puzzling from equilibrium glassy perspective where increased
$\chi_{4}^{P}$ implies higher $\tau_{\alpha}$, are entirely consistent in the
light of our theory.

Finally, we want to highlight that this work for the first time demonstrates
that the dynamical heterogeneity in active glassy systems is inherently
different from their equilibrium behaviour. Via a combination of theory and
simulation, we have shown that whereas the relaxation dynamics can be
described as an equilibrium-like behavior at a certain $T_{eff}$, the DH has
completely different behavior. This also highlights the decoupling of
relaxation dynamics and DH, an effect that has also been stressed in some
equilibrium systems Karmakar  _et al._ (2009). Studying the active glass,
thus, in general promises deeper understanding of the role of DH in glassy
dynamics.

## Materials and Methods

We study two different glass-forming models, one is Kob-Andersen binary
mixture with 80% A-type and 20% B-type particles interacting via the Lennard-
Jones pair potential

$\Phi\left(\frac{r_{ij}}{\sigma_{ij}}\right)=4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]$
(12)

where $r$ is the distance between two particles and the indices $i$ and $j$
can be A or B. The values of $\sigma_{ij}$ and $\epsilon_{ij}$ are chosen to
be: $\sigma_{AB}=0.8\sigma_{AA}$, $\sigma_{BB}=0.88\sigma_{AA}$,
$\epsilon_{AB}=1.5\epsilon{AA}$, and $\epsilon_{AB}=0.5\epsilon_{AA}$. We set
a cutoff in the potential at $r_{ij}=2.5\sigma_{ij}$ and shift it accordingly.
We set the unit of length and energy as $\sigma_{AA}=1$ and $\epsilon_{AA}=1$
and fix the overall density $\rho$ at 1.2. We use a quadratic polynomial to
make the potential and its first two derivatives smooth at the cutoff
distance. We chose another model that interpolates between finite-temperature
glasses and hard-sphere glasses and has been studied extensively in the
context of jamming physics. This is a $50:50$ binary mixture with diameter
ratio of $1.4$. Details of the potential for this model is given in the SM.

We have done all our simulations in constant number of particle ($N$), volume
($V$) and temperature ($T$) (NVT) ensemble. To keep the temperature constant
we use Gaussian operator-splitting thermostat Zhang (1997) throughout our
simulation. Following Ref. Mandal  _et al._ (2016), we introduce activity by
making 10% of the particles active (see Sm for details). In this work, we keep
$\tau_{p}=1.0$ fixed and study the effect of activity as a function of $f_{0}$
alone. The opposite trends of activity, discussed in Ref. Nandi  _et al._
(2018), is not important in our study since both these models have same trends
as a function of $f_{0}$.

## Dynamical quantities

To characterise the dynamics of glassy systems under various active forcing,
$f_{0}$, we measure the following dynamical quantities.

Mean-squared displacement: The mean-squared displacement (MSD) $\left<\Delta
r(t)^{2}\right>$ is defined as

$\left<\Delta
r(t)^{2}\right>=\left<\frac{1}{N}\sum_{i=1}^{N}\mid{\textbf{r}_{i}(t+t_{0})-\textbf{r}_{i}(t_{0})}\mid^{2}\right>.$
(13)

Overlap correlation function: The two point overlap correlation function
$Q(t)$ is defined as

$Q(t)=\langle\tilde{Q}(t)\rangle=\left<\frac{1}{N}\sum_{i=1}^{N}w(\mid{\textbf{r}_{i}(t+t_{0})-\textbf{r}_{i}(t_{0})}\mid)\right>.$
(14)

where

$w(r)=\begin{cases}1,\quad\text{if $r\leq a$}\\\
0,\quad\text{otherwise}\end{cases}$ (15)

$\left<...\right>$ represents the average over time origin $t_{0}$ as well as
10 independently generated configurations, $N$ is the total number of
particles in the system and the parameter $a$ is associated with the typical
vibrational amplitude of the caged particles. Throughout our analysis, we have
used the parameter $a=0.3$ and we have verified that our results do not change
with moderate variation in $a$.

Four-point dynamical susceptibility: The four-point dynamical susceptibility
$\chi_{4}(t)$ is defined in terms of the fluctuations of the two point overlap
correlation function as

$\chi_{4}(t)=N[\left<\tilde{Q}^{2}(t)\right>-\left<\tilde{Q}(t)\right>^{2}]$
(16)

Excess displacement correlation function: The excess displacement correlation
$\Gamma(r,\Delta t)$ defined as

$\Gamma(r,\Delta t)=\frac{g_{uu}(r,\Delta t)}{g(r)}-1,$ (17)

where $g(r)$ is the radial pair correlation function, defined as

$g(r)=\frac{1}{\rho N}\left<\sum_{i,j=1,j\neq
i}^{N}\delta(\textbf{r}+\textbf{r}_{i}(0)-\textbf{r}_{j}(0))\right>,$ (18)

and the spatial correlation of the particle displacements, $g_{uu}(r,\Delta
t)$, defined as

$g_{uu}(r,\Delta t)=\frac{\left<\sum_{i,j=1,j\neq i}^{N}u_{i}(t,\Delta
t)u_{j}(t,\Delta t)\delta(r-\mid\textbf{r}_{ij}(t)\mid)\right>}{4\pi
r^{2}\Delta rN\rho\left<u(\Delta t)\right>^{2}},$ (19)

where $u_{i}(t,\Delta t)=\mid\textbf{r}_{i}(t+\Delta t)-\textbf{r}_{i}(t)\mid$
is the scalar displacement of the particle between time $t$ and $t+\Delta t$.

van Hove correlation function: The self part of the van Hove correlation
function is defined as

$G_{s}(x,\tau)=\left<\delta[x-(x_{i}(\tau)-x_{i}(0))]\right>$ (20)

Author Contributions: SK designed the research project. KP performed the
research. SKN developed the active MCT formalism. KP, SKN and SK analysed the
data. SK and SKN wrote the paper.

Acknowledgments: We thank R. Cerbino, F. GIavazzi, G. Scita, Indrajit Tah and
Rituparno Mandal for fruitful discussions. We acknowledge support of the
Department of Atomic Energy, Government of India, under Project Identification
No. RTI 4007. SK acknowledges support from Core Research Grant CRG/2019/005373
from Science and Engineering Research Board (SERB) as well as Swarna Jayanti
Fellowship grants DST/SJF/PSA-01/2018-19 and SB/SFJ/2019-20/05. KP
acknowledges financial support from SB/SFJ/2019-20/05.

## References

  * Berthier and Biroli (2011) Ludovic Berthier and Giulio Biroli, “Theoretical perspective on the glass transition and amorphous materials,” Rev. Mod. Phys. 83, 587 (2011).
  * Angelini  _et al._ (2011) T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg, and D. A Weitz, “Glass-like dynamics of collective cell migration,” Proc. Natl. Acad. Sci. (USA) 108, 4717 (2011).
  * Garcia  _et al._ (2015) Simon Garcia, Edouard Hannezo, Jens Elgeti, J. F. Joanny, Pascal Silberzan, and Nir S. Gov, “Physics of active jamming during collective cellular motion in a monolayer,” Proc. Natl. Acad. Sci. (USA) 112, 15314 (2015).
  * Parry  _et al._ (2014) B. R. Parry, I. V. Surovtsev, M. T. Cabenn, C. S. O’Hern, E. R. Dufresne, and C. Jacobs-Wagner, “The bacterial cytoplasm has glass-like properties and is fluidized by metabolic activity,” Cell 156, 183 (2014).
  * Sadati  _et al._ (2014) Monirosadat Sadati, Amir Nourhani, Jeffrey J. Fredberg, and Nader Taheri Qazvini, “Glass-like dynamics in the cell and in cellular collectives,” WIREs Syst. Biol. Med. 6, 137 (2014).
  * Zhou  _et al._ (2009) E. H. Zhou, X. Trepat, C. Y. Park, G. Lenormand, M. N. Oliver, S. M. Mijailovich, C. Hardin, D. A. Weitz, J. P. Butler, and J. J. Fredberg, “Universal behavior of the osmotically compressed cell and its analogy to the colloidal glass transition,” Proc. Natl. Acad. Sci. (USA) 106, 10632 (2009).
  * Nishizawa  _et al._ (2017) Kenji Nishizawa, Kei Fujiwara, Masahiro Ikenaga, Nobushige Nakajo, Miho Yanagisawa, and Daisuke Mizuno, “Universal glass-forming behavior of in vitro and living cytoplasm,” Sci. Rep. 7, 15143 (2017).
  * Takatori and Mandadapu (2020) Sho C. Takatori and Kranthi K. Mandadapu, “Motility-induced buckling and glassy dynamics regulate three-dimensional transitions of bacterial monolayers,” arXiv: 2003:05618 (2020).
  * Gravish  _et al._ (2013) Nick Gravish, Daria Monaenkova, Michael A. D. Goodisman, and Daniel I. Goldman, “Climbing, falling, and jamming during ant locomotion in confined environments,” Proc. Natl. Acad. of Sci. (USA) 110, 9746 (2013).
  * Gravish  _et al._ (2015) Nick Gravish, Gregory Gold, Andrew Zangwill, Michael A. D. Goodisman, and Daniel I. Goldman, “Glass-like dynamics in confined and congested ant traffic,” Soft Matter 11, 6552 (2015).
  * Dijksman  _et al._ (2011) Joshua A. Dijksman, Geert H. Wortel, Louwrens T. H. van Dellen, Olivier Dauchot, and Martin van Hecke, “Jamming, yielding, and rheology of weakly vibrated granular media,” Phys. Rev. Lett. 107, 108303 (2011).
  * Dauchot  _et al._ (2005) O. Dauchot, G. Marty, and G. Biroli, “Dynamical heterogeneity close to the jamming transition in a sheared granular material,” Phys. Rev. Lett. 95, 265701 (2005).
  * J.Deseigne  _et al._ (2010) J.Deseigne, O. Dauchot, and H. Chaté, “Collective motion of vibrated polar disks,” Phys. Rev. Lett. 105, 135702 (2010).
  * Dreyfus  _et al._ (2005) R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, “Microscopic artificial swimmers,” Nature 437, 862 (2005).
  * Palacci  _et al._ (2010) J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, “Sedimentation and effective temperature of active colloidal suspensions,” Phys. Rev. Lett. 105, 088304 (2010).
  * Klongvessa  _et al._ (2019) Natsuda Klongvessa, F. Ginot, Christophe Ybert, C. Cottin-Bizonne, and Mathieu Leocmach, “Active glass: Ergodicity breaking dramatically affects response to self-propulsion,” Phys. Rev. Lett. 123, 248004 (2019).
  * Ramaswamy (2010) S. Ramaswamy, “The mechanics and statistics of active matter,” Annu. Rev. Condens. Matter Phys. 1, 323 (2010).
  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, “Hydrodynamics of soft active matter,” Rev. Mod. Phys. 85, 1143 (2013).
  * Ni  _et al._ (2013) R. Ni, M. A. Cohen Stuart, and M. Dijkstra, “Pushing the glass transition towards random close packing using self-propelled hard spheres,” Nat. Commun 4, 2704 (2013).
  * Berthier (2014) L. Berthier, “Nonequilibrium glassy dynamics of self-propelled hard disks,” Phys. Rev. Lett. 112, 220602 (2014).
  * Mandal  _et al._ (2016) Rituparno Mandal, Pranab Jyoti Bhuyan, Madan Rao, and Chandan Dasgupta, “Active fluidization in dense glassy systems,” Soft Matter 12, 6268–6276 (2016).
  * Flenner  _et al._ (2016) Elijah Flenner, Grzegorz Szamel, and Ludovic Berthier, “The nonequilibrium glassy dynamics of self-propelled particles,” Soft Matter 12, 7136 (2016).
  * Mandal  _et al._ (2017) Rituparno Mandal, Pranab Jyoti Bhuyan, Pinaki Chaudhuri, Madan Rao, and Chandan Dasgupta, “Glassy swirls of active dumbbells,” Phys. Rev. E 96, 042605 (2017).
  * Mandal  _et al._ (2020) Rituparno Mandal, Pranab Jyoti Bhuyan, Pinaki Chaudhuri, Chandan Dasgupta, and Madan Rao, “Extreme active matter at high densities,” Nat. Comm. 11, 2581 (2020).
  * Berthier and Kurchan (2013) Ludovic Berthier and Jorge Kurchan, “Non-equilibrium glass transitions in driven and active matter,” Nat. Phys. 9, 310 (2013).
  * Nandi and Gov (2017) Saroj Kumar Nandi and Nir S. Gov, “Nonequilibrium mode-coupling theory for dense active systems of self-propelled particles,” Soft Matter 13, 7609 (2017).
  * Liluashvili  _et al._ (2017) Alexander Liluashvili, Jonathan Ónody, and Thomas Voigtmann, “Mode-coupling theory for active brownian particles,” Phys. Rev. E 96, 062608 (2017).
  * Feng and Hou (2017) M. Feng and Z. Hou, “Mode coupling theory for nonequilibrium glassy dynamics of thermal self-propelled particles,” Soft Matter 13, 4464 (2017).
  * Nandi (2018) Saroj Kumar Nandi, “Activity-dependent self-regulation of viscous length scales in biological systems,” Phys. Rev. E 97, 052404 (2018).
  * Szamel (2016) G. Szamel, “Theory for the dynamics of dense systems of athermal self-propelled particles,” Phys. Rev. E 93, 012603 (2016).
  * Chaki and Chakrabarty (2020) S. Chaki and R. Chakrabarty, “Escape of a passive particle from an activity-induced energy landscape: emergence of slow and fast effective diffusion,” Soft Matter 16, 7103 (2020).
  * Merrigan  _et al._ (2020) Carl Merrigan, Kabir Ramola, Rakesh Chatterjee, Nimrod Segall, Yair Shokef, and Bulbul Chakraborty, “Arrested states in persistent active matter: Gelation without attraction,” Phys. Rev. Res. 2, 013260 (2020).
  * Bi  _et al._ (2016) Dapeng Bi, Xingbo Yang, M. Cristina Marchetti, and M. Lisa Manning, “Motility-driven glass and jamming transitions in biological tissues,” Phys. Rev. X 6, 021011 (2016).
  * Caprini  _et al._ (2020) Lorenzo Caprini, U. Marini Bettolo Marconi, C. Maggi, M. Paoluzzi, and Andrea Puglisi, “Hidden velocity ordering in dense suspensions of self-propelled disks,” Phys. Rev. Res. 2, 023321 (2020).
  * Vicsek  _et al._ (1995) Tamás Vicsek, András Czirók, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet, “Novel type of phase transition in a system of self-driven particles,” Phys. Rev. Lett. 75, 1226 (1995).
  * Narayan  _et al._ (2007) Vijay Narayan, Sriram Ramaswamy, and Narayanan Menon, “Long-lived giant number fluctuations in a swarming granular nematic,” Science 317, 105 (2007).
  * Ramaswamy  _et al._ (2003) S. Ramaswamy, R. Aditi Simha, and J. Toner, “Active nematics on a substrate: Giant number fluctuations and long-time tails,” Europhys. Lett. 62, 196 (2003).
  * Ediger (2000) M. D. Ediger, “Spatially heterogeneous dynamics in supercooled liquids,” Annu. Rev. Phys. Chem. 51, 99 (2000).
  * Berthier (2011) Ludovic Berthier, “Dynamic heterogeneity in amorphous materials,” Physics 4, 42 (2011).
  * Yamamoto and Onuki (1998) R. Yamamoto and Akira Onuki, “Heterogeneous diffusion in highly supercooled liquids,” Phys. Rev. Lett. 81, 4915 (1998).
  * Franz and Parisi (2000) S. Franz and G. Parisi, “On non-linear susceptibility in supercooled liquids,” J. Phys.: Condens. Matter 12, 6335 (2000).
  * Karmakar  _et al._ (2009) Smarajit Karmakar, Chandan Dasgupta, and Srikanth Sastry, “Growing length and time scales in glass-forming liquids,” Proc. Natl. Acad. Sci. (USA) 106, 3675 (2009).
  * Weeks  _et al._ (2000) Eric R. Weeks, J. C. Crocker, Andrew C. Levitt, Andrew Schofield, and D. A. Weitz, “Three-dimensional direct imaging of structural relaxation near the colloidal glass transition,” Science 287, 627 (2000).
  * Avila  _et al._ (2014) Karina E. Avila, Horacio E. Castillo, Andrea Fiege, Katharina Vollmayr-Lee, and Annette Zippelius, “Strong dynamical heterogeneity and universal scaling in driven granular fluids,” Phys. Rev. Lett. 113, 025701 (2014).
  * Malinverno  _et al._ (2017) Chiara Malinverno, Salvatore Corallino, Fabio Giavazzi, Martin Bergert, Qingsen Li, Marco Leoni, Andrea Disanza, Emanuela Frittoli, Amanda Oldani, Emanuele Martini, Tobias Lendenmann, Gianluca Deflorian, Galina V. Beznoussenko, Dimos Poulikakos, Kok Haur Ong, Marina Uroz, Xavier Trepat, Dario Parazzoli, Paolo Maiuri, Weimiao Yu, Aldo Ferrari, Roberto Cerbino, and Giorgio Scita, “Endocytic reawakening of motility in jammed epithelia,” Nat. Mat. 16, 587 (2017).
  * Cerbino  _et al._ (2021) Roberto Cerbino, Stefano Villa, Andrea Palamidessi, Emanuela Frittoli, Giorgio Scita, and Fabio Giavazzi, “Disentangling collective motion and local rearrangements in 2d and 3d cell assemblies,” Soft Matter , (Advanced Article) (2021).
  * Biroli  _et al._ (2006) Giulio Biroli, Jean-Philippe Bouchaud, Kunimasa Miyazaki, and David R. Reichman, “Inhomogeneous mode-coupling theory and growing dynamic length in supercooled liquids,” Phys. Rev. Lett. 97, 195701 (2006).
  * Kob and Andersen (1995) Walter Kob and Hans C. Andersen, “Testing mode-coupling theory for a supercooled binary lennard-jones mixture i: The van hove correlation function,” Phys. Rev. E 51, 4626 (1995).
  * S and Sastry (2003) Ashwin S S and Srikanth Sastry, “Low-temperature behaviour of the kob–andersen binary mixture,” J. Phys.: Condens. Matter 15, S1253 (2003).
  * Giavazzi  _et al._ (2017) Fabio Giavazzi, Chiara Malinverno, Salvatore Corallino, Francesco Ginelli, Giorgio Scita, and Roberto Cerbino, “Giant fluctuations and structural effects in a flocking epithelium,” J. Phys. D: Appl. Phys. 50, 384003 (2017).
  * Karmakar  _et al._ (2015) Smarajit Karmakar, Chandan Dasgupta, and Srikanth Sastry, “Length scales in glass-forming liquids and related systems: a review,” Rep. Prog. Phys. 79, 016601 (2015).
  * Ghoshal and Joy (2020) Dipanwita Ghoshal and Ashwin Joy, “Connecting relaxation time to a dynamical length scale in athermal active glass formers,” Phys. Rev. E 102, 062605 (2020).
  * Chakrabarty  _et al._ (2017) Saurish Chakrabarty, Indrajit Tah, Smarajit Karmakar, and Chandan Dasgupta, “Block analysis for the calculation of dynamic and static length scales in glass-forming liquids,” Phys. Rev. Lett. 119, 205502 (2017).
  * Bhowmik  _et al._ (2018) Bhanu Prasad Bhowmik, Indrajit Tah, and Smarajit Karmakar, “Non-gaussianity of the van hove function and dynamic-heterogeneity length scale,” Phys. Rev. E 98, 022122 (2018).
  * Flenner and Szamel (2010) Elijah Flenner and Grzegorz Szamel, “Dynamic heterogeneity in a glass forming fluid: Susceptibility, structure factor, and correlation length,” Phys. Rev. Lett. 105, 217801 (2010).
  * Tah and Karmakar (2020) Indrajit Tah and Smarajit Karmakar, “Signature of dynamical heterogeneity in spatial correlations of particle displacement and its temporal evolution in supercooled liquids,” Phys. Rev. Research 2, 022067 (2020).
  * Royall and Kob (2017) C Patrick Royall and Walter Kob, “Locally favoured structures and dynamic length scales in a simple glass-former,” J. Stat. Mech. , 024001 (2017).
  * Sadhukan and Nandi (2021) Souvik Sadhukan and Saroj Kumar Nandi, “Theory and simulation for equilibrium glassy dynamics in cellular potts model of confluent biological tissue,” Phys. Rev. E , (accepted) (2021).
  * Palamidessi  _et al._ (2019) Andrea Palamidessi, Chiara Malinverno, Emanuela Frittoli, Salvatore Corallino, Elisa Barbieri, Sara Sigismund, Galina V. Beznoussenko, Emanuele Martini, Massimiliano Garre, Ines Ferrara, Claudio Tripodo, Flora Ascione, Elisabetta A. Cavalcanti-Adam, Qingsen Li, Pier Paolo Di Fiore, Dario Parazzoli, Fabio Giavazzi, Roberto Cerbino, and Giorgio Scita, “Unjamming overcomes kinetic and proliferation arrest in terminally differentiated cells and promotes collective motility of carcinoma,” Nat. Mat. 18, 1252 (2019).
  * Zhang (1997) Fei Zhang, “Operator-splitting integrators for constant-temperature molecular dynamics,” J. Chem. Phys. 106, 6102 (1997).
  * Nandi  _et al._ (2018) Saroj Kumar Nandi, Rituparno Mandal, Pranab Jyoti Bhuyan, Chandan Dasgupta, Madan Rao, and Nir S. Gov, “A random first-order transition theory for an active glass,” Proc. Natl. Acad. Sci. (USA) 115, 7688–7693 (2018).