# Rheological response of a glass-forming liquid having large bidispersity

Vinay Vaibhav  The Institute of Mathematical Sciences, CIT Campus, Taramani,
Chennai 600113, India  Homi Bhabha National Institute, Anushaktinagar, Mumbai
400094, India  Jürgen Horbach  Institut für Theoretische Physik II, Heinrich-
Heine-Universität Düsseldorf, Universitätsstraße 1, 40225 Düsseldorf, Germany
Pinaki Chaudhuri  The Institute of Mathematical Sciences, CIT Campus,
Taramani, Chennai 600113, India  Homi Bhabha National Institute,
Anushaktinagar, Mumbai 400094, India

###### Abstract

Using extensive numerical simulations, we investigate the flow behaviour of a
model glass-forming binary mixture whose constituent particles have a large
size ratio. The rheological response to applied shear is studied in the regime
where the larger species are spatially predominant. We demonstrate that the
macroscopic rigidity that emerges with increasing density occurs in the regime
where the larger species undergo a glass transition while the smaller species
continue to be highly diffusive. We analyse the interplay between the
timescale imposed by the shear and the quiescent relaxation dynamics of the
two species to provide a microscopic insight into the observed rheological
response. Finally, by tuning the composition of the mixture, we illustrate
that the systematic insertion of the smaller particles affects the rheology by
lowering of viscosity of the system.

##  I I. Introduction

The rheological properties of soft glassy materials (colloids, gels, foams,
emulsions, granular matter etc.) larson ; yodhreview ; ludovicRMP2017 ;
joshipetekidis ; coussot2014yield ; nicolas2018deformation ;
rodney2011modeling  are very important for a wide range of applications in our
daily lives as well as in industries, and also lead to various natural
phenomena. A characteristic feature of soft glassy materials is the existence
of a yield stress, i.e. a threshold stress value which needs to be overcome to
make the material flow. Once flowing, soft glasses typically behave as shear-
thinning materials. Understanding the processes that lead to yielding and flow
from a microscopic perspective is fundamental to our physical knowledge of the
properties of soft glasses and could be utilized for the development of soft
glasses that are functionalized towards specific applications.

Many soft glass systems consist of particles having a wide range of sizes.
This is an important issue with respect to their mechanical properties.
However, in experiments pusey2002 ; besseling2007 ; schall2007 ; chikkadi12
and simulations yamamoto98 ; ludovic2002 ; varnik2004 ; gauravPRE15 ;
gauravJRheo15 ; anshul17 ; falkLanger98 ; sollich2012 ; sollich2013  on the
rheology of soft glasses, the focus has been on model systems with a moderate
size dispersity. Here, the size dispersity has essentially been introduced to
prevent the system from crystallization. However, recent studies poon2015 ;
laurati17  have revealed that the degree of polydispersity has a strong impact
on the glassy dynamics. In systems with a large polydispersity, there is a
time-scale separation between the structural relaxation of large and small
particles that leads to strong dynamical heterogeneities. How such a
dispersity in particle size affects the rheological response is not well
understood and needs to be further explored.

To address this issue, it is useful to consider disparate-sized binary
mixtures. In the quiescent state, such systems are known to exhibit a complex
phase behaviour dhont95  and a wealth of fascinating phenomena with respect to
their glassy dynamics colmenero06 ; colmenero06pre ; juergen2009 ; laurati15 ;
laurati19 . In the regime where the large particles predominantly occupy
volume, there is a distinct difference in the dynamics of large and small
particles as the density or packing fraction is increased. The larger species
are observed to undergo a glass transition first, while the smaller species
remain highly diffusive. Subsequently, at much higher density, the smaller
species reach dynamical arrest, with their motion reminiscent to the
localization dynamics of tracer particles in an amorphous matrix of frozen-in
soft spheres skinner2013 ; schnyder2017 ; schnyder2018 .

Only recently, there has been systematic experimental investigations of the
transient and steady-state rheological response of binary colloidal mixtures
having large size disparity between the constituent particles sentjabrskaja13
; egelhaaf13 ; sentjabrskaja14 ; sentjabrskaja18 ; sentjabrskaja19 . These
investigations, typically done at a fixed volume fraction and using different
shear protocols, reveal that the shear response changes varying the
concentration of the smaller species. As long as the larger particles occupy
more volume one observes a softening of the system, while it hardens when the
volume occupied by the smaller species becomes predominant. This changing
rheological response has been attributed to the change in local packing
structures as the composition of the mixture is varied.

In this work, we report an extensive numerical study of the rheological
behaviour of a model glass-forming binary mixture having large bidispersity,
considering the regime where the larger particles predominantly occupy volume.
We consider a model system where the two constituent species undergo separate
dynamical arrests at vastly different densities juergen2009 . We probe how the
timescale introduced by the external shear interacts with the relaxation
dynamics of the constituent species, which is then manifested in the flow
behaviour. The main finding is that the measured rheological curves reveal the
onset of macroscopic rigidity in the density regime where the larger particles
undergo a mode coupling transition in the quiescent state. Importantly, when
this rigidity is observed at macro-scale, the smaller species continue to
express a fluid-like dynamics that is not affected by the external shear.
Thus, there is a macroscopic solid-like response even though internally there
still exists a fast dynamic species. Only at much larger densities, the
external shear is seen to affect the dynamics of the smaller particles.
However, these fast particles do influence the rheology near the regime where
rigidity sets in; we observe that the systematic insertion of the smaller
particles into the mixture leads to a lowering of the viscosity, i.e. to the
softening of the material, consistent with earlier experimental observations.

The paper is organised as follows. After the introductory discussion in
Section I, the model system and the numerical methods are discussed in Section
II. In Section III, we discuss our findings, with subsections on (i)
recapitulation of equilibrium results that have been reported previously, (ii)
elaboration of macroscopic and microscopic analysis of the response to applied
shear, and (iii) illustration of how the smaller species influence the
rheological response. Finally, in Section IV, we have a concluding discussion.

##  II II. Model and computer simulation

The model glass-former that we consider for our study is a $50:50$ binary
mixture of repulsive particles whose quiescent dynamics has been well studied
juergen2009 . The bigger particles (labelled A) have a slight polydispersity
in size to avoid crystallization; their diameters $d_{\rm A}$ are sampled from
a uniform distribution, i.e., $d_{\rm A}\in[0.85,1.15]$. The diameter of the
smaller particles (labelled B) are fixed to $d_{\rm B}=0.35$. The average size
ratio of A and B particles is $\langle d_{\rm A}\rangle/d_{\rm B}\approx 2.85$
where $\langle d_{\rm A}\rangle/\approx 1$. Details about the model
interactions between these constituent particles, as well as the units for
length, energy, and time can be found in Refs. juergen2009 ; vinay .

We study the rheological behaviour of this model system using molecular
dynamics (MD) simulation via LAMMPS lammps . We consider a system consisting
of $N=2000$ particles in a three-dimensional box with periodic boundary
conditions at constant temperature $T=2/3$, which is maintained via a
dissipative particle dynamics (DPD) thermostat dpd . The time step used for
the simulations is $dt=0.00075$. Initial states at different densities are
sampled from our study on the quiescent behaviour of this system; see Ref.
vinay . Apart from the 50:50 composition, we also study a couple of other
compositions, viz. 100:0 and 75:25, in order to understand how the inclusion
of B species changes the system’s response to an external shear.

Figure 1: Snapshot of the sheared binary mixture. The color panel represents
the diameter of the polydisperse larger species (A). Smaller particles (B) are
shown in black. The shear is applied in the $x$-direction along the $xy$
plane.

To investigate the rheological response of the system, we impose shear along
the $xy$ plane in the $x$-direction (see Fig. 1) using different shear rates
($\dot{\gamma}$) ranging between $1.5{\times}10^{-6}$ and $10^{-3}$. Lees-
Edwards boundary conditions are utilised during the shearing. The rheological
behaviour is studied for a range of densities varying from very small, viz.
$\rho_{\rm A}=1.05$ to very large $\rho_{\rm A}=1.75$, where $\rho_{\rm A}$ is
the partial density of A species. Note that all reported densities will be in
terms of $\rho_{A}$, since we will be eventually probing how changing
$\rho_{\rm B}$ influences the rheology at fixed $\rho_{\rm A}$, as discussed
above. A snapshot of the sheared system is shown in Fig. 1. As is evident, the
larger species indeed populate most of the volume and the smaller species dot
the intervening space. In this study, we inquire how the microscopic dynamics
of these two species are manifested in the macro-rheology during the shear
process.

##  III III. Results

###  III.1 Quiescent dynamics

Figure 2: Quiescent dynamics. Mean-squared displacement of bigger (upper
panel) and smaller (lower panel) species of the unperturbed system i.e., in
the absence of shear for different partial densities of bigger species.
Densities are labelled in terms of partial density of A species $\rho_{\rm
A}$.

As an introductory step, we summarize the dynamical behaviour observed by
scanning across densities and monitoring the time evolution of the mean
squared displacement (MSD) of the two species; see Fig. 2. As reported earlier
in Ref. juergen2009 , the dynamics of the A species falls out of equilibrium
around the critical density of mode coupling theory, $\rho_{\rm A}^{\rm
MCT}=1.115$). In the same density regime, the B species remain diffusive and
undergo a dynamical arrest at a much higher density juergen2009 . Or in other
words, we observe two separate onsets of glassiness in the sub-populations,
due to the large size ratio of the two species. Thus, at any density, the
dynamics of the B species are faster than that of the A species, i.e., the
relaxation timescales are relatively much slower. This large separation of
timescales plays an important role in the rheological response, when a shear
is introduced by driving the system at different shear-rates. The shear
imposes a new timescale into the system and its interplay with the intrinsic
timescale leads to the diversity in rheological response.

We also note here that for the A particles the MSD is measured in their
center-of-mass reference frame. With this, we avoid finite-size effects in the
high density regime, i.e. for $\rho>\rho_{\rm A}^{\rm MCT}$, occurring due to
a relatively large center-of-mass diffusion coefficient of the A species for
the considered relatively small system sizes. As a consequence, the A
particles move collectively with their center of mass, despite their dynamical
arrest with respect to particle rearrangements. For further discussion on
this, see Ref. vinay .

Figure 3: Rheological response. Variation of steady-state shear stress
$\sigma_{xy}$ (top) and corresponding viscosity $\eta$ (bottom) with applied
shear-rate $\dot{\gamma}$, at different partial densities ($\rho_{\rm A}$) of
bigger species as marked.  Figure 4: Herschel-Bulkley fits. Flow curves, i.e.,
shear stress vs. strain-rate at different partial densities ($\rho_{\rm A}$)
of bigger species as marked, shown with corresponding Herschel-Bulkley fit
($\sigma=\sigma_{0}+K\dot{\gamma}^{n}$) using dotted line. Fit estimates are
$\sigma_{0}=1.06,3.52,5.54$, $K=8.33,13.92,18.63$ and $n=0.36,0.25,0.09$ for
$\rho_{\rm A}=1.25,1.50,1.75$. (Inset) Variation of estimated value of yield
stress $\sigma_{0}$ with $\rho_{\rm A}$. Dotted vertical line marks $\rho_{\rm
A}^{\rm MCT}$.

###  III.2 Steady-state macroscopic shear response

Next, we discuss the rheological response of the system over a range of
densities and shear rates. A snapshot of the system during shear is shown in
Fig. 1, with the large size disparity of the two species distinctly visible.

Since the shear is imposed along the $xy$ plane, we measure the corresponding
shear stress $\sigma_{xy}$ using the following Irving-Kirkwood expression:
$\sigma_{xy}=\langle\frac{1}{V}\sum_{\alpha\beta}f^{x}_{\alpha\beta}r^{y}\rangle$,
where $f^{x}_{\alpha\beta}$ is the $x$-component of the force with respect to
a pair of particles and $r^{y}$ is the $y$-component of the distance vector
between two particles labelled $\alpha$ and $\beta$ which could belong to
either of the species A and B. $V$ is the total volume of the simulated
system. $\langle\cdot\rangle$ corresponds to averaging over states sampled in
steady state.

At each density and imposed shear-rate, we obtain the steady state by shearing
the system to large strains and then ensuring that the measured shear stress
$\sigma_{xy}$ fluctuates around a steady mean value. The data gathered for the
variation of shear-stress with imposed shear rate, i.e., the flow curve, for
the wide range of densities is shown in the top panel of Fig. 3. And, the
corresponding viscosity defined by
$\eta(\dot{\gamma})=\sigma_{xy}(\dot{\gamma})/\dot{\gamma}$ is shown in the
bottom panel of Fig. 3.

The main findings from the rheological data are the following. At small
densities ($\rho_{\rm A}<1.09$), a Newtonian regime, i.e.,
$\eta(\dot{\gamma})$ is constant, is observed at small shear-rates. At larger
shear rates, shear thinning is observed, i.e., the viscosity decreases with
increasing $\dot{\gamma}$. The Newtonian regime shrinks with increasing
density, i.e., pushed to smaller and smaller shear rates. These are distinct
characteristics of complex liquids. As the density is increased, an apparent
yield stress becomes visible at $\rho_{\rm A}\approx 1.125$, with
$\sigma_{xy}$ becoming flat with decreasing $\dot{\gamma}$. This implies that
there is an onset of macroscopic rigidity for the binary mixture. Around this
density regime, there is also a change in behaviour of $\eta(\dot{\gamma})$ –
within the window of small shear rates, the viscosity continues to be an
increasing function, i.e., there is no tendency to change curvature towards
plateauing out. With further increase of density, the apparent yield stress of
the system, viz. $\sigma_{xy}(\dot{\gamma}\rightarrow 0)$, increases steadily,
which is characteristic to amorphous solids. One can fit the variation of
stress with shear rate in this large density regime using the Herschel-Bulkley
(HB) function, $\sigma=\sigma_{0}+K\dot{\gamma}^{n}$, where $\sigma_{0}$ is
the estimated yield stress, $K$ is a constant and $n$ is the HB exponent. The
fits to the measured data for different $\rho_{\rm A}$ are shown in Fig. 4 and
the values of the fit parameters are listed in the caption. In the inset of
Fig. 4, we display the emergence and increase of the estimated yield stress
$\sigma_{0}$ with density, beyond $\rho_{\rm A}^{\rm MCT}$ .

Overall, such variation of rheological flow curves, viz. the emergence of an
yield stress with the variation in a control parameter (e.g. temperature or
packing fraction or density) are very similar to the ones obtained in other
soft-sphere glassy systems ludovic2002 ; sollich2012 ; sollich2013 ;
golkia2020 .

Figure 5: Microscopic dynamics under applied shear: Mean squared displacement
of bigger (a-c) and smaller (d-f) species in the presence (with solid lines)
and in the absence (with dotted lines) of shear, measured in the vorticity
direction. Corresponding self part of van Hove function, $G_{s}(y,t)$ of
bigger species (g-i) and smaller (j-l) species in the presence (with symbols)
and absence (dotted lines) of shear, at $t=8075.4$, also measured in the
vorticity direction. In both cases, first column represents system at partial
density $\rho_{\rm A}=1.0500$, middle at $\rho_{\rm A}=1.1125$ and last at
$\rho_{\rm A}=1.7500$.

###  III.3 Steady-state microscopic dynamics

In order to further analyse the rheological response, we study the microscopic
dynamics of the system in steady state. In the vorticity direction, i.e.,
normal to the shear plane, we compute the MSD for the A and B species,
separately, as well as the self part of the van Hove function of which MSD is
the second moment. The van Hove Function, $G_{s}(r,t)$ is measured at a late
time ($t=8075.4$) in each case. The corresponding aggregated plots showing the
variation of the dynamics with changing density are shown in Fig. 5, for MSD
and $G_{s}(r,t)$, for a range of imposed shear rates. All measurements are
done in the center-of-mass reference of respective species, to remove the
effects resulting from finite center-of-mass diffusion in the glassy regime of
A species.

At the smaller density ($\rho_{\rm A}=1.0500$), we do not observe any
variation in the dynamics with imposed shear rate, for both A and B species,
and the MSD curves as well as the $G_{s}(r,t)$ data are indistinguishable
between the equilibrium and non-equilibrium regimes. This corresponds to the
linear response regime, which is characterised by the Newtonian viscosity, see
Fig. 5 (a), (d) for the MSD data and the corresponding data for $G_{s}(r,t)$
in Fig. 5 (g), (j).

If we now consider the higher density, $\rho_{\rm A}=1.1125$, which is in the
vicinity of $\rho_{\rm A}^{\rm MCT}$, we observe that for the B species, the
indistinguishability between the equilibrium and non-equilibrium data for MSD
and $G_{s}(r,t)$ continue, i.e., the external shear is having no effect; see
Fig. 5 (e), (k). However, for the A species, the shear introduces deviation
from equilibrium behaviour, with an enhancement in long-time diffusive
dynamics with increasing shear rate, as is visible in the MSD data (Fig. 5
(b)) which is also evident in the distinct change in the shape of the
corresponding $G_{s}(r,t)$ (Fig. 5(h)). Thus, at this density, the imposed
shear rate corresponds to a timescale which is faster than the equilibrium
relaxation timescale of the A species, which is not the case for B species.
Note that, for this density, an apparent macroscopic yield stress is visible.
Thus, we can conclude that the apparent emergence of macroscopic rigidity and
subsequent yielding under shear is a consequence of the interplay between the
external drive and the structure formed by the A species, while the B species
remain oblivious to this process and thus seems inconsequential to the overall
rheological response, i.e., variation of viscosity with shear rate in this
density regime.

However, if we go to even larger densities, $\rho_{\rm A}=1.7500$, the imposed
shear affects the dynamics of both the A and B species; see Fig. 5 (c), (f),
(i), (l); for both species, enhanced diffusion is observed with increasing
shear rate. Thus, the imposed timescales via the shear are faster compared to
the equilibrium relaxation timescale for both species, even though there is
still a large separation in dynamical timescales in the quiescent state.

We can extract a diffusion coefficient from the long-time diffusive dynamics
observed in the MSD data, for the various imposed shear rates at different
densities. The shear-rate dependence of the diffusion constant, for the range
of densities explored, is shown in Fig. 6, for both A and B species.

At small densities, the diffusion coefficient for both A and B species has no
shear-rate dependence, which is just a reflection of the fact that the
dynamics remains unaffected by the shear, as discussed above. This behaviour
continues for B species up to very large densities where the shear starts to
have an effect. Also, in this density regime, the diffusion coefficient for A
species has a power-law dependence on shear rate and appears to vanish at
$\dot{\gamma}\rightarrow 0$, which is typical to the shear response of
amorphous systems. In the intermediate density regime, the diffusion
coefficient also decreases with decreasing shear rate, but deviates from a
power-law behaviour and seems to reach a constant at vanishing shear rates.
This would imply that there is a diffusive behaviour in the absence of shear
at these densities. This behaviour is observed around the density where the
onset of glassiness occurs for the A species, viz. in the vicinity of
$\rho_{\rm A}^{\rm MCT}$. We also note that for B species a similar behaviour
is observed at the largest density that we have explored, implying that in the
absence of shear, the smaller particles would be undergoing equilibrium
diffusive dynamics in that density regime.

To summarise, from the rheological response of the binary mixture, we infer
that beyond $\rho_{\rm MCT}$, an apparent yield stress emerges which is linked
to the glassiness undergone by the A species. On contrary, the B species
continue to be very mobile inside the matrix formed by the A species and their
dynamics does not couple to the shear, and thus their contribution to the
system’s shear response is minimal. Only at much higher densities, the
yielding process involves the collective interplay of both A and B species.

Figure 6: Diffusion coefficient in the presence of shear. Variation of
diffusion coefficient of bigger (top) and smaller (bottom) species with shear
rate, measured in the vorticity direction, at different partial densities of
bigger species as marked.

Figure 7: Effect of varying composition. (Top) Quiescent behavior. MSD
measured for the A species, for different composition of the mixture (circle:
100-0, square: 75-25, triangle: 50-50).(Inset) Variation of larger particle
diffusion coefficient with partial density $\rho_{\rm A}$ for different
composition of the system. (Bottom) Rheological response. Shear stress and
viscosity as a function of shear-rate for different composition of the
mixture. And, diffusion coefficient versus shear rate for the same
compositions.

###  III.4 Role of smaller species: varying the composition

We now address the question of what role do the B species play in the
rheological response. We address this questions by varying the composition of
the mixture, and even considering the situation where the B species are
completely absent.

We begin by first examining how the equilibrium behaviour of the system is
influenced by the presence or absence of the B species. In the top panel of
Fig. 7, we show the MSD data for the A species for varying A:B compositions,
viz. 100:0, 75:25 and 50:50 for three different $\rho_{\rm A}$. We note that
at small densities, there is no variation in the dynamics, but at relatively
larger densities ($\rho_{\rm A}=1.113$) there is distinct variation with the
long-time dynamics becoming faster with increased presence of B species. This
is also reflected in the measured diffusion coefficients as shown in the inset
of the top panel of Fig. 7. From this data, it would imply that $\rho_{\rm
A}^{\rm MCT}$ shifts to larger densities with the increasing inclusion of B
species. It is interesting to note that, the insertion of the smaller B
species has a very different effect than insertion of bigger particles in a
system of smaller particles. In the latter case, the bigger particles act as
pinning sites and slow down the dynamics.

Next, we study how the rheological response changes, in the presence or
absence of the small particles. In the bottom panel of Fig. 7, we do a
comparison of the rheological curves for the different compositions listed
above. The effect of enhancement in equilibrium dynamics due to the presence
of B species is reflected in the rheology data, especially in the intermediate
density regime (as the glassy regime of A species is approached), where the
viscosity of the system distinctly decreases at small shear rates with
increasing presence of B species. This effect is absent at small densities and
also at very large densities where the entire system becomes glassy. Thus, in
the intermediate density regime, the mobile B population does have a softening
effect on the overall shear response, an effect which can be exploited in
various applications.

Our findings are consistent with experiments involving colloidal glasses with
constituents having large size ratio, where a similar softening of the
material was observed with increasing insertion of small particles till a
50:50 composition was reached sentjabrskaja18 ; sentjabrskaja19 .

##  IV IV. Conclusions and perspective

In conclusion, we have studied the shear response of a glass forming model
binary mixture having a large size ratio between the constituent species. The
mixture having a composition of 50:50 corresponds to the case where the larger
species occupy most of the volume. For such a system, it had earlier been
demonstrated that as the overall density is increased, the larger particles
undergo dynamical arrest first while the small particles continue to be
diffusive within the matrix formed by the large ones. In our study, we probe
how the external drive applied via a fixed shear rate, of varying magnitude,
influences this microscopic dynamics and how that manifests in the observed
steady-state rheological response.

The main observation comes via the measured shear stress, in steady state, as
a function of the imposed shear rate. At small densities, the low shear-rate
regime corresponds to that of a Newtonian fluid. However, in the density
regime where the larger species within the mixture undergo a mode coupling
dynamical transition in the quiescent state, the shape of the flow curves
change at vanishing shear rates and an effective yield stress can clearly be
estimated. Thus, there is a macroscopic rigidity. If we now focus on the
microscopic dynamics, we observe that while the motion of the larger species
indeed get influenced by the external shear, the dynamics of the smaller
species remain completely unchanged vis-a-vis their motion in the absence of
shear, within the explored range of shear rates. Thus, the macroscopic
rigidity at vanishing shear rate is occurring while there is a highly mobile
population of the smaller species within the system.

With increasing density, the yield stress increases as expected and the
diffusion coefficient of the larger species exhibit a distinct power-law
behaviour as a function of applied shear rate which is a distinctive feature
of glassy systems. Only at large densities, we start to observe that the
external shear rate is able to influence the motion of the smaller species,
i.e., the timescale imposed by the applied shear rate is large enough to
compete with the relaxation timescale of the smaller species.

Finally, we probe the role of the smaller species in the rheological response
of the mixture, since the macroscopic rigidity emerges in the system even
though these particles remain highly mobile. We observe that if we start with
a system of just the large particles and then start inserting the smaller
species, i.e., systematically change the mixture’s composition, the
rheological response changes in the small shear-rate regime, viz. the
viscosity of the system decreases. Thus, the presence of the smaller particles
softens the mixture’s mechanical response, which is consistent with
experimental observation in similar colloidal mixtures. Note that similar
plasticizing effects have also been reported for polymeric systems
zaccarelli2005 ; kalathiprl2012  where the inclusion of very small colloidal
particles led to softening. Hence, this indicates to a very generic mechanism
at play in how such inclusions can influence rheology. We also note here that
this behaviour is in contrast to the situation where the inclusion of large
particles in a matrix of small particles can lead to hardening of the system,
i.e., the rheological behaviour of the mixture changes as the composition is
changed.

The next step should be to study the flow behaviour of glass forming liquids
where particles have large polydispersity in size laurati17 ; ranjini2017 , to
investigate how our findings on the connection between microscopic dynamics
and macroscopic mechanical response get translated to more complex mixtures.
Also, there is a need to study how thermal fluctuations influence rheology in
such mixtures, by comparing the response in the athermal limit where there
have been some initial studies on probing the jamming behaviour, the
plasticity and also the shear response in asymmetric binary mixtures ching10 ;
atsushi21 ; luding20 ; pednekar2018 ; eric2015 .

##  V Acknowledgments

We acknowledge the use of HPC facility at IMSc Chennai for our computational
work.

## References

  * (1) R. G. Larson, The structure and rheology of complex fluids (Oxford Univ. Press, New York, 1999).
  * (2) D. T. Chen, Q. Wen, P. A. Janmey, J. C. Crocker, and A. G. Yodh, Annu. Rev. Condens. Matter Phys. 1, 301 (2010).
  * (3) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Rev. Mod. Phys. 89, 035005 (2017).
  * (4) Y. Joshi and G. Petekidis, Rheol. Acta 57, 521 (2018).
  * (5) P. Coussot, J. Non-Newton. Fluid Mech. 211, 31 (2014).
  * (6) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Rev. Mod. Phys. 90, 045006 (2018).
  * (7) D. Rodney, A. Tanguy, and D. Vandembroucq, Model. Simul. Mater. Sci. Eng. 19, 083001 (2011).
  * (8) G. Petekidis, A. Moussaid, and P. N. Pusey, Phys. Rev. E 66, 051402 (2002).
  * (9) R. Besseling, E. R. Weeks, A. B. Schofield, and W. C. K. Poon, Phys. Rev. Lett. 99, 028301 (2007).
  * (10) P. Schall, D. A. Weitz, and F. Spaepen, Science 318, 1895 (2007).
  * (11) V. Chikkadi and P. Schall, Phys. Rev. E 85, 031402 (2012).
  * (12) R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 (1998).
  * (13) L. Berthier and J. L. Barrat, J. Chem. Phys. 116, 6228 (2002).
  * (14) F. Varnik, L. Bocquet, and J. L. Barrat, J. Chem. Phys. 120, 2788 (2004).
  * (15) G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 94, 042605 (2016).
  * (16) G. P. Shrivastav, P. Chaudhuri, and J. Horbach, J. Rheol. 60, 835 (2016).
  * (17) P. Leishangthem, A. D.S. Parmar, and S. Sastry, Nat. Commun. 8, 14653 (2017).
  * (18) M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
  * (19) A. Ikeda, L. Berthier, and P. Sollich, Phys. Rev. Lett. 109, 018301 (2012).
  * (20) A. Ikeda, L. Berthier, and P. Sollich, Soft Matter 9, 7669 (2013).
  * (21) E. Zaccarelli, S. M. Liddle, and W. C. K. Poon, Soft Matter 11, 324 (2015).
  * (22) D. Heckendorf, K. J. Mutch, S. U. Egelhaaf, and M. Laurati, Phys. Rev. Lett. 119, 048003 (2017).
  * (23) A. Imhof and J. K. G. Dhont, Phys. Rev. Lett. 75, 1662 (1995).
  * (24) A. J. Moreno and J. Colmenero, J. Chem. Phys. 125, 164507 (2006).
  * (25) A. J. Moreno and J. Colmenero, Phys. Rev. E 74, 021409 (2006).
  * (26) J. Horbach and T. Voigtmann, Phys. Rev. Lett. 109, 205901 (2009).
  * (27) J. Hendricks, R. Capellmann, A. B. Schofield, S. U. Egelhaaf, and M. Laurati, Phys. Rev. E 91, 032308 (2015).
  * (28) E. Lázaro-Lázaro, J. A. Perera-Burgos, P. Laermann, T. Sentjabrskaja, G. Pérez-Ángel, M. Laurati, S. U. Egelhaaf, M. Medina-Noyola, T. Voigtmann, R. Castañeda-Priego, and L. F. Elizondo-Aguilera, Phys. Rev. E 99, 042603 (2019).
  * (29) T. O. E. Skinner, S. K. Schnyder, D. G. A. L. Aarts, J. Horbach, and R. P. A. Dullens, Phys. Rev. Lett. 111, 128301 (2013).
  * (30) S. K. Schnyder, T. O. E. Skinner, A. L. Thorneywork, D. G. A. L. Aarts, J. Horbach, and R. P. A. Dullens, Phys. Rev. E 95, 032602 (2017).
  * (31) S. K. Schnyder and J. Horbach, Phys. Rev. Lett. 120, 078001 (2018).
  * (32) T. Sentjabrskaja, D. Guu, M. P. Lettinga, S. U. Egelhaaf, and M. Laurati, AIP Conf. Proc. 1518, 206 (2013).
  * (33) T. Sentjabrskaja, E. Babaliari, J. Hendricks, M. Laurati, G. Petekidis, and S. U. Egelhaaf, Soft Matter 9, 4524 (2013).
  * (34) T. Sentjabrskaja, M. Hermes, W. C. K. Poon, C. D. Estrada, R. Castaneda-Priego, S. U. Egelhaaf, and M. Laurati, Soft Matter 10, 6546 (2014).
  * (35) T. Sentjabrskaja, J. Hendricks, A. R. Jacob, G. Petekidis, S. U. Egelhaaf, and M. Laurati, J. Rheol. 62, 149 (2018).
  * (36) T. Sentjabrskaja, A. R. Jacob, S. U. Egelhaaf, G. Petekidis, T. Voigtmann, and M. Laurati, Soft Matter 15, 2232 (2019).
  * (37) V. Vaibhav, J. Horbach, and P. Chaudhuri, arxiv:2202.12189.
  * (38) S. Plimpton, J. Comp. Phys. 117, 1 (1995).
  * (39) T. Soddemann, B. Dünweg, and K. Kremer, Phys. Rev. E 68, 046702 (2003).
  * (40) M. Golkia, G. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 102, 023002 (2020).
  * (41) E. Zaccarelli, C. Mayer, A. Asteriadi, C. N. Likos, F. Sciortino, J. Roovers, H. Iatrou, N. Hadjichristidis, P. Tartaglia, H. Löwen, and D. Vlassopoulos, Phys. Rev. Lett. 95, 268301 (2005).
  * (42) J. T. Kalathi, G. S. Grest, and S. K. Kumar, Phys. Rev. Lett. 109, 198301 (2012).
  * (43) S. K. Behera, D. Saha, P. Gadige, and R. Bandyopadhyay, Phys. Rev. Mater. 1, 055603 (2017).
  * (44) S. Pednekar, J. Chun, and J. F. Morris, J. Rheol. 62, 513 (2018).
  * (45) Y. Hara, H. Mizuno, and A. Ikeda, Phys. Rev. Res. 3, 023091 (2021).
  * (46) J. C. Petit, N. Kumar, S. Luding, and M. Sperl, Phys. Rev. Lett. 125, 215501 (2020).
  * (47) N. Xu and E. S. C. Ching, Soft Matter 6, 2944 (2010).
  * (48) J. Clara-Rahola, T. A. Brzinski, D. Semwogerere, K. Feitosa, J. C. Crocker, J. Sato, V. Breedveld, and E. R. Weeks, Phys. Rev. E 91, 010301(R) (2015).


[a]L. Ya. Glozman

# Symmetries of temporal correlators and the nature of hot QCD

Y. Aoki  S. Hashimoto  C. Rohrhofer

###### Abstract

The temperature of the chiral restoration phase transition at $\sim 130$ MeV
as well as the temperature of the center symmetry ("deconfinement") phase
transition in a pure glue theory at $\sim 300$ MeV are two independent
temperatures and their interplay determines a structure of different regimes
of hot QCD. Given a chiral spin symmetry of the color charge and of the
chromoelectric interaction we can conclude from observed symmetries of spatial
and temporal correlators of $N_{F}=2$ QCD with domain wall Dirac operator at
physical quark masses that above the chiral symmetry restoration crossover
around $T_{pc}$ but below rougly $3T_{pc}$ there should exist an intermediate
regime (the stringy fluid) of hot QCD that is characterized by approximate
chiral spin symmetry and where degrees of freedom are chirally symmetric
quarks bound into color singlet objects by the chromoelectric field. Above
this intermediate regime the color charge and the chromoelectric field are
Debye screened and one observes a transition to QGP with magnetic confinement.

##  1 Introduction

Before the RHIC era there was a belief that at some critical temperature
$T_{c}$ there should exist a deconfinement phase transition from the hadron
gas to the quark gluon plasma (QGP) with free deconfined quarks and gluons.
There is no dynamical breaking of chiral symmetry in a system of free quarks
and consequently the spontaneously broken chiral symmetry below $T_{c}$ should
be restored above $T_{c}$. So the critical temperature $T_{c}$ should be a
common temperature of both deconfinement and of chiral restoration. At the
RHIC era this view changed since both on the lattice and experimentally no
phase transition was seen and instead there was a smooth but rather fast
crossover observed on the lattice first in Wuppertal [1]. The quark condensate
drops from its vacuum value at $T\sim 100$ MeV to practically zero at $T\sim
200$ MeV with a pseudocritical temperature around $T_{pc}\sim 155$ MeV. The
Polyakov loop showed an inflection point only slightly above this temperature
so after RHIC a new lore arised that in QCD there is a fast common
deconfinement - chiral restoration crossover from hadron gas to QGP around the
pseudocritical temperature $T_{pc}\sim 155$ MeV. This picture was also
confirmed by the Bielefeld lattice group.

However, the inflection point of the nonrenormalized Polyakov loop was used to
establish a "deconfinement temperature". If the renormalized Polyakov loop is
studied, which is physical, see Fig. 1 [2], then one clearly sees that there
is no hint of a deconfinement transition between 100 and 200 MeV, since the
deconfinement transition should be accompanied by Polyakov loop evolution from
0 to 1. The renormalized Polyakov loop shows its evolution from small values
to 1 in a broad temperature interval with inflection point around 300 MeV.
This temperature coincides with the temperature of the center symmetry
("deconfinement") phase transition in a pure glue theory. Hence the evolution
of the renormalized Polyakov loop suggests that in full QCD with light quarks
a dramatic rearrangement in gluodynamics happens still around 300 MeV, which
is however strongly smeared out by dynamical light quarks. This suggests that
properties of hot QCD should be influenced by two independent temperatures,
the chiral phase transition around $T_{c}\sim 130$ MeV [3] and by center
symmetry phase transition at $T_{d}\sim 300$ MeV. There should be no unique
pseudocritical temperature of combined chiral restoration and deconfinement
crossover.

Figure 1:  Renormalized Polyakov loop at physical quark masses. From Ref. [2]

From the Polyakov loop correlators one can extract an effective potential
between static sources (which is however dependent from an assumption whether
such a potential is real or complex). Such a potential below and above
$T_{pc}$ at physical and smaller than physical quark masses is shown in Fig.
2. [4]. We clearly see essentially the same potential at both temperatures. In
both cases there is a flattening of the linear part of the potential that is
related to the string breaking. At the same time there is not yet Debye
screening of the color charge. The Debye screening is by definition a
screening of the negative Coulomb potential. Deconfinement could be associated
with the Debye screening of the color charge [5]. Hence this potential tells
that both below and above $T_{pc}$ QCD is in the confining regime.

Figure 2:  Effective potential between the static sources at physical and
below than physical quark masses extracted from the correlators of
renormalized Polyakov loops. Left: $T=141$ MeV; right: $T=166$ MeV. From Ref.
[4]

Another hint that in QCD there is no a simple fast crossover from hadron gas
to QGP is suggested by experiments. The hadron gas is a gas. The QGP is also a
gas with a large $\eta/s$. Experimentalists tell us that at RHIC and LHC
temperatures between $\sim 150$ and $\sim 400$ MeV they observe a liquid with
a very small $\eta/s$. Hence what is seen in experiments cannot be a QGP.

One obviously needs an objective information about degrees of freedom above
$T_{pc}$ to decide what physics and regimes take place there. The objective
and model independent information can be delivered by symmetries that can be
extracted from correlators on the lattice.

##  2 Symmetries of the color charge and chromoelectric interaction in QCD

The chromoelectric interaction is defined via the color charge (Lorentz-
invariant)

$Q^{a}=\int d^{3}x\psi^{\dagger}(x)\frac{t^{a}}{2}\psi(x).$ (1)

The color charge has both $SU(N_{F})_{L}\times SU(N_{F})_{R}$ and $U(1)_{A}$
chiral symmetries which are also symmetries of the QCD Lagrangian. However it
has an additional symmetry: the chiral spin $SU(2)_{CS}$ symmetry [6,
7].111This symmetry was reconstructed from a large degeneracy of hadrons
observed on the lattice upon artifical truncation of the near zero modes of
the Dirac operator at zero temperature [8, 9].

The $SU(2)_{CS}$ chiral spin transformations are defined as follows:

$\psi\rightarrow\psi^{\prime}=\exp\left(i\frac{\varepsilon^{n}\Sigma^{n}}{2}\right)\psi\;$
(2)

with generators constructed from Euclidean gamma matrices
$\gamma_{k};~{}~{}~{}k=1,2,3,4$

$\Sigma=\\{\gamma_{k},-i\gamma_{5}\gamma_{k},\gamma_{5}\\}.$ (3)

The direct product of $SU(2)_{CS}$ and of $SU(N_{F})$ can be embedded into
$SU(2N_{F})$ group that contains as subgroups both $SU(N_{F})_{L}\times
SU(N_{F})_{R}$ and $U(1)_{A}$ symmetries. The $SU(2N_{F})$ is also a symmetry
of the color charge.

In Minkowski space in a given reference frame the quark-gluon interaction part
of the QCD Lagrangian can be split into temporal and spatial parts:

$\overline{\psi}\gamma^{\mu}D_{\mu}\psi=\overline{\psi}\gamma^{0}D_{0}\psi+\overline{\psi}\gamma^{i}D_{i}\psi.$
(4)

The covariant derivative $D_{\mu}$ includes interaction of the quark field
$\psi$ with the gluon field $A^{a}_{\mu}$,

$D_{\mu}\psi=(\partial_{\mu}-ig\frac{t^{a}A^{a}_{\mu}}{2})\psi.$ (5)

The temporal term contains interaction of the color-octet charge density

$\bar{\psi}(x)\gamma^{0}\frac{t^{a}}{2}\psi(x)=\psi(x)^{\dagger}\frac{t^{a}}{2}\psi(x)$
(6)

with the chromoelectric part of the gluonic field. It is a singlet under the
$SU(2)_{CS}$ and $SU(2N_{F})$ transformations. The spatial part consists of a
quark kinetic term and interaction with the chromomagnetic part of the gluonic
field. It breaks $SU(2)_{CS}$ and $SU(2N_{F})$. Hence interaction of quarks
with the electric and magnetic components of the gluonic field can be
distinguished by symmetry. Obviously one needs to fix a reference frame to
make a discussion of electric and magnetic components sensible.

Emergence of approximate chiral spin and $SU(2N_{F})$ symmetries in hadrons at
zero temperature upon truncation of the near-zero modes of the Dirac operator
observed in [8, 9] tells that while the magnetic interaction is located
predominantly in the near-zero modes of the Dirac operator, a confining
electric interaction is distributed among all modes of the Dirac operator. The
quark condensate of the vacuum is associated only with the near-zero modes.
Hence one can conclude from those results that confinement and chiral symmetry
breaking in QCD are not directly related phenomena. Given this observation it
was predicted that above $T_{pc}$, where the chiral symmetry is restored and
the near-zero modes of the Dirac operator are suppressed by temperature
effects, there should emerge the chiral spin and $SU(2N_{F})$ symmetries and
hence QCD should be still in the confining regime [10], for an overview see
[11]. In order to understand physics and degrees of freedom in QCD above
$T_{pc}$ one should study symmetry properties of correlators with respect to
chiral spin and $SU(2N_{F})$ transformations.

##  3 Emergence of the chiral spin and $SU(2N_{F})$ symmetries in spatial
correlators.

In this section we shortly overview the results on symmetries of spatial
correlators obtained with JLQCD domain wall $N_{F}=2$ configurations at
physical quark masses in Refs. [12, 13] that have been reported at two
previous lattice conferences.

In Fig. 3 we show a complete set of $J=0,1$ correlators at different
temperatures. One observes three distinct multiplets above the chiral symmetry
restoration crossover and below roughly $T\sim 500$ MeV. The multiplet $E_{1}$
represents degenerate scalar and pseudoscalar correlators. This signals
effective restoration of $U(1)_{A}$ symmetry that was previously observed in
Ref. [14]. The approximately degenerate multiplets $E_{2}$ and $E_{3}$ arise
from emerged approximate $SU(2)_{CS}$ and $SU(4)$ symmetries. These symmetries
suggest that degrees of freedom in QCD above the chiral crossover but below
$T\sim 500$ MeV are chirally symmetric quarks connected by the chromoelectric
field into color singlet objects. This regime of QCD is referred to as a
stringy fluid. At $T\sim 500$ MeV the chiral spin and $SU(4)$ symmetries
disappear and only chiral symmetries remain. It suggests that the color charge
and electric field get Debye screened and one observes a smooth transition to
QGP with quasiquarks and quasigluons being the effective degrees of freedom.

Figure 3: A complete set of spatial J=0,1 correlators at different
temperatures. From Ref. [13]

##  4 Temporal correlators above $T_{pc}$

If above the chiral crossover the QCD action is approximately chiral spin and
$SU(4)$ symmetric then it should also be observed in temporal correlators that
are connected to observable spectral density.

The t-direction correlators in the hadron rest frame are

$C_{\Gamma}(t)=\sum\limits_{x,y,z}<\mathcal{O}_{\Gamma}(x,y,z,t)\mathcal{O}_{\Gamma}(\mathbf{0},0)^{\dagger}>,$

where $\mathcal{O}_{\Gamma}=\bar{q}\Gamma\frac{\vec{\tau}}{2}q$ are all
possible $J=0$ and $J=1$ local operators. The transformation properties of the
$J=1$ operators with respect to $U(1)_{A}$, $SU(2)_{L}\times SU(2)_{R}$,
$SU(2)_{CS}$ and $SU(4)$ are shown in Fig. 4 [7]. Emergence of these
symmetries should be seen as a degeneracy of the corresponding correlators.

Figure 4: Transformations between interpolating vector operators, $i=1,2,3$.
The left columns indicate the chiral representation for each operator. Red and
blue arrows connect operators that transform into each other under
$SU(2)_{L}\times SU(2)_{R}$ and $U(1)_{A}$, respectively. Green arrows connect
operators that form triplets of $SU(2)_{CS},k=4$. The $f_{1}$ and $a_{1}$
operators are the $SU(2)_{CS},k=4$ – singlets. Purple arrows show the 15-plet
of $SU(4)$. The $f_{1}$ operator is a $SU(4)$-singlet. From Ref. [7]

These correlators with the domain wall Dirac operator at physical quark masses
with $N_{F}=2$ QCD (JLQCD configurations) normalized at $n_{t}=1$ calculated
on $48^{3}\times 12$ lattices at $220$ MeV [15] are shown in Fig. 5. Note that
one needs a sufficiently large size $N_{t}$ along the time direction to
observe a real evolution of the temporal correlation functions. E.g. with
$N_{t}=2$ by construction only chiral symmetries can be obtained since
everything is fixed by the free quark fields at source and sink. Consequently
at the moment we are limited only to the temperature $T=220$ MeV.

On the left side of Fig. 5 we show the correlators calculated with free,
noninteracting quarks. Dynamics of free quarks are governed by the Dirac
equation and only chiral symmetries exist. Indeed only degeneracies due to
$U(1)_{A}$ and $SU(2)_{L}\times SU(2)_{R}$ symmetries are seen in meson
correlators calculated for free quark gas. This pattern reflects correlators
at very high temperatures since due to the asymptotic freedom at a very high
$T$ the quark-gluon interactions can be neglected.

On the right side of Fig. 5 we show the correlators calculated in full QCD at
$T=220$ MeV. The pattern in full QCD is qualitatively different as compared to
the free quark gas. In full QCD one clearly sees multiplets of all symmetries
under discussion: $U(1)_{A}$, $SU(2)_{L}\times SU(2)_{R}$, $SU(2)_{CS}$ and
$SU(4)$. This confirms results obtained with the spatial correlators.

Figure 5: Temporal correlation functions for $48^{3}\times 12$ lattices. The
l.h.s. shows correlators calculated with free noninteracting quarks on the
same lattice, and features a symmetry pattern expected from chiral symmetry.
The r.h.s. presents full QCD data at a temperature of $T=220$ MeV, which shows
multiplets of all $U(1)_{A}$, $SU(2)_{L}\times SU(2)_{R}$, $SU(2)_{CS}$ and
$SU(4)$ groups. From Ref. [15]

##  5 Conclusions

Our results on symmetries of correlators of $N_{F}=2$ QCD at physical quarks
masses allow to distinguish three different regimes of QCD at nonzero
temperatures.

Figure 6: Three regimes of QCD. From Ref. [13]

.

At temperatures $0-T_{pc}$ we have a hadron gas consisting of confined mesons
with broken chiral symmetry. From the hadron gas there is a smooth chiral
symmetry restoration crossover and one arrives at the stringy fluid that is
characterized by chiral, $SU(2)_{CS}$ and $SU(4)$ symmetries. The electric
confinement still persists and degrees of freedom are chirally symmetric
quarks bound by the electric field into color singlet compounds. Around $T\sim
3T_{pc}$ the color charge gets Debye screened, the $SU(2)_{CS}$ and $SU(4)$
symmetries disappear (only chiral symmetries remain) and one observes a very
smooth transition to QGP with quasiquarks and quasigluons being effective
degrees of freedom.

## References

  * [1] Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643 (2006), 46-54 doi:10.1016/j.physletb.2006.10.021 [arXiv:hep-lat/0609068 [hep-lat]].
  * [2] P. Petreczky and H. P. Schadler, Phys. Rev. D 92 (2015) no.9, 094517 doi:10.1103/PhysRevD.92.094517 [arXiv:1509.07874 [hep-lat]].
  * [3] H. T. Ding et al. [HotQCD], Phys. Rev. Lett. 123 (2019) no.6, 062002 doi:10.1103/PhysRevLett.123.062002 [arXiv:1903.04801 [hep-lat]].
  * [4] D. A. Clarke, O. Kaczmarek, F. Karsch and A. Lahiri, PoS LATTICE2019 (2020), 194 doi:10.22323/1.363.0194 [arXiv:1911.07668 [hep-lat]].
  * [5] E. V. Shuryak, Sov. Phys. JETP 47 (1978), 212-219 IYF-77-34.
  * [6] L. Y. Glozman, Eur. Phys. J. A 51 (2015) no.3, 27 doi:10.1140/epja/i2015-15027-x [arXiv:1407.2798 [hep-ph]].
  * [7] L. Y. Glozman and M. Pak, Phys. Rev. D 92 (2015) no.1, 016001 doi:10.1103/PhysRevD.92.016001 [arXiv:1504.02323 [hep-lat]].
  * [8] M. Denissenya, L. Y. Glozman and C. B. Lang, Phys. Rev. D 89 (2014) no.7, 077502 doi:10.1103/PhysRevD.89.077502 [arXiv:1402.1887 [hep-lat]].
  * [9] M. Denissenya, L. Y. Glozman and C. B. Lang, Phys. Rev. D 91 (2015) no.3, 034505 doi:10.1103/PhysRevD.91.034505 [arXiv:1410.8751 [hep-lat]].
  * [10] L. Y. Glozman, Acta Phys. Polon. Supp. 10 (2017) 583 doi:10.5506/APhysPolBSupp.10.583 [arXiv:1610.00275 [hep-lat]].
  * [11] L. Y. Glozman, Int. J. Mod. Phys. A https://doi.org/10.1142/S0217751X20440315 [arXiv:1907.01820 [hep-ph]].
  * [12] C. Rohrhofer, Y. Aoki, G. Cossu, H. Fukaya, L. Y. Glozman, S. Hashimoto, C. B. Lang and S. Prelovsek, Phys. Rev. D 96 (2017) no.9, 094501 Erratum: [Phys. Rev. D 99 (2019) no.3, 039901] doi:10.1103/PhysRevD.96.094501, 10.1103/PhysRevD.99.039901 [arXiv:1707.01881 [hep-lat]].
  * [13] C. Rohrhofer et al., Phys. Rev. D 100 (2019) no.1, 014502 doi:10.1103/PhysRevD.100.014502 [arXiv:1902.03191 [hep-lat]].
  * [14] A. Tomiya, G. Cossu, S. Aoki, H. Fukaya, S. Hashimoto, T. Kaneko and J. Noaki, Phys. Rev. D 96 (2017) no.3, 034509 Addendum: [Phys. Rev. D 96 (2017) no.7, 079902] doi:10.1103/PhysRevD.96.034509, 10.1103/PhysRevD.96.079902 [arXiv:1612.01908 [hep-lat]].
  * [15] C. Rohrhofer, Y. Aoki, L. Y. Glozman and S. Hashimoto, Phys. Lett. B 802 (2020), 135245 doi:10.1016/j.physletb.2020.135245 [arXiv:1909.00927 [hep-lat]].