# Anti-site disorder and Berry curvature driven anomalous Hall effect in spin
gapless semiconducting Mn2CoAl Heusler compound

Nisha Shahi  School of Materials Science and Technology, Indian Institute of
Technology (Banaras Hindu University), Varanasi 221005, India  Ajit K. Jena
Indo-Korea Science and Technology Center (IKST), Bangalore 560065, India
Gaurav K. Shukla  School of Materials Science and Technology, Indian Institute
of Technology (Banaras Hindu University), Varanasi 221005, India  Vishal Kumar
School of Materials Science and Technology, Indian Institute of Technology
(Banaras Hindu University), Varanasi 221005, India  Shivani Rastogi  School of
Materials Science and Technology, Indian Institute of Technology (Banaras
Hindu University), Varanasi 221005, India  K. K. Dubey  School of Materials
Science and Technology, Indian Institute of Technology (Banaras Hindu
University), Varanasi 221005, India  Indu Rajput  UGC-DAE Consortium for
Scientific Research, Indore, India  Sonali Baral  UGC-DAE Consortium for
Scientific Research, Indore, India  Archana Lakhani  UGC-DAE Consortium for
Scientific Research, Indore, India  Seung-Cheol Lee  Indo-Korea Science and
Technology Center (IKST), Bangalore 560065, India  Satadeep Bhattacharjee
Indo-Korea Science and Technology Center (IKST), Bangalore 560065, India
Sanjay Singh  School of Materials Science and Technology, Indian Institute of
Technology (Banaras Hindu University), Varanasi 221005, India

###### Abstract

Spin gapless semiconductors exhibit a finite band gap for one spin channel and
closed gap for other spin channel, emerged as a new state of magnetic
materials with a great potential for spintronic applications. The first
experimental evidence for the spin gapless semiconducting behavior was
observed in an inverse Heusler compound Mn2CoAl. Here, we report a detailed
investigation of the crystal structure and anomalous Hall effect in the
Mn2CoAl using experimental and theoretical studies. The analysis of the high-
resolution synchrotron x-ray diffraction data shows anti-site disorder between
Mn and Al atoms within the inverse Heusler structure. The temperature-
dependent resistivity shows semiconducting behavior and follows Mooij’s
criteria for disordered metal. Scaling behavior of the anomalous Hall
resistivity suggests that the anomalous Hall effect in the Mn2CoAl is
primarily governed by intrinsic mechanism due to the Berry curvature in
momentum space. The experimental intrinsic anomalous Hall conductivity (AHC)
is found to be $\sim$35 S/cm, which is considerably larger than the
theoretically predicted value for ordered Mn2CoAl. Our first-principle
calculations conclude that the anti-site disorder between Mn and Al atoms
enhances the Berry curvature and hence the value of intrinsic AHC, which is in
a very well agreement with the experiment.

##  I Introduction

The Hall effect is defined as the generation of a transverse voltage across a
current-carrying conductor, when an external magnetic field is applied
perpendicular to the direction of the current flow [1, 2, 3, 4]. An additional
contribution in the transverse voltage was observed in materials with broken
time-reversal symmetry (TRS), known as anomalous Hall effect (AHE) that arises
due to the mutual interaction of magnetization and spin-orbit coupling (SOC)
[2, 5]. The AHE was discovered more than a century ago, attained a vast
interest in current years due to its role in understanding the fundamental
physics [2, 6, 7, 8, 9, 10, 11] and potential for applications in spintronics-
based data storage devices and Hall sensors [12, 13, 14, 15, 16, 17, 18, 19].
Theoretical studies suggest that the AHE could originate from both extrinsic
and intrinsic mechanisms [6, 20, 99]. The extrinsic mechanism receives
contribution from skew scattering and side jump mechanisms that are related to
the asymmetric scattering and transverse shift of propagation direction of
spin-polarized charge carriers, respectively [22, 23]. On the other hand, the
intrinsic mechanism is related to the momentum space Berry curvature
associated with the electronic band structure [20, 26, 24, 25]. The momentum
space Berry curvature acts as a fictitious magnetic field and introduces an
anomalous velocity in addition to the group velocity of the electronic wave,
perpendicular to the electric field direction, which creates intrinsic
anomalous Hall conductivity (AHC) in the system [2, 20, 24].

The large AHC due to non-vanishing momentum space Berry curvature has been
observed in various kinds of systems such as Nd2Mo2O7 [27], (In, Mn)As [28],
(Ga, Mn)As [28], Mn3Sn [29], Mn3Ge [8], Fe3Sn2 [30], Co3Sn2S2 [31], CoNb3S6
[32], and SrIrO3 [33]. Besides these materials, Heusler alloys emerged as
promising candidate for the realization of large Berry curvature originated
from their peculiar band structure due to integrated effect of crystal
symmetry and SOC [34, 35, 36, 37, 38]. It has been observed that any small
perturbation in the electronic band structure for example a local disorder may
affect the momentum space Berry curvature and therefore the intrinsic AHC in
the system [40, 41, 42, 43, 39]. Recently, an enhanced intrinsic AHC has been
observed in Fe2-based inverse Heusler compounds due to the anti-site disorder
[44]. Also, a Berry curvature driven enhanced intrinsic AHC has been observed
in Co2FeAl Heusler compound, originated via anti-site disorder between Fe and
Al atoms [45].

The discovery of the spin gapless semiconducting (SGS) behavior in the Mn2CoAl
Heusler compound puts this material forward as an important candidate for
technological application in the field of spintronics [47, 48, 13, 16].
Experimental studies on thin film as well as bulk systems of Mn2CoAl reflect
that most of the systems crystallize with compositional and/or anti-site
disorder and the reported value of AHC is not in an agreement with the theory
[47, 49, 50, 51, 80, 58, 52, 53, 54, 55, 56, 59]. A theoretical investigation
considering anti-site disorder suggests an enhancement in AHC, necessitates a
detailed experimental and theoretical investigation to understand the origin
of AHC in Mn2CoAl Heusler compound [60].

In the present manuscript, we investigate AHE in the Mn2CoAl Heusler compound
experimentally as well as theoretically. The high-resolution synchrotron x-ray
diffraction (SXRD) analysis exhibits 25% anti-site disorder between Mn and Al
atoms. The scaling behavior of anomalous Hall resistivity reveals the major
contribution of intrinsic mechanism in the observed AHE. Experimentally, we
found a higher value of intrinsic AHC as compared to the value reported for an
ordered Mn2CoAl from theory [60, 47]. Our theoretical calculations show that
the anti-site disorder between Mn and Al atoms enhances the Berry curvature
and hence the value of intrinsic AHC, which is in well agreement with the
experiment.

##  II Methods

Figure 1: (a) Rietveld profile fitting of the room temperature (RT)
synchrotron x-ray diffraction (SXRD) pattern of Mn2CoAl compound by
considering the XA-type ordered structure. The observed (Iobs), calculated
(Ical) and difference between observed and calculated profiles (Iobs-Ical) are
shown by black circle, red continuous line and green continuous line,
respectively. The blue tick bars and the $\chi^{2}$ value indicate the Bragg
peak positions and the goodness of fit, respectively. (b) Rietveld profile
fitting of the RT SXRD pattern of Mn2CoAl compound by considering the L21 type
ordered structure. Insets in (a) and (b) show an enlarged view of (111) and
(200) peaks. (c) Simulated SXRD patterns considering the different percentage
of anti-site disorder between Mn4c and Al atoms. Inset shows the variation in
(111) peak intensity with the different percentage of anti-site disorder. (d)
Rietveld profile fitting of the RT SXRD pattern of Mn2CoAl by considering 25%
anti-site disorder between Mn4c and Al atoms. Inset shows an enlarged view
around (111) and (200) peaks.

Polycrystalline Mn2CoAl Heusler compound is prepared by the arc-melting method
[61, 62, 63]. The constituent elements (purity higher than 99.9%) of the
intermetallic system are melted in a water-cooled copper hearth under argon
atmosphere (purity better than 99.999%). The sample is remelted 5 times to
ensure the homogeneous mixing of the constituents. The chemical composition is
verified using the energy dispersive analysis of x-ray technique. The average
composition is found to be Mn1.99Co1.00Al0.99, which corresponds to Mn2CoAl. A
small piece of the sample is powdered and room temperature SXRD is performed
in PETRA III DESY, Germany (wavelength $\sim$0.207
$\mathrm{\SIUnitSymbolAngstrom}$) for structural analysis. The direct current
(DC) magnetization measurements are performed using a 9 Tesla (T) physical
property measurement system (PPMS) of quantum design(QD). Transport
measurements are carried out using alternating current transport (ACT) option
of 9 T PPMS of QD. Resistivity and magneto-resistance measurements are
performed by the four-probe method, whereas the five-probe method is used for
the Hall measurement. To remove the longitudinal resistivity contribution in
the Hall data due to voltage probe misalignment, we have anti-symmetrized the
Hall resistivity data by using the formula $\rho$H = [$\rho$H(+H) -
$\rho$H(-H)]/2.

The spin-polarized Kohn-Sham Hamiltonian was solved within the framework of
pseudo-potentials (PP) and plane waves as implemented in Quantum ESPRESSO (QE)
density-functional theory (DFT) package [64]. Exchange-correlation part of the
above Hamiltonian is approximated by generalized gradient approach [65]
through ONCV Vanderbilt PPs [66]. The kinetic energy cutoff of 80 Ry is used
to fix the number of plane waves. The same Gaussian smearing value (0.01 Ry)
is used both for the self-consistent (SC) and non-self-consistent (NSC)
calculations to carry out electronic integration over the Brillouin zone (BZ).
A tight energy threshold ($10^{-8}$ Ry) is considered for the SC energy
calculations. We have used the WANNIER90 tool (implemented within QE) in order
to compute the Wannier interpolated bands, Berry curavture and AHC [64, 67,
68, 69]. We note that spin-orbit coupling (SOC) is introduced in all the AHC
related calculations. Collinear spin-polarized calculations are also performed
to compare some of our results with earlier reports and are mentioned at the
appropriate places. The Monkhorst-Pack k-grid of $8\times 8\times 8$ of the BZ
is considered in all the calculations, viz.: SC, NSC, and WANNIER90. We found
that the use of transition metal-$d$ orbitals as the projections in the
WANNIER90 calculations provide the good interpolation. A denser BZ k-grid of
$75\times 75\times 75$ is taken to calculate the intrinsic AHC. Through the
adaptive refinement technique, a further fine mesh of $5\times 5\times 5$ is
added around the points wherever the mode of the Berry curvature
($\absolutevalue{\Omega(\textbf{k})}$) exceeds 100 bohr2. The calculations are
carried out using experimental refined lattice constant of 5.858 Å.

##  III Results and discussion

###  a Structural characterization

To investigate the crystal structure of the Mn2CoAl Heusler compound, the SXRD
pattern was collected at room temperature. Mn2YZ Heusler compounds generally
exhibit XA-type (prototype Hg2CuTi, space group F$\overline{4}$3m) crystal
structure, if atomic number of Y is higher than the Mn like Mn2CoAl [70].
Therefore, in the primary step, we performed Rietveld refinement of SXRD data
by assuming the XA-type crystal structure (Fig. 1(a)), using FULLPROF software
package [71]. In the refinement, the Mn atoms were considered at 4a (0, 0, 0)
and 4c (0.25, 0.25, 0.25) Wyckoff positions and will be denoted as Mn4a and
Mn4c, respectively. The Co and Al atoms were considered at 4b (0.50, 0.50,
0.50) and 4d (0.75, 0.75, 0.75) Wyckoff positions, respectively. It can be
noticed from Fig. 1(a) that the observed and calculated peak profiles do not
match with each other as the calculated intensity of (111) superlattice peak
is much larger than the observed intensity (depicted in the inset of Fig.
1(a)). This mismatch between the calculated and observed intensities of Bragg
peaks indicates that Mn2CoAl does not crystallize in an ordered XA-type
structure. In the next step, we tried the Rietveld refinement with the ordered
L21 type (Fm$\overline{3}$m) crystal structure (Fig. 1(b)) [63, 70], which
also fails to account the intensities of both the (111) and (200) superlattice
reflections as depicted in the Fig. 1(b) and its inset.

In literature it has been suggested that the most of the inverse Heusler
compounds crystallize with anti-site disorder [72, 73, 70, 44, 63]. The anti-
site disorder in the Heusler compounds can be generally captured by analyzing
the intensities of ordering dependent (111) and (200) superlattice reflections
[74]. For Mn2CoAl compound, the structure factor for (111), (200) and (220)
reflections can be written as-

$F_{111}=4[(f\textsubscript{Mn\textsubscript{4a}}-f\textsubscript{Co})-i(f\textsubscript{Mn\textsubscript{4c}}-f\textsubscript{Al})]$
(1)

$F_{200}=4[(f\textsubscript{Mn\textsubscript{4a}}+f\textsubscript{Co})-(f\textsubscript{Mn\textsubscript{4c}}+f\textsubscript{Al})]$
(2)

$F_{220}=4[(f\textsubscript{Mn\textsubscript{4a}}+f\textsubscript{Co})+(f\textsubscript{Mn\textsubscript{4c}}+f\textsubscript{Al})]$
(3)

Eq. (1) suggests that the intensity of (111) superlattice reflection is due to
the difference between the atomic scattering factor (SF) of Mn4c and Al atoms
as the difference between the atomic SF of Mn4a and Co atoms is negligible
(since both are neighbouring elements in the periodic table) and also the
intensity of (111) peak may change, if the system has the anti-site disorder
between the Mn4c and Al atoms. From Eq. (2) and Eq. (3), it is obvious that
the intensities of (200) and (220) reflection will be unaffected from the
anti-site disorder between the Mn4c and Al atoms. Hence, the mismatch between
the observed and calculated intensities of (111) peak in the Rietveld fitted
SXRD pattern with XA-type crystal structure suggests the presence of anti-site
disorder between Mn4c and Al atoms.

Now, to find out the amount of anti-site disorder, we simulated the SXRD
pattern using VESTA software [75] by incorporating the different percentage of
anti-site disorder between Mn4c and Al atoms within XA-type crystal structure
as shown in Fig. 1(c). The inset of Fig. 1(c) represents an enlarged view of
the change in intensity of (111) superlattice reflection with anti-site
disorder. We found that for 25% anti-site disorder between Mn4c and Al atoms,
the intensity ratio of (111) superlattice reflection to the (220) fundamental
reflection is in good agreement with the experiment. It is worthwhile to
mention here that SXRD patterns simulated using considering other types of
disorder like disorder between Mn4a and Mn4c, between Mn4a and Al and between
Co and Al do not provide intensity ratio of Bragg reflections as observed in
experimental SXRD pattern. Therefore, based on the information obtained from
the Eq. (1-3) and simulated SXRD patterns, in the final step we performed the
Rietveld refinement considering the 25% anti-site disorder between the Mn4c
and Al atoms. A very good match between the observed and calculated peak
profile was observed as depicted in Fig 1(d). The inset of Fig. 1(d) shows an
enlarged view of (111) and (200) peaks. The refined lattice parameter was
obtained $\sim$5.858 $\mathrm{\SIUnitSymbolAngstrom}$, which is in well
agreement with the literature [51].

###  b Magnetization and resistivity measurements

The temperature variation of magnetization measured in a temperature range of
2-400 K under the magnetic field of 1 T, shown in Fig. 2(a). The magnetization
decreases with increasing temperature, as expected for ferrimagnetic system
[76]. The observed behavior is similar to the earlier reports [51, 47]. Fig.
2(b) presents the field-dependent magnetic isotherms (M(H) curves) at 2 K and
300 K. The negligible hysteresis is due to the soft magnetic nature of the
Mn2CoAl compound. The M(H) curves show similar behavior in the entire
temperature range of 2-300 K. To calculate the Curie temperature (Tc) of
present compound, the saturation magnetization (Ms) vs temperature (Ms(T))
data obtained from the M(H) curves (black spheres in the inset of Fig. 2(b))
were fitted by the empirical law Ms= M0[1-(T/Tc)2]1/2 as shown by the red
curve in the inset of Fig. 2(b) [80]. From this fitting, the Tc was found to
be 720 K, which is in well agreement with the value reported in the literature
[80, 51, 47]. The saturation magnetization of present compound is found to be
$\sim$1.9 $\mu$B/f.u. at 2 K, which is in well agreement with the literature
[47, 80].

The monotonous decrease in the longitudinal resistivity ($\rho$xx) with
temperature (Fig. 2(c)) illustrates the semiconducting behavior of the present
compound. A slight decrease in the value of $\rho$xx with increasing magnetic
field indicates the presence of small negative magnetoresistance in the
system. The resistivity at 2 K (that corresponds to residual resistivity) is
$\sim$387 $\mu$$\mathrm{\SIUnitSymbolOhm}$-cm, which is comparable to the
reported one [51]. The residual resistivity ratio [$\rho$xx(300 K)/$\rho$xx(2
K)] was found $\sim$ 0.92, which is comparable to the value in the literature
[51, 47] and also suggest the presence of disorder in the Mn2CoAl compound.

Figure 2: (a) Temperature-dependent magnetization curve at 1 T magnetic field.
(b) The isothermal field-dependent magnetization at different temperatures.
Inset shows the fitting of saturation magnetization (Ms) vs temperature (T)
with empirical relation Ms=M0[1-(T/Tc)2]1/2. The black dots and the red
continuous line in the inset represent the Ms vs T data and fitted curve,
respectively. (c) Resistivity ($\rho$xx) vs temperature (T) plot at different
magnetic fields. Insets (i) and (ii) show the fitting of magneto-conductivity
with field by Hikami-Larkin-Nagaoka equation and the variation of temperature
coefficient of resistivity ($\alpha$) with temperature (T), respectively. (d)
Fitting of conductivity ($\sigma$xx) vs temperature (T) data in the
temperature range below 150 K and in higher temperature range above 150 K at
zero field.

A theoretical calculation considering anti-site disorder into account predicts
a half-metallic character of Mn2CoAl Heusler compound in contrast to the SGS
behavior [60, 77]. It is interesting to note that the Mn2CoAl compound shows
the semiconducting behavior (Fig. 2(c)) although a large (25%) anti-site
disorder is observed. It has been suggested in literature that if the system
has weak localization (WL) then inelastic scattering increases with increasing
the temperature, which may destroy the phase coherence, and hence the
resistivity decreases with increasing temperature [78]. To check the
possibility of WL in Mn2CoAl, the magneto-conductivity (MC) i.e. difference of
conductivity between finite magnetic field and zero field at room temperature
is fitted by the Hikami-Larkin-Nagaoka (HLN) equation [79]. Although the HLN
model explains the WL effect in two-dimensional systems but it has been
observed that the several three dimensional (3D) systems also follow this
model [80, 81, 82, 83, 84, 85]. The HLN equation is given as-

$$\Delta$G(H)=\frac{$\gamma$e^{2}}{$\pi$h}\left[$\psi$\left(\frac{h}{4eHL\textsubscript{$\phi$}^{2}}+\frac{1}{2}\right)-ln\left(\frac{h}{4eHL\textsubscript{$\phi$}^{2}}\right)\right]$
(4)

where $\Delta$G is MC, $\psi$ is the digamma function, L$\phi$ is the phase
coherence length, e is the electronic charge, h is the Planck’s constant, H is
the magnetic field, and $\gamma$ is a prefactor. The sign of prefactor
$\gamma$ decides whether there is WL or weak anti-localization in the system.
The positive MC (black dots in the inset (i) of Fig. 2(c)), fitted (red solid
line Fig. 2(c)) by above equation with positive $\gamma$ and a large
logarithmic increment for lower magnetic field region indicates the presence
of WL in Mn2CoAl [86]. Moreover, we may see that the magnetoconductivity in
the low field region follows H0.5 power law variation, which is well
established for WL in 3D systems [87, 88, 89, 90, 91], also validate the
presence of WL in bulk Mn2CoAl. Further, Mooij established a correlation
between temperature coefficient of resistivity and residual resistivity for
disordered metal which describes that a disordered metal may exhibit
semiconducting behavior, if the residual resistivity exceeded from 150
$\mu$$\mathrm{\SIUnitSymbolOhm}$-cm [92]. For the Mn2CoAl, a larger value of
the residual resistivity ($\sim$387 $\mu$$\mathrm{\SIUnitSymbolOhm}$-cm) is
observed and follows the Mooij’s criteria.

The temperature coefficient of resistivity ($\alpha$) for Mn2CoAl changes from
-0.9$\times$10-9 to -1.9$\times$10-9 $ohm{\mathrm{\SIUnitSymbolOhm}}$-m/K in
the temperature range of 2-350 K with a slope change around 150 K (shown in
the inset (ii) of Fig. 2(c)). To understand the origin of slope change in
$\alpha$, we analyzed the longitudinal conductivity ($\sigma$xx = 1/$\rho$xx)
vs temperature data in the absence of magnetic field, and we found that
$\sigma$xx is proportional to the T1/2 below the temperature 150 K as shown in
Fig. 2(d). This conductivity behavior follows the interaction theory at low-
temperature according to Kaveh and Mott [93], which might be due to anti-site
disorder as supported by our SXRD analysis. Above the 150 K, the conductivity
variation with temperature follows the semiconducting-like character and is
well approximated by the relation ($\sigma$xx = $\sigma$0 \+ $\sigma$a
exp${}^{-E_{g}/K_{B}T}$) as shown in Fig. 2(d) [94]. The calculated band gap
(Eg) from the best fitting was found to be $\sim$50 meV, that is of the same
order as reported for the thin film of Mn2CoAl [80]. Recently, the appreciable
variation in Seebeck effect around 150 K has been observed in the polycrystal
Mn2CoAl [47]. This shows that our prepared compound is similar to the reported
one and the disorder present in Mn2CoAl compound is inherent and does not
change with the sample preparation conditions.

###  c Anomalous Hall measurement

Figure 3: (a) Field-dependent Hall resistivity curves at different
temperatures and the inset shows the fitting of Hall data in higher field
region $>$ 0.6 T at 2 K. (b) Variation of ordinary Hall coefficient with
temperature and the inset shows carrier concentration vs temperature plot. (c)
Field-dependent Hall conductivity curves at different temperatures. The inset
indicates zero field extrapolation of high field Hall conductivity data at 2
K, shown by red dotted line. (d) Fitting of ratio of anomalous Hall
resistivity and longitudinal resistivity ($\rho$AH/$\rho$xx) vs $\rho$xx data.
(e) Different contributions in $\rho$AH with temperature are plotted on the
same scale. (f) Linear fitting of log($\rho$AH) vs log($\rho$xx) data. (g)
Variation of anomalous Hall scaling coefficient (SH) with temperature. (h)
Variation of AHC with temperature. (i) Variation of AHC with longitudinal
conductivity.

Hall measurements were carried out in the temperature range of 2-350 K to
investigate the anomalous transport behavior of Mn2CoAl compound. In general
the total Hall resistivity ($\rho$H) consists of two parts namely ordinary
Hall and anomalous Hall and can be written as [1, 95]-

$$\rho$\textsubscript{H}={R_{0}}H+{R_{s}}{M_{s}}$ (5)

where R0 and Rs are the ordinary and anomalous Hall coefficients,
respectively. Ms corresponds to the saturation magnetization and RsMs
represents the magnitude of anomalous Hall resistivity ($\rho$AH).The field-
dependent $\rho$H data were recorded at different temperatures upto magnetic
field of 1.5 T as shown in the Fig. 3(a). The $\rho$H increases steeply upto
$\sim$ 0.6 T field, which is observed due to AHE. At the higher field region
($>$0.6 T), $\rho$H changes linearly and shows positive slope with magnetic
field, which is due to the ordinary Hall effect [37]. To separate out the
ordinary and anomalous Hall contributions, we performed fitting of the $\rho$H
vs H data by using the Eq. (5) in the higher field region ($>$0.6 T). The
fitting of the $\rho$H vs H data at 2 K, shown by the continuous red line in
the inset of Fig. 3(a), provides the value of R0 and RsMs that correspond to
the slope and intercept on the y axis of the fitted line. A similar analysis
(i.e. fitting) is performed at various temperatures in the temperature range
of 10-350 K and Fig. 3 (b) shows a temperature variation of R0 obtained from
the fitting. The positive value of R0 in whole temperature range (2-350 K)
reveals that holes are the dominating charge carriers in the transport. The
carrier concentration (n) determined by the expression n = $\frac{1}{eR_{0}}$,
was found to be $\sim$1.5$\times$1020cm-3 and $\sim$2$\times$1020cm-3 at 2 K
and 350 K, respectively. Thus, the value of n is nearly temperature
independent. The temperature variation of n for temperature range 2-350 K is
depicted in the inset of Fig. 3(b), which is similar to as reported in the
literature [47]. The Hall conductivity ($\sigma$H) has been extracted by using
the equation [41, 35]-

$$\sigma$\textsubscript{\tiny{H}}=\frac{$\rho$\textsubscript{\tiny{H}}}{$\rho$\textsubscript{\tiny{H}}^{2}+$\rho$\textsubscript{xx}^{2}}$
(6)

The field-dependent Hall conductivity at different temperatures is shown in
Fig. 3(c). The value of AHC at 2 K is determined by zero field extrapolation
of the higher field Hall conductivity curve as elucidated by red dotted line
in the inset of the Fig. 3(c). The obtained value of AHC is 27 S/cm, which is
comparable to the experimentally found value in the literature [47, 51]. From
the Berry curvature calculations, the AHC value was reported about 3 S/cm
[47], which is an order of magnitude lower than the value obtained from the
experiment ($\sim$ 27 S/cm). To understand this discrepancy, it is necessary
to address whether there is contribution in the AHC from an extrinsic
mechanisms, such as skew scattering and side jumps, or from intrinsic
mechanism due to the momentum space Berry curvature associated with electronic
band structure.

To calculate the separate contributions of extrinsic and intrinsic mechanism
in the total AHC, the $\rho$AH/$\rho$xx vs $\rho$xx data (black dots in Fig.
3(d)) was fitted (red line in Fig. 3(d)) using the following relation [96, 97,
6]-

$$\rho$\textsubscript{\tiny{AH}}/$\rho$\textsubscript{xx}=a+b$\rho$\textsubscript{xx}$
(7)

where the parameters a and b contain information about extrinsic skew
scattering and the combined effect of extrinsic side jump and intrinsic
contribution, respectively. From the $\rho$AH/$\rho$xx vs $\rho$xx data
fitting with above equation, we obtained the value of a $\sim$ -0.003 and b
$\sim$ 35 ($$\mathrm{\SIUnitSymbolOhm}$.cm)^{-1}$. The negative value of
coefficient a indicates that extrinsic skew scattering contribution is
opposite to both side jump and intrinsic contribution due to momentum space
Berry curvature. The value of parameter b contains the contributions in AHC
due to both side jump and momentum space Berry curvature. By using the
coefficients a and b, we calculated skew scattering term (a$\rho$xx) and
intrinsic plus side jump term ($b$\rho$\textsubscript{xx}{{}^{2}}$) and
plotted on the same scale as shown in Fig. 3(e). We can clearly see that side
jump together with intrinsic contribution dominates over the skew scattering
contribution in the overall AHC in the temperature range 2-350 K.

The dominating mechanism in the AHE can be alternatively evaluated by the
exponent $\beta$ using the scaling relation $\rho$AH$\propto$$\rho$xx$\beta$
[2, 45, 46]. If $\beta$ = 1, the skew scattering mechanism will dominant in
AHE and if $\beta$ = 2, the combination of side jump and intrinsic mechanism
will largely contribute to the AHE. The exponent $\beta$ determined by linear
fitting of log($\rho$AH) vs log($\rho$xx) is turned out to be 2, which also
supports that the AHC is dominated by the side jump and intrinsic mechanism
[45, 46]. It is not possible practically to separate out the side jump and
intrinsic contributions because both the contributions show quadratic
dependency on $\rho$xx [2]. However, the AHC due to side jump mechanism can be
approximated using an expression (e2/(ha)(Eso/EF), where Eso is the spin–orbit
interaction energy and EF is Fermi energy [98, 99]. The physical quantities e,
h and a are the electronic charge, Planck’s constant and lattice parameter,
respectively. Here, for Mn2CoAl the $\sigma$xx is $\sim$2.6$\crossproduct$103
S/cm. Therefore, we may conclude that Mn2CoAl is in the moderately dirty
regime for which Eso $\sim$$\hbar$/$\tau$, where $\hbar$ is reduced Planck’s
constant i.e. h/2$\pi$ and $\tau$ is the scattering time [2, 100]. The
calculated value of the ratio Eso/EF for Mn2CoAl is of the order of $10^{-3}$,
and hence in AHC, side jump contribution is negligible as compared to the
Berry curvature induced intrinsic contribution. It has been suggested that the
intrinsic AHC is nearly proportional to the magnetization, thus the scaling
coefficient SH = $\frac{$\sigma$_{\textsubscript{\tiny{AH}}}}{M_{s}}$ should
show temperature independent behavior, if there is dominant intrinsic
contribution to AHE [30, 46, 101]. Fig. 3(g) depicts that the SH for Mn2CoAl
is indeed independent of temperature, which indicates that intrinsic
contribution due to momentum space Berry curvature dominates over the
extrinsic scattering contributions in the AHE [30, 46, 101]. The AHC shows a
temperature independent behavior (Fig. 3(h)) and does not vary with $\sigma$xx
as well (Fig. 3(i)). This robust behavior of the AHC with the temperature and
$\sigma$xx further confirms that AHE in Mn2CoAl is mainly originated from the
intrinsic mechanism, and is thus dominated by Berry curvature in the momentum
space [102].

###  d First-principle calculations

Figure 4: Spin-polarized band structure of ordered Mn2CoAl. Down and up arrows
represent the spin-minority and majority electrons, respectively. Figure 5:
Top: Full electronic band structure (blue continuous line) and Wannier
interpolated band structure (red dashed line) of ordered Mn2CoAl. The Fermi
energy is set to 0 eV. Bottom: Distribution of the Berry curvature along high
symmetry path in the Brillouin zone. Figure 6: Top: Full electronic band
structure (blue continuous line) and Wannier interpolated band structure (red
dashed line) of disordered Mn2CoAl. The Fermi energy is set to 0 eV. Bottom:
Distribution of the Berry curvature along high symmetry path in the Brillouin
zone. Figure 7: (a) Density of states (DOS) and variation of intrinsic AHC of
ordered Mn2CoAl as function of the Fermi energy. The Fermi energy is set to 0
$eV$. (b) DOS and variation of intrinsic AHC for the disordered Mn2CoAl as
function of Fermi energy. The Fermi energy is set to 0 $eV$. Insets in (a) and
(b) show an enlarged view of intrinsic AHC around Fermi energy.

In literature, the value of AHC obtained from theory is very less as compared
to the experiment [47]. Here, from the experiment, we found that the Mn2CoAl
has 25% disorder between Mn4c and Al atoms. It is possible that this disorder
may play an important role in the experimentally observed AHC, therefore we
investigated the AHC in ordered and disordered Mn2CoAl using theoretical
calculations.

Before the discussion of AHC calculation, we first present our electronic
structure results. The spin-polarized band structure of ordered Mn2CoAl is
shown in Fig. 4. The spin-minority channel has a finite gap and the majority
one has zero gap, which suggests SGS feature of the Mn2CoAl compound. Spin
gapless semiconductors offer intriguing transport properties because both the
electrons and holes can be $100\%$ spin-polarized, so that spin can be
controlled via only a small applied external energy. The band structure for up
and down spins shown in Fig. 4 matches very well with the literature [47]. Our
calculated magnetic moment for the stoichiometric Mn2CoAl is found to be 2
$\mu$B/f.u. with anti-parallel coupling between Mn4a and Mn4c sites, which is
also in agreement with the literature [47].

Now we will discuss the theoretical results on anomalous Hall transport
obtained from Wannier90. We compared the full DFT band structure of ordered
Mn2CoAl with the Wannier interpolated one in Fig. 5 (top) and obtained a very
good interpolation. The distribution of calculated Berry curvature along the
high symmetry path is shown in Fig. 5 (bottom). An efficient first-principle
approach has been used in which the maximally localized Wannier functions are
first constructed from the Bloch states, on a relatively coarse $k$-grid. Then
the quantities of interest e.g. Berry curvature is interpolated onto a dense
$k$-mesh in calculating the intrinsic AHC as a Brillouin zone summation of the
Berry curvature over all occupied states [69, 103]-

$\sigma_{xy}(AHC)=\frac{e^{2}}{\hbar}\frac{1}{NV}\sum_{\textbf{k}\in(BZ)}(-1)\Omega_{xy}(\textbf{k})f(\textbf{k}),$
(8)

where the indices $x$ and $y$ are the Cartesian coordinates. $f(\textbf{k})$
is the Fermi distribution function, $\Omega_{xy}(\textbf{k})$ denotes the
Berry curvature for the wave vector k, $N$ is the number of electrons in the
crystal and $V$ is the cell volume.

Our theoretical calculation gives intrinsic AHC value $\sim$8.88 S/cm for the
ordered Mn2CoAl compound, which is slightly higher than the reported
theoretical intrinsic AHC value of 3 S/cm [47]. This difference could be due
to two different approximations to compute AHC. Thus, the theoretically
calculated AHC for the ordered structure of the Mn2CoAl compound is an order
of magnitude smaller than the experimental value of intrinsic AHC $\sim$35
S/cm. Therefore, to understand this difference, we performed the theoretical
calculations by incorporating the amount of anti-site disorder obtained from
the SXRD experiment to compute the AHC. The magnetic moment obtained in this
(disordered Mn2CoAl) case also found to be $\sim$2 $\mu$B/f.u., which
indicates that the anti-site disorder does not affect the magnetization and
therefore, electronic structure is most likely deciding the observed AHC in
disordered Mn2CoAl. In Fig. 6 (top), we have compared the full DFT band
structure with the Wannier interpolated one for disordered structure. The
distribution of calculated Berry curvature along the high symmetry path is
shown in Fig. 6 (bottom). Interestingly, the calculated intrinsic AHC in the
disordered Mn2CoAl increased to the 26.30 S/cm, which is close to our
experimentally obtained intrinsic AHC. In the present case, the disordered
sample still possesses certain space group symmetry, which suggests that the
disorder in the atomic positions will also be repeated in the periodic units
and the theoretical calculation also considers the same. It is true that the
value in the real samples could be slightly different from the theoretical
results. The possible source of significantly less discrepancy in the value of
experimental and theoretical intrinsic AHC for the disordered Mn2CoAl could be
the involved approximations in all steps of the theoretical calculations
(starting from obtaining ground state density to electronic and magnetic
properties to Berry curvature computation to AHC). This suggests that the
anti-site disorder modifies the Berry curvature and hence, enhances the
intrinsic AHC in the Mn2CoAl Heusler compound.

Recently, a comparative study of AHC between SGS Mn2CoGa and half-metallic
Co2VGa compounds has been performed, which reports that the Mn2CoGa compound
exhibits almost zero AHC, whereas half-metallic Co2VGa (which has the similar
magnetization and number of valence electrons as Mn2CoGa) shows a larger AHC
about 140 S/cm [35]. The origin of this behavior was explained based on the
electronic structure of these two compounds. SGS compounds are characterized
by a finite gap in the minority spin channel and a closed gap in the majority
spin channel, which results in a few number of majority spin states around
Fermi level. In the SGS compounds like Mn2CoGa, the Berry curvature of
majority spin states is reported to cancel out by the minority spin states
resulting into a small/almost zero intrinsic AHC [35]. As compared to the SGS
compounds, the half-metallic systems have a finite gap in minority spin
channel and a large number of states in majority spin channel around the Fermi
level [104]. In the half-metallic Co2VGa compound, the large intrinsic AHC has
been observed as compared to the Mn2CoGa compound due to the uncompensated
Berry curvature associated with a large number of majority spin states [35].
Our study follow this behavior as the ordered Mn2CoAl exhibits the small
intrinsic AHC due to its SGS nature (depicted in Fig, 7(a)), while with
introducing the 25% disorder the system seems to be half-metallic due to the
formation of new majority spin states around the Fermi level, provides larger
Berry curvature and leads to enhanced intrinsic AHC (depicted in Fig, 7(b)),
which is in good agreement with the experiment. The inset of the Fig. 7(a) and
(b) show an enlarged view of the intrinsic AHC value around the Fermi level
for ordered and disordered Mn2CoAl compounds, respectively. Thus, our combined
experimental theoretical results show that atomic disorder enhances the
intrinsic AHC due to the modification in the Berry curvature linked with the
electronic structure in the Mn2CoAl compound.

##  IV Conclusion

To conclude, we have presented here the evidence of anti-site disorder
enhanced intrinsic AHC in the Mn2CoAl Heusler compound by comprehensive
analysis of the crystal structure and anomalous Hall effect using experimental
and theoretical tools. The high resolution SXRD data reveals 25% anti-site
disorder between Mn4c and Al atoms within inverse Heusler structure. The
temperature-dependent resistivity shows semiconducting behavior and follows
Mooij’s criteria for disordered metal. Scaling behavior suggests that the
intrinsic mechanism dominates over the extrinsic mechanism in the AHE. The
experimental intrinsic AHC is found to be larger than the theoretically
reported value for the ordered Mn2CoAl. The first-principle calculations
conclude that the anti-site disorder enhances the Berry curvature induced
intrinsic AHC, which is in well agreement with the experimentally found
intrinsic AHC.

## Acknowledgment

We thankfully acknowledge UGC-DAE CSR, Indore for experimental support. SS
thanks Science and Engineering Research Board of India for financial support
through the award of Core Research Grant (grant no: CRG/2021/003256). NS
thanks to DST INSPIRE scheme for financial support. Parts of this research
were carried out at PETRA III of DESY, a member of the Helmholtz Association.
Financial support from the Department of Science and Technology, Government of
India within the framework of the India@DESY is thankfully acknowledged. We
would like to show our gratitude towards the beamline scientist Dr. Martin
Etter for his help in setting up the experiments.

## References

  * [1] C. M. Hurd, The Hall Effect in Metals and Alloys (Plenum, New York, 1972).
  * [2] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).
  * [3] L. Šmejkal, A. H. MacDonald, J. Sinova, S. Nakatsuji, and T. Jungwirth, Nat. Rev. Mater. (2022).
  * [4] S. Friedemann, M. Brando, W. J. Duncan, A. Neubauer, C. Pfleiderer, and F M. Grosche, Phys. Rev. B 87, 024410 (2013).
  * [5] T. Golod, A. Rydh, P. Svedlindh, and V. M. Krasnov, Phys. Rev. B 87, 104407 (2013).
  * [6] Y. Tian, L. Ye, and X. Jin, Phys. Rev. Lett. 103, 087206 (2009).
  * [7] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Zhong Fang, Science 329, 61-64 (2010).
  * [8] A. K. Nayak, J. E. Fischer, Y. Sun, B. Yan, J. Karel, A. C. Komarek, C. Shekhar, N. Kumar, W. Schnelle, J. Kübler, C. Felser, and S. S. P. Parkin, Sci. Adv. 2, e1501870 (2016).
  * [9] C.-Z. Chang, J. Zhang, X. Feng, Jie Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, Science 340, 167-170 (2013).
  * [10] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. G.-Hernandez, L. $\check{S}$mejkal, Y. Chen, S. S. P. Parkin, S. D. Wilson, E. S. Toberer, T. McQueen, and M. N. Ali, Sci. Adv. 6, eabb6003 (2020).
  * [11] T. Taniguchi, J. Grollier, and M. D. Stiles, Phys. Rev. Appl. 3, 044001 (2015).
  * [12] G. E. Bauer, E. Saitoh , and B. J. v. Wees, Nat. Mater. 11(5), 391-9 (2012).
  * [13] S. K. Kim, G. S. D. Beach, K. J. Lee, T. Ono, T. Rasing, and H. Yang, Nat. Mater. 21, 24–34 (2022).
  * [14] T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich Nat. Nanotechnol. 11, 231–241 (2016).
  * [15] L. Šmejkal, Y. Mokrousov, B. Yan, and A. H. MacDonald Nat. Phys. 14, 242–251 (2018).
  * [16] S.D. Bader, S.S.P. Parkin, Annu. Rev. Condens. Matter. Phys. 1, 71–88 (2010).
  * [17] R. J. Gambino and T. Mcguire, IBM Tech. Disclosure Bull. 18, 4214 (1976).
  * [18] A. Gerber, J. Magn. Magn. Mater. 310, 2749-2751 (2007).
  * [19] J. Moritz, B. Rodmacq, S. Auffret, and B. Dieny, J. Phys. D: Appl. Phys. 41, 135001 (2008).
  * [20] R. Karplus and J. M. Luttinger, Phys. Rev. 95, 1145 (1954).
  * [21] S. Onoda, N. Sugimoto, and N. Nagaosa, Phys. Rev. Lett. 97, 126602 (2006).
  * [22] J. Smit and J. Volger, Physica (Amsterdam) 24, 39 (1958).
  * [23] L. Berger, Phys. Rev. B 2, 4559-4566 (1970).
  * [24] M. V. Berry, Proc. R. Soc. London 392, 45-57 (1984).
  * [25] X. Wang, D. Vanderbilt, J. R. Yates, and I. Souza, , Phys. Rev. B 76, 195109 (2007).
  * [26] J. M. Luttinger, Phys. Rev. 112, 739-751 (1958).
  * [27] Y. Taguchi, Y. Oohara, H. Yoshizawa, N. Nagaosa, and Y. Tokura, Science 291, 2573-2576 (2001).
  * [28] T. Jungwirth, Q. Niu, and A. H. MacDonald, Phys. rev. lett., 88, 207208 (2002).
  * [29] S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212–215 (2015).
  * [30] Q. Wang, S. Sun, X. Zhang, F. Pang, and H. Lei, Phys. Rev. B 94, 075135 (2016).
  * [31] Q. Wang, Y. Xu, R. Lou, Nat. Commun. 9, 3681 (2018).
  * [32] N.J. Ghimire, A. S. Botana, J. S. Jiang, J. Zhang, Y. -S. Chen, and J. F. Mitchell, Nat. Commun. 9, 3280 (2018).
  * [33] M. W. Yoo, J. Tornos, A. Sander, L. F. Lin, N. Mohanta, A. Peralta, D. Sanchez-Manzano, F. Gallego, D. Haskel, J. W. Freel, D. J. Keavney, Y. Choi, J. Strempfer, X. Wang, M. Cabero, H. B. Vasili, M. Valvidares, G. Sanchez-Santolino, J. M. Gonzalez-Calbet, A. Rivera, C. Leon, S. Rosenkranz, M. Bibes, A. Barthelemy, A. Anane, E. Dagotto, S. Okamoto, S. G. E. Te Velthuis, J. Santamaria, J. E. Villegas, Nat. Commun. 12, 3283 (2021).
  * [34] B. Ernst, R. Sahoo, Y. Sun, J. Nayak, L. Müchler, A. K. Nayak, N. Kumar, J. Gayles, A. Markou, G. H. Fecher, and C. Felser, Phys. Rev. B 100, 054445 (2019).
  * [35] K. Manna, L. Muechler, T. H. Kao, R. Stinshoff, Y. Zhang, J. Gooth, N. Kumar, G. Kreiner, K. Koepernik, and R. Car, J. Kübler, G. H. Fecher, C. Shekhar, Y. Sun, and C. Felser, Phys. Rev. X 8, 041045 (2018).
  * [36] P. Li, J. Koo, W. Ning, J. Li, L. Miao, L. Min, Y. Zhu, Y. Wang, N. Alem, C. X. Liu, Z. Mao, and B. Yan, Nat. Commun. 11, 3476 (2020).
  * [37] G. K. Shukla, J. Sau, N. Shahi , A. K. Singh, M. Kumar, and S. Singh, Phys. Rev. B 104, 195108 (2021).
  * [38] J. Noky, Y. Zhang, J. Gooth, C. Felser, and Y. Sun, Npj Comput. Mater. 6, 77 (2020).
  * [39] J. Shen, Q. Yao, Q. Zeng, H. Sun, X. Xi, G. Wu, W. Wang, B. Shen, Q. Liu, and E. Liu, Phys. Rev. Lett. 125, 086602 (2020).
  * [40] E. V. Vidal, H. Schneider, and G. Jakob, Phy. Rev. B 83, 174410 (2011).
  * [41] B. K. Hazra, M. M. Raja, R. Rawat, A. Lakhani, S. Kaul, and S. Srinath, J. Magn. Magn. Mater. 448, 371 (2018).
  * [42] Z. Wang, Q. Liu, J. W. Luo, and A. Zunger, Mater. Horiz. 6, 2124 (2019).
  * [43] S. Kastbjerg, N. Bindzus, M. Søndergaard, S. Johnsen, N. Lock, M. Christensen, M. Takata, M. A. Spackman, and B. B. Iversen, Adv. Funct. Mater. 23, 5477–5483 (2013).
  * [44] F. Mende, J. Noky, S. N. Guin, G. H. Fecher, K. Manna,P. Adler, W. Schnelle, Y. Sun, C. Fu, and C. Felser, Adv. Sci. 8, 2100782 (2021).
  * [45] G. K. Shukla, A. K. Jena, N. Shahi, K. K. Dubey, I. Rajput, S. Baral, K. Yadav, K. Mukherjee, A. Lakhani, K. Carva, S.-C. Lee, S. Bhattacharjee, S. Singh, Phys. Rev. B 105, 035124 (2022).
  * [46] S. Roy , R. Singha, A. Ghosh, A. Pariari, and P. Mandal, Phys. Rev. B. 102, 085147 (2020).
  * [47] S. Ouardi, G. Fecher, C. Felser, and J. Kübler, Phys. Rev. Lett. 110, 100401 (2013).
  * [48] X. L. Wang, Phys. Rev. Lett., 100, 156404 (2008).
  * [49] K. Arima, F. Kuroda, S. Yamada, T. Fukushima, T. Oguchi, and K. Hamaya, Phys. Rev. B 97, 054427 (2018).
  * [50] R. G. Buckley, T. Butler, C. Pot, N. M. Strickland and S. Granville, Mater. Res. Express 6, 106113 (2019).
  * [51] X. D. Xu, Z. X. Chen, Y. Sakuraba, A. Perumal, K. Ma-suda, L. S. R. Kumara, H. Tajiri, T. Nakatani, J. Wang, W. Zhou, Y. Miura, T. Ohkubo, and K. Hono, Acta Mater. 176, 33-42 (2019).
  * [52] P. Chen, C. Gao, G. Chen, K. Mi, M. Liu, P. Zhang, and D. Xue, Appl. Phys. Lett. 113, 122402 (2018).
  * [53] Y. Feng, C.-l. Tian, H.-k. Yuan, A.-l. Kuang and H. Chen, J. Phys. D: Appl. Phys. 48, 445003 (2015).
  * [54] N. Y. Sun, Y. Q. Zhang, H. R. Fu, W. R. Che, C. Y. You, and R. Shan, AIP Adv. 6, 015006 (2016).
  * [55] Y. Xin, H. Hao, Y.Ma, H. Luo, F. Meng, H. Liu, E. Liu, and G. Wu, Intermetallics 80, 10-15 (2017).
  * [56] Y. J. Zhang, G. J. Li, E. K. Liu, J. L. Chen, W. H. Wang, and G. H. Wu, J. Appl. Phys. 113, 123901 (2013).
  * [57] G. Z. Xu, Y. Du, X. M. Zhang, H. G. Zhang, E. K. Liu, W. H. Wang, and G. H. Wu, Appl. Phys. Lett. 104, 242408 (2014).
  * [58] M. E. Jamer, B. A. Assaf, T. Devakul, and D. Heiman, App. Phys. Lett. 103, 142403 (2013).
  * [59] K. Ueda, S. Hirose, and H. Asano, App. Phys. Lett. 110, 202405 (2017).
  * [60] J. Kudrnovský, V. Drchal, and I. Turek, Phys. Rev. B 88, 014422 (2013).
  * [61] B. Holt, J. Diaz, J. Huber, and C. A. Luengo, Rev. Bras. Ensino Fis. 8, 155-163 (1978).
  * [62] S. Singh, B. Dutta, S. W. D’Souza, M. G. Zavareh, P. Devi, A. S. Gibbs, T. Hickel, S. Chadov, C. Felser, and D. Pandey, Nat Commun 8, 1006 (2017).
  * [63] S. Singh, R. Rawat, S. E. Muthu, S. W. D’Souza, E. Suard, A. Senyshyn, S. Banik, P. Rajput, S. Bhardwaj, A. M. Awasthi, R. Ranjan, S. Arumugam, D. L. Schlagel, T. A. Lograsso, A. Chakrabarti, and S. R. Barman, Phys. Rev. Lett. 109, 246601 (2012).
  * [64] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Cond. Matt. 21, 395502 (2009).
  * [65] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * [66] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * [67] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847(1997).
  * [68] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
  * [69] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. I.-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, J. Phys. Cond. Matt. 32, 165902 (2020).
  * [70] T. Graf, C. Felser, and S. S. P. Parkin, Prog. Solid. State Ch. 39, 1-50 (2011).
  * [71] J. Rodríguez-Carvajal, FULLPROF, a Rietveld and pattern matching and analysis programs version 2016, Laboratoire Leon Brillouin, CEA-CNRS, France, http://www.ill.eu/sites/ fullprof/
  * [72] A. Nayak, M. Nicklas, S. Chadov, P. Khuntia, C. Shekhar, A. Kalache, M. Baenitz, Y. Skourski, V. K. Guduru, A. Puri, U. Zeitler, J. M. D. Coey, and C. Felser, Nat. Mat. 14, 679–684 (2015).
  * [73] A. K. Nayak, M. Nicklas, S. Chadov, C. Shekhar, Y. Skourski, J. Winterlik, and C. Felser, Phys. Rev. Lett. 110, 127204 (2013).
  * [74] G. D. Liu, X. F. Dai, H. Y. Liu, J. L. Chen, Y. X. Li, G. Xiao, and G. H. Wu, Phys. Rev. B 77, 014424 (2008).
  * [75] K. Momma and F. Izumi, “VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data”, J. Appl. Crystallogr. 44, 1272-1276 (2011).
  * [76] L. N. Bezmaternykh, E. M. Kolesnikova, E. V. Eremin, S. N. Sofronova, N.V. Volkov, and M.S. Molokeev, J. Magn. Magn 364, 55-59 (2014).
  * [77] I. Galanakis, K. Özdoan, E. Sa sioglu, and S. Blügel, J. Appl. Phys. 115, 093908 (2014).
  * [78] S. Dugdale, The Electrical Properties of Disordered Metals (Cambridge University Press, Cambridge) (1995).
  * [79] S. Hikami, A. I. Larkin, and Y. Nagaoko, Prog. Theor. Phys. 63, 707 (1980).
  * [80] G. Xu, W. Wang, X. Zhang, Y. Du, E. Liu, S. Wang, G. Wu, Z. Liu, and X. X. Zhang, Sci. Rep. 4, 5709 (2014).
  * [81] M. Singh, L. Ghosh, V. K. Gangwar, Y. Kumar, D. Pal, P. Shahi, S. Kumar, S. Mukherjee, K. Shimada, and S. Chatterjee, Appl. Phys. Lett. 121, 032403 (2022).
  * [82] J. R. Chamorro, A. Topp, Y. Fang, M. J. Winiarski, C. R. Ast, M. Krivenkov, A. Varykhalov, B. J. Ramshaw, L. M. Schoop, and T. M. McQueen, APL Mater. 7, 121108 (2019).
  * [83] S. Sasmal, R. Mondal, R. Kulkarni, A. Thamizhavel, and B. Singh, J. Phys. Condens. Matter 32, 335701 (2020).
  * [84] A. Laha, S. Malick, R. Singha, P. Mandal, P. Rambabu, V. Kanchana, and Z. Hossain, Phys. Rev. B 99, 241102 (2019).
  * [85] J. Chen, H. Li, B. Ding, Z. Hou, E. Liu, X. Xi, G. Wu, and W. Wang, Appl. Phys. Lett. 116, 101902 (2020).
  * [86] R. P. Dulal, B. R.Dahal, A. Forbes, N. Bhattarai, I. L. Pegg, and J. Phili, Sci Rep 9, 3342 (2019).
  * [87] A. Kawabata, J. Phys. Soc. Jpn. 49, 628-637 (1980)
  * [88] D. V. Baxter, R. Richter, M. L. Trudeau, R. W. Cochrane, and J. O. Strom-Olsen, J. Phys. France 50, 1673-1688 (1989).
  * [89] , W. Chen, H.-Z. Lu, and O. Zilberberg, Phys. Rev. Lett. 122, 196603 (2019).
  * [90] S. Malick, A. Ghosh, C. K. Barman, A. Alam, Z. Hossain, P. Mandal, and J. Nayak, Phys. Rev. B 105, 165105 (2022).
  * [91] H. Afzal, S. Bera, A.K. Mishra, M. Krishnan, M. M. Patidar, R. Venkatesh, and V. Ganesan, J. Supercond. Nov. Magn., 1-8 (2020).
  * [92] J. H. Mooij, Phys. Stat. Sol. (a) 17, 521 (1973).
  * [93] M. Kaveh, and N. F. Mott, J. Phys. C: Solid State Phys. 15, L707 (1982).
  * [94] K. Gofryk, D. Kaczorowski, and T. Plackowski, Phys. Rev. B 72, 094409 (2005).
  * [95] T. Kida, L. A. Fenner, A. A. Dee, I. Terasaki, M. Hagi-wara1, and A. S. Wills, J. Phys.: Condens. Matter 23, 112205 (2011).
  * [96] P. He, L. Ma, Z. Shi, G. Y. Guo, and J.-G. Zheng, and Y. Xin, and S. M. Zhou, Phys. Rev. Lett. 109, 066402 (2012).
  * [97] J. Kötzler, and W. Gil, Phys. Rev. B, 72,060412 (2005).
  * [98] P. Nozieres and C. Lewiner, Journal de Physique 34, 901 (1973).
  * [99] S. Onoda, N. Sugimoto, and N. Nagaosa, Phys. Rev. Lett. 97, 126602 (2006).
  * [100] K. Kim, J. Seo, E. Lee, K.-T. Ko, B. S. Kim, B. G. Jang, J. M. Ok, J. Lee, Y. J. Jo, W. Kang, J. H. Shim, C. Kim, H. W. Yeom, B. I. Min, B.-J. Yang, and J. S. Kim, Nat. Mater. 17, 794–799 (2018).
  * [101] C. Zeng, Y. Yao, Q. Niu, and H. H. Weitering, Phys. Rev. Lett. 96, 037204 (2006).
  * [102] E. Liu, Y. Sun, N. Kumar, L. Muechler, A. Sun, L. Jiao,S. Y. Yang, D. Liu, A. Liang, Q. Xu, J. Kroder, V. Süß, H. Borrmann, C. Shekhar, Z. Wang, C. Xi, W. Wang, W.Schnelle, S. Wirth, Y. Chen, S. T. B. Goennenwein, and C. Felser, Nat. Phys. 14, 1125–1131 (2018).
  * [103] J. Kübler and C. Felser, Phys. Rev. B 85, 012405 (2012).
  * [104] D. Rani, L. Bainsla, A. Alam, K. G. Suresh, Journal of Applied Physics 128, 220902 (2020).


# Supplementary Information – On the origin of viscosity saturation at high
densities during Zebrafish morphogenesis

Rajsekhar Das1 Sumit Sinha2 Xin Li1 T. R. Kirkpatrick3 D. Thirumalai1,2 1
Department of Chemistry, University of Texas at Austin, Austin, Texas 78712,
USA  2 Department of Physics, University of Texas at Austin, Austin, Texas
78712, USA  3 Institute for Physical Science and Technology, University of
Maryland, College Park, Maryland 20742, USA

##  I Cell polydispersity is needed to account for viscosity saturation

To explain the observed saturation in viscosity (Fig. 1C in the main text)
when cell area fraction, $\phi$, exceeds $\phi_{S}$ ($\approx 0.90$), we first
simulated a monodisperse cell system using the model described in the Method
section. The monodisperse cell system ($R=8.5\mu m$) crystallizes (see Fig. S1
(A) and (B)), which excludes it from being a viable model for explaining the
experimental findings. We also find that a system consisting of 50:50 binary
mixture of cells can not account for the experimental data even though
crystallization is avoided (see the next section). These findings forced us to
take polydispersity into account.

In order to develop an agent based model that accounts for polydispersity
effects, we first extracted the distribution, $P(A_{i})$, of the area
($A_{i}$) in the Zebrafish cells from the experimental data PETRIDOU20211914 .
Figure S1 (D), shows that $P(A_{i})$ is broad, implying that cell sizes are
highly heterogeneous. Based on this finding, we sought a model of the non-
confluent tissue that approximately mimics polydispersity found in
experiments. In other words, $\tfrac{\Delta A}{\left\langle A\right\rangle}$
($\Delta A$ is the dispersion in $P(A_{i})$ and $\left\langle A\right\rangle$
is the mean value) should be similar to the data in Fig. S1 (D). The radii
($R_{i}$) of the cells in the simulations are sampled from a Gaussian
distribution $\sim\exp\left(-(R_{i}-\left\langle R\right\rangle)^{2}/2\Delta
R^{2}\right)$ with $\left\langle R\right\rangle=8.5\mu m$ and $\Delta R=4.5\mu
m$. The resulting $P(A_{i})$ for one of the realizations (see Fig. S1 (C)) is
shown in Fig. S1 (E). The value of $\tfrac{\Delta A}{\left\langle
A\right\rangle}$ is $\sim 0.5$, which compares favourably with the
experimental estimate ($\sim 0.4$). The unusual dependence of the viscosity
($\eta$) as a function of packing fraction ($\phi$) cannot be reproduced in
the absence of polydispersity.

Figure S1: Area distribution of the cells: (A) Simulation snapshot for
monodisperse cell system. The number of cells in the two-dimensional periodic
box is $N=500$. (B) Pair correlation function, $g(r)$, as a function of $r$.
There is clear evidence of order, as reflected in the sharp peaks at regular
intervals, which reflects the packing in (A). (C) A schematic picture of
polydisperse cell system from the simulations. Color bar on the right shows
the scale of radii in $\mu m$. There is no discernible order. (D) Distribution
of cell area extracted from experiment during morphogenesis of Zebrafish
blastoderm (extracted from Fig. S2 (A)) PETRIDOU20211914 . (E) Same as (D)
except, $P(A_{i})$, used in a typical simulation. Cell radii vary from $2\mu
m$ to $15\mu m$.  Figure S2: Structure and relaxation behavior for a binary
mixture of cells: (A) A typical simulation snapshot for binary mixture of
cells at $\phi=0.93$. (B) The corresponding pair correlation function, $g(r)$,
between all the cells. The vertical dashed line is at the first peak position
($r_{max}$). (C) $F_{s}(q,t)$, with $q=\tfrac{2\pi}{r_{max}}$, where $r_{max}$
is the location of the first peak in the $g(r)$, as function of time at
various $\phi$ values. (D) The logarithm of the relaxation time,
$\tau_{\alpha}$, as a function of $\phi$. Over the entire range of $\phi$, the
increase in $\tau_{\alpha}$ is well fit by the Vogel-Fucher-Tamman (VFT)
relation. Most importantly, the relaxation time does not saturate, which means
the evolving tissue cannot be modeled using a 50:50 binary mixture. (E)
Effective shear viscosity $\bar{\eta}$ as a function of $\phi$ reflects the
behavior of $\tau$ as a function of $\phi$ in (D).

##  II Relaxation time in the binary mixture of cells does not saturate at
high $\phi$

In the main text, we showed that the saturation in the relaxation time above a
critical area fraction, $\phi_{S}$, is related to two factors. One is that
there ought to be substantial dispersion in the cell sizes (Fig. S1). The
second criterion is that the cells should be soft, allowing them to overlap at
high area fraction. In other words, the cell diameters are explicitly non-
additive. The Hertz potential captures the squishy nature of the cells.

In order to reveal the importance of polydispersity, we first investigated if
a binary mixture of cells (see Fig. S2 (A)) would reproduce the observed
dependence of $\tau_{\alpha}$ on $\phi$. We created a $50:50$ mixture with a
cell size ratio $\sim 1.4$ ($R_{B}=8.5\mu m$ and $R_{S}=6.1\mu m$). All other
parameters are the same as in the polydisperse system (see Table 1). The
radial distribution function, $g(r)$, calculated by considering both the cell
types, exhibits intermediate but not long range translational order (Fig. S2
(B)). A peak in $g(r)$ appears at $r_{max}=14.0\mu m$.

In order to determine the dependence of the relaxation time, $\tau_{\alpha}$,
on $\phi$, we first calculated $F_{s}(q,t)$ at $q=2\pi/r_{max}$ (Fig. S2 (C)).
The decay of $F_{s}(q,t)$ is similar to what one finds in typical glass
forming systems. As $\phi$ increases, the decay of $F_{s}(q,t)$ slows down
dramatically. When $\phi$ exceeds $0.93$ there is a visible plateau at
intermediate times followed by a slow decay. By fitting the long time decay of
$F_{s}(q,t)$ to $\exp(-(t/\tau_{\alpha})^{\beta})$ ($\beta$ is the stretching
exponent), we find that $\tau_{\alpha}(\phi)$ as a function of $\phi$ is well
characterized by the VFT relation (Fig. S2 (D)). There is no evidence of
saturation in $\tau_{\alpha}(\phi)$ at high $\phi$. The $\phi$ dependence of
shear viscosity $\bar{\eta}$ (Fig. S2 (E)), calculated using equation (6),
also shows no sign of saturation. The VFT behavior, with $\phi_{0}\approx 1$
and $D\approx 1.2$, shows that the two component cell system behaves as a
“fragile” glass ANGELL199113 .

A comment regarding the binary system of cells is in order. The variables that
characterize this system are, $\lambda=R_{B}/R_{S}$ ($R_{B}$ ($R_{S}$) is the
radius of the big (small) cells), $\Phi_{B}=N_{B}/(N_{A}+N_{B})$ ($N_{B}$
($N_{A}$) is the number of big (small) cells) with $\Phi_{A}=1-\Phi_{B}$, and
the packing fraction, $\phi=\pi(N_{B}R_{B}^{2}+N_{A}R_{A}^{2})/A_{S}$ where
$A_{S}$ is the area of the sample. The results in Fig. S2 were obtained using
$\lambda=1.4$, $\Phi_{B}=0.5$. With this choice the value of polydispersity is
$(\lambda-1)/(\lambda+1)\approx 0.17$, which is smaller than the experimental
value. It is possible that by thoroughly exploring the parameter space,
$\lambda$ and $\phi$, one could find regions in which the $\tau_{\alpha}$ in
the two component cell system would saturate beyond $\phi_{S}$. However, in
light of the results in Fig. S1 (D)), we choose to simulate systems that also
have high degree of polydispersity (Fig. S1 (E)).

Figure S3:  Free area decreases monotonically for the binary mixture of cells:
(A) Mean voronoi cell size, $\left\langle A\right\rangle$, as a function of
$\phi$ for the 50:50 binary system. (B) The free area fraction, $\phi_{free}$,
as a function of $\phi$ shows that $\phi_{free}$ decreases monotonically as
$\phi$ increases.

##  III Absence of saturation in the free area in binary mixture of cells

In the main text, we established that the viscosity ($\eta$) and relaxation
time ($\tau_{\alpha}$) in the polydisperse cell system saturates beyond
$\phi_{S}$. The dynamics saturates beyond $\phi_{S}$ because the available
area per cells (quantified as $\phi_{free}$) is roughly a constant in this
region ( see Fig. 5 (E) in main text). In the binary system, however, we do
not find any saturation in the $\tau_{\alpha}$. Therefore, if $\phi_{free}$
controls the dynamics, one would expect that $\phi_{free}$ should decrease
monotonically with $\phi$ for the binary cells system. We calculated the
average Vornoi cell size $\left\langle A\right\rangle$ (Fig. S3 (A)) and
$\phi_{free}$ (Fig. S3 (B)) for the binary system. Both $\left\langle
A\right\rangle$ and $\phi_{free}$ decrease with $\phi$, which is consistent
with the free volume picture proposed in the context of glasses.

Figure S4:  Relaxation in the polydisperse cell system is independent of the
waiting time: (A) $F_{s}(q,t)$ for $\phi=0.92$ at different waiting times
($\tau_{\omega}=10^{6}(s))$. Regardless of the value of $\tau_{\omega}$, all
the $F_{s}(q,t)$ curves collapse onto a master curve. (B) Relaxation time,
$\ln(\tau_{\alpha})$, as a function of $\tau_{\omega}$. Over a $3$ orders of
magnitude change in $t_{\omega}$, the variation in relaxation times is less
than the sample to sample fluctuations, as shown by the error bar.

##  IV Absence of aging

Our main result, which we explain by adopting the free volume theory developed
in the context of glasses freeV1 ; freeV2 ; freeV3 ; freeV4 , is that above a
critical packing fraction $\phi_{S}\sim 0.90$ the viscosity saturates.
Relaxation time, $\tau_{\alpha}$, measured using dynamic light scattering, in
nearly monodisperse microgel poly(N- isopropylacrylamide) (PNiPAM)
LudoSaturation  was found to depend only weakly on the volume fraction (3D),
if $\phi_{V}$, exceeds a critical value. It was suggested that the near
saturation of $\tau_{\alpha}$ at high $\phi_{V}$ is due to aging, which is a
non-equilibrium effect. We wondered if aging could explain the observed
saturation of $\eta$ in the embryonic tissue above $\phi_{S}$. If aging causes
the plateau in the tissue dynamics, then $\eta$ or $\tau_{\alpha}$ should be
an increasing function of the waiting time, $\tau_{\omega}$. To test the
effect of $\tau_{\omega}$ on $\tau_{\alpha}$, we calculated the self-
intermediate scattering function $F_{s}(q,t+\tau_{\omega})$ as a function of
$t$ by varying $\tau_{\omega}$ over three orders of magnitude at $\phi=0.92$
(Fig. S4 (A)). There is literally no change in $F_{s}(q,t+\tau_{\omega})$ over
the entire range of $t_{\omega}$. We conclude that, $\tau_{\alpha}$, extracted
from $F_{s}(q,t+\tau_{\omega})$ is independent of $\tau_{\omega}$. The
variations in $\tau_{\alpha}$ (Fig. S4 (B)), with respect to $\tau_{\omega}$,
is significantly smaller than the errors in the simulation. Thus, the
saturation in $\eta$ or $\tau_{\alpha}$ when $\phi>\phi_{S}$ is not a
consequence of aging.

There are two implications related to the absence of aging in the dynamics of
the non-confluent embryonic tissues. (i) Although active forces drive the
dynamics of the cells, as they presumably do in reality, the cell collectives
can be treated as being near equilibrium, justifying the use of Green-Kubo
relation to calculate $\eta$. (ii) Parenthetically, we note that the absence
of significant non-equilibrium effects, even though Zebrafish is a living
system, further justifies the use of equilibrium rigidity percolation theory
to analyze the experimental data PETRIDOU20211914 .

Figure S5:  Measure of ergodicity: (A) Ergodic measure $\Omega(t)$ scaled by
the value at $t=0$ ($\Omega(0)$) as a function of $t$ for $\phi=0.85$. (B) and
(C) Similar plots for $\phi=0.90$ and $\phi=0.92$ respectively. At long time
$\Omega(t)/\Omega(0)$ reach $0.01$, $0.016$ and $0.026$ for $\phi=0.85$,
$\phi=0.90$ and $\phi=0.92$ respectively. (D) $\Omega(0)/\Omega(t)$ as a
function of $t$ for $\phi=0.90$. The dashed line shows a linear fit.

##  V Absence of broken ergodicity

We have shown in the earlier section that, the saturation of $\tau_{\alpha}$
above $\phi_{S}$ is not a consequence of aging. In the context of glasses and
supercooled liquids ergodicity , it has been shown that if ergodicity is
broken then a variety of observables would depend on initial conditions. In
order to test if ergodicity is broken in our model of non-confluent tissues,
we define a measure $\Omega(t)$ given by,

$\Omega(t)=\left[\omega^{k}(t)-\omega^{l}(t)\right]^{2},$ (1)

where, $k$ and $l$ represent systems with different initial conditions, and
$\omega^{k}(t)=\frac{1}{t}\int_{0}^{t}P_{\mu\nu}^{k}(s)ds$ with
$P_{\mu\nu}(s)$ being the value of stress (see Eq. (2) in the main text) at
time $s$. If the system is ergodic, implying the system has explored the whole
phase space on the simulation time scale, the values of $\omega^{k}(t)$ and
$\omega^{l}(t)$ would be independent of $k$ and $l$. Therefore, $\Omega(t)$
should vanish at very long time for ergodic systems. On the other hand, if
ergodicity is broken then $\Omega(t)$ would be a constant whose value would
depend on $k$ and $l$.

The long time values of $\Omega(t)$, normalized by $\Omega(0)$ (at $t=0$), are
$0.01$, $0.016$ and $0.026$ for $\phi=0.85,0.90$ and $0.92$ respectively (Fig.
S5 (A), (B) & (C)). Because these values are sufficiently small, we surmise
that effectively ergodicity is established. Therefore, our conclusion in the
earlier section that the polydisperse cell system is in near equilibrium is
justified, and also explains the absence of aging. Furthermore, it was also
predicted previously ergodicity  that in long time the ergodic measure
($\Omega(t)$ in our case) should decay as $\approx 1/t$. Fig. S5 (D) shows
that this is indeed the case - at long time $\Omega(t)/\Omega(0)$ decays
approximately as $1/t$.

##  VI Calculation of viscosity

Figure S6: Fit of the stress-stress correlation functions to stretched
exponential functions:(A) The stress-stress correlation function $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ divided by the value at $t=0$
$\left\langle P_{\mu\nu}(0)^{2}\right\rangle$ as a function of $t$ for
$\phi\in(0.75-0.87)$. (B) Similar plot for $\phi\in(0.89-0.93)$. (C) The long
time decay of $\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ is fit to
$C_{s}\exp\left[-(\tfrac{t}{\tau_{\eta}})^{\beta}\right]$, as shown by the
dashed lines. The inset shows the dependence of $\beta$ on $\phi$. (D) The
data that is fit using the stretched exponential function (black dashed line)
is combined with the short time data (blue solid line), which is fit using the
cubic spline function. The resulting fits produces a smooth curve$\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{combined}}$, as shown in the
inset.

We calculated the effective viscosity ($\bar{\eta}$) for various values of
$\phi$ by integrating the off-diagonal part of the stress-stress correlation
function $\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$, using the
Green-Kubo relation Book  (without the pre-factor $\tfrac{A}{k_{B}T}$ ),

$\bar{\eta}=\int_{0}^{\infty}dt\sum_{(\mu\nu)^{\prime}}\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle,$ (2)

where $\mu$ and $\nu$ denote Cartesian components ($x$ and $y$) of the stress
tensor $P_{\mu\nu}(t)$ (see main text for the definition of $P_{\mu\nu}(t)$ ).
The time dependence of $\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$,
normalized by $\left\langle P_{\mu\nu}(0)^{2}\right\rangle$, for different
values of $\phi$ (Fig. S6 (A) & (B)) shows that the stress relaxation is
clearly non-exponential, decaying to zero in two steps. After an initial rapid
decay followed by a plateau at intermediate times (clearly visible for
$\phi\geq 0.91$) the normalized $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ decays to zero as a stretched
exponential. The black dashed lines in Fig. S6 (C)) show that a stretched
exponential function,
$C_{s}\exp\left[-\left(\tfrac{t}{\tau_{\eta}}\right)^{\beta}\right]$, where
$\tau_{\eta}$ is the characteristic time in which stress relax and $\beta$ is
the stretching exponent, provides an excellent fit to the long time decay of
$\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ (from the plateau
region to zero) as a function of $t$. Therefore, we utilized the fit function,
$C_{s}\exp\left[-\left(\tfrac{t}{\tau_{\eta}}\right)^{\beta}\right]$, to
replace the noisy long time part of $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ by a smooth fit data before
evaluating the integral in Eqn. (2). The details of the procedure to compute
$\bar{\eta}$ is described below.

We divided $\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle$ in two
parts.

(i) The short time part ($\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{short}}$) – the smooth initial
rapid decay until the plateau is reached (for example see the blue circles in
Fig. S6(D) for $\phi=0.93$). For the $n$ data points at short times,
$\left(t_{1},\left\langle
P_{\mu\nu}(t_{1})P_{\mu\nu}(0)\right\rangle_{\text{short}}\right)$,…,
$\left(t_{n},\left\langle
P_{\mu\nu}(t_{n})P_{\mu\nu}(0)\right\rangle_{\text{short}}\right)$, we
constructed a spline $S(t)$ using a set of cubic polynomials:

$\displaystyle S_{1}(t)=$ $\displaystyle\ \left\langle
P_{\mu\nu}(t_{1})P_{\mu\nu}(0)\right\rangle_{\text{short}}+b_{1}(t-t_{1})+c_{1}(t-t_{1})^{2}+d_{1}(t-t_{1})^{3}$
(3) $\displaystyle S_{2}(t)=$ $\displaystyle\ \left\langle
P_{\mu\nu}(t_{2})P_{\mu\nu}(0)\right\rangle_{\text{short}}+b_{2}(t-t_{2})+c_{2}(t-t_{2})^{2}+d_{2}(t-t_{2})^{3}$
(4) $\displaystyle S_{n-1}(t)=$ $\displaystyle\ \left\langle
P_{\mu\nu}(t_{n-1})P_{\mu\nu}(0)\right\rangle_{\text{short}}+b_{n-1}(t-t_{n-1})+c_{n-1}(t-t_{n-1})^{2}+d_{n-1}(t-t_{n-1})^{3}.$
(5)

The polynomials satisfy the following properties. (a)
$S_{i}(t_{i})=\left\langle
P_{\mu\nu}(t_{i})P_{\mu\nu}(0)\right\rangle_{\text{short}}$ and
$S_{i}(t_{i+1})=\left\langle
P_{\mu\nu}(t_{i+1})P_{\mu\nu}(0)\right\rangle_{\text{short}}$ for $i=1,..,n-1$
which guarantees that the spline function $S(t)$ interpolates between the data
points. (b) $S_{i-1}^{\prime}(t)=S_{i}^{\prime}(t)$ for $i=2,..,n-1$ so that
$S^{\prime}(t)$ is continuous in the interval $\left[t_{1},t_{n}\right]$. (c)
$S_{i-1}^{\prime\prime}(t)=S_{i}^{\prime\prime}(t)$ for $i=2,..,n-1$ so that
$S^{\prime\prime}(t)$ is continuous in the interval
$\left[t_{1},t_{n}\right]$. By solving for the unknown parameters,
$b_{i},c_{i}$ and $d_{i}$ using the above mentioned properties, we constructed
the function S(t). We used $S(t)$ to fit $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{short}}$ to get an evenly
spaced ($\delta t=10s$) smooth data (solid blue line in Fig. S6 (D)). The
fitting was done using the software “Xmgrace”.

(ii) The long time part ($\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{long}}$) – from the plateau
until it decays to zero- is shown by the red circles in Fig. S6 (D). The long
time part was fit using the analytical function
$C_{s}\exp\left[-\left(\tfrac{t}{\tau_{\eta}}\right)^{\beta}\right]$ (black
dashed line in Fig. S6 (D)). We refer to the fit data ($\delta t=10s$) as
$\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{long}}^{\text{fit}}$. We then
combined $\left\langle P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{short}}$
and $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{long}}^{\text{fitted}}$ to
obtain $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{combined}}$ (see inset of Fig.
S6 (D)). Finally, we calculated $\bar{\eta}$ using the equation,

$\displaystyle\bar{\eta}$ $\displaystyle=\lim_{\delta t\to
0}\sum_{i=0}^{T}\delta t\sum_{(\mu\nu)^{\prime}}\left\langle
P_{\mu\nu}(i\delta t)P_{\mu\nu}(0)\right\rangle_{\text{combined}}$
$\displaystyle=\lim_{\delta t\to 0}\sum_{i=0}^{t_{1}}\delta
t\sum_{(\mu\nu)^{\prime}}\left\langle P_{\mu\nu}(i\delta
t)P_{\mu\nu}(0)\right\rangle_{\text{short}}+\lim_{\delta t\to
0}\sum_{i=t_{1}}^{T}\delta t\sum_{(\mu\nu)^{\prime}}\left\langle
P_{\mu\nu}(i\delta t)P_{\mu\nu}(0)\right\rangle_{\text{long}}^{\text{fit}},$
(6)

where $t_{1}\delta t$ is the end point of $\left\langle
P_{\mu\nu}(t)P_{\mu\nu}(0)\right\rangle_{\text{short}}$ and $T\delta t$ is the
end point of $\left\langle P_{\mu\nu}(i\delta
t)P_{\mu\nu}(0)\right\rangle_{\text{combined}}$.

Figure S7: Cell size dependent structures and dynamics: (A) Radial
distribution function $g_{SS}(r)$ between small sized cells ($R_{S}\leq 4.5\mu
m$) at $\phi=0.905$ (blue) and $0.92$ (red). These values are greater than
$\phi_{s}\approx 0.90$ (B) Same as (A) except the results are for $g_{BB}(r)$
between large cells ($R_{B}\geq 12.0\mu m$). (C) $F_{s}(q,t)$ for cells with
$R_{S}\leq 4.5\mu m$ at $\phi=0.905$ and $\phi=0.92$. Note that even at these
dense packings, the mobility of the smaller sized cells is substantial, which
is reflected in the time dependence of $F_{s}(q,t)$. (D) $F_{s}(q,t)$ for
cells with $R_{B}\geq 12.0\mu m$ at $\phi=0.905$ and $\phi=0.92$. The black
dashed lines are fits to stretched exponential functions,
$F_{s}(q,t)\sim\exp(-(\tfrac{t}{\tau_{\alpha}})^{\beta})$ where
$\tau_{\alpha}$ is the relaxation time and $\beta$ is the stretching exponent.
The dotted lines correspond to the value $F_{s}(q,t)=\tfrac{1}{e}$.

##  VII Dynamics of small and large cells are dramatically different

The structural and the dynamical behavior of the small ($R_{S}\leq 4.5\mu m$)
and large ($R_{B}\geq 12.0\mu m$) cells are dramatically different in the non-
confluent tissue. The pair correlation functions between small cells
($g_{SS}(r)$) and between large cells ($g_{BB}(r)$) (Fig. S7 (A) & (B)) for
$\phi=0.905$ and $\phi=0.92$ show that, both small and large cells exhibit
liquid like disordered structures. However, it is important to note that,
$g_{SS}(r)$ has only one prominent peak and a modest second peak. In contrast,
$g_{BB}(r)$ has three prominent peaks. Thus, the smaller sized cells exhibit
liquid-like behavior whereas the large cells are jammed. This structural
feature is reflected in the decay of $F_{s}(q,t)$ with
$q=\tfrac{2\pi}{r_{max}}$, where $r_{max}$ is the position of the first peak
in the $g(r)$ (Fig. S7 (C) & (D)). There is a clear difference in the decay of
$F_{s}(q,t)$, spanning nearly four orders of magnitude, between small and
large cells (compare Figs. S7 (C) and (D)). At the highly jammed packing
fraction ($\phi=0.92$), small sized cells have the characteristics of fluid-
like behavior.

Figure S8: Simulation snapshot and trajectories for a few smaller and bigger
sized cells: (A) Cells ($\phi=0.91$) are coloured according to their sizes
(gray colors). A few small sized cells are shown in different colors (pink,
blue, orange, purple, cyan, light purple and yellow). (B) The corresponding
trajectories are shown over the entire simulation time. (C) Similar plot as
(A) but for a few bigger sized cells shown in purple, yellow, light green,
red, cyan and green colors. (D) Same as (B) expect the trajectories of the
large sized cells are highlighted. Clearly, the large cells are jammed.

A typical snapshot from one of the simulations ($\phi=0.91$) and the
trajectories for a few small and big sized cells are displayed in Fig. S8 (A),
(B), (C) and (D). The figures reflect the decay in $F_{s}(q,t)$. Not only are
the mobilities heterogeneous it is also clear that the the Mean Square
Displacements (MSDs) of the small cells are greater than the large cells.

Figure S9: Dynamical rearrangement of jammed cells: The changing local
environment of a randomly selected cell(black) over time. (Top Panels: From
left to right $t=9.41\tau_{\alpha},10.01\tau_{\alpha}$ and
$25.39\tau_{\alpha}$). The black coloured cell is completely jammed by other
cells. (Bottom Panels: Form left to right
$t=10.97\tau_{\alpha},25.44\tau_{\alpha}$ and $27.49\tau_{\alpha}$). Dynamical
facilitation, resulting in collective rearrangement of the cells surrounding
the black cell, enables it to move in the dynamically created free volume.

##  VIII Dynamical changes in local packing fraction cause jammed cells to
move

The mobilities of all the cells, even under highly jammed conditions
($\phi\geq\phi_{S}$), is only possible if the local area fraction changes
dynamically. This is a collective effect, which is difficult to quantify
because it would involve multi-cell correlation function. As explained in the
main text, the movement of a jammed cell can only occur if several neighboring
cells move to create space. This picture is not that dissimilar to kinetic
facilitation in glass-forming materials ScienceChandler ; BiroliJCP . However
in glass like systems facilitation is due to thermal excitation but in the
active system, it is self-propulsion that causes the cells to move.

The creation of free space may be visualized by tracking the positions of the
nearest neighbour cells. In Fig. S9, we display the local free area of a black
coloured cell at different times. The top panels show the configurations where
the black cell is completely jammed by other cells. The cell, coloured in
black, can move if the neighboring cells rearrange (caging effect in glass
forming systems) in order to increase the available free space. The bottom
panels show that upon rearrangement of the cells surrounding the black cell
its mobility increases. Such rearrangement occurs continuously, which
qualitatively explains the saturation in viscosity in the multicomponent cell
system.

Figure S10: Finite size effects: $F_{s}(q,t)$ for $N=200$ (A) and $N=750$ (B).
Logarithm of $\tau_{\alpha}$ as a function of $\phi$ for $N=200$ (C) and for
$N=750$ (D). The dashed lines are the VFT fits.

##  IX Finite system size effects

In the main text, we report results for $N=500$. To asses if the unusual
dynamics is not an effect of finite system size, we performed additional
simulations with $N=200$ and $N=750$. As shown in Fig. S10 (A) & (B),
$F_{s}(q,t)$ saturates at $\phi\geq\phi_{S}$, which is reflected in the
logarithm of $\tau_{\alpha}$ as a function of $\phi$ (Fig. S10 (C) & (D)). The
saturation value $\phi_{S}\sim 0.90$ is independent of the system size. The
value of $\phi_{0}$ ($\approx 0.95$) is also nearly independent of system
size. Therefore, the observed dynamics, reflected in the plateau in the
viscosity at high $\phi$, is likely not an effect of finite system size.

Figure S11: Mean coordination number and cell area fraction: (A), (B) and (C)
shows the distribution of coordination number $P(N_{c})$ for $\phi$ =0.85,
$0.90$ and $0.93$ respectively. The orange lines are Gaussian fits to the
histograms. (D) Shows mean $\left\langle N_{c}\right\rangle$ as a function of
$\phi$. The dashed line shows the linear relationship between them.

##  X Dependence of viscosity on average coordination number and average
connectivity

Figure S12: Viscosity and coordination number: (A) Shows $\left\langle
C\right\rangle$ as a function of $\left\langle N_{c}\right\rangle$. Clearly
they are linearly related as shown by the dashed line. Viscosity $\bar{\eta}$
as a function of $\left\langle N_{c}\right\rangle$ (B) and $\left\langle
C\right\rangle$ (C).

In the Zebrafish blastoderm experiment PETRIDOU20211914 , the change in
viscosity ($\eta$) was shown as a function of mean connectivity $\left\langle
C\right\rangle$. To test whether our two dimensional tissue simulations also
exhibit similar behavior, we calculated viscosity as a function of mean
coordination number $\left\langle N_{c}\right\rangle$ (defined below), which
is equivalent to $\left\langle C\right\rangle$. We define the coordination
number, $N_{c}$, as the number of cells that are in contact with a given cell.
Two cells with indices $i$ and $j$ are in contact if
$h_{ij}=R_{i}+R_{j}-r_{ij}>0$. We calculate $N_{c}$ for all the cells for each
$\phi$ and calculate the histogram $P(N_{c})$. The distributions of
$P(N_{c})$s are well fit by a Gaussian distribution function
$A\exp\left[-(\frac{x-\mu}{\sigma})^{2}\right]$ (Figs. S11 (A), (B) & (C)).
The calculated mean, $\left\langle N_{c}\right\rangle$, from the fit is
linearly related to the cell area fraction $\phi$ (Fig. S11 (D)).

We also calculate the average connectivity $\left\langle C\right\rangle$
defined in the following way. Each cell is defined as a node and an edge is
defined as the line connecting two nodes. If a snapshot has $n$ nodes and $m$
edges then the connectivity is defined as $C=\frac{2m}{n}$ PETRIDOU20211914 .
We calculate $C$ for all the snapshots for each $\phi$ and estimated its mean
value $\left\langle C\right\rangle$. We find that $\left\langle
C\right\rangle$ and $\left\langle N_{c}\right\rangle$ are of similar values
(Fig. S12 (A)), and the dependence of viscosity $\bar{\eta}$ on $\left\langle
N_{c}\right\rangle$ and $\left\langle C\right\rangle$ is similar to
experimental results (Fig. S12 (B) & (C)) PETRIDOU20211914 .

##  XI Connectivity Map

Figure S13: Connectivity profile: Connectivity map for
$\phi=0.80,0.85,0.89,0.90,0.92$ and $0.93$ are shown in (A), (B), (C), (D),
(E) and (F) respectively. For $\phi\geq 0.89$ there is a path that connects
the cells in the entire sample. The percolation transition occurs over a very
narrow range of $\phi$ (roughly at $\phi\approx 0.89$ -orange map), which also
coincides with the sharp increase in $\eta$, thus linking equilibrium
transition to geometric connectivity.

To pictorially observe the percolation transition in the simulations, we plot
the connectivity map for various values $\phi$ in Fig. S13. Note that for
smaller $\phi\leq 0.89$ the map shows that cells are loosely connected,
suggesting a fluid-like behavior. For $\phi\geq 0.89$ the connectivity between
the cells spans the entire system, which was noted and analysed elsewhere
thoroughly PETRIDOU20211914 . Cells at one side of simulation box are
connected to cells at the side. The cell connectivity extends throughout the
sample. The transition from a non-percolated state to a percolated state
occurs over a very narrow range of $\phi$, which corresponds to the onset of
rigidity percolation transition occurs PETRIDOU20211914 . It is gratifying
that the simple simulations reproduce the experimental observations.

##  XII Simulation parameters

The parameters for the simulation model described in the main text are given
in Table 1. In the present study, we do not consider the growth and division
of cells. Therefore, we investigate the steady-state dynamics of the tissue.
We have run each of the simulations at least ($5-10$)$\tau_{\alpha}$ before
storing the data.

Table 1: The parameters used in the simulation.

Parameters |  Values |  References
---|---|---
Timestep ($\Delta t$) |  10$\mathrm{s}$ |  This paper
Self-propulsion ($\mu$) |  0.045 $\mathrm{\mu m/\sqrt{s}}$ |  This paper
Friction coefficient ($\gamma_{o}$) |  0.1 $\mathrm{kg/(\mu m~{}s)}$ |  This paper
Mean Cell Elastic Modulus ($E_{i})$ |  $10^{-3}\mathrm{MPa}$ |  galle2005modeling ; malmi2018cell
Mean Cell Poisson Ratio ($\nu_{i}$) |  0.5 |  schaller2005multicellular ; malmi2018cell

## References

  * [1] N. I. Petridou, B. Corominas-Murtra, C.-P. Heisenberg, and E. Hannezo.  Rigidity percolation uncovers a structural basis for embryonic tissue phase transitions.  Cell, 184(7):1914–1928.e19, 2021.
  * [2] C. Angell.  Relaxation in liquids, polymers and plastic crystals — strong/fragile patterns and problems.  Journal of Non-Crystalline Solids, 131-133:13–31, 1991.  Proceedings of the International Discussion Meeting on Relaxations in Complex Systems.
  * [3] M. H. Cohen and D. Turnbull.  Molecular transport in liquids and glasses.  The Journal of Chemical Physics, 31(5):1164–1169, 1959.
  * [4] D. Turnbull and M. H. Cohen.  Free‐volume model of the amorphous phase: Glass transition.  The Journal of Chemical Physics, 34(1):120–125, 1961.
  * [5] D. Turnbull and M. H. Cohen.  On the free‐volume model of the liquid‐glass transition.  The Journal of Chemical Physics, 52(6):3038–3041, 1970.
  * [6] K. Falk, D. Savio, and M. Moseler.  Nonempirical free volume viscosity model for alkane lubricants under severe pressures.  Physical Review Letters, 124:105501, Mar 2020.
  * [7] A.-M. Philippe, D. Truzzolillo, J. Galvan-Myoshi, P. Dieudonné-George, V. Trappe, L. Berthier, and L. Cipelletti.  Glass transition of soft colloids.  Physical Review E, 97:040601, Apr 2018.
  * [8] D. Thirumalai, R. D. Mountain, and T. R. Kirkpatrick.  Ergodic behavior in supercooled liquids and in glasses.  Physical Review A, 39:3563–3574, Apr 1989.
  * [9] J.-P. Hansen and I. R. McDonald.  Theory of Simple Liquids: With Applications to Soft Matter.  Academic press, 2013.
  * [10] L. O. Hedges, R. L. Jack, J. P. Garrahan, and D. Chandler.  Dynamic order-disorder in atomistic models of structural glass formers.  Science, 323(5919):1309–1313, 2009.
  * [11] G. Biroli and J. P. Garrahan.  Perspective: The glass transition.  The Journal of Chemical Physics, 138(12):12A301, 2013.
  * [12] J. Galle, M. Loeffler, and D. Drasdo.  Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in vitro.  Biophysical Journal, 88(1):62–75, 2005.
  * [13] A. N. Malmi-Kakkada, X. Li, H. S. Samanta, S. Sinha, and D. Thirumalai.  Cell growth rate dictates the onset of glass to fluidlike transition and long time superdiffusion in an evolving cell colony.  Physical Review X, 8(2):021025, 2018.
  * [14] G. Schaller and M. Meyer-Hermann.  Multicellular tumor spheroid in an off-lattice voronoi-delaunay cell model.  Physical Review E, 71(5):051910, 2005.