# Yielding behaviour of active particles in bulk and in confinement

Yagyik Goswami  Laboratory of Nanoscale Biology, Paul Scherrer Institut,
Villigen, Aargau, Switzerland  Theoretical Sciences Unit and School of
Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research,
Rachenahalli Lake Road, Bengaluru-560064, India  G. V. Shivashankar
Laboratory of Nanoscale Biology, Paul Scherrer Institut, Villigen, Aargau,
Switzerland  Department of Health Sciences and Technology, ETH Zürich, Zürich,
Switzerland  Srikanth Sastry<EMAIL_ADDRESS>Theoretical Sciences Unit
and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced
Scientific Research, Rachenahalli Lake Road, Bengaluru-560064, India

###### Abstract

The investigation of collective behaviour in dense assemblies of self-
propelled active particles has been motivated by a wide range of biological
phenomena. Of particular interest are dynamical transitions of cellular and
sub-cellular biological assemblies, including the cytoskeleton and the cell
nucleus. Motivated by observations of mechanically induced changes in the
dynamics of such systems, and the apparent role of confinement geometry, we
study the transition between jammed and fluidized states of active particles
assemblies, as a function of the strength and temporal persistence of active
forces, and in different confinement geometries. Our results show that the
fluidization transition broadly resembles yielding in amorphous solids,
consistently with recent suggestions. More specifically, however, we find that
a detailed analogy holds with the yielding transition under cyclic shear
deformation, for finite persistence times. The fluidization transition is
accompanied by driving induced annealing, strong dependence on the initial
state of the system, a divergence of time scales to reach steady states, and a
discontinuous onset of diffusive motion. We also observe a striking dependence
of the transition on persistence times and on the nature of the confinement.
Collectively, our results have implications in epigenetic cell state
transitions induced by alterations in confinement geometry.

###  I Introduction

Active matter, composed of interacting self-propelled particles or entities,
has been investigated intensely in recent years [1, 2], and displays novel
aspects of dynamics and organization not realised in conventional condensed
matter. Although avenues for synthesizing active materials have been well
explored, many obvious examples of active matter arise in biological systems
across length scales, from bird flocks and animal herds[3], bacterial
suspensions[4], tissues [5, 6, 7], to the cytoskeleton[8], which have been
investigated theoretically and computationally [9] to comprehend experimental
findings. In several contexts, the systems of interest may be characterised as
dense, disordered assemblies, displaying features characteristic of glassy
systems or phenomenology associated with packings that undergo jamming.

In particular, novel behaviour arising in models of active glasses [10, 11,
12, 13, 14, 5, 6, 7, 15, 16, 17, 18, 35, 32, 19] have been explored, seeking
similarities with and departures from glass and jamming phenomenology in the
absence of active forces or self-propulsion. In addition to investigations of
the nature of relaxation dynamics, and approach to structurally arrested
states with a decrease in thermal fluctuations and active forces and/or
increase in density [13, 5, 6, 7], some recent works [15, 14, 17, 18, 19] have
addressed the converse problem of loss of rigidity or fluidization, with an
increase in the magnitude of active forces. Such a fluidization transition has
an appealing analogy with the phenomenon of yielding in amorphous solids, upon
an increase of external driving through the application of, e. g., shear
stress or deformation. The active forces on individual particles, in the
simple case where their orientations remain fixed, may be viewed as a
spatially randomised generalisation of shear induced forces or displacements
that arise when a global stress or strain is applied. Such an analogy has
indeed been fruitfully pursued for the case when the orientations of the
active forces is fixed.

However, if one considers such active glass models in specific experimental
contexts, including biological contexts where the active forces are generated
by energy-consumption, one must consider a situation where the orientation of
the active forces is not fixed but may evolve in time, in a manner
characterised by a persistence time. For active glasses with finite
persistence times, we find that the fluidization transition bears striking
similarities to the yielding transition under cyclic shear deformation,
wherein the shear deformation is applied cyclically (typically) around the
zero strain state. The response of amorphous solids to cyclic shear is
characterised [20, 21] by (i) annealing effects at small deformations, (ii)
strong dependence on the initial annealing state of the amorphous solids,
(iii) a discontinuous transition from an absorbing state with zero diffusivity
of particles to a diffusive state with finite diffusivities, and (iv)
diverging times to reach steady states as the yielding transition is
approached from either side. As we show, active glasses with finite (but
large) persistence times show a transition from arrested to fluidized states
display these features remarkably.

In addition to the persistence time, we also study the role of confinement on
the transition to the fluidized state. The presence of a confining boundary
has been found to induce non-trivial behaviour in active systems[22, 23, 24].
Such effects can be expected to play an important role in determining the
properties of many biological assemblies that experience strong confinement
conditions. A particular example of this kind is the organization of chromatin
within the nucleus, where an interplay of active forces and the effect of the
confinement boundary are of potential interest to understand in mechanically
induced cell-state changes[25, 26, 27]. While we do not attempt to model the
specific details of such biological assemblies in this work, the results we
obtain should be of relevance in several such contexts.

Instead, we aim to consider an idealised system, studying which we aim to
comprehend the the role of perturbing factors which may be biologically
relevant. We thus consider the arrested-to-fluidised, or yielding, transition
in a dense assembly of particles at very low temperatures, subjected to active
forces which have a finite persistence time. We investigate the yielding
behaviour of local energy minimum configurations, referred to as glasses or
amorphous solids, which are prepared with different degrees of annealing (see
Materials and Methods and SI Appendix, Fig. S1 for details), followed by local
energy minimization, so that their initial energies vary over a wide range. We
simulate a two dimensional $65:35$ (Kob-Andersen) binary mixture of particles
interacting with the Lennard-Jones potential, whose dynamics is described by
the Langevin equation, with forces from interparticle interactions, and active
forces whose orientation changes diffusively with a specified persistence
time. Further details of the model and simulations methods are given under
Materials and Methods.

We first investigate the yielding behaviour as a function of the strength of
the activity and persistence time, and demonstrate that the yielding behaviour
indeed strikingly resembles the yielding transition under cyclic shear. Since
we also expect the effects of confinement to be significant, we consider
different confinement geometries, where the interacting particles experience
short range attractive interactions and strong repulsion at smaller distances
to the confining boundary.

###  II Results

We first consider a system of $N=1000$ particles in two dimensions with a
number density of $\rho=1.2$, with a small reduced temperature of $T=10^{-3}$,
with periodic boundary conditions. Each of these cases is simulated using
Langevin dynamics with forces arising from interparticle interactions and an
active force, with a time step of $dt=0.01$, over a range of values of the
active force magnitudes $f$, which reorients on a time scale given by the
persistence time $\tau_{p}$ (see Materials and Methods). We thus present our
results for a fixed value of $\tau_{p}$ as a function of $f$, considering
amorphous solids with a range of initial energies as described earlier (see SI
Appendix, Fig. S1). In all cases, the systems evolve to reach steady states
and we report the properties of the system in the steady state (unless
otherwise stated).

####  II.1 Yielding transition

In cyclic shear simulations, it is found that for high initial energies (or,
poorly annealed amorphous solids) the final (non-diffusive or absorbing)
steady states exhibit progressively larger degrees of mechanically induced
annealing (decrease in the internal energy) with increasing amplitude of
strain (the control parameter, like the strength of active force $f$ in the
present case) till the yielding or fluidization transition, beyond which there
is a drop in stress and a discontinuous onset of diffusive behavior. Well
annealed (low initial energy) amorphous solids exhibit little annealing
before, and strongly discontinuous change at, the yield point. We thus
consider the corresponding behaviour of the steady state energies and stresses
(Fig. 1) and diffusion (Fig. 2) for the active glasses we now investigate.
Fig. 1 (A) shows, as a function of the active force strength $f$, the final
(steady state) energies reached, for a number of initial energies, for
$\tau_{p}=2.31\times 10^{2}$. The resulting curves bear striking similarities
with the cyclic shear yielding diagram [21] wherein, as a function of the
strain amplitude (i) poorly annealed (high energy) glasses show progressive
mechanical annealing before yielding, and (ii) well annealed glasses show
strong discontinuous yielding. In Fig. S2 in the SI Appendix, we show
corresponding steady state energies at two other, higher, values of
$\tau_{p}$.

The stress can be computed from the first derivative of the energy with
respect to strain, and corresponding procedures for computing the active
stress (discussed, e.g in [16, 18]) requires identifying a strain step
$\Delta\gamma$ associated with displacements in the presence of active forces.
Following [18], we first consider that the alignment of particle velocities
with the direction of active force produces a net strain rate (see Eq. 3 in
Materials and Methods and SI Appendix), which can be used to obtain the strain
step $\Delta\gamma$ over a time interval (see Materials and Methods). We
analyse the parametric dependence of $\Delta E$ on $\Delta\gamma$ to compute
the active stress. The resultant average stress, $\langle\sigma_{act}\rangle$,
defined in Eq. 4 is found to increase linearly with $f$ until the yield point
where it attains a maximal value, as shown in Fig. 1 (B). The net velocity
alignment, from which we determine an effective strain rate, remains
negligible in magnitude below the yield point, with a subsequent increase for
active forces larger than the yield value (see SI Appendix, Fig. S3) as also
noted in similar contexts [28, 29].

For comparison, we show the corresponding data for cyclic shear, in Fig. 1 (C)
and (D), which indeed confirms the similarity of the results in these two
different driving protocols. Similar results have also been reported recently
for active ‘run and tumble’ driving [30] at considerably smaller persistence
times, in the thermal limit, as discussed, e. g., in [13] and below. In Fig. 1
(D), we compare the shear stress at maximal strain, $\gamma_{max}$, in the
cyclically sheared system, obtained from the virial expression with that
obtained from the procedure outlined above (see also SI Appendix) and employed
for the active case (Fig. 1 (B)). The close agreement validates the latter
procedure, although improvements to the procedure merit further investigation.
Comparing results in Fig. 1 (B) and Fig. 1 (D) we note that while the computed
stresses clearly mark the fluidization transition for the active case, the
stress overshoot, the presence of a constant flow stress, as well as the
difference between the peak stresses of the well-annealed samples and the
poorly annealed samples, is less pronounced.

(a)  (b)  (c)  (d)

Figure 1: Yielding diagram under active driving. (A) Amorphous solids subject
to active forces display annealing to lower energies for high energy initial
configurations for small active forces, with the degree of annealing
diminished or absent for lower energy initial states. The steady state
energies are shown. A common, ergodic, fluidized state is observed above a
critical value of active forces. These observations are analogous to those
observed for yielding under cyclic shear deformation. (B) The yielding
transition is accompanied by a saturation of the active stress, computed as
explained in Materials and Methods. Corresponding results for (C) energy at
zero strain, and (D) stress at strain amplitude $\gamma_{max}$, for cyclic
shear are shown for comparison. Also shown in (D) is the comparison of
stresses obtained from the procedure outlined in Materials and Methods (open
symbols) and the conventional virial stress (filled symbols).

####  II.2 Fluidisation transition under active forcing and relaxation to the
steady state

We next consider the mean squared displacements (MSD) of particles as a
function of elapsed time, in the steady state, for different values of $f$, as
shown in Fig. 2, for two different initial conditions (Fig. 2(A) and 2(B)). It
can be seen that the MSD curves exhibit a saturation below a critical value of
$f$ (which can be deduced from data in Fig. 1) and a linear dependence on time
for higher $f$ values. The corresponding diffusion coefficients jump
discontinuously at the yield point, analogously with cyclic shear yielding
(see SI Appendix, Fig. S4 for supporting data).

(a)  (b)  (c)  (d)

Figure 2: Diffusivity and time to reach steady state. (A) Mean squared
displacement for different magnitudes of the active force showing a transition
from an absorbing (non-diffusive) to a diffusive state for the poorly annealed
case, and (B) the well-annealed case. Lines are shown with linear (navy blue)
and no time dependence (dashed maroon) for reference. (C) Relaxation curves
for energy vs. active force, obtained by averaging over $8$ independent
trajectories in logarithmically spaced time intervals (D) Relaxation times
exhibiting divergence at yielding transition. Dashed orange and turquoise
lines are guides to the eye, indicating the diverging time taken to reach the
steady state.

Another characteristic of the cyclic shear yielding transition is the time
taken to reach the steady state for different strain amplitudes. To obtain the
corresponding data, we consider the relaxation of the energy (shown in Fig. 2
(C) for $T_{p}=0.35$) for different values of $f$. A stretched exponential fit
to the average energy over time yields a relaxation timescale, $\langle
t_{ss}\rangle$ (see SI Appendix, Fig. S5, S6 for the fits) which is shown in
Fig. 2 (D), indicating that the time to approach the steady state becomes
progressively longer as the yielding transition is approached from either
side, with an apparent divergence at the transition. Data for force magnitudes
smaller than the critical force are shown for a poorly annealed sample,
prepared at $T_{p}=1.00$, while those for larger forces are shown for
$T_{p}=0.35$. The estimated times to steady state obtained from the procedure
described above compare well with corresponding estimates from the average
first passage time to reach the steady state energy, as shown in the SI
Appendix in Fig. S7.

The above results clearly demonstrate that yielding in a system with active
forces with a finite persistence time has the characteristics of yielding
under cyclic shear. Intuitively, a way to understand this is to note that in
addition to having no macroscopically well-defined direction of driving as in,
say, simple shear deformation, with finite persistence time, the direction of
driving also changes over time, thus bearing a broad similarity to the ‘back
and forth’ driving under cyclic shear.

In the thermal limit of low persistence times $\tau_{p}$, the active driving
we consider is equivalent to a higher effective temperature, with the
contribution from active driving increasing quadratically with $f$ (see SI
Appendix, Fig. S8-S9 for comparison with kinetic energies in the low
persistence limit and passive dynamics). Extensions of such a relationship to
finite and large persistence times have also been discussed [31, 13, 32]. The
persistence times we consider are larger than those in the thermal regime
[13]. Indeed, they span a range that includes intermediate (up to $\sim
10^{3}$ and large persistence times ($10^{4}-10^{7}$) for which one observes
clear deviations from the expectations in [31] and the yielding active force
value (see below) increases with active force $f$, rather than saturate to a
constant value.

We next focus on the dependence of this yielding transition on parameters that
may be of relevance not only in the biological context of our interest, but
more generically. We thus consider the effect of the change in persistence
times, and the nature of the confinement, on the yielding transition.

####  II.3 Role of persistence time in the yield point

A key control factor of relevance in various biological contexts is the
persistence time, which determines the time scale over which the active force
reorients. Changes to the persistence time have strong implications for the
efficacy of exploration of phase space, producing different dynamical
regimes[29, 13, 33, 34, 35, 32], motivating us to study changes in yielding
behaviour arising from changes in the persistence time. We study the effect of
changing the persistence time on mechanical annealing below, and on the
location of, the yield point, for both well annealed (Fig. 3 (A)) and poorly
annealed samples (Fig. 3 (B)). In both cases, the yield point shifts to larger
values of active force $f$ as the persistence time increases. For the poorly
annealed case, the mechanical annealing effects diminish as the persistence
time increases. In both cases, the apparently continuous increase of energies
in the yielded branch have to do with the difficulty in reaching the steady-
state close to the transition, as we explain further in the SI Appendix. Thus,
with an increase in the persistence time, the particles are more and more
readily ‘jammed’ in kinetically arrested higher energy states, reminiscent of
phenomenology explored in [13, 35, 32]. Unlike the small persistence time
limit, yielding occurs at higher values of the active force with increasing
persistence time, consistently with data and discussion in [13, 32]. We also
show the results for the limit of zero persistence time ($\tau_{p}=dt$,
mentioned above), which reveal negligible change in the range of $f$ values
shown (which correspond to very small increments in the bath temperature, as
elaborated in the SI Appendix, Fig. S8).

In Fig. 3 (C), we show the change in $t_{ss}$ as the persistence time is
increased, where $t_{ss}$ below the yield point are obtained for poorly
annealed initial samples, and for well annealed samples above the yield point,
as before. The point of divergence shifts to larger values of $f$ with larger
$\tau_{p}$, consistently with results in Fig. 3 (A) and (B). Fig. 3 (D) shows
the active stress, $\sigma_{act}$, which increases up to the yield point (as
seen earlier (Fig. 1 (B)), which are slightly different for the two cases
shown. The schematic in Fig. 3 (E) illustrates how the efficiency of
exploration of phase space, and in turn the activity induced annealing and the
yield point, may be altered by the persistence time. Such a picture, while
intuitive for the range of $\tau_{p}$ that we have explored, is not expected
to hold for small $\tau_{p}$ values, where the system begins to resemble a
thermal system as the $\tau_{p}\rightarrow 0$ limit is approached. Preliminary
data also indicate that the yield values increase non-monotically upon
lowering $\tau_{p}$, consistently with results in [13].

(a)  (b)  (c)  (d)  (e)

Figure 3: Dependence of yielding behaviour on persistence time. The yielding
transition shifts to larger active force magnitudes $f$ as the persistence
time increases, for (A) a well annealed and (B) poorly annealed initial
configurations. The mechanical annealing in the latter case diminishes with an
increase in the persistence time. (C) The divergence of the time to reach the
steady state at different persistence times. The dashed lines are guides to
the eye. (D) Stress as a function of active force, with a shift in the value
of force ($f$) where deviation from linearity first occurs. The value of $f$
corresponding to the peak stress is indicated by a vertical magenta line for
$\tau_{p}=2.31\times 10^{2}~{}\tau$ and a vertical orange for
$\tau_{p}=1.01\times 10^{4}\tau$. (E) schematic description of the role of
persistence time. Small persistence times facilitate rapid exploration of
routes to escape the cages constituted by their nearest neighbours, leading to
yielding at small active driving magnitudes. As the persistence time
increases, this capacity for exploration decreases and particles instead
“break through" their cages, which requires larger magnitude active forces.

####  II.4 Simulations in confinement

We next consider the influence on the yielding behaviour of confinement
geometry, a factor that may be relevant in several, including biological,
contexts. For example, experiments studying the dynamics and spatial
organisation of chromatin within the nucleus have pointed to the significant
role of confinement geometry, with possible implications on cell-state
transitions[36, 37]. Active particles in confinement have previously been
noted to display rich behaviour such as wall-aggregation and the formation of
cavities[22, 23]. We consider here a case of high density and investigate
yielding behaviour as the strength of active driving is increased. We compare
the yielding behaviour for two confinement geometries: (i) a nearly circular
confinement with aspect ratio (ratio of major axis to minor axis) $e=1.2$, and
(ii) a strongly elliptical confinement with aspect ratio $e=1.8$. Particles
experience attractive interactions of the Lennard-Jones form with the wall,
with parameters as given in Materials and Methods. Initial configurations with
different annealing are prepared by a procedure analogous to the bulk case as
described in Materials and Methods and the SI Appendix, Fig. S12. These
configurations are then subjected to active forces, at a persistence time of
$\tau_{p}=200\tau$ and the potential energy is tracked to identify a
transition from an absorbing to a diffusive state (see SI Appendix, Fig. S13).

The geometries considered, and the resulting yielding behaviour is shown in
Fig. 4. The results indicate that in going from a more circular to a more
elliptical geometry, the yielding point is pushed to larger values of the
active force. In both geometries, the post-yield states are ergodic, as
observed from the dependence of the potential energy per particle on $f$ in
Fig. 4 (B). It has been observed that confined active systems exhibit global
rotations, which should be subtracted in computing diffusive displacements of
particles. This is accomplished, following [38], by computing the cage-
relative mean-squared displacements, $\langle|r_{CR}|^{2}(t)\rangle$(see
Materials and Methods), which are shown in Fig. 4 (C) for the more circular
case, and in Fig. 4 (D) for the more elliptical case starting from well-
annealed samples. The superposition of global rotations and local
rearrangements results in a less reliable measure of diffusive motion despite
the correction employed. Nevertheless, one observes that the system remains in
a non-diffusive absorbing state in the pre-yield regime, whereas diffusive
motion sets in beyond the yielding $f$ value in each of the two geometries,
$f_{yield}\approx 0.3$ for $e=1.2$ and $f_{yield}\approx 0.6$ for $e=1.8$.
Such a large change in the yield value of the active force is not readily
explained by the perimeter length, of about 6%, and point to a sensitive
dependence of the fluidization transition on the confinement geometry. Such
dependence, including the role of wall interactions, needs to be explored in
greater detail.

These results show that at a given magnitude of the active forces, a change to
a more symmetric geometry can induce yielding and fluidization. They highlight
the role of confinement geometry in state transitions in assemblies subject to
active forces, including the biological assemblies mentioned in the
introduction, as discussed briefly below.

(a)  (b)  (c)  (d)

Figure 4: Dependence of yielding behaviour on confinement. Simulations are
performed in the two confinement geometries shown in (A), with eccentricities
of $e=1.2$ (more isotropic) and $e=1.8$, more anisotropic, respectively. Red
(blue) discs represent A(B) type particles. (B) The yielding diagram with the
steady state potential energy per particle at different values of the active
force in the two geometries, comparing the behaviour for a well-annealed case
($T_{p}=0.4$) with that of a poorly annealed case ($T_{p}=1.0$). The cage-
relative mean squared displacements, $r_{CR}^{2}(t)$, are shown for (C)
($e=1.2$ and (D) $e=1.8$, for $T_{p}=0.4$, the well annealed case. Note that
in the absorbing regime of low force, $r_{CR}^{2}(t)$, remains below $0.01$,
indicating an average displacement of $\approx~{}10\%\sigma_{AA}$ for each
particle with respect to its caging neighbours. The dashed teal lines in
panels (C) and (D) indicate a lineaar dependence characteristic of the
diffusive regime.

###  III Discussion

Our results demonstrate the existence of a strong analogy between cyclically
sheared amorphous assemblies and actively driven assemblies. Both modes of
drive result in a yielding transition from an absorbing state to a fluidised
state beyond a threshold. We find a strong correspondence in the yielding
phase diagram, with regimes of annealing and yielding to an “ergodic” state,
similar dependence on the initial state of the amorphous solids, the onset of
a diffusive state for large perturbations and the divergence in the time to
reach the steady state.

Active driving with a finite persistence time, and in confinement, as
implemented in this work, is motivated by our interest in biological systems
that display a transition from rigid to flowing states when subjected to
analogous mechanical conditions. Our results demonstrate that the persistence
time is indeed a key parameter, with an increase in the persistence time
leading to an increase in the critical perturbation strength (active force) at
which the system fluidizes. We also find that changing the confinement
geometry leads to a change in the yielding behaviour with the transition to
fluidised states occurring at much lower values of the active force for the
more isotropic confinement geometry. That confinement geometry alone affects
yielding behaviour in such a manner is, in our view, an interesting effect
with numerous possible implications in various contexts. While the model
system we investigate does not serve as a model for specific biological
assemblies, we believe that the results are of interest to contexts such as
the organization and dynamics of chromatin in the nucleus. In such contexts,
the interplay and complex causal relationships of the nature of active
driving, confinement and the fluidization state are important to elucidate in
order to develop faithful theoretical models. The analysis presented here, we
believe, forms an important first step in that direction.

###  IV Methods

#### Model - 2D binary mixture

We simulate a bidisperse Lennard-Jones mixture at constant $N,V$, with
periodic boundary conditions (wraparound), at high density, $\rho=1.2$ in 2D.
The interaction is truncated at a cut-off of
$r_{c}^{\alpha\beta}=\sigma_{\alpha\beta}$, where $\alpha$ and $\beta$ denote
the type of the interacting particles $i$ and $j$, either $A$ or $B$, and is
given by:

$\displaystyle\phi_{LJ}(r_{ij})$
$\displaystyle=\begin{cases}4\epsilon_{\alpha\beta}\left[\left({\sigma_{\alpha\beta}\over
r_{ij}}\right)^{12}-\left({\sigma_{\alpha\beta}\over
r_{ij}}\right)^{6}\right]\quad\text{if}\quad r_{ij}\leq r_{c}^{\alpha\beta}\\\
=0\quad\text{otherwise}\end{cases}$

We choose $N=1000$ particles in a $65:35$ mixture with $N_{A}=650$. Reduced
units are used throughout, in terms of $\sigma_{AA}$, $\epsilon_{AA}$ and the
particle mass $M$, with the unit of time
$\tau=\sqrt{M\sigma_{AA}^{2}/\epsilon_{AA}}$,
${\sigma_{AB}\over\sigma_{AA}}=0.8$, ${\sigma_{BB}\over\sigma_{AA}}=0.88$,
${\epsilon_{AB}\over\epsilon{AA}}=1.5$ and
${\epsilon_{BB}\over\epsilon{AA}}=0.5$.

#### Initial sample preparation

We perform constant temperature ($NVT$) molecular dynamics simulations using
the Nosé-Hoover thermostat with a time step of $dt=0.01$ in the LAMMPS
software suite[39] at different preparation temperatures, $T_{p}$, to generate
equilibrated liquid configurations. Energy minimum configurations (inherent
structures) are produced by performing conjugate gradient minimisation using
LAMMPS. The inherent structure energies show a $T_{p}$ dependence as shown in
the SI Appendix.

#### Shear simulations

Cyclic shear simulations are performed at a finite shear rate and temperature,
propagating the SLLOD equations of motion using the LAMMPS software suite[39].
Simulations are performed with a timestep of $dt=0.01$, with the strain being
varied sinusoidally, for a $100$ cycles for each strain amplitude
$\gamma_{max}$. The strain rate, computed at zero strain, is
$\dot{\gamma}=0.001$.

#### Equations of motion for active dynamics

Active dynamics are performed starting from the energy minimised
configurations by integrating the following stochastic equations of motion.

$\displaystyle\mathbf{\dot{p}}_{i}$
$\displaystyle=-\gamma\mathbf{p}_{i}+\sum\limits_{j\neq
i=1}^{N}\mathbf{f_{ij}}+f\mathbf{\hat{n}}_{i}+\mathbf{\xi_{t}}^{i}$
$\displaystyle\mathbf{\dot{x}}_{i}$ $\displaystyle=\mathbf{v}_{i}\equiv
M^{-1}\mathbf{p}_{i}$ $\displaystyle\mathbf{\dot{\theta}}_{i}$
$\displaystyle=\mathbf{\xi_{r}}^{i}$ (1)

The translational noise $\xi_{t}$ has zero mean and variance given by
$<~{}\xi_{t}^{i\alpha}(t)~{}\xi_{t}^{j\beta}(t^{{}^{\prime}})~{}>=2k_{B}TM\gamma\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{{}^{\prime}})$.
The rotational diffusion of the direction of active forcing likewise has zero
mean and variance given by
$<\xi_{r}^{i}(t)~{}\xi_{r}^{j}(t^{{}^{\prime}})>=2\tau_{p}^{-1}\delta_{ij}\delta(t-t^{{}^{\prime}})$.
The persistence time, $\tau_{p}$, determines the timescale of change in
$\mathbf{\hat{n}_{i}}$ ($(\cos(\theta_{i}),\sin(\theta_{i}))$ in 2D), the
direction of the active force with strength $f$. $\gamma$ is the dissipative
friction term, and $\mathbf{f_{ij}}=-\nabla\phi(\mathbf{\\{r_{ij}\\}})$ are
the interparticle forces. We integrate the equations of motion with a time
step of $dt=0.01$ following a BAOAB update scheme described in [40, 41].
Details are included in the SI Appendix.

#### Strain, stress and strain rate resulting from active forcing

The alignment of velocities with the active force direction produces a
deformation field with which we identify a strain rate, $\gamma_{act}$[15,
18].

$v_{par}=\left\langle\sum\limits_{i=1}^{N}\mathbf{v_{i}(t)}\cdot\mathbf{n_{i}(t)}\right\rangle$
(2)

Following the approach in [16], we consider the case of cyclic shear at finite
strain rate, where the affine component of the velocity of particle $i$ is of
the form $\dot{\gamma}\mathbf{r}_{i}\cdot\hat{y}$, and the magnitude of the
resultant field is given by $\dot{\gamma}[\sum(y_{i})^{2}]^{1/2}$. Discounting
centre of mass motion and assuming a uniform distribution of particles along
the $y$ direction one obtains a relationship of
$v_{par}^{cs}=\dot{\gamma}L_{y}\sqrt{N/12}$. Retaining the pre-factor as an
approximation for active dynamics, one can write the effective strain rate in
active systems in terms of the deformation field as:

$\dot{\gamma}_{act}=\frac{\sqrt{12}}{\sqrt{N}L}v_{par}^{act}$ (3)

One can then define an instantaneous strain step,
$\Delta\gamma=\dot{\gamma}_{act}\Delta t$. The derivative of the potential
energy with respect to the strain, $dPE/d\gamma$, yields the stress[16]
$\sigma_{act}$, which we estimate from the parametric dependence of $\Delta
PE$ with respect to $\Delta\gamma$ over a given timestep, $t\rightarrow
t+\Delta t$(see SI Appendix):

$\sigma_{act}=\frac{1}{L^{2}}\frac{dE}{d\gamma},$ (4)

while $L$ is the box length. We employ a similar procedure in cyclically
sheared simulations in order to compare the stress measured using this
procedure, with that obtained directly from the virial stress tensor at the
maximal strain value, $\gamma_{max}$.

#### Boundary interaction in confinement

We perform simulations of the same $65:35$ Lennard-Jones mixture of $N=1000$
particles in elliptical confinement of constant shape and volume. Each
particle interacts with the boundary wall with a truncated Lennard-Jones
interaction where the wall-specific parameters are given by
$\epsilon_{w}=0.5\epsilon_{\alpha\alpha}$,
$\sigma_{w}=0.5\sigma_{\alpha\alpha}$ and $r_{cw}=3\sigma_{AA}$ ($\alpha$
denotes the particle type, either $A$ or $B$).

The instantaneous distance of each particle from the confinement boundary is
obtained by implementing an iterative procedure is described in [42] (see SI
Appendix for details). Heatmaps showing the distance from an ellipse boundary
calculated using this scheme are shown in the SI Appendix. Visualisations of
the particle assembly in confinement are prepared using the OVITO software
suite[43].

#### Initial sample preparation in confinement

The initial configurations are prepared by first equilibrating a high
temperature liquid at $T=2.5$. Uncorrelated configurations are then thermally
annealed at different preparation temperatures $T_{p}$ to obtain samples in
confinement (see SI Appendix for time evolution during thermal annealing). The
annealed samples are then subjected to an instantaneous thermal quench to a
very low temperature, $T=10^{-4}$, to approximate an energy minimisation
procedure.

#### The cage-relative mean squared displacement

The cage-relative mean squared displacement[38] is given by:

$\displaystyle\Delta\mathbf{r}_{i,CR}(t)=\Delta\mathbf{r}_{i}-\frac{1}{N_{nn}}\sum\limits_{j=1}^{N_{nn}}\Delta\mathbf{r}_{j}(t)$
$\displaystyle\langle\|\Delta
r_{CR}(t)\|^{2}\rangle=\left\langle\frac{1}{N}\sum\limits_{i=1}^{N}\|\Delta\mathbf{r}_{i,CR}(t)\|^{2}\right\rangle$
(5)

where $\Delta\mathbf{r}_{i}(t)=\mathbf{r}_{i}(t)-\mathbf{r}_{i}(0)$, $N_{nn}$
is the number of nearest neighbours within a shell of
$1.5\sigma_{\alpha_{i}\alpha_{j}}$ from particle $i$ at time $t=0$. Further,
in order to correct for motion due to global rotations without rearrangement,
we identify the optimally rotated configuration,
$\\{\mathbf{R}(\theta_{min})\mathbf{r_{i}}(t)\\}$, where $\mathbf{R(\theta)}$
is a matrix that performs a rotation by $\theta$ in 2D. The angle,
$\theta_{min}$, is determined by an optimisation procedure that rotates each
particle along a scaled ellipse of identical aspect ratio as the confinement
geometry (but with scaled axes).

$\theta_{min}=\text{argmin}_{\theta}\sqrt{\langle\left[\mathbf{r_{i}(t+\Delta
t)}-\mathbf{R}(\theta)\mathbf{r_{i}(t)}\right]^{2}\rangle}.$ (6)

The accumulated displacement vector between $\\{\mathbf{r}(t+\Delta t)\\}$ and
the optimally rotated $\\{\mathbf{R}(\theta_{min})\mathbf{r_{i}}(t)\\}$, is
tracked and it’s magnitude gives the instantaneous cage relative mean squared
displacement.

$\displaystyle\Delta\mathbf{r}_{i,CR}(t)$
$\displaystyle=\Delta\mathbf{r}_{i,CR}(t-\Delta
t)+\Delta\mathbf{r}_{i}^{\theta_{min}}(t)-\frac{1}{N_{nn}}\sum\limits_{j=1}^{N_{nn}}\Delta\mathbf{r}_{j}(t)^{\theta_{min}}$
$\displaystyle\Delta\mathbf{r}_{i}^{\theta_{min}}(t)$
$\displaystyle=\mathbf{r_{i}}(t)-\mathbf{R}(\theta_{min})\mathbf{r_{i}}(t-\Delta
t)$ (7)

Note that while the initial centre of mass force and angular momentum is set
to zero by appropriately adjusting the active force orientation vectors, the
subsequent free diffusion of the orientation vectors does not prevent global
motion.

### References

  * Ramaswamy [2010] S. Ramaswamy, The mechanics and statistics of active matter, Annual Review of Condensed Matter Physics 1, 323 (2010), https://doi.org/10.1146/annurev-conmatphys-070909-104101 .
  * 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).
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, Bird flocks as condensed matter, Annual Review of Condensed Matter Physics 5, 183 (2014), https://doi.org/10.1146/annurev-conmatphys-031113-133834 .
  * Hatwalne  _et al._ [2004] Y. Hatwalne, S. Ramaswamy, M. Rao, and R. A. Simha, Rheology of active-particle suspensions, Phys. Rev. Lett. 92, 118101 (2004).
  * Bi  _et al._ [2015] D. Bi, J. H. Lopez, J. M. Schwarz, and L. M. Manning, A density-independent rigidity transition in biological tissues, Nature Physics 11, 1074–1079 (2015).
  * Bi  _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, Phys. Rev. X 6, 021011 (2016).
  * Das  _et al._ [2021] A. Das, S. Sastry, and D. Bi, Controlled neighbor exchanges drive glassy behavior, intermittency, and cell streaming in epithelial tissues, Phys. Rev. X 11, 041037 (2021).
  * Jülicher  _et al._ [2007] F. Jülicher, K. Kruse, J. Prost, and J.-F. Joanny, Active behavior of the cytoskeleton, Physics Reports 449, 3 (2007), nonequilibrium physics: From complex fluids to biological systems III. Living systems.
  * Shaebani  _et al._ [2020] M. Shaebani, A. Wysocki, R. Winkler, G. Gompper, and H. Rieger, Computational models for active matter, Nature Reviews Physics 2, 181–199 (2020).
  * Janssen [2019] L. M. C. Janssen, Active glasses, Journal of Physics: Condensed Matter 31, 503002 (2019).
  * Henkes  _et al._ [2011] S. Henkes, Y. Fily, and M. C. Marchetti, Active jamming: Self-propelled soft particles at high density, Phys. Rev. E 84, 040301 (2011).
  * Mandal  _et al._ [2016] R. Mandal, P. J. Bhuyan, M. Rao, and C. Dasgupta, Active fluidization in dense glassy systems, Soft Matter 12, 6268 (2016).
  * Mandal  _et al._ [2020] R. Mandal, P. J. Bhuyan, P. Chaudhuri, C. Dasgupta, and M. Rao, Extreme active matter at high densities, Nature communications 11, 2581 (2020).
  * Mo  _et al._ [2020] R. Mo, Q. Liao, and N. Xu, Rheological similarities between dense self-propelled and sheared particulate systems, Soft Matter 16, 3642 (2020).
  * Liao and Xu [2018] Q. Liao and N. Xu, Criticality of the zero-temperature jamming transition probed by self-propelled particles, Soft Matter 14, 853 (2018).
  * Morse  _et al._ [2021] P. K. Morse, S. Roy, E. Agoritsas, and M. L. Manning, A direct link between active matter and sheared granular systems, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2019909118 (2021).
  * Agoritsas [2021] E. Agoritsas, Mean-field dynamics of infinite-dimensional particle systems: global shear versus random local forcing, Journal of Statistical Mechanics: Theory and Experiment 2021, 033501 (2021).
  * Villarroel and Düring [2021] C. Villarroel and G. Düring, Critical yielding rheology: from externally deformed glasses to active systems, Soft Matter 17, 9944 (2021).
  * Amiri  _et al._ [2022] A. Amiri, C. Duclut, F. Jülicher, and M. Popović, Random traction yielding transition in epithelial tissues (2022), arXiv:2211.02159 [cond-mat.soft] .
  * Leishangthem  _et al._ [2017] P. Leishangthem, A. D. S. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nature Communications 8, 14653 (2017).
  * Bhaumik  _et al._ [2021] H. Bhaumik, G. Foffi, and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2100227118 (2021).
  * Yang  _et al._ [2014] X. Yang, M. L. Manning, and M. C. Marchetti, Aggregation and segregation of confined active particles, Soft matter 10, 6477 (2014).
  * Fily  _et al._ [2014] Y. Fily, A. Baskaran, and M. F. Hagan, Dynamics of self-propelled particles under strong confinement, Soft matter 10, 5609 (2014).
  * Vishen  _et al._ [2018] A. S. Vishen, J.-F. Rupprecht, G. Shivashankar, J. Prost, and M. Rao, Soft inclusion in a confined fluctuating active gel, Physical Review E 97, 032602 (2018).
  * Talwar  _et al._ [2013] S. Talwar, A. Kumar, M. Rao, G. I. Menon, and G. Shivashankar, Correlated spatio-temporal fluctuations in chromatin compaction states characterize stem cells, Biophysical journal 104, 553 (2013).
  * Makhija  _et al._ [2016] E. Makhija, D. Jokhun, and G. Shivashankar, Nuclear deformability and telomere dynamics are regulated by cell geometric constraints, Proceedings of the National Academy of Sciences 113, E32 (2016).
  * Uhler and Shivashankar [2017] C. Uhler and G. V. Shivashankar, Regulation of genome organization and gene expression by nuclear mechanotransduction, Nature Reviews Molecular Cell Biology 18, 717–727 (2017).
  * Xu  _et al._ [2021] H. Xu, J. C. Andresen, and I. Regev, Yielding in an amorphous solid subject to constant stress at finite temperatures, Physical Review E 103, 052604 (2021).
  * Hopkins  _et al._ [2022] A. Hopkins, M. Chiang, B. Loewe, D. Marenduzzo, and M. C. Marchetti, Local yield and compliance in active cell monolayers, Physical Review Letters 129, 148101 (2022).
  * Sharma and Karmakar [2023] R. Sharma and S. Karmakar, Activity-induced annealing leads to ductile-to-brittle transition in amorphous solids, arXiv preprint arXiv:2305.17545 (2023), arXiv:2305.17545 [cond-mat.soft] .
  * Nandi  _et al._ [2018] S. K. Nandi, R. Mandal, P. J. Bhuyan, C. Dasgupta, M. Rao, and N. S. Gov, A random first-order transition theory for an active glass, Proceedings of the National Academy of Sciences 115, 7688 (2018).
  * Mandal  _et al._ [2022] R. Mandal, S. K. Nandi, C. Dasgupta, P. Sollich, and N. S. Gov, The random first-order transition theory of active glass in the high-activity regime, Journal of Physics Communications 6, 115001 (2022).
  * Szamel  _et al._ [2015] G. Szamel, E. Flenner, and L. Berthier, Glassy dynamics of athermal self-propelled particles: Computer simulations and a nonequilibrium microscopic theory, Physical Review E 91, 062304 (2015).
  * Keta  _et al._ [2023] Y.-E. Keta, R. Mandal, P. Sollich, R. L. Jack, and L. Berthier, Intermittent relaxation and avalanches in extremely persistent active matter, Soft Matter 19, 3871 (2023).
  * Mandal and Sollich [2021] R. Mandal and P. Sollich, How to study a persistent active glassy system, Journal of Physics: Condensed Matter 33, 184001 (2021).
  * Roy  _et al._ [2017] B. Roy, S. Venkatachalapathy, P. Ratna, Y. Wang, D. S. Jokhun, M. Nagarajana, and G. V. Shivashankar, Laterally confined growth of cells induces nuclear reprogramming in the absence of exogenous biochemical factors, Proceedings of the National Academy of Sciences 115, E4741–E4750 (2017).
  * Roy  _et al._ [2020] B. Roy, L. Yuan, Y. Lee, and G. V. Shivashankr, Fibroblast rejuvenation by mechanical reprogramming and redifferentiation, Proceedings of the National Academy of Sciences 117, 10131 (2020).
  * Shiba  _et al._ [2016] H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, Unveiling dimensionality dependence of glassy dynamics: 2d infinite fluctuation eclipses inherent structural relaxation, Physical review letters 117, 245701 (2016).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computational physics 117, 1 (1995).
  * Leimkuhler and Matthews [2013] B. Leimkuhler and C. Matthews, Rational construction of stochastic numerical methods for molecular sampling, Applied Mathematics Research eXpress 2013, 34 (2013).
  * Allen and Tildesley [2017] M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford university press, 2017).
  * Chatfield [2017] C. Chatfield, A simple method for distance to ellipse, https://chatfield.io/simple-method-for-distance-to-ellipse/. (2017).
  * Stukowski [2009] A. Stukowski, Visualization and analysis of atomistic simulation data with ovito–the open visualization tool, Modelling and simulation in materials science and engineering 18, 015012 (2009).
  * Maisonobe [2006] L. Maisonobe, Quick computation of the distance between a point and an ellipse (2006).

##

### Supplementary Information

#### Integration of equations of motions for active dynamics

The discretized time integration is done following the prescription in [40,
41] for the BAOAB operator splitting to include the active force along the
orientation vector, $\hat{n}=(cos(\theta),sin(\theta))$.

Given the equations of motion:

$\displaystyle\mathbf{\dot{p}}_{i}$
$\displaystyle=-\gamma\mathbf{p}_{i}+\sum\limits_{j\neq
i=1}^{N}\left(\mathbf{f_{ij}}\right)+f\mathbf{\hat{n}}_{i}+\mathbf{\xi_{t}}^{i}$
$\displaystyle\mathbf{\dot{x}}_{i}$ $\displaystyle=\mathbf{v}_{i}\equiv
M^{-1}\mathbf{p}_{i}$ $\displaystyle\mathbf{\dot{\theta}}_{i}$
$\displaystyle=\mathbf{\xi_{r}}^{i}$ (S1)

we write

${\rm d}\left[\\!\\!\begin{array}[]{c}\mathbf{x}\\\ \mathbf{p}\\\
\theta\end{array}\\!\\!\right]=\underbrace{\left[\\!\\!\begin{array}[]{c}M^{-1}\mathbf{p}\\\
0\\\ 0\end{array}\\!\\!\right]{\rm d}t}_{\rm
A}+\underbrace{\left[\\!\\!\begin{array}[]{c}0\\\ \mathbf{-\nabla
U}+f\mathbf{{\hat{n}}}\\\ 0\end{array}\\!\\!\right]{\rm d}t}_{\rm
B}+\underbrace{\left[\\!\\!\begin{array}[]{c}0\\\ -\gamma\mathbf{p}{\rm
d}t+\sqrt{2\gamma Mk_{B}T}{\rm d}\mathbf{w}\\\ \sqrt{2/\tau_{p}}{\rm
d}\mathbf{w}\end{array}\\!\\!\right],}_{\rm O}$ (S2)

The stochastic component is modelled on the Weiner process $\mathbf{w}$, which
follows:

${\rm d}\mathbf{w}=\mathbf{w}(t+{\rm d}t)-w(t)=\sqrt{{\rm
d}t}\mathbf{GR(0,1;t)}$ (S3)

where $\mathbf{GR(0,1;t)}$ is a vector of independent, Gaussian distributed
random variables with means $0$ and variance $1$. For the active particle
system, one needs to include the active forces in the velocity update of the
BAOAB integrator, and as well perform the time integration for the
orientation, $\theta(t)$.

Defining constants,

$c_{1}=e^{-\gamma\delta t},\quad c_{2}=\gamma^{-1}(1-c_{1}),$

$c_{3}=\sqrt{k_{B}T(1-c_{1}^{2})},\quad c_{4}=\sqrt{2~{}dt/\tau_{p}},$

we have the update steps:

$\displaystyle\mathbf{p}_{i}\left(t+\frac{{\rm d}t}{2}\right)$
$\displaystyle=\mathbf{p}_{i}(t)+\frac{{\rm d}t}{2}\left[\,\mathbf{-\nabla
U(x}_{i}(t))+f\mathbf{{\hat{n}_{i}(t)}}\right];$ (S4a)
$\displaystyle\mathbf{x}_{i}\left(t+\frac{{\rm d}t}{2}\right)$
$\displaystyle=\mathbf{x}_{i}(t)+\frac{{\rm
d}t}{2}\,\mathbf{M^{-1}p}_{i}\left(t+\frac{{\rm d}t}{2}\right);$ (S4b)
$\displaystyle\mathbf{\hat{p}}_{i}\left(t+\frac{{\rm d}t}{2}\right)$
$\displaystyle=c_{1}\mathbf{p}_{i}\left(t+\frac{{\rm
d}t}{2}\right)+c_{3}\mathbf{M^{1/2}GR}_{i}(0,1;t+dt);$ (S4c)
$\displaystyle\theta_{i}(t+dt)$
$\displaystyle=\theta_{i}(t)+c_{4}\mathbf{GR}^{{}^{\prime}}_{i}(0,1;t+dt);$
(S4d) $\displaystyle\mathbf{x}_{i}(t+{\rm d}t)$
$\displaystyle=\mathbf{x}_{i}\left(t+\frac{{\rm d}t}{2}\right)+\frac{{\rm
d}t}{2}\,\mathbf{M^{-1}\hat{p}}_{i}\left(t+\frac{{\rm d}t}{2}\right);$ (S4e)
$\displaystyle\mathbf{p}_{i}(t+{\rm d}t)$
$\displaystyle=\mathbf{\hat{p}}_{i}\left(t+\frac{{\rm
d}t}{2}\right)+\frac{{\rm d}t}{2}\left[\,\mathbf{-\nabla U(x}_{i}(t+{\rm
d}t))+f\mathbf{{\hat{n}}}_{i}(t+{\rm d}t)\right].$ (S4f)

#### Sample preparation

We show the inherent structure potential energies, $E_{IS}$, as a function of
time as well as the dependence of the steady state inherent structure
energies, $E_{IS}$, as a function of the temperature at which the simulation
is performed, $T_{p}$, in Fig. S1.

Figure S1: Average inherent structure energy per particle as a function of
time for $N=1000$ particles of the binary Lennard-Jones mixture in 2D,
simulated in the NVT ensemble at different temperatures. The steady state
average inherent structure energy per particle, $E_{IS}$, is shown as a
function of simulation temperature for the cases at which a steady state is
reached, i.e., $T~{}\geq~{}0.45$.

#### Yielding phase diagram at different $\mathbf{\tau_{p}}$

In Fig. S2 we show the full yielding phase diagram with active dynamics at
$\tau_{p}=1.01\times 10^{4}\tau$ and $\tau_{p}=1.01\times 10^{5}\tau$.

Figure S2: Yielding phase diagram for $N=1000$ particles of the $65:35$ binary
Lennard-Jones mixture in 2D subjected to active dynamics at
$\tau_{p}=1.01\times 10^{4}\tau$ and $\tau_{p}=1.01\times 10^{5}\tau$ (right).

#### Parameteric dependence of change in potential energy, $\Delta E$, on
strain step $\Delta\gamma$

The alignment of instantaneous particle velocities with the respective
directions of active forcing yield a measure of the instantaneous strain rate
as in Eq. S5 below [18]. In systems driven by active forces or using a shear
stress, the strain rate response changes drastically across the yield point
[28, 18]. We define the instantaneous strain rate from the alignment of
velocities with the active force direction as

$\displaystyle v_{par}^{act}$
$\displaystyle=\frac{1}{N}\left\langle\sum\limits_{i=1}^{N}\mathbf{v_{i}(t)}\cdot\mathbf{n_{i}(t)}\right\rangle_{t}$
$\displaystyle\dot{\gamma}_{act}$
$\displaystyle=\frac{12\sqrt{N}}{L}v_{par}^{act}.$ (S5)

From the instantaneous strain rate, $\dot{\gamma}_{act}(t)$, we define an
instantaneous step $\Delta\gamma=\dot{\gamma}_{act}\Delta t$ and consider the
change in potential energy with respect to it. We assume a linear parametric
dependence of $\Delta E$ on $\Delta\gamma$ of the form $y=\beta x+\xi$, where
$\xi$ is a noise term, and identify the slope $\beta$, taken as the Pearson
correlation times the ratio of the standard deviations,
$\sigma_{\Delta\gamma}/\sigma_{\Delta E}$. As mentioned in the main text, this
is an approximate procedure that provides estimates of the stress using the
response of the change in potential energy to changes in the measured strain
step. As shown in Fig. 1 (d) of the main manuscript, and discussed below, such
an approximation is found to give estimates in reasonable agreement close to
the maximal strain $\gamma_{max}$ when used in the context of cyclic shear
simulations at finite rate.

#### Development of non-zero strain rate in active yielding

In Fig. S3, we show the average measured strain rate in the steady state
condition for differently annealed samples subjected to active driving at two
values of the persistence time, $\tau_{p}$. We find that the strain rate is
approximately zero below the yield value of the active force and increases
monotonically beyond. The point at which this change in behaviour occurs
shifts with the yield point as $\tau_{p}$ is changed. In both cases, the
increase in strain rate occurs for $f$ at the yielding value identified from
Fig. 1 in the main manuscript, instead occurring at slightly larger $f$.

Figure S3: We compute the resultant strain rate at different values of the
active force, $f$, using the definition in Eq. S5. The strain rate fluctuates
around $0$ at all pre-yield values of force, subsequently rising non-linearly
beyond the yield point (Top). The onset of non-zero strain rate shifts with
the yield point at higher persistence time, $\tau_{p}$, (Bottom).

#### Calculating stress for cyclic shear from the parametric dependence of
$\Delta E$ on $\Delta\gamma$

In Fig. 1 (d) of the main text, we show the dependence of the maximum stress
on $\gamma_{max}$ in the case of cyclic shear, where the stress is computed by
finding the slope of the best linear fit to the scatter of $\Delta E$ vs
$\Delta\gamma$ from cyclic shear simulations at constant shear rate. The
velocity stress component, $\sum\limits_{i=1}^{N}\delta v_{x}^{i}\delta
v_{y}^{i}$ (the $\delta$ represents change with respect to any global flow),
is subsequently added. This stress value is compared to that obtained directly
from the virial stress tensor for $\gamma$ near $\gamma_{max}$. More
specifically, we first consider a reference $\gamma$ close to $\gamma_{max}$
($\gamma_{max}-10^{-4}$ and note the energy $E_{ref}$. We then consider the
energy of configurations within a neighbourhood of $\Delta\gamma=10^{-4}$ of
$\gamma_{ref}$ (identified from the strain rate multiplied with time windows
$\Delta t$) and store the corresponding $\Delta\gamma$ and $\Delta
E=E-E_{ref}$. The parametric dependence of $\Delta E$ and $\Delta\gamma$ is
used as in active systems to determine the slope. Such a procedure is limited
in this case to conditions where the strain step is very small, such that
plastic rearrangements and other non-affine motion does not obscure the
dependence of $\Delta E$ on the strain step. The data of $\Delta E$ vs
$\Delta\gamma$ is aggregated over $10$ steady state cycles from $16$
independent simulations.

#### Diffusivity beyond the yield point

In Fig. S4, we show the diffusivities extracted from the mean squared
displacements shown in Fig. 2 of the main manuscript. Data is shown for
trajectories simulated with large active forces displaying a diffusive regime
where the scaling of the MSD with time is close to, or exactly linear. For
both the well-annealed and poorly annealed cases, a free fit of the MSD yields
an exponent less than $0.1$ in the pre-yield regime where the MSD is nearly
flat. For the case of $T_{p}=1.00$ and $f=0.7$, the measured exponent is close
to $0.8$, which we interpret as reflecting the difficulty of reaching the
steady state in this case. For the case of $T_{p}=1.00$ and $f=0.7$ (not
shown), the measured exponent is close to $0.8$, which we interpret as
reflecting the difficulty of reaching the steady state in this case. A best
fit through the aggregated data in Fig. S4 reveals that the diffusivity
vanishes at $f<0.5$. However, the data in Fig. 1 and Fig. 2 in the main text
show that the system remains in an absorbing state up to a higher value of
$f=0.6$, suggesting a discontinuity in the range between $f=0.6$ and $0.8$.

Figure S4: Diffusivity obtained from the mean squared displacement data in
Fig. 2 of the main manuscript shown for cases where the MSD shows at or near
linear scaling with time.

#### Energy relaxation and timescale to steady state

We show the time series of average potential energy per particle over time for
$8$ independent trajectories simulated with different values of the active
force. The energy vs time curves are fit to stretched exponential fits as
shown in Fig. S5 for configurations prepared at $T_{p}=0.35$ and in Fig. S6
for configurations prepared at $T_{p}=1.00$. The fitting function has the
following form:

$\langle E(t)\rangle=E_{0}+\Delta E\exp(-(t/\tau)^{\beta}),$ (S6)

controlling for $E_{0}$ and $\Delta E$ based on the bounds in the trajectory.
In each case, the region used for the fit includes the time interval
immediately preceding the sharp change that signals the transition and the
steady state reached after the transition (see for example the case of
$T_{p}=1.00$, $\tau_{p}=1.01\times 10^{5}\tau$ and $f\geq 1.2$, where the
system undergoes an initial relaxation before transitioning to the fluidised
state). Simulations performed with $f$ values close to the transition point
($f=0.8$ in Fig. S5 (A), $f=0.9$ in Fig. S5 (B) and $f=1.0$ in Fig. S5 (C))
reach the steady state very slowly, as a result of which the observed final
values of $\langle PE/N\rangle$, reported in Fig. 3 of the main manuscript,
are apparently intermediate to the values in the pre-yield and the post-yield
branches.

Figure S5: Potential energy vs time, averaged over $8$ independent
simulations, starting from a well-annealed sample at $T_{p}=0.35$ at $3$
values of the persistence time, $\tau_{p}$. Data are obtained from active
dynamics simulations performed with a timestep of $dt=0.01$. Dashed lines are
stretched exponential fits to the data.

Figure S6: Potential energy vs time, averaged over $8$ independent
simulations, starting from a poorly-annealed sample at $T_{p}=1.00$ at $3$
values of the persistence time, $\tau_{p}$. Data are obtained from active
dynamics simulations performed with a timestep of $dt=0.01$. Dashed lines are
stretched exponential fits to the data.

#### Comparison of timescale to steady state from stretched exponential fits
and from first passage time identification

In Fig. S7, we show the average time to reach state obtained from two
different procedures. In the first we consider the time series, $\langle
E\rangle$ vs time, averaged over $8$ independent simulations and fitted a
stretched exponential to the curve, to identify the relaxation timescale
$\tau$ from Eq. S6. Next, we consider the average first passage time to the
steady state energy $<E_{ss}^{i}>_{t_{ss}}$ in each trajectory $i$. Here, the
average is taken over a window of time $t_{ss}$ where the system is in its
final steady state. The average of these first passage times is then taken
over $8$ independent simulations.

We note that for the poorly annealed case, trajectories undergo further slow
annealing at small active forces, reminiscent of creep, after they enter into
the neighbourhood of its final steady state. Likewise, for trajectories close
to the transition point, the approach to steady state is extremely slow. This
leads to the appearance of a shorter first passage time to steady state than
that obtained from estimates of the relaxation timescale in a stretched
exponential process, as seen in Fig. S7. While this may be ameliorated by
improved estimates of the steady state energy, $\langle
E_{ss}^{i}\rangle_{t_{ss}}$, reached by each trajectory, we opt for the
stretched exponential procedure in our work in order to identify an asymptotic
steady state energy. Trajectories subjected to active force significantly
larger than the critical force yield similar estimates for $\langle
t_{ss}\rangle$ from these two approaches.

Figure S7: Timescale to steady state shown at different values of $\tau_{p}$
for well-annealed samples at $T_{p}=0.35$ subjected to active dynamics at
forces greater than the yielding force. The hollow symbols with dashed lines
are obtained from stretched exponential fits to the aggregate $\langle
E\rangle$ vs time curves where the energy at each time point is the average
over $8$ independent simulations. The solid symbols with dotted lines are
obtained from identifying the average of the first passage time to the steady-
state value of the potential energy $E_{ss}$ in each trajectory, where the
average is over $8$ independent simulations.

#### Comparison with passive dynamics

We address the question of whether purely passive dynamics at comparable
temperatures can lead to fluidisation in the same way that active dynamics
does. In order to do this, we consider the average kinetic energy per particle
across different cases. In Fig. S8, we show the average potential energy per
particle for active dynamics where the active direction re-orients every
timestep, $dt$, which we use to represent the limit of zero persistence where
the direction of active force is delta-correlated in time. In this limit of
small persistence, the active temperature is given by the expression
$T_{a}=f^{2}\tau_{p}/4$ in 2D and the effective temperature is given by
$T_{eff}=T_{bath}+T_{a}$. This relation holds exactly in the small persistence
time limit, as we see in Fig. S8. We reiterate here that the kinetic energies
at finite persistence are considerably lower than that for the small
persistence time case. One does not, however, expect the same expression for
$T_{a}$ to hold at large $\tau_{p}$. It has been argued in the literature that
$T_{a}$ goes as $Af^{2}\tau_{p}/(1+B\tau_{p})$[31, 32], but results reported
in [13, 32] suggest that for the large persistence times we consider in the
present manuscript, the above expression may not be valid, with the arrest
line exhibiting non-monotonic dependence on $f$.

Figure S8: No discontinuity in energy for zero persistence active force with a
quadratic dependence of potential energy indicative of increased thermal
motion. Quadratic dependence of kinetic energy on active force.

In Fig. S9 we show the average kinetic energy for simulations at different
active force at finite and small $\tau_{p}$. The average kinetic energy per
particle remains low in these cases, and this is underscored by the fact that
it remains below the kinetic energy for the small persistence case, where the
kinetic energy can be compared with an exact expression of the effective
temperature. The average kinetic energy per particle is higher in this low
persistence case compared to that at finite persistence time.

Figure S9: Effective temperature from the kinetic energy per particle for
differently annealed samples subjected to active forces at either finite
persistence time, or at zero persistence for $T_{p}=0.35$.

In order to determine whether these, or comparable, temperatures are
sufficient to induce fluidisation with passive dynamics, we subject
differently annealed samples, i.e., different parent temperature $T_{p}$, to
$4\times 10^{7}$ MD steps of Langevin dynamics at constant volume at a range
of simulation temperatures, $T_{sim}$, the highest of which ($T_{sim}=0.1$),
exceeds the highest observed average kinetic energy ($\approx~{}0.01$) by
almost an order of magnitude. In Fig. S10 A, we show the potential energy per
particle averaged over $4$ independent samples at different thermostat
temperatures, $T_{sim}$. We find that a steady state is reached in each case
with a linear dependence of the potential energy per particle on simulation
temperature (inset of Fig. S10 A). Examining the inherent structure energies
per particle in each case shows a constant energy at even $T_{sim}=0.1$ for
$T_{p}=0.35$, and thermal annealing for $T_{p}=1.00$, as shown in Fig. S10 B.
Note that the highest simulation temperature, $T_{sim}$, exceeds the highest
$KE~{}/~{}N$ in active simulations by a factor of $\sim 25$. Finally, we
consider trajectories at the highest simulation temperature, $T_{sim}=0.1$,
and observe the mean squared displacements for the two samples, finding that
the trajectories are in an absorbing state, as shown in Fig. S10 C. Based on
this, we conclude that the fluidisation observed in the actively driven
simulations at finite persistence is not correlated with the temperature of
the system reflected in the average kinetic energy per particle. The
nontrivial scaling of the kinetic energy with active force $f$ and $\tau_{p}$
has been discussed in [13] and more generally, the departure from equilibrium-
like behaviour for large persistence times in [8].

Figure S10: Bi-disperse assemblies heated through the temperatures
corresponding to the average kinetic energy of the actively driven systems. We
perform NVT MD of the bidisperse assembly (see methods) and heat the
differently annealed initial samples directly to target temperatures,
$T_{sim}$, in the desired range. We find a linear dependence between potential
energy and temperature (left). A steady state is reached in each case with the
inherent structure energies being either constant, or showing annealing
depending on whether we start from a high initial energy or a low initial
energy (middle). The MSD obtained after ignoring a transient of $2\times
10^{7}$ MD steps for the raw trajectories shows no diffusive behaviour at even
the highest target temperature (right).

#### Confinement boundary implementation

In order to compute the wall interaction, we need to identify the vector from
the centre of any particle, $p$, to the closest point, $q$, on the elliptical
confinement boundary. We follow an iterative procedure with the following
steps [44, 42]:

  1. 1.

Pick a random initial $a$ on the ellipse.

  2. 2.

Compute the distance $r_{pa}$ between $p$ and $a$.

  3. 3.

The circle of radius $r_{pa}$ centred at $p$ will intersect the ellipse at
another point $b$.

  4. 4.

Identify the midpoint $a^{\prime}$ between $a$ and $b$.

  5. 5.

Calculate $r_{pa^{\prime}}$ and go to step 3.

In order to follow this procedure, we describe the relevant details from [42]
for completeness:

$\displaystyle x$ $\displaystyle=a\cos(t)$ $\displaystyle y$
$\displaystyle=b\sin(t)$ (S7)

We can approximate the curvature of the ellipse at $a$ to a circle centred at
$ev_{x},ev_{y}$ with radius $|\omega_{a}|$. The centres of curvature of the
ellipse are obtained using

$\displaystyle ev_{x}$ $\displaystyle=\frac{(a^{2}-b^{2})}{a}\cos^{3}(t)$
$\displaystyle ev_{y}$ $\displaystyle=\frac{(b^{2}-a^{2})}{b}\sin^{3}(t)$ (S8)

Vectors from the centre of curvature to $a$ and to $a^{\prime}$, $\omega_{a}$
and $\omega_{a^{\prime}}$ are separated by an arclength $\Delta c$. We can
write:

$sin\left(\frac{\Delta
c}{|\omega_{a}|}\right)=\frac{\omega_{a}\times\omega_{a^{\prime}}}{|\omega_{a}||\omega_{a^{\prime}}|}$
(S9)

One can then write

$\displaystyle\frac{dc}{dt}$
$\displaystyle=\sqrt{\left(\frac{dx}{dt}\right)^{2}-\left(\frac{dy}{dt}\right)^{2}}$
$\displaystyle\frac{\Delta c}{\Delta t}$
$\displaystyle=\sqrt{a^{2}\sin^{2}(t)+b^{2}\cos^{2}(t)}$ (S10)

We can then identify $\Delta t$ as

$\Delta t=\frac{\Delta c}{\sqrt{a^{2}+b^{2}-x^{2}-y^{2}}}$ (S11)

Finally, one can find the coordinates of $a^{\prime}$ as

$\displaystyle x^{\prime}$ $\displaystyle=a\cos(t+\Delta t)$ $\displaystyle
y^{\prime}$ $\displaystyle=b\sin(t+\Delta t)$ (S12)

Once $a^{\prime}$ has been identified, we can iterate the procedure till
convergence is achieved. In Fig. S11, we show heatmaps of the distance
computed from ellipses of different aspect ratio, including the symmetric
circular case.

Figure S11: Heat map of distance from the boundary of the ellipse, as computed
by the optimisation procedure to identify the closest point on a given ellipse
from a chosen query point.

#### Sample preparation in confinement

In order to prepare initial samples in confinement, we first thermally anneal
the binary mixture at different simulation temperatures, $T$, shown in Fig.
S12.

Figure S12: Evolution of potential energy per particle over time for a system
in two confinement geometries, $e=1.2$ (left), and $e=1.8$ (right) at
different preparation temperatures, $T_{p}$.

Steady state configurations from the annealed trajectories are then subjected
to an instantaneous thermal quench to $T_{q}=0.0001$ and simulated for
$5\times 10^{4}$ MD steps to approximate a mechanically stable initial
configuration.

#### Time evolution of energy in confinement

In Fig. S13 we show the time evolution of the average energy per particle for
the two confinement geometries for both the well-annealed sample ($T_{p}=0.4$)
and the poorly annealed one ($T_{p}=1.0$). Well-annealed samples show a
constant $\langle PE/N\rangle$ with increasing $f$ for small $f$, followed by
a steep rise as $f$ is increased beyond the yield point. For the poorly
annealed sample, the system undergoes further annealing with increasing $f$
upto the yield point, with a corresponding monotonic decrease in $\langle
PE/N\rangle$. The time taken to reach the steady state is longer close to the
yield point, which we interpret will introduce difficulty in obtaining
reliable estimates of the final steady state value of the potential energy, as
can be seen for the case of $f=0.6$ in Fig. S13 (a) and (c) and for the case
of $f=0.8$ in Fig. S13 (b) and (d).

(a)  (b)  (c)  (d)

Figure S13: Time evolution of average energy per particle at the values of $e$
and $T_{p}$ shown in the plots, averaged in logarithmically spaced bins.
Curves reach a steady state for all values of active driving force, $f$,
except those close to the respective transition values.

### References

  * [1] Benedict Leimkuhler and Charles Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Applied Mathematics Research eXpress, 2013(1):34–56, 2013.
  * [2] Michael P Allen and Dominic J Tildesley.  Computer simulation of liquids.  Oxford university press, 2017.
  * [3] Carlos Villarroel and Gustavo Düring.  Critical yielding rheology: from externally deformed glasses to active systems.  Soft Matter, 17:9944–9949, 2021.
  * [4] Haiyan Xu, Juan Carlos Andresen, and Ido Regev.  Yielding in an amorphous solid subject to constant stress at finite temperatures.  Physical Review E, 103(5):052604, 2021.
  * [5] 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.  Proceedings of the National Academy of Sciences, 115(30):7688–7693, 2018.
  * [6] Rituparno Mandal, Saroj Kumar Nandi, Chandan Dasgupta, Peter Sollich, and Nir S Gov.  The random first-order transition theory of active glass in the high-activity regime.  Journal of Physics Communications, 6(11):115001, 2022.
  * [7] Rituparno Mandal, Pranab Jyoti Bhuyan, Pinaki Chaudhuri, Chandan Dasgupta, and Madan Rao.  Extreme active matter at high densities.  Nature communications, 11(1):2581, 2020.
  * [8] Étienne Fodor, Cesare Nardini, Michael E. Cates, Julien Tailleur, Paolo Visco, and Frédéric van Wijland.  How far from equilibrium is active matter?  Phys. Rev. Lett., 117:038103, Jul 2016.
  * [9] L Maisonobe.  Quick computation of the distance between a point and an ellipse, 2006\.
  * [10] Carl Chatfield.  A simple method for distance to ellipse.  https://chatfield.io/simple-method-for-distance-to-ellipse/., 2017\.


# Scalar Mesons and Axial-vector Mesons Via Dyson-Schwinger Equation and
Bethe-Salpeter Equation Approach 111This work is supported by: the National
Natural Science Foundation of China under contracts No. 11947108 and No.
12005060.

Muyang Chen 1<EMAIL_ADDRESS>1 Department of Physics, Hunan Normal
University, Changsha 410081, China

###### Abstract

The spectrum of scalar mesons and axial-vector mesons are systematically
studied via Dyson-Schwinger equation and Bethe-Salpeter equation approach in
the rainbow-ladder approximation. An interaction model with a repulsive term
added to the one used for the pseudoscalar and vector mesons is proposed. The
Dyson-Schwinger equation results are consistent with the experiment values and
other model results, which shows that this interaction model is effective for
all the heavy, heavy-light, and light scalar and axial-vector mesons. This is
the first time a systematic study of the spectrum of the scalar mesons and
axial-vector mesons in the Dyson-Schwinger equation approach is achieved.

###### keywords:

Dyson-Schwinger Equation, Bethe-Salpeter Equation, Scalar Meson, Axial Vector
Meson

pacs

##  1 Introduction

A systematic study of the hadrons in a Quantum Chromdynamics (QCD) related way
is important for the hadron phenomenon. The Dyson-Schwinger equations (DSEs)
[1, 2] provide such a way. This approach preserves important properties of QCD
such as dynamical chiral symmetry breaking (DCSB), and many phenomenological
successfulness have been gained [3, 4].

Studying the meson spectrum via Dyson-Schwinger equation and Bethe-Salpeter
equation (DSBSE) approach [1, 2, 5] using realistic interaction begun 20 years
ago. A gaussian interaction model in the rainbow ladder (RL) approximation is
quite successful for the light pseudoscalar and vector mesons [6, 7, 8, 9].
This model was extended to study the heavy mesons [10, 11] and heavy-light
mesons [12] in recent years. Combining with the Faddeev equation, this
approach can also study the baryons [3, 13].

The DSBSE approach can be briefly stated as following. The quark propagator is
solved by the Gap equation [1, 2, 4] 222We work in the Euclidean space, where
the inner product of the four vector is defined by $a\cdot
b=\delta_{\mu\nu}a_{\mu}b_{\nu}=\Sigma_{i=1}^{4}a_{i}b_{i}$, with
$\delta_{\mu\nu}$ being the Kronecker delta. The Dirac matrices satisfy the
algebra $\\{\gamma_{\mu},\gamma_{\nu}\\}=2\delta_{\mu\nu}$, and
$\gamma_{5}=-\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}$.,

$\displaystyle S_{f}^{-1}(k)$ $\displaystyle=$ $\displaystyle
Z_{2}(i\gamma\cdot k+Z_{m}m_{f})$ (1)
$\displaystyle+\frac{4}{3}\bar{g}^{2}Z_{1}\int^{\Lambda}_{dq}D_{\mu\nu}(l)\gamma_{\mu}S_{f}(q)\Gamma^{f}_{\nu}(k,q),$

where $f=\\{u,d,s,c,b,t\\}$ represents the quark flavor. $S_{f}(k)$ is the
quark propagator, which can be decomposed as
$S_{f}(k)=\frac{Z_{f}(k^{2})}{i\not{k}+M_{f}(k^{2})}$. $Z_{f}(k^{2})$ is the
quark dressing function and $M_{f}(k^{2})$ is the quark mass function.
$l=k-q$. $m_{f}$ is the current quark mass. $\Gamma^{f}_{\nu}$ is the quark-
gluon-vertex. $D_{\mu\nu}$ is the gluon propagator. $\bar{g}$ is the coupling
constant. $Z_{1}$, $Z_{2}$, $Z_{m}$ are the renormalisation constants of the
quark-gluon-vertex, the quark field and the quark mass respectively.
$\int^{\Lambda}_{dq}=\int^{\Lambda}d^{4}q/(2\pi)^{4}$ stands for a
Poincar$\acute{\text{e}}$ invariant regularized integration, and $\Lambda$ is
the regularization scale.

The meson Bethe-Salpeter amplitude (BSA), $\Gamma^{fg}(k;P)$, is solved by the
Bethe-Salpeter equation (BSE). $k$ and $P$ are the relative and total momentum
of the meson. $P^{2}=-M^{2}_{fg}$ and $M_{fg}$ is the mass of the meson with
valence quark flavor $(f,g)$. The BSE expresses as [5, 4]

$\big{[}\Gamma^{fg}(k;P)\big{]}^{\alpha}_{\beta}=\int^{\Lambda}_{dq}\big{[}K^{fg}(k,q;P)\big{]}^{\alpha\delta}_{\sigma\beta}\big{[}\chi^{fg}(q;P)\big{]}^{\sigma}_{\delta},$
(2)

where $K^{fg}(k,q;P)$ is the quark-antiquark scattering kernel, and $\alpha$,
$\beta$, $\sigma$ and $\delta$ are the Dirac indexes.
$\chi^{fg}(q;P)=S_{f}(q_{+})\Gamma^{fg}(q;P)S_{g}(q_{-})$ is the wave
function. $q_{+}=q+\eta P/2$, $q_{-}=q-(1-\eta)P/2$. $\eta$ is the
partitioning parameter describing the momentum partition between the quark and
antiquark and varying its value dosen’t affect the physical observables.

The RL approximation is replacing the self energy in Eq. (1) and the
scattering kernel in Eq. (2) by

$\displaystyle\bar{g}^{2}Z_{1}D_{\mu\nu}(l)\Gamma^{f}_{\nu}(k,q)\to[Z_{2}]^{2}\tilde{D}^{f}_{\mu\nu}(l)\gamma_{\nu},$
(3)
$\displaystyle\big{[}K^{fg}(k,q;P)\big{]}^{\alpha\delta}_{\sigma\beta}\to-\frac{4}{3}[Z_{2}]^{2}\tilde{D}^{fg}_{\mu\nu}(l)[\gamma_{\mu}]^{\alpha}_{\sigma}[\gamma_{\nu}]^{\delta}_{\beta},$
(4)

where
$\tilde{D}^{f}_{\mu\nu}(l)=\left(\delta_{\mu\nu}-\frac{l_{\mu}l_{\nu}}{l^{2}}\right)\mathcal{G}^{f}(l^{2})$
and
$\tilde{D}^{fg}_{\mu\nu}(l)=\left(\delta_{\mu\nu}-\frac{l_{\mu}l_{\nu}}{l^{2}}\right)\mathcal{G}^{fg}(l^{2})$
are the flavor dependent effective interaction between the quark and
antiquark. In Ref. [12] the dressing function $\mathcal{G}^{f}(l^{2})$ and
$\mathcal{G}^{fg}(l^{2})$ are modeled as

$\displaystyle\mathcal{G}^{f}(s)$ $\displaystyle=$
$\displaystyle\mathcal{G}^{f}_{IR}(s)+\mathcal{G}_{UV}(s),$ (5)
$\displaystyle\mathcal{G}^{f}_{IR}(s)$ $\displaystyle=$ $\displaystyle
8\pi^{2}\frac{D_{f}^{2}}{\omega_{f}^{4}}e^{-s/\omega_{f}^{2}},$ (6)
$\displaystyle\mathcal{G}^{fg}(s)$ $\displaystyle=$
$\displaystyle\mathcal{G}^{fg}_{IR}(s)+\mathcal{G}_{UV}(s),$ (7)
$\displaystyle\mathcal{G}^{fg}_{IR}(s)$ $\displaystyle=$ $\displaystyle
8\pi^{2}\frac{D_{f}}{\omega_{f}^{2}}\frac{D_{g}}{\omega_{g}^{2}}e^{-s/(\omega_{f}\omega_{g})},$
(8) $\displaystyle\mathcal{G}_{UV}(s)$ $\displaystyle=$
$\displaystyle\frac{8\pi^{2}\gamma_{m}\mathcal{F}(s)}{\text{ln}[\tau+(1+s/\Lambda^{2}_{QCD})^{2}]},$
(9)

where $s=l^{2}$. $\mathcal{G}^{f}_{IR}(s)$ and $\mathcal{G}^{fg}_{IR}(s)$ are
the infrared interaction responsible for DCSB. $1/\omega_{f}$ is the
interaction width and $D_{f}$ expresses the infrared interaction strength.
$\mathcal{G}_{UV}(s)$ keeps the one-loop perturbative QCD limit in the
ultraviolet. $\mathcal{F}(s)=[1-\exp(-s^{2}/[4m_{t}^{4}])]/s$,
$\gamma_{m}=12/(33-2N_{f})$, with $m_{t}=1.0\textmd{ GeV}\,$, $\tau=e^{10}-1$,
$N_{f}=5$, and $\Lambda_{\text{QCD}}=0.21\textmd{ GeV}\,$. This model turned
out to be successful for all the ground state pseudoscalar and vector mesons,
from heavy, heavy-light to light mass scale.

However, the orbital and radial excited states are less studied. The
interaction models successful for pseudoscalar and vector mesons usually don’t
produce a good spectrum for the orbital and radial excited states [9, 14].
Even in the heavy quark region, the masses of the excited states are lower
than the experiment data systematically [10, 11, 15, 16, 17] because the
repulsive interaction from higer order corrections of the quark-gluon vertex
is lacking [18]. One way to solve this problem is using beyond RL
approximations. The effects of beyond RL approximations have been researched
for years [19, 20, 21]. However, no systematic study of the hadron spectrum
was achieved due to the technique difficulties. In this background,
phenomenological method is of practical significance. In this work, I follow
the latter way. The previous model, Eq. (5) $\sim$ (9), is extended to study
the scalar mesons and axial-vector mesons. In section 2, I introduce the
modified interaction model. In section 3, I show the results and give some
discussions. In section 4, I give a summary.

##  2 The modified model

To make the interaction more reliable for scalar mesons and axial-vector
mesons, Eq. (8) is modified into

$\mathcal{G}^{fg}_{IR}(s)=8\pi^{2}(\frac{D_{f}}{\omega_{f}^{2}}\frac{D_{g}}{\omega_{g}^{2}}-\Delta_{fg})e^{-s/(\omega_{f}\omega_{g})},$
(10)

where $\Delta_{fg}>0$ represents a repulsive interaction that is missing in
the RL approximation. An interaction universial for different mass scales
should depend on the quark masses. The following form is proposed,

$\Delta_{fg}=\frac{\delta}{[M_{f}(0)M_{g}(0)]^{\alpha/2}(\omega_{f}\omega_{g})^{\beta/2}},$
(11)

where $M_{f}(0)$ is the value of the quark mass function at $k^{2}=0$.
$\Delta_{fg}$ in Eq. (11) depends on the quark masses and the interaction
width explicitly. $\delta$, $\alpha$ and $\beta$ are 3 parameters to be fitted
by the observables.

In this paper, I consider 5 quark flavors and work in the isospin symmetry
limit. So $f,g=\\{n,s,c,b\\}$, with $n=u\,\textmd{or}\,d$. Herein I choose the
masses of scalar ($J^{PC}=0^{++}$) charmonium ($M_{c\bar{c}}$), bottomium
($M_{b\bar{b}}$) and the scalar ($J^{P}=0^{+}$) D meson ($M_{n\bar{c}}$) to
fit the three parameters. The masses of the other ground state scalar mesons
and all the ground state axial-vector mesons are outputs of this model.

Before discussing the results, there is one more ingredient to be stated. The
quark propagators in Eq. (2) are functions of the complex momenta
$q^{2}_{\pm}$ which lie in parabolic regions. As the dressed quark propagators
possess time-like complex conjugate mass poles [22], there is an upper bound
of the bound state mass obtainable directly, i.e. $P^{2}>-M_{max}^{2}$, where
$M_{max}^{2}$ defines the contour border of the parabolic region. See the
appendix of Ref. [11] for this problem.

We adopt the extrapolation scheme to overcome this difficulty. This method has
been used to study the radial excited $B_{c}$ mesons [17] and a brief
description is given below. Eq. (2) is solved as a $P^{2}$-dependent
eigenvalue problem,

$\lambda^{fg}(P^{2})\big{[}\Gamma^{fg}(k;P)\big{]}^{\alpha}_{\beta}=\int^{\Lambda}_{dq}\big{[}K(k,q;P)\big{]}^{\alpha\delta}_{\sigma\beta}\big{[}\chi^{fg}(q;P)\big{]}^{\sigma}_{\delta}.$
(12)

The meson mass is determined by $\lambda^{fg}(P^{2}=-M^{2}_{fg})=1$. If the
meson mass is beyond the contour border, i.e. $M_{fg}>M_{max}$, then the
following form is used to fit $\lambda^{fg}(P^{2})$,

$\frac{1}{\lambda^{fg}(P^{2})}=\frac{1+\sum^{N_{o}}_{n=1}\,a_{n}(P^{2}+s_{0})^{n}}{1+\sum^{N_{o}}_{n=1}\,b_{n}(P^{2}+s_{0})^{n}}\,.$
(13)

In Eq. (13) $N_{o}$ is the order of the series, $s_{0}$, $a_{n}$ and $b_{n}$
are the parameters to be determined by the least square method.

(color online) The $P^{2}$ dependence of $1/\lambda^{n\bar{n}}$ for
$J^{PC}=0^{++}$ state using the middle group of parameters in Table. 5. The
vertical dot-dashed orange line is the contour border on the right of which
the direct calculation can be applied. The black points are the calculated
values. The dashed black line, dot-dashed green line and dotted blue line are
the extrapolation functions, $1/\lambda^{n\bar{n}}(P^{2})$, with
$N_{o}=1,\,2,\,3$ in Eq. (13) respectively. The red stars present the location
of the meson masses.

The extrapolation results are illustrated by the lightest scalar meson, shown
in figure 2. As the mass pole is not too far away from the contour border, the
series in Eq. (13) converge very well. So the error from extrapolation is
controllable. In practice $N_{o}=1,2,3$ is considered, and the maximal
difference between them is estimated as the error from extrapolation.

##  3 Results and Discussion

In Ref. [12], three group of parameters are used for
$\\{\bar{m}_{f}^{\zeta=2\textmd{GeV}},\omega_{f},D_{f}\\}$, $f=\\{n,s,c,b\\}$,
$n=u\,\textmd{or}\,d$. Herein I also consider these three cases. In each case,
$\delta$, $\alpha$ and $\beta$ in Eq.(11) are fitted by three observables: the
masses of scalar ($J^{PC}=0^{++}$) charmonium ($M_{c\bar{c}}$), bottomium
($M_{b\bar{b}}$) and the scalar ($J^{P}=0^{+}$) D meson ($M_{n\bar{c}}$). The
parameters are listed in table 5 in the appendix. The masses of the other
ground state scalar mesons and all the ground state axial-vector mesons are
the outputs, listed in table 3.1, 3.2 and 3.3. In these tables, two errors of
the DSE results are given. The first error is from the extrapolation, the
second one is from varying the parameters. Next, let’s discuss the results of
heavy mesons, heavy-light mesons and light mesons respectively.

###  3.1 The heavy mesons

The results of the heavy mesons are listed in table 3.1. While the scalar
($J^{PC}=0^{++}$) charmonium and bottomium masses are used to fit the
parameters, the axial vector ($J^{PC}=1^{++}\textmd{ or }1^{+-}$) counterparts
are the outputs. Comparing the DSE results with the experiment values, the
largest deviation comes from the $J^{PC}=1^{+-}$ charmonium, being about
$20\textmd{ MeV}$. The other three mass deviations are within $10\textmd{
MeV}$. The scalar ($J^{P}=0^{+}$) and lowest two axial-vector
($J^{P}=1^{+}\textmd{ or }1^{-}$) $B_{c}$ meson masses are also output
quantities. Note that the DSE results from the modified model
(Eq.(5),(6),(7),(10),(11),(9)) are almost the same as those from the original
model (Eq.(5),(6),(7),(8),(9)) after subtracting the errors carefully [16]. As
there is no experiment data for these $B_{c}$ meson masses yet, the DSE
results are compared with the lattice QCD [23] and quark model predictions
[24]. The deviations are all within $20\textmd{ MeV}$. The systematic error of
the modified model for the heavy meson mass is about $20\textmd{ MeV}$.

Masses of ground state heavy scalar and axial-vector mesons (in GeV).
$M^{\textmd{DSE}}_{f\bar{g}}$ is the Dyson-Schwinger equation results of
mesons with quark flavor content $f\bar{g}$. The first error is from the
extrapolation, and the second is from varying the parameters.
$M^{\textmd{expt.}}_{f\bar{g}}$ is the corresponding experiment results [25].
The underlined values in this table and TABLE 3.2 are used to fit the
parameters in Eq. (11). $\bar{M}^{\textmd{RL}}_{c\bar{b}}$ are the results by
estimating the errors from the original RL approximation calculation carefully
[16]. $M^{\textmd{LQCD}}_{c\bar{b}}$ are the lattice QCD results [23].
$M^{\textmd{QM}}_{c\bar{b}}$ are the quark model results [24].
$J^{\textmd{PC}}$ $0^{++}$ $1^{+-}$ $1^{++}$ $M^{\textmd{DSE}}_{c\bar{c}}$
3.415(0)(0) 3.502(4)(10) 3.507(4)(5) $M^{\textmd{expt.}}_{c\bar{c}}$ 3.415
3.525 3.511 $M^{\textmd{DSE}}_{b\bar{b}}$ 9.859(0)(0) 9.905(1)(3) 9.897(1)(2)
$M^{\textmd{expt.}}_{b\bar{b}}$ 9.859 9.899 9.893 $J^{\textmd{P}}$ $0^{+}$
$1^{+}_{1}$ $1^{+}_{2}$ $M^{\textmd{DSE}}_{c\bar{b}}$ 6.700(1)(4) 6.743(2)(3)
6.780(3)(3) $\bar{M}^{\textmd{RL}}_{c\bar{b}}$ 6703 6745 6781
$M^{\textmd{LQCD}}_{c\bar{b}}$ 6712 6736 – $M^{\textmd{QM}}_{c\bar{b}}$ 6714
6757 6776

###  3.2 The heavy-light mesons

The results of the heavy-light meson masses are listed in table 3.2. Except
the scalar ($J^{P}=0^{+}$) D meson mass ($M_{n\bar{c}}$) being used to fit the
parameters, all the other masses are output quantities. Comparing the DSE
results with the well established experiment value, i.e. the masses of the
lowest axial-vector D meson ($M^{\textmd{expt.}}_{n\bar{c}}$), B meson
($M^{\textmd{expt.}}_{n\bar{b}}$) and $B_{s}$ meson
($M^{\textmd{expt.}}_{s\bar{b}}$), the systematic error of the modified model
for the heavy-light mesons is conservatively estimated to be $40\textmd{
MeV}$. Now look at the fourth row of table 3.2. The DSE results of the $D_{s}$
meson masses are unnaturally larger than the experiment value, from around
$60\textmd{ MeV}$ to $120\textmd{ MeV}$. This supports that the ground state
scalar and axial vector $D_{s}$ mesons may be not pure $q\bar{q}^{\prime}$
states [26, 27, 28]. The other experiment values are still missing at present
and my calculations make a prediction for them.

Masses of ground state heavy-light scalar and axial-vector mesons (in GeV).
The notations are the same as TABLE 3.1. For isospin multiplets, the average
value of the experiment data is cited. The experiment data with a ${\dagger}$
on the shoulder indicate that they may not be pure $q\bar{q}^{\prime}$ states
[26, 27, 28].  $J^{\textmd{P}}$ $0^{+}$ $1^{+}_{1}$ $1^{+}_{2}$
$M^{\textmd{DSE}}_{n\bar{c}}$ 2.333(6)(0) 2.443(8)(9) 2.569(10)(21)
$M^{\textmd{expt.}}_{n\bar{c}}$ 2.333 2.421 – $M^{\textmd{DSE}}_{n\bar{b}}$
5.636(13)(5) 5.693(10)(2) 5.836(10)(6) $M^{\textmd{expt.}}_{n\bar{b}}$ – 5.726
– $M^{\textmd{DSE}}_{s\bar{c}}$ 2.434(0)(10) 2.517(6)(17) 2.596(11)(18)
$M^{\textmd{expt.}}_{s\bar{c}}$ $2.318^{\dagger}$ $2.460^{\dagger}$
$2.535^{\dagger}$ $M^{\textmd{DSE}}_{s\bar{b}}$ 5.755(7)(8) 5.792(9)(12)
5.878(12)(13) $M^{\textmd{expt.}}_{s\bar{b}}$ – 5.829 –

###  3.3 The light mesons

Masses of ground state light scalar and axial-vector $q\bar{q}$-states (in
GeV). The notations are the same as TABLE 3.1. For isospin multiplets, the
average value of the experiment data is cited. The experiment data with a
${\dagger}$ on the shoulder indicate that they may not a single
$q\bar{q}^{\prime}$ pole. The experiment value for the light scalar and
isosinglet axial-vector meson masses are leaved empty, see the text.
$J^{\textmd{PC}}$ $0^{++}$ $1^{+-}$ $1^{++}$ $M^{\textmd{DSE}}_{n\bar{n}}$
1.098(6)(8) 1.220(10)(40) 1.252(13)(31) $M^{\textmd{expt.}}_{n\bar{n}}$ –
1.230 1.230 $M^{\textmd{DSE}}_{s\bar{s}}$ 1.442(1)(57) 1.533(5)(76)
1.558(6)(62) $M^{\textmd{expt.}}_{s\bar{s}}$ – – – $J^{\textmd{P}}$ $0^{+}$
$1^{+}_{1}$ $1^{+}_{2}$ $M^{\textmd{DSE}}_{n\bar{s}}$ 1.331(4)(51)
1.433(11)(68) 1.502(11)(79) $M^{\textmd{expt.}}_{n\bar{s}}$ –
$1.253^{\dagger}$ $1.403^{\dagger}$

The light scalar and axial-vector mesons are much more illusive. The
$a_{0}(980)$ and $f_{0}(980)$ may be $K\bar{K}$ molecule states due to the
color hyperfine interaction, and the actual $q\bar{q}^{\prime}$ scalar mesons
may become too light to be detected due to the same reason [29]. For the
isosinglet light axial-vector mesons, flavor mixing should be taken into
account [30]. It’s beyond the scope of this paper to explore those questions.
However, we can still draw some conclusions from present DSE results. All the
DSE results for the light mesons are the output quantities, and are listed in
table 3.3. The masses of the lightest two axial vector mesons, i.e.
$M^{\textmd{DSE}}_{n\bar{n}}$ with $J^{\textmd{PC}}=1^{++}$ and $1^{+-}$ are
consistent with the experiment values, the masses of $a_{1}$ and $b_{1}$. So
Eq. (10) and (11) is still effective for the light mesons.

Note that the parameter dependency of the light meson masses is much larger
than those of heavy and heavy-light mesons. The overall systematic error of
the light meson masses is also expected larger. However, the mass splittings
of the two lowest axial-vector mesons of the same content should be reliable.
Then look at the K mesons. From the DSE results,
$M^{\textmd{DSE}}_{n\bar{s}}(1^{+}_{2})-M^{\textmd{DSE}}_{n\bar{s}}(1^{+}_{1})\approx
70\textmd{ MeV}$. The experiment value
$M^{\textmd{expt.}}_{n\bar{s}}(1^{+}_{2})-M^{\textmd{expt.}}_{n\bar{s}}(1^{+}_{1})\approx
150\textmd{ MeV}$ seems unnaturally large. Some studies indicate that the
experiment measured $K(1)(1270)$ may be actually two poles [31, 32, 33]. This
is a possible reason of the deviation herein.

##  4 Summary

In this work, I studied the spectrum of the scalar mesons and axial-vector
mesons via Dyson-Schwinger equation and Bethe-Salpeter equation approach in
the rainbow-ladder approximation. The interaction model is modified, i.e.
adding a repulsive term to the original one that are effective for the
pseudoscalar and vector mesons. The added term is flavor dependent and has 3
parameters. With the three parameters fitted by the masses of scalar
($J^{PC}=0^{++}$) charmonium ($M_{c\bar{c}}$), bottomium ($M_{b\bar{b}}$) and
the scalar ($J^{P}=0^{+}$) D meson ($M_{n\bar{c}}$), the other outputted
masses are all consistent with those experiment value that are well
established. So this modified model is effective for all the heavy, heavy-
light and light scalar and axial-vector mesons. This is the first time a
systematic study of the spectrum of the scalar mesons and axial-vector mesons
in the DSE approach is achieved. The systematic errors of the DSE results are
estimated carefully. Comparing the DSE results with the experiment value, we
find that the scalar and axial-vector $D_{s}$ mesons may not be pure
$q\bar{q}^{\prime}$ states and the experiment measured $K(1)(1270)$ may
actually be two poles.

###### Acknowledgements.

Acknowledgments: thanks Professor Xianhui Zhong for usefull discussion.

##  5 Appendix A

The three groups of parameters correspond to
$\omega_{u}=0.45,0.50,0.55\textmd{ GeV}$ are listed in Tab.5. The quark mass
$\bar{m}_{f}^{\zeta}$ is defined by

$\displaystyle\bar{m}_{f}^{\zeta}$ $\displaystyle=$
$\displaystyle\hat{m}_{f}/\left[\frac{1}{2}\textmd{Ln}\frac{\zeta^{2}}{\Lambda^{2}_{\textmd{QCD}}}\right]^{\gamma_{m}},$
(14) $\displaystyle\hat{m}_{f}$ $\displaystyle=$
$\displaystyle\lim_{p^{2}\to\infty}\left[\frac{1}{2}\textmd{Ln}\frac{p^{2}}{\Lambda^{2}_{\textmd{QCD}}}\right]^{\gamma_{m}}M_{f}(p^{2}),$
(15)

with $\zeta$ the renormalisation scale, $\hat{m}_{f}$ the renormalisation-
group invariant current-quark mass and $M_{f}(p^{2})$ the quark mass function.

Fitted parameters correspond to $\omega_{u/d}=0.45,0.50,0.55\textmd{ GeV}$.
$\bar{m}_{f}^{\zeta=2\textmd{GeV}}$, $\omega_{f}$ and $D_{f}$ are measured in
GeV. $\alpha$ and $\beta$ are of unit 1. $\delta$ is measured in
$\textmd{GeV}^{\alpha+\beta-2}$.  Parameter-1 Parameter-2 Parameter-3 flavor
$\bar{m}_{f}^{\zeta=2\textmd{GeV}}$ $w_{f}$ $D_{f}^{2}$ $w_{f}$ $D_{f}^{2}$
$w_{f}$ $D_{f}^{2}$ $u/d$ 0.0049 0.450 1.133 0.500 1.060 0.550 1.014 $s$ 0.112
0.490 1.090 0.530 1.040 0.570 0.998 $c$ 1.17 0.690 0.645 0.730 0.599 0.760
0.570 $b$ 4.97 0.722 0.258 0.766 0.241 0.792 0.231 $\delta$ 0.02385 0.02210
0.01959 $\alpha$ 0.2691 0.2530 0.2902 $\beta$ 7.893 8.632 9.826

## References

  * Dyson [1949] F. J. Dyson, Phys. Rev. 75, 1736 (1949).
  * Schwinger [1951] J. S. Schwinger, Proc. Nat. Acad. Sci. 37, 452 (1951).
  * Eichmann et al. [2016] G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, Prog. Part. Nucl. Phys. 91, 1 (2016), 1606.09602.
  * Roberts [2017] C. D. Roberts, Few Body Syst. 58, 5 (2017), 1606.03909.
  * Salpeter and Bethe [1951] E. E. Salpeter and H. A. Bethe, Physical Review 84, 1232 (1951).
  * Maris and Roberts [1997] P. Maris and C. D. Roberts, Phys. Rev. C56, 3369 (1997), nucl-th/9708029.
  * Maris et al. [1998] P. Maris, C. D. Roberts, and P. C. Tandy, Phys. Lett. B420, 267 (1998), nucl-th/9707003.
  * Maris and Tandy [1999] P. Maris and P. C. Tandy, Phys. Rev. C60, 055214 (1999), nucl-th/9905056.
  * Qin et al. [2011] S.-x. Qin, L. Chang, Y.-x. Liu, C. D. Roberts, and D. J. Wilson, Phys. Rev. C84, 042202 (2011), 1108.0603.
  * Fischer et al. [2015] C. S. Fischer, S. Kubrak, and R. Williams, Eur. Phys. J. A51, 10 (2015), 1409.5076.
  * Hilger et al. [2015] T. Hilger, C. Popovici, M. Gomez-Rocha, and A. Krassnigg, Phys. Rev. D 91, 034013 (2015), 1409.3205.
  * Chen and Chang [2019] M. Chen and L. Chang, Chinese Physics C43, 114103 (2019), 1903.07808.
  * Qin and Roberts [2020] S.-x. Qin and C. D. Roberts, Chin. Phys. Lett. 37, 121201 (2020), 2008.07629.
  * Krassnigg and Blank [2011] A. Krassnigg and M. Blank, Phys.Rev. D83, 096006 (2011), 1011.6650.
  * Chen et al. [2017] J. Chen, M. Ding, L. Chang, and Y.-x. Liu, Phys. Rev. D95, 016010 (2017), 1611.05960.
  * Chen et al. [2020] M. Chen, L. Chang, and Y.-x. Liu, Phys. Rev. D 101, 056002 (2020), 2001.00161.
  * Chang et al. [2020] L. Chang, M. Chen, and Y.-x. Liu, Phys. Rev. D 102, 074010 (2020), 1904.00399.
  * Bender et al. [1996] A. Bender, C. D. Roberts, and L. Smekal, Physics Letters B 380, 7 (1996).
  * Chang and Roberts [2009] L. Chang and C. D. Roberts, Phys. Rev. Lett. 103, 081601 (2009), 0903.5461.
  * Qin et al. [2014] S.-X. Qin, C. D. Roberts, and S. M. Schmidt, Phys. Lett. B 733, 202 (2014), 1402.1176.
  * Binosi et al. [2016] D. Binosi, L. Chang, J. Papavassiliou, S.-X. Qin, and C. D. Roberts, Phys. Rev. D 93, 096010 (2016), 1601.05441.
  * Bhagwat et al. [2003] M. Bhagwat, M. A. Pichowsky, and P. C. Tandy, Phys. Rev. D67, 054019 (2003), hep-ph/0212276.
  * Mathur et al. [2018] N. Mathur, M. Padmanath, and S. Mondal, Phys. Rev. Lett. 121, 202002 (2018), 1806.04151.
  * Li et al. [2019] Q. Li, M.-S. Liu, L.-S. Lu, Q.-F. Lü, L.-C. Gui, and X.-H. Zhong, Phys. Rev. D99, 096020 (2019), 1903.11927.
  * Zyla et al. [2020] P. Zyla et al. (Particle Data Group), Progress of Theoretical and Experimental Physics 8, 083C01 (2020).
  * Chen and Li [2004] Y.-Q. Chen and X.-Q. Li, Phys. Rev. Lett. 93, 232001 (2004), hep-ph/0407062.
  * Faessler et al. [2007a] A. Faessler, T. Gutsche, V. E. Lyubovitskij, and Y.-L. Ma, Phys. Rev. D 76, 014005 (2007a), 0705.0254.
  * Faessler et al. [2007b] A. Faessler, T. Gutsche, V. E. Lyubovitskij, and Y.-L. Ma, Phys. Rev. D 76, 114008 (2007b), 0709.3946.
  * Godfrey and Isgur [1985] S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985).
  * Chen et al. [2015] K. Chen, C.-Q. Pang, X. Liu, and T. Matsuki, Phys. Rev. D 91, 074025 (2015), 1501.07766.
  * Geng et al. [2007] L. S. Geng, E. Oset, L. Roca, and J. A. Oller, Phys. Rev. D75, 014017 (2007), hep-ph/0610217.
  * Meißner [2020] U.-G. Meißner, Symmetry 12, 981 (2020), 2005.06909.
  * Wang et al. [2020] G.-Y. Wang, L. Roca, E. Wang, W.-H. Liang, and E. Oset, Eur. Phys. J. C80, 388 (2020), 2002.07610.