# Brittle yielding in supercooled liquids below the critical temperature of
mode coupling theory

Konstantin Lamp, Niklas Küchler, and Jürgen Horbach  Institut für Theoretische
Physik II: Weiche Materie, Heinrich-Heine-Universität Düsseldorf,
Universitätsstraße 1, 40225 Düsseldorf, Germany

###### Abstract

Molecular Dynamics (MD) computer simulations of a polydisperse soft-sphere
model under shear are presented. Starting point for these simulations are
deeply supercooled samples far below the critical temperature, $T_{c}$, of
mode coupling theory. These samples are fully equilibrated with the aid of the
swap Monte Carlo technique. For states below $T_{c}$, we identify a life time
$\tau_{\rm lt}$ that measures the time scale on which the system can be
considered as an amorphous solid. The temperature dependence of $\tau_{\rm
lt}$ can be well described by an Arrhenius law. The existence of transient
amorphous solid states below $T_{c}$ is associated with the possibility of
brittle yielding, as manifested by a sharp stress drop in the stress-strain
relation and shear banding. We show that brittle yielding requires on the one
hand low shear rates and on the other hand, the time scale corresponding to
the inverse shear rate has to be smaller or of the order of $\tau_{\rm lt}$.
Both conditions can be only met for large life time $\tau_{\rm lt}$, i.e. for
states far below $T_{c}$.

##  I Introduction

Glassforming liquids exhibit a dramatic slowing down of their dynamics with
decreasing temperature $T$. Important insight on the origin of this slowing
down has been given by the mode coupling theory (MCT) of the glass transition
goetze2009 . This theory predicts a divergence of the structural relaxation
time of the liquid when decreasing $T$ towards a critical temperature $T_{c}$.
At $T_{c}$, a transition from an ergodic liquid state to a non-ergodic
amorphous solid state occurs. The order parameter of this transition is
associated with the localization of each particle in the cage that is formed
by neighboring particles. Thus, in the framework of MCT, the glass transition
can be seen as a localization transition where, approaching the transition
from temperatures $T<T_{c}$, i.e. from below, the critical temperature $T_{c}$
marks the stability limit of the amorphous solid. At $T_{c}$, the length scale
$\xi$, that measures the localization of the particles in their cages, reaches
a critical value such that the amorphous solid state cannot be stable anymore
(note the analogy with the Lindemann criterion for crystalline solids
solyom2007 ).

In real glassforming systems, a transition, as predicted by MCT, is not
observed. However, using the predictions of MCT, a critical temperature
$T_{c}$ can be identified around which the dynamics of the supercooled liquid
gradually changes from a liquid-like to a solid-like dynamics cavagna2009 . As
a consequence, far below $T_{c}$, the supercooled liquid can be found in the
state of an amorphous solid, albeit this state has only a finite life time
$\tau_{\rm lt}$ and there is a diffusional time scale $\tau_{D}\gg\tau_{\rm
lt}$ where the ergodicity of the system is restored via structural
rearrangements of the particles. Below $T_{c}$, the decrease of the
localization length $\xi$ with decreasing temperature is accompanied by a
rapid increase of $\tau_{D}$ and therefore also with an increase of the life
time $\tau_{\rm lt}$ of the amorphous solid state, such that at sufficiently
low temperatures below $T_{c}$, the life time $\tau_{\rm lt}$ may reach
macroscopic time scales.

One may expect that the response of a supercooled liquid to an external
mechanical load such as a shear field is qualitatively different far below
$T_{c}$ from the response above and around $T_{c}$. This is due to the solid-
like behavior over a large time scale $\tau_{\rm lt}$ in the former case. A
system in an ideal amorphous solid state (i.e. with $\tau_{\rm lt}=\infty$) is
associated with a broken continuous translation symmetry which implies its
rigidity and the presence of long-range density correlations szamel2011  as
well as a far-field decay of frozen-in stress fluctuations maier2017 . When
shearing a three-dimensional ideal amorphous solid with a constant strain rate
$\dot{\gamma}$ in a planar Couette flow geometry, in the steady state, a
flowing fluid state with a constant shear stress $\sigma_{\rm ss}$ is
obtained. In the limit $\dot{\gamma}\to 0$, the stress $\sigma_{\rm ss}$ is
non-zero and reaches the yield stress $\sigma_{\rm yield}$. Note that
extensions of MCT to glassforming liquids under shear have been proposed
fuchs2002 ; miyazaki2002 ; miyazaki2004 ; szamel2004 . In the framework of the
MCT by Fuchs and Cates fuchs2002 , a yield stress is predicted for systems
below $T_{c}$ amann2013 ; amann2015 .

Thus, in an ideal amorphous solid, due to the broken translation symmetry, the
shear viscosity $\eta$ is infinitely large and one does not obtain a Newtonian
behavior with $\sigma_{\rm ss}=\eta\dot{\gamma}$ in the limit $\dot{\gamma}\to
0$. However, this is certainly different in a supercooled liquid far below
$T_{c}$ that is associated with a large but finite value of the time scale
$\tau_{\rm lt}$ on which it can be considered to be in an amorphous solid
state. In such a system, one expects on the one hand a Newtonian behavior for
$\dot{\gamma}^{-1}>\tau_{D}$ and on the other hand a solid-like response for
shear rates with $\dot{\gamma}^{-1}<\tau_{\rm lt}<\tau_{D}$. In the latter
case, shear rates have to be sufficiently small such that the resulting
steady-state stress $\sigma_{\rm ss}$ is only slightly larger than an apparent
yield stress that can be obtained via extrapolation to the limit
$\dot{\gamma}\to 0$ (see below).

In this work, the latter regime is studied for a model glassformer using non-
equilibrium molecular dynamics (NEMD) computer simulation. The model under
consideration is a polydisperse soft-sphere system that has been recently
proposed by Ninarello et al. ninarello2017 . It allows the application of the
swap Monte Carlo technique grigera2001  in combination with MD simulation from
which we obtain equilibrated samples far below $T_{c}$, that we use as
starting configurations for NEMD simulations under shear. At sufficiently low
shear rates, the simulations of the sheared samples far below $T_{c}$ show
features that, in computer simulations, have been encountered so far only for
out-of-equilibrium glass states at very low or zero temperature. In
particular, we observe the occurrence of brittle yielding schuh2007 , as
manifested by a sharp stress drop in the stress-strain relation at a strain of
the order of 0.1 ozawa2018 ; popovic2018 ; barlow2020 . Thereby, we
demonstrate that, for an appropriate choice of the shear rate and temperature
$T<T_{c}$, brittle yielding and shear banding can be seen in a supercooled
liquid state, provided that this state exhibits transient elasticity over a
significant time scale $\tau_{\rm lt}$.

Our investigations are complementary to a recent study by Ozawa et al.
ozawa2018  where, for the same model glassformer, first fully equilibrated
samples at different initial temperatures $T_{\rm ini}$ above, around and far
below $T_{c}$ were generated, followed by a quench to zero temperature and
subsequent shear simulations using the athermal quasi-static shear (aqs)
protocol. As we shall see below, our findings are similar to those of Ozawa et
al. when comparing the stress-strain relation of our shear simulations at a
given temperature $T$ and finite shear rate with their aqs calculations for
the corresponding temperature $T_{\rm ini}=T$. As in our case, they observe
brittle yielding for “well-annealed” samples at $T_{\rm ini}\ll T_{c}$ while
for temperatures $T_{\rm ini}$ around and above $T_{c}$ a more ductile
response is seen. The similar response in the aqs calculations and our shear
simulations is remarkable, keeping in mind that, in our simulations, we shear
supercooled liquids at a finite shear rate. In the limit $\dot{\gamma}\to 0$,
i.e. in the “quasi-static” limit, these supercooled liquid states always show
the ductile mechanical response of a Newtonian liquid. This is also true for
temperatures below $T_{c}$ where elasticity has to be considered as a
transient phenomenon, albeit over a very long time scale $\tau_{\rm lt}$ for
temperatures far below $T_{c}$. The fact that the aqs simulations do not show
a Newtonian response for initial temperatures $T_{\rm ini}<T_{c}$ indicates
that for well-annealed samples processes that would lead to a Newtonian
response are suppressed in the framework of the aqs scheme and one obtains the
response of a solid with a finite yield stress.

The occurrence of brittle yielding is associated with the formation of shear
bands. Shear banding is a ubiquitous phenomenon in glasses under mechanical
load schuh2007 ; ozawa2018 ; besseling2010 ; divoux2010 ; chikkadi2011 ;
divoux2016 ; maass2015 ; bokeloh2011 ; binkowski2016 ; hubek2020 ; varnik2003
; bailey2006 ; shi2006 ; shi2007 ; ritter2011 ; sopu2011 ; chaudhuri2012 ;
dasgupta2012 ; dasgupta2013 ; albe2013 ; shiva2016_2 ; golkia2020 ; singh2020
; parmar2019 . Especially in metallic glasses, shear bands lead to
inhomogeneities in the microstructure and can cause a catastrophic failure of
the material schuh2007 ; maass2015 ; hubek2020 . In aqs simulations of a
glassforming binary Lennard-Jones mixture, Parmar et al. parmar2019  have
demonstrated that shear-banded states can be stabilized by applying
oscillatory shear with an appropriate strain amplitude, thereby obtaining
states where a fluidized band coexists with a stress-released amorphous solid.
This indicates that at a given strain above the yield strain, shear-banded
states minimize the energy of the system.

Unlike previous studies, in this work, we observe brittle yielding and shear
banding in transient amorphous solids under equilibrium conditions. We find
two types of shear-banded states right after the yielding transition, namely
states with horizontal and states with vertical shear bands. The formation of
both types of shear bands is an efficient way of releasing stresses, i.e. the
magnitude of the stress drops is similar in both cases. However, in the case
of the vertical bands, the stress shows an increase with strain up to a second
maximum and a second, albeit smaller, stress drop which is associated with the
formation of a horizontal shear band in addition to the vertical one. The
formation of shear bands is also associated with a drop of the potential
energy such that, after the drop, the potential energy is monotonously
increasing towards the steady state value. Recently, the occurrence of
horizontal and vertical shear bands has been also observed in sheared low-
temperature glass states of a binary Lennard-Jones mixture golkia2020 ;
however, in the present study, we find these features in equilibrated systems.

The rest of the paper is organized as follows: In the next section (Sec. II),
the details of the model potential, the simulation techniques, and the
simulation protocols are reported. Section 3 presents results on the
equilibrium dynamics of supercooled liquids, focussing on the change of the
dynamics around the MCT critical temperature. Section 4 is devoted to the
analysis of supercooled liquids under shear. Here, we address the question
under which conditions brittle yielding and shear banding occur. Finally, Sec.
5 summarizes the results and draws conclusions.

##  II Model and details of the simulation

We consider a model of polydisperse non-additive soft spheres that has been
recently proposed by Ninarello et al. ninarello2017 . In this model,
interactions between particles are pair-wise additive. To each particle $i$, a
diameter $\sigma_{i}$ is assigned according to a probability distribution
$P(\sigma)=A\sigma^{-3}$ with $A=2/(\sigma_{\rm min}^{-2}-\sigma_{\rm
max}^{-2})$. We have chosen $\sigma_{\rm min}=0.725\,\bar{\sigma}$ and
$\sigma_{\rm max}=\sigma_{\rm min}/(2\sigma_{\rm min}-1)\approx
1.611\,\bar{\sigma}$. This choice of $\sigma_{\rm min}$ and $\sigma_{\rm max}$
provides that the first moment of $P(\sigma)$ is equal to $\bar{\sigma)}$;
$\bar{\sigma}=1.0$ is used as the length unit in the following. The
interactions between pairs of particles depend on the variable
$x_{ij}=r_{ij}/\sigma_{ij}$ where $r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$ is the
distance between particle $i$ at position $\vec{r}_{i}$ and particle $j$ at
position $\vec{r}_{j}$ and
$\sigma_{ij}=0.5(\sigma_{i}+\sigma_{j})(1-0.2|\sigma_{i}-\sigma_{j}|)$
introduces the non-additivity of the particle diameters. Note that the non-
additivity is essential to avoid any crystallization when the swap Monte Carlo
method is applied (see below).

The interaction potential between a pair of particles is defined by

$V(x)=\begin{cases}V_{0}(x^{-12}+c_{0}+c_{2}x^{2}+c_{4}x^{4})&\mathrm{for}\quad
x<x_{c}\\\ 0&\mathrm{for}\quad x\geq x_{c}\,,\end{cases}$ (1)

where the cut-off $x_{c}=1.25$ is chosen. The terms with the parameters
$c_{0}=-28/x_{c}^{12}$, $c_{2}=48/x_{c}^{14}$, $c_{4}=-21/x_{c}^{16}$ ensure
the smoothness of the function $V(x)$ at $x=x_{c}$. The parameter $V_{0}=1.0$
sets the unit of energy in the following.

The simulations at constant particle number $N$, constant volume $V$, and
constant temperature $T$ are performed with the LAMMPS package plimpton1995 .
The number density is fixed at $\varrho=N/V=1.0$. The masses of the particles
are set to $m=1.0$. In the molecular dynamics (MD) simulations, Newton’s
equations of motion are integrated by the velocity Verlet algorithm allenbook
, using a time step of $\delta t=0.01\,\tau_{\rm MD}$ (with $\tau_{\rm
MD}=(\bar{\sigma}^{2}m/V_{0})^{1/2}$). The temperature is kept fixed by a DPD
thermostat allenbook ; soddemann2003 , using a similar implementation as in
Ref. golkia2020  with the friction coefficient $\zeta=1.0$ and the weight
function $\omega(r)=1.0$ for $r\leq 1.3x_{c}$ and $\omega(r)=0$ otherwise [cf.
Eqs. (2)-(6) in Ref. golkia2020 ]. The DPD thermostat locally conserves the
momentum and is Galilean invariant. This is especially advantageous for the
non-equilibrium MD simulations under shear, because the Galilean-invariant
thermostat does not introduce any bias with respect to the direction of the
velocity flow.

To obtain fully equilibrated samples at very low temperatures, a combination
of MD simulation and the swap-Monte-Carlo (SMC) technique grigera2001  is
used. In a “trial SMC move”, one randomly selects a pair of particles and
exchanges their diameters. Then, this move is accepted or rejected according
to a Metropolis criterion allenbook . In our hybrid scheme, every 25 MD steps
$N$ trial SMC moves are performed. In the considered temperature range,
$0.01\leq T\leq 0.3$, the acceptance rate for the SMC moves varies between 10
and 22% (with a decreasing acceptance rate with decreasing temperature). The
longest equilibration runs with the hybrid MD-SMC method were over $10^{7}$
time steps which allowed to fully equilibrate samples with $N=1372$, 2048,
6000, and 10000 particles at the temperature $T=0.06$, corresponding to the
glass transition temperature $T_{g}$ in our study.

Non-equilibrium MD simulations are employed to shear the samples in a planar
Couette flow geometry. The shear is imposed via Lees-Edwards boundary
conditions lees1972  along the $xz$ plane in the direction of $x$. For the
simulations under shear, we have integrated the equations of motion with the
time step $\delta t=0.001\,\tau_{\rm MD}$. Most of the data shown below
correspond to the temperatures $T=0.15$, 0.11, 0.09, 0.07, and 0.06 for a
system of $N=10000$ particles. At each temperature, 30 runs were performed,
starting from statistically independent samples that were fully equilibrated
via the MD-SMC method. The considered shear rates range from
$\dot{\gamma}=10^{-6}$ to $\dot{\gamma}=10^{-3}$. For the calculation of the
stress-strain relations, we have performed a running average over strain
windows of width $\delta\gamma=10^{-4}$.

##  III From liquid to amorphous solid: Equilibrium dynamics

The dynamics of supercooled liquids is associated with the cage effect. At
sufficiently low temperatures, the particles are trapped in cages formed by
the surrounding particles and the breaking of cages requires collective
particle rearrangements that slow down with decreasing temperature. As we
shall see below, around the critical temperature of mode coupling theory
(MCT), $T_{c}$, the system gradually transforms from a liquid-like state to a
state that can be characterized as an amorphous solid. This transition is due
to the localization of the particles in their cages and, as we shall see in
the next section, the response to an external shear changes drastically from
the liquid-like state above $T_{c}$ to the amorphous solid well below $T_{c}$,
especially with respect to the yielding behavior. In this section, we first
present the “equation of state” of our system, i.e. the temperature dependence
of the potential energy per particle, and then study the one-particle dynamics
in terms of the mean-squared displacement (MSD) of a tagged particle. From the
MSD, a localization length is determined that indicates the transition from
liquid to solid-like behavior around $T_{c}$. Furthermore, we estimate the
life time $\tau_{\rm lt}$ of the amorphous solid as a function of temperature.
We have computed the equation of state from fully equilibrated configurations
that we have obtained via hybrid MD-SMC simulations at constant temperature.
For the calculation of the MSD, we have used such fully equilibrated samples
as starting configurations for microcanonical runs where we have switched off
the SMC and the coupling to the thermostat.

Figure 1: Potential energy per particle, $e_{\rm pot}$, as a function of
temperature $T$. The solid line is a fit with the function
$f(T)=A+BT^{\gamma}$ with $A=0.0120017$, $B=1.4999$, and $\gamma=0.637135$.
Indicated are the locations of the critical mode coupling temperature,
$T_{c}=0.104$, and the glass transition temperature, $T_{g}=0.06$.

Figure 1 shows the potential energy per particle, $e_{\rm pot}$, as a function
of temperature. In this plot, the critical MCT temperature at $T_{c}=0.104$ as
well as the glass transition temperature at $T_{g}=0.06$ are indicated. The
MCT temperature $T_{c}$ was determined from fits to dynamic quantities such as
the mean-square displacement (see below). Below $T_{g}$, the hybrid MD-SMC
runs on the time scale of $10^{5}\,\tau_{\rm MD}$ are no longer sufficient to
fully equilibrate the system. The data for $T\geq T_{g}$ can be well described
by the function (solid line in Fig. 1)

$f(T)=A+BT^{\gamma}$ (2)

with $A$, $B$, and $\gamma$ being fit parameters. While the density functional
theory of Rosenfeld and Tarazona rosenfeld1998  predicts the exponent
$\gamma=0.6$ for simple high-density soft-sphere fluids, we find the exponent
$\gamma\approx 0.64$ which is very close to this prediction. Note that Eq. (2)
with a value of $\gamma$ around 0.6 also provides a good approximation for
other glassforming liquids with a $1/r^{n}$-type interactions at low
temperature (for a detailed discussion see Ref. ingebrigtsen2013 ).

Figure 2: a) MSD as a function of time for temperatures $T<T_{g}$ (dashed blue
lines), $T_{g}\leq T<T_{c}$ (solid blue lines), and $T>T_{c}$ (red solid
lines). b) MSDs at $T=0.11$, $T=0.09$, $T=0.07$, and $T=0.06$. The dashed
lines are fits to Eq. (4) and the filled circles mark the location of the life
time $\tau_{\rm lt}$ for the different temperatures (see text). The inset is a
schematic plot of the MSD that illustrates the definition of $\tau_{\rm lt}$.
c) Life time $\tau_{\rm lt}$ as a function of inverse temperature. The solid
line is a fit with an Arrhenius law (see text). The MSDs in a) and b)
correspond to systems with $N=2048$ particles.

Now we come to the one-particle dynamics of the system and investigate the MSD
of a tagged particle, defined by

$\left\langle\delta
r^{2}(t)\right\rangle=\frac{1}{N}\sum_{i=1}^{N}\left\langle|\vec{r}_{i}(t)-\vec{r}_{i}(0)|^{2}\right\rangle$
(3)

with $\vec{r}_{i}(t)$ the position of particle $i$ at time $t$. The brackets
$\langle\dots\rangle$ represent an ensemble as well as a time average over the
different samples. Note, however, that for states below $T_{g}$ we have only
applied an ensemble average. The MSDs are calculated from microcanonical MD
simulation for a system of $N=2048$ particles, using as initial configurations
60 independent samples from the MD-SMC simulations.

In Fig. 2a, MSDs are plotted double-logarithmically for different
temperatures. Here, we have marked the different temperature regimes. The red
solid lines correspond to temperatures above $T_{c}$ at $T=0.105$, 0.11,
0.115, 0.12, 0.125, 0.13, 0.14, 0.15, 0.16, 0.18, 0.20, 0.22, 0.25, and 0.3.
At the highest temperature, $T=0.3$, the MSD displays a ballistic regime
$\propto t^{2}$ at very short times, an emerging shoulder at intermediate
times, and a diffusive regime $\propto t$ in the long-time limit. With
decreasing temperature, the diffusive regime shifts to longer times and the
intermediate time regime evolves into a plateau. The blue solid lines show the
MSDs for temperatures $T_{g}<T<T_{c}$ at $T=0.06$, 0.065, 0.0675, 0.07, 0.075,
0.0775, 0.08, 0.0825, 0.085, 0.0875, 0.09, 0.0925, 0.095, 0.0975, 0.10, and
0.1025. Here, the initial configurations are fully equilibrated samples from
the MD-SMC simulations. However, the microcanonical MD runs over a time scale
of $10^{5}\,\tau_{\rm MD}$ are not long enough to reach a diffusive regime far
below $T_{c}$. So at $T=0.06$, we hardly see deviations from the plateau at
long times. The MSDs below $T_{g}$ in Fig. 2a (blue dashed lines) correspond
to the temperatures $T=0.01$, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045,
0.05, and 0.055. Here, the MSDs display a plateau for $1\leq t\leq 10^{5}$,
the height of which decreases with decreasing temperature. Note that the small
overshoot in the low-temperature MSDs around $t\approx 0.8$ is associated with
the microscopic dynamics horbach1996 ; horbach2001 . This feature disappears
for larger system sizes (e.g. for our model it cannot be seen anymore for
systems with $N=10000$ particles).

The emergence of a shoulder that evolves into a plateau at low temperature
manifests the caging of the particles. MCT provides detailed predictions about
the behavior of the MSD around the plateau (as well as corresponding
predictions for the plateau-like regions in intermediate scattering functions
goetze2009 ). One of them describes the initial increase of the MSD from the
plateau and is given by goetze2009

$\phi(t)=\delta r^{2}_{\rm plateau}+ht^{b}+h_{2}t^{2b}\,.$ (4)

This equation corresponds to a von Schweidler law, extended by a correction
term $\propto t^{2b}$. $\delta r^{2}_{\rm plateau}$ quantifies the height of
the (emerging) plateau in the MSD, $h$ and $h_{2}$ are temperature-dependent
amplitudes, and the exponent $b$ is expected to be universal for a given
system (but it may vary for different systems in the range $0<b\leq 1$).
Figure 2b shows the MSDs at $T=0.11$, $T=0.09$, $T=0.07$, and $T=0.06$
together with fits to Eq. (4). These fits and also the fits to the MSDs at the
other temperatures were performed with the constant exponent value $b=0.59$.
Note, however, that the values for $\delta r^{2}_{\rm plateau}$, as obtained
from the fit to Eq. (4), are not very sensitive with respect to the choice of
the exponent $b$.

Using the fits to Eq. (4), we can now introduce a definition of the life time
$\tau_{\rm lt}$ of the transient amorphous solid state for the different
temperatures. To this end, we define $\tau_{\rm lt}$ as the time for which
$\left\langle\delta r^{2}(\tau_{\rm lt})\right\rangle/\delta r^{2}_{\rm
plateau}=1.3$ (see the inset of Fig. 2b for an illustration of this
definition). The locations of $\tau_{\rm lt}$ for the MSDs in Fig. 2b are
marked as filled circles.

Figure 2c shows the logarithm of the time scale $\tau_{\rm lt}$ as a function
of inverse temperature. For $T\lesssim 0.09$, the data can be well fitted by
an Arrhenius law $f(T)=\tau_{0}\,\exp\left(E_{\rm A}/T\right)$, which is
represented by the bold solid line in the figure. The values of the fit
parameters are $\tau_{0}=3.3\times 10^{-6}$ and $E_{\rm A}=1.43641$. Here, the
energy $E_{\rm A}$ can be interpreted as an activation energy. The application
of the Arrhenius law and thus the interpretation of a kinetic process as an
activated one are only sensible if the ratio of the activation energy to the
thermal energy, $E_{\rm A}/T$, is much larger than unity riskenbook . In our
case, this ratio varies between about 16 at $T=0.09$ and about 24 at $T=0.06$
which is consistent with the condition $E_{\rm A}/T\gg 1$. At temperatures
$T\gtrsim T_{c}$, we observe significant deviations from the Arrhenius
behavior and $\tau_{\rm lt}$ is close to the microscopic time scale $\tau_{\rm
MD}$. From the temperature dependence of $\tau_{\rm lt}$ we can conclude that
around $T_{c}$ there is a gradual crossover towards an activated dynamics with
decreasing temperature.

Figure 3: Localization length divided by the mean nearest-neighbor distance,
$\xi/\bar{d}$, as a function of temperature. The dashed red line marks the
critical value of the reduced localization length, $\xi_{c}/\bar{d}\approx
0.077$.

In the framework of the Gaussian approximation hansenbook ; thorneywork2016 ;
fuchs1998 , one can relate $\delta r^{2}_{\rm plateau}$ to a localization
length $\xi$ as

$\xi^{2}=\frac{1}{6}\delta r^{2}_{\rm plateau}\,.$ (5)

Figure 3 shows the temperature dependence of $\xi$, scaled with the average
nearest-neighbor distance $\bar{d}\approx 1.07$ (we have estimated $\bar{d}$
from the location of the first peak of the radial distribution function at
$T=0.06$). At $T=0.01$, i.e. far below $T_{g}$, the reduced localization
length is $\xi/\bar{d}\approx 0.02$. It increases with increasing temperature.
At $T_{g}$, $\xi/\bar{d}$ slightly changes slope and then increases roughly
linearly up to $T_{c}$ where it reaches the constant $\xi_{c}/\bar{d}\approx
0.077$. The critical value, $\xi_{c}$, of the localization length marks the
stability limit of the amorphous solid, i.e. for $T>T_{c}$ the system is in a
liquid state. In analogy to crystalline solids, the critical value
$\xi_{c}/\bar{d}$ can be interpreted as a Lindemann criterion for the
stability of an amorphous solid goetze2009 . Note that Fuchs et al. fuchs1998
have obtained $\xi_{c}/\bar{d}\approx 0.0746$ in a calculation for a hard
sphere system in the framework of MCT, thus a value that is very close to our
finding.

The behavior of both $\tau_{\rm lt}$ and $\xi$ indicate a gradual change of
the dynamics around $T_{c}$. Below $T_{c}$, the localization of particles in
their cages, as quantified by $\xi_{c}/\bar{d}$, is below the stability limit,
given by $\xi_{c}/\bar{d}\approx 0.077$. As a consequence, there is the
emergence of transient amorphous solid state for $T<T_{c}$, the life time
$\tau_{\rm lt}$ of which follows an Arrhenius law with an activation energy of
about 1.44. The gradual change from liquid-like to solid-like dynamics is also
associated with a qualitative change of the system’s response to an external
shear. As we shall see in the next section, brittle yielding and the formation
of shear bands can be observed in the supercooled liquid below $T_{c}$. These
features are typical for the response of low-temperature glasses to a
mechanical load. In the following, we shall analyze the conditions for the
occurrence of brittle yielding and shear banding in deeply supercooled
liquids. An important parameter in this context is the time scale $\tau_{\rm
lt}$. For example, for $T=0.06$, the life time $\tau_{\rm lt}$ is of the order
of $10^{5}$ (Fig. 2c). Therefore, for $\dot{\gamma}\gtrsim 10^{-5}$, the
product $\dot{\gamma}\tau_{\rm lt}$ is lower equal unity and one may expect
the shear response of an amorphous solid.

##  IV Supercooled liquids under shear

Now we analyze the results for equilibrated supercooled liquids under shear.
Our focus is on the temperature range $0.06\leq T\leq 0.15$ to study the
response to the external shear from liquid-like states slightly above $T_{c}$
to the solid states far below $T_{c}$. As we have seen in the previous
section, the latter states can be characterized via the localization length
$\xi$ being significantly lower than the critical value $\xi_{c}$.

Figure 4: Stress-strain relation for the temperature $T=0.06$ and the shear
rate $\dot{\gamma}=10^{-3}$. Indicated are the maximum $\sigma_{\rm max}$, the
steady-state stress $\sigma_{\rm ss}$, and the definition of the stress drop
$\Delta\sigma$.

A typical stress-strain relation, indicating a non-Newtonian response of the
supercooled liquid, is shown in Fig. 4 for the temperature $T=0.06$ and the
shear rate $\dot{\gamma}=10^{-3}$. While the strain is given by
$\gamma=\dot{\gamma}t$, the stress $\sigma$ was computed from the virial
equation, as described in Ref. golkia2020 . Different regimes can be
identified in the figure. First the stress increases almost linearly up to a
maximum value $\sigma_{\rm max}$ which is reached at a strain $\gamma\approx
0.135$ in this case. The maximum in the stress marks the transition from an
elastic deformation of the “solid” to the onset of plastic flow. During the
plastic deformation, the stress drops from $\sigma_{\rm max}$ towards the
steady-state stress $\sigma_{\rm ss}$ which can be quantified by
$\Delta\sigma=\sigma_{\rm max}-\sigma_{\rm ss}$. In the steady state, the
system can be described by a flowing homogeneous liquid. Heterogeneous flow
patterns are observed between the onset of plastic flow and the steady state.
The morphology of these flow patterns, especially with respect to the
dependence on temperature and shear rate, shall be elucidated in the
following.

Figure 5: Flow curves for the temperatures $T=0.06$, 0.09, 0.11, and 0.15. The
solid line is a fit with a Herschel-Bulkley law to the data for $T=0.06$ (see
text). The dashed line indicates a linear behavior, $\sigma_{\rm
ss}\propto\dot{\gamma}$.

In Fig. 5, the steady-state stress $\sigma_{\rm ss}$ as a function of the
shear rate $\dot{\gamma}$ (i.e. the flow curve) is plotted double-
logarithmically for different temperatures above and below $T_{c}$. For
sufficiently low shear rates, one expects that the system behaves like a
Newtonian fluid with a linear increase of the stress as a function of the
shear rate, $\sigma_{\rm ss}=\eta\dot{\gamma}$ (with $\eta$ the shear
viscosity). At $T=0.15$, we can still identify a Newtonian regime (dashed
line), followed by sublinear shear-thinning regime for $\dot{\gamma}>5\times
10^{-4}$. At $T=0.11$, the Newtonian regime is not anymore in the window of
considered shear rates $\dot{\gamma}\geq 10^{-6}$. Here, we observe an
emerging plateau around $\dot{\gamma}=10^{-4}$ that becomes more pronounced at
$T=0.09$, and eventually, at $T=0.06$, the data can be well fitted by a
Herschel-Bulkley law herschel1926  (solid line), $\sigma_{\rm ss}=\sigma_{\rm
yield}+A\dot{\gamma}^{\alpha}$ with the yield stress $\sigma_{\rm
yield}=0.0900265$, the amplitude $A=1.57446$, and the exponent
$\alpha=0.477008$. Note that at $T=0.07$ the flow curve is very similar to
that at $T=0.06$. So at the lowest considered temperatures where we are able
to obtain a fully equilibrated state, our system can be seen as a yield stress
material in equilibrium (although we also expect at these temperatures the
occurrence of a Newtonian regime at extremely low shear rates).

Figure 6: Reduced stress $\sigma^{\star}$ as a function of $\gamma-\gamma_{\rm
max}$ at the different temperatures $T=0.06$, 0.07, and 0.09 for the shear
rates a) $\dot{\gamma}=10^{-3}$, b) $\dot{\gamma}=10^{-4}$, and c)
$\dot{\gamma}=10^{-5}$. The data corresponds to systems with $N=10000$
particles. In c), the dotted lines correspond to data for $N=12000$ and
$N=96000$, as adapted from Ref. ozawa2018 . The dashed and the dotted line in
a) are fits with with a compressed exponential and a logarithm, respectively
(see text).

Having characterized the steady-state behavior of our system under shear, we
now investigate the relaxation of the stress from the onset of plastic flow
(marked by the maximum stress $\sigma_{\rm max}$ at a given temperature) to
the steady-state stress. To this end, we define the reduced stress

$\sigma^{\star}=\frac{\sigma-\sigma_{\rm ss}}{\sigma_{\rm max}-\sigma_{\rm
ss}}$ (6)

which is displayed in Fig. 6 for three shear rates and three temperatures
below $T_{c}$ as a function of $\gamma-\gamma_{\rm max}$ (with $\gamma_{\rm
max}$ the strain corresponding to $\sigma_{\rm max}$). At
$\dot{\gamma}=10^{-3}$ (Fig. 6a), the decay of $\sigma^{\star}$ for $T=0.09$
can be described by the compressed exponential $\exp[-((\gamma-\gamma_{\rm
max})/\delta\gamma^{\star})^{a_{\rm ce}}]$ with
$\delta\gamma^{\star}=0.115946$ and $a_{\rm ce}=1.33297$. Also for the two
lower temperatures $T=0.06$ and 0.07, the reduced stress decays on a strain
scale $\delta\gamma^{\star}\approx 0.1$, but the functional form of its decay
changes around $\gamma-\gamma_{\rm max}\approx 0.08$ in that the compressed-
exponential-like decay is followed by a logarithmic one
$\propto\ln[(\gamma-\gamma_{\rm max})/1.07353]$ (dotted line in Fig. 6a,
fitted to the “tail” of the $T=0.06$ curve). The difference in the decay of
$\sigma^{\star}$ with respect to temperature becomes more pronounced at the
lower shear rates $10^{-4}$ (Fig. 6b) and $10^{-5}$ (Fig. 6c). While at
$T=0.09$, the reduced stress still decays essentially with a compressed
exponential on the strain scale $\delta\gamma^{\star}\approx 0.1$, at the two
lower temperatures the initial decay is significantly faster and
$\sigma^{\star}$ exhibits a local maximum around $\gamma-\gamma_{\rm
max}\approx 0.2$. The strain scale of the initial decay decreases with
decreasing temperature and shear rate. At $T=0.06$ and $\dot{\gamma}=10^{-5}$,
the reduced stress $\sigma^{\star}$ decays on the strain scale
$\delta\gamma^{\star}\approx 0.01$.

Also included in Fig. 6c are data for $N=12000$ and $N=96000$, as adapted from
the simulation study of Ozawa et al. ozawa2018  using an aqs protocol. It is
remarkable that the reduced stress for $N=12000$ from Ozawa et al. agrees well
with our data for a comparable system of $N=10000$ particles, although we
consider a system at a finite temperature as well as a finite strain rate and,
moreover, our system is in an equilibrated supercooled liquid state (note,
however, that $\sigma_{\rm max}$ is significantly larger in the athermal
case).

Figure 7: Reduced stress $\sigma^{\star}$ as a function of $\gamma-\gamma_{\rm
max}$ at the shear rate $\dot{\gamma}=10^{-5}$ and temperatures a) $T=0.09$,
b) $T=0.07$, and c) $T=0.06$ for different system sizes, as indicated. The
dashed lines are fits with compressed exponentials (see text). Also included
in b) is the reduced stress for $\dot{\gamma}=10^{-6}$ and $N=10000$.

It is tempting to interpret the rapid drop of the stress as a first-order
phase transition, as proposed by Ozawa et al. ozawa2018 . The interpretation
of the stress drop as a phase transition would be appropriate in the limit of
zero shear rate, $\dot{\gamma}\to 0$. So we have to take this limit in some
sensible manner, keeping in mind that the expected true behavior of the system
in the zero shear-rate limit is that of a Newtonian fluid for which
$\sigma_{\rm ss}\propto\dot{\gamma}$ and the absence of any stress drop in the
stress-strain relation. However, the fluid curves for $T\leq 0.07$ suggest
that the systems can be considered as a yield stress fluid also at very low
shear rates and one obtains $\sigma_{\rm yield}$ by extrapolation via the
Herschel-Bulkley law. Below, we perform a similar extrapolation to obtain the
initial strain scale $\delta\gamma^{\star}$ with which the stress decays from
$\sigma_{\rm max}$ to $\sigma_{\rm ss}$ in the limit $\dot{\gamma}\to 0$.

The existence of a first-order transition in the thermodynamic limit implies
that it is rounded for finite systems and it becomes sharper with increasing
system size. Figure 7 shows the decay of the reduced stress $\sigma^{\star}$
for different system sizes at the temperatures $T=0.09$, $T=0.07$, and
$T=0.06$ in panels a), b), and c), respectively. For all three temperatures,
the shear rate is $\dot{\gamma}=10^{-5}$. While for $T=0.09$ there is almost
no dependence of $\sigma^{\star}$ on system size, for the two lower
temperatures the decay of $\sigma^{\star}$ becomes significantly sharper with
increasing system size. For $T=0.07$ (Fig. 7b), we have also included the
reduced stress for the lower shear rate $\dot{\gamma}=10^{-6}$ and $N=10000$
which exhibits a less rapid decay than the corresponding result for
$\dot{\gamma}=10^{-5}$. This can be explained in terms of the life time
$\tau_{\rm lt}$ in relation to the shear rate $\dot{\gamma}$. Above we have
estimated $\tau_{\rm lt}\approx 3300$ for $T=0.07$ (cf. Fig. 2c) and thus the
time scale $\dot{\gamma}^{-1}$ is much larger than $\tau_{\rm lt}$ for both
shear rates $10^{-5}$ and $10^{-6}$. The yielding of the system interferes
with structural relaxation processes in this case and this certainly in a more
pronounced manner for $\dot{\gamma}=10^{-6}$ than for $\dot{\gamma}=10^{-5}$.
Therefore, the reduced stress decays faster for the higher shear rate of
$\dot{\gamma}=10^{-5}$.

The dashed lines in Fig. 7 are fits with compressed exponentials. In these
fits, the exponent $a_{\rm ce}$ is around 3.0 and the strain scale changes
from $\delta\gamma^{\star}\approx 0.027$ for $N=1000$ to
$\delta\gamma^{\star}\approx 0.012$ for $N=10000$. While the initial decay of
$\sigma^{\star}$ strongly depends on $N$, the second feature in
$\sigma^{\star}$, the appearance of a local maximum at $\gamma-\gamma_{\rm
max}\approx 0.2$, does not show significant finite-size effects.

Figure 8: Strain scale $\delta\tilde{\gamma}^{\star}$ as a function of shear
rate $\dot{\gamma}$ at the temperatures $T=0.07$ [a)] and $T=0.06$ [b)] for
different system sizes. In b), the solid lines are fits with the function
$g(\dot{\gamma})=\delta\tilde{\gamma}^{\star}_{0}+A_{\gamma}\dot{\gamma}^{c_{\gamma}}$
(for details see text). The inset shows $\delta\tilde{\gamma}^{\star}_{0}$ as
a function of $N^{-1}$.

In the following, we do not use the strain scale, as directly obtained from
the fits to the compressed exponentials. In lieu thereof, we use the value
where the reduced stress, as described by the compressed exponential, has
decayed to 0.2. We denote this quantity by $\delta\tilde{\gamma}^{\star}$.
Figure 8 displays the shear-rate dependence of $\delta\tilde{\gamma}^{\star}$
for $T=0.07$ in a) and $T=0.06$ in b) and different system sizes. In the case
of $T=0.07$, the transition becomes significantly sharper for all shear rates
when changing the particle number from $N=1500$ to $N=6000$. However, only
small changes are observed when going from $N=6000$ to $N=10000$ (at
$\dot{\gamma}=10^{-6}$, the values for $\delta\tilde{\gamma}^{\star}$ are
essentially equal for the two system sizes). This is due to the fact that the
life time $\tau_{\rm lt}$ is smaller than the time scale $\dot{\gamma}^{-1}$
for $\dot{\gamma}\leq 10^{-4}$ and thus the yielding transition interferes
with relaxation processes in the liquid.

In Fig. 8b for $T=0.06$, the solid lines correspond to the fit function
$g(\dot{\gamma})=\delta\tilde{\gamma}^{\star}_{0}+A_{\gamma}\dot{\gamma}^{c_{\gamma}}$,
with $\delta\tilde{\gamma}^{\star}_{0}$ the estimate of
$\delta\tilde{\gamma}^{\star}$ at zero shear rate, $A_{\gamma}$ an amplitude,
and $c_{\gamma}$ an exponent that has value of about $0.62$ in the fits of
Fig. 8b. Thus, the zero shear-rate values of the strain scale can be well
estimated via a “Herschel-Bulkley-like” law. The inset of Fig. 8b shows
$\delta\tilde{\gamma}^{\star}_{0}$ as a function of $N^{-1}$. In the
considered range of system sizes, we observe a weak dependence of
$\delta\tilde{\gamma}^{\star}_{0}$ on system size. The data suggests that
there might be a regime $\propto 1/N$ for large $N$, similar to what one
expects for a first-order phase transition. However, our data is not
conclusive to support this interpretation.

But what happens at this yielding transition? And what is the meaning of the
second feature in $\sigma^{\star}$, i.e. the increase of $\sigma^{\star}$ up
to a maximum around $\gamma-\gamma_{\rm max}=0.2$? Below we show that both
features are connected to the formation of shear bands. The rapid initial
decay of $\sigma^{\star}$ is a manifestation of brittle yielding.

Figure 9: Stress-strain relations for a) horizontal shear bands and b)
vertical shear bands at the temperature $T=0.06$ and the shear rate
$\dot{\gamma}=10^{-5}$ in the main figures and for different shear rates in
the insets. The grey lines in the main plots are the stress-strain relations
for the individual runs and the black lines correspond to the average over
these runs.

To analyze the behavior of the system around yielding, we now consider
individual runs at the temperature $T=0.06$ and the shear rate
$\dot{\gamma}=10^{-5}$. Among the 30 independent runs for the systems with
$N=10000$ particles, we find two types of stress-strain relations. In both
cases, we observe an initial sharp drop, indicating brittle yielding. However,
while we see in the first type only the initial stress drop (Fig. 9a), in the
second one there is the additional increase after the first drop up to
$\gamma\approx 0.32$, followed by a second drop of the stress (Fig. 9b). The
first type of stress-strain relation corresponds to the formation of a
horizontal shear band, i.e. the occurrence of a thin melted layer with an
orientation parallel to the flow direction. The second type of stress-strain
relation corresponds to the initial formation of a vertical shear band where
the melted thin layer is oriented perpendicular to the flow direction. Note
that among the 30 runs, we have observed horizontal and vertical shear bands
in 18 and in 12 cases, respectively. In Fig. 9, the grey lines correspond to
the individual runs and black ones to the average over these runs in each
case. The insets of Fig. 9 show the averaged stress-strain relations for
different shear rates. For $\dot{\gamma}\leq 10^{-4}$, a qualitatively similar
behavior is seen with essentially the initial stress drop getting slightly
sharper with decreasing shear rate (cf. Fig. 8). At $\dot{\gamma}=10^{-3}$,
however, the yielding transition is washed out and, instead of the second
maximum in the stress-strain relation for the vertical shear bands, there is a
logarithmic decay around $\dot{\gamma}=0.3$ (cf. Fig. 6a).

Figure 10: Potential energy per particle as a function of strain for a)
horizontal shear bands and b) vertical shear bands at the temperature $T=0.06$
and the shear rate $\dot{\gamma}=10^{-5}$ in the main figures and for
different shear rates in the insets. The grey lines in the main plots are the
stress-strain relations for the individual runs and the black lines correspond
to the average over these runs. The horizontal dashed lines in both panels
mark the average steady-state value, $\bar{e}_{\rm pot}=0.305$.

The brittle yielding is also reflected in the behavior of the potential energy
per particle, $e_{\rm pot}$, as a function of the strain $\gamma$. In the case
of the horizontal shear bands at $T=0.06$ and $\dot{\gamma}=10^{-5}$ (Fig.
10a), there is first a drop of $e_{\rm pot}$ at the yield point, followed by a
slow increase towards the steady-state value which is at about $\bar{e}_{\rm
pot}=0.305$ (horizontal dashed line). As can be inferred from the figure,
there is a large scatter in the values of $e_{\rm pot}$ from sample to sample.
This is due to the polydispersity of the samples. However, the shape of the
curves for the different samples is very similar and they are essentially
shifted with respect to each other. This is also true for the behavior of
$e_{\rm pot}$ vs. $\gamma$ for the case of the vertical shear bands (Fig.
10b). Here, after the first drop of the energy, it increases to a value which
is close to the steady-state value, then, around $\gamma\approx 0.3$
(corresponding to the local maximum in the stress-strain relation), it
slightly decreases before it increases towards the steady-state value. For the
case of the vertical shear bands, the system reaches the steady state much
faster than in the case of the horizontal shear bands. The behavior of $e_{\rm
pot}(\gamma)$ for the different shear rates (see insets of Fig. 10) is similar
to that of the corresponding stress-strain relations.

Figure 11:  Mobility maps at $T=0.06$ and $\dot{\gamma}=10^{-5}$ for a sample
with the formation of a horizontal shear band for a) $\gamma=0.1024$, b)
$\gamma=0.1472$, and c) $\gamma=0.495$ and a sample with the formation of a
vertical shear band for the same values of $\gamma$, d) - f).

To visualize the shear bands, mobility color maps shriva2016_1  are computed.
To this end, we determine, for each particle $i$, the non-averaged MSDs
$\delta y_{i}^{2}(t)$ and $\delta z_{i}^{2}(t)$ in the neutral $y$ direction
and the shear-gradient $z$ direction, respectively. From this, we obtain the
“mobility displacement” $\Delta_{i}(t)=\sqrt{\delta y_{i}^{2}(t)+\delta
z_{i}^{2}}$ and we assign a color to the magnitude of $\Delta_{i}$. The time
origin for the calculation of $\Delta_{i}$, $t=0$, corresponds to the time
where the external shear is switched on. For the snapshots in Fig. 11 at
different values of the strain, we have selected a sample with a horizontal
shear band [a)-c)] and one with a vertical shear band [d)-f)], both samples
are at $T=0.06$ and $\dot{\gamma}=10^{-5}$. At $\gamma=0.1024$, i.e. just
before the onset of plastic flow, the system is in a homogeneously deformed
state and therefore the mobility of the particles is close to zero, as
represented by the blue color. At the strain $\gamma=0.1472$ a horizontal
shear band has formed in the first sample (Fig. 11b) and a vertical one in the
second sample (Fig. 11e). In both cases, the fluidized regions along the band
are represented by particles, colored in green. The horizontal shear band
exhibits a slow growth as a function of strain, as reflected, e.g., in a slow
increase of the potential energy of the system (cf. Fig. 10a). At
$\gamma=0.495$, the thickness of the horizontal shear band corresponds to
about 5-6 $\sigma$, i.e. a few liquefied layers (Fig. 11c). The behavior is
different in the case of the vertical shear band. Here, the thickness of the
vertical band first increases which is accompanied by an increase of the
stress with increasing strain (cf. Fig. 9b). The stress drop at $\gamma\approx
0.3$ is associated with the formation of an additional horizontal shear band
which grows with increasing strain (cf. Fig. 11f).

##  V Summary and conclusions

In summary, we have investigated the yielding behavior of a glassforming soft-
sphere model under shear. Using molecular dynamics (MD) simulation in
combination with the swap Monte Carlo (SMC) technique, fully equilibrated
supercooled liquid samples around and far below the critical mode coupling
temperature $T_{c}$ were obtained.

First, these samples served as starting configurations for simulations in the
microcanonical ensemble to study how the dynamics of the supercooled liquid
changes when decreasing the temperature from above to far below $T_{c}$. In
qualitative agreement with mode coupling theory (MCT), we have seen that the
reduced localization length $\xi/\bar{d}$, as extracted from the mean squared
displacement, shows a kink at $T_{c}$, changing from $\xi_{c}/\bar{d}=0.077$
for $T>T_{c}$ to a roughly linearly decreasing function for decreasing
temperature below $T_{c}$. Here, the critical value $\xi_{c}/\bar{d}$ marks
the stability limit of the amorphous solid. In fact, the decrease of $\xi$
with decreasing temperature is accompanied by an exponential increase of a
time scale $\tau_{\rm lt}$ that measures the life time of the amorphous solid
state. The Arrhenius law that we find for the temperature dependence of
$\tau_{\rm lt}$ is consistent with the interpretation of an activated dynamics
for structural relaxation processes below $T_{c}$.

The gradual change of structural relaxation from a liquid-like to a solid-like
dynamics around $T_{c}$ is associated with a change of the system’s response
to a mechanical load, in particular with respect to the yielding of the
system. In this work, we have studied sheared supercooled liquids in a planar
Couette flow geometry, applying a constant shear rate $\dot{\gamma}$. We have
shown that the emergence of a transient amorphous solid state implies the
possibility of brittle yielding which is characterized a sharp stress drop in
the stress-strain relation. This means that around a strain of the order of
0.1, the stress shows a sudden decrease on a strain scale
$\delta\tilde{\gamma}^{\star}$ much less than 0.1 (this value is found for the
stress decay at yielding for temperatures above and around $T_{c}$). For
example, at a temperature $T=0.06$ and a shear rate $\dot{\gamma}=10^{-5}$, we
find $\delta\tilde{\gamma}^{\star}=0.014$. While at low temperatures, $T\ll
T_{c}$, $\delta\tilde{\gamma}^{\star}$ significantly decreases with increasing
system size, our data is not conclusive with respect to the question whether
brittle yielding can be interpreted in terms of an underlying kinetic first-
order transition in the limit $\dot{\gamma}\to 0$. Anyway, at finite
temperature, such an interpretation has to be taken with a grain of salt. On
the one hand, the signatures of a first-order transition can be only seen on
the time scale $\tau_{\rm lt}$ and thus for shear rates $\dot{\gamma}$ with
$\dot{\gamma}\tau_{\rm lt}\gtrsim 1$ (note that for $\dot{\gamma}\tau_{\rm
lt}\gg 1$, one expects Newtonian behavior). On the other hand, at a given
temperature $T<T_{c}$, the shear rate has to be small enough that the steady-
state stress is close to the apparent yield stress, as obtained from the
extrapolation to $\dot{\gamma}\to 0$ in terms of a Herschel-Bulkley law. Thus,
the time scale $\tau_{\rm lt}$ has to be very large in order to see the
signatures of a first-order transition and this is the case for temperatures
far below $T_{c}$. As similar interplay of time scales has been recently found
by Shrivastav and Kahl shrivastav2021 , studying the yielding in a cluster
crystal.

Brittle yielding is associated with horizontal or vertical shear bands. Both
types of shear bands are equally efficient to release the stress at the yield
strain. The mechanism for the formation of such shear bands at finite
temperatures and shear rates is still not well understood, but for the
transient amorphous solid states under equilibrium conditions, as studied in
this work, techniques and theoretical frameworks can be adapted that have been
previously mainly used for athermal systems, such as an analysis of soft modes
determined from the dynamical matrix jaiswal2016 ; procaccia2017 ; parisi2017
. Another promising framework to investigate the yielding transition and shear
banding in amorphous solids is the analysis in terms of non-affine
displacements falk1998 ; lemaitre2006 ; ganguly2013 ; ganguly2015 ;
baggioli2021 , as recently applied to elucidate plasticity and yielding in
crystalline solids ganguly2017 ; nath2018 ; reddy2020 . Work in this direction
is in progress.

###### Acknowledgements.

J.H. gratefully acknowledges useful discussions with Parswa Nath and Surajit
Sengupta.

## References

  * (1) W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory (Oxford University Press, Oxford, 2009).
  * (2) J. Solyom, Fundamentals of the Physics of Solids, Volume 1 Structure and Dynamics (Springer, Berlin, 2007).
  * (3) A. Cavagna, Phys. Rep. 476, 51 (2009).
  * (4) G. Szamel and E. Flenner, Phys. Rev. Lett. 107, 105505 (2011).
  * (5) M. Maier, A. Zippelius, and M. Fuchs, Phys. Rev. Lett. 119, 265701 (2017).
  * (6) M. Fuchs and M. E. Cates, Phys. Rev. Lett. 89, 248304 (2002).
  * (7) K. Miyazaki and D. R. Reichman, Phys. Rev. E 66, 050501 (2002).
  * (8) K. Miyazaki, D. R. Reichman, and R. Yamamoto, Phys. Rev. E 70, 011501 (2004).
  * (9) G. Szamel, Phys. Rev. Lett. 93, 178301 (2004).
  * (10) C. P. Amann, M. Siebenbürger, M. Krüger, F. Weysser, M. Ballauff, and M. Fuchs, J. Rheol. 57, 149 (2013).
  * (11) C. P. Amann, M. Siebenbürger, M. Ballauff, and M. Fuchs, J. Phys.: Condens. Matter 27, 194121 (2015).
  * (12) A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X 7, 021039 (2017).
  * (13) T. S. Grigera and G. Parisi, Phys. Rev. E 63, 045102 (2001).
  * (14) C. A. Schuh, T. C. Hufnagel, and U. Ramamurty, Acta Mater. 55, 4067 (2007).
  * (15) M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Proc. Natl. Acad. Sci. USA 115, 6656 (2018).
  * (16) M. Popović, T. W. J. de Geus, and M. Wyart, Phys. Rev. E 98, 040901(R) (2018).
  * (17) H. J. Barlow, J. O. Cochran, and S. M. Fielding, Phys. Rev. Lett. 125, 168003 (2020).
  * (18) R. Besseling, L. Isa, P. Ballesta, G. Petekidis, M. E. Cates, and W. C. K. Poon, Phys. Rev. Lett. 105, 268301 (2010).
  * (19) T. Divoux, D. Tamarii, C. Barentin, and S. Manneville, Phys. Rev. Lett. 104, 208301 (2010).
  * (20) V. Chikkadi, G. Wegdam, D. Bonn, B. Nienhuis, and P. Schall, Phys. Rev. Lett. 107, 198303 (2011).
  * (21) T. Divoux, M. A. Fardin, S. Manneville, and S. Lerouge, Ann. Rev. Fluid Mech. 48, 81 (2016).
  * (22) R. Maaß and J. F. Löffler, Adv. Funct. Mater. 25, 2353 (2015).
  * (23) J. Bokeloh, S. V. Divinski, G. Reglitz, and G. Wilde, Phys. Rev. Lett. 107, 235503 (2011).
  * (24) I. Binkowski, G. P. Shrivastav, J. Horbach, S. V. Divinski, and G. Wilde, Acta Mater. 109, 330 (2016).
  * (25) R. Hubek, S. Hilke, F. A. Davani, M. Golkia, G. P. Shrivastav, S. V. Divinski, H. Rösner, J. Horbach, and G. Wilde, Front. Mater. 7, 144 (2020).
  * (26) F. Varnik, L. Bocquet, J.-L. Barrat, and L. Berthier, Phys. Rev. Lett. 90, 095702 (2003).
  * (27) N. P. Bailey, J. Schiotz, and K. W. Jacobsen, Phys. Rev. B 73, 064108 (2006).
  * (28) Y. Shi and M. L. Falk, Phys. Rev. B 73, 214201 (2006).
  * (29) Y. Shi, M. B. Katz, H. Li, and M. L. Falk, Phys. Rev. Lett. 98, 185505 (2007).
  * (30) Y. Ritter and K. Albe, Acta Mater. 59, 7082 (2011).
  * (31) D. Sopu, Y. Ritter, H. Gleiter, and K. Albe, Phys. Rev. B 83, 100202(R) (2011).
  * (32) P. Chaudhuri, L. Berthier, and L. Bocquet, Phys. Rev. E 85, 021503 (2012).
  * (33) R. Dasgupta, H. G. E. Hentschel, and I. Procaccia, Phys. Rev. Lett. 109, 255502 (2012).
  * (34) R. Dasgupta, O. Gendelman, P. Mishra, I. Procaccia, and C. A. B. Z. Shor, Phys. Rev. E 88, 032401 (2013).
  * (35) K. Albe, Y. Ritter, and D. Sopu, Mech. Mater. 67, 94 (2013).
  * (36) G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 94, 042605 (2016).
  * (37) M. Golkia, G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 102, 023002 (2020).
  * (38) M. Singh, M. Ozawa, and L. Berthier, Phys. Rev. Mater. 4, 025603 (2020).
  * (39) A. D. S. Parmar, S. Kumar, and S. Sastry, Phys. Rev. X 9, 021018 (2019).
  * (40) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * (41) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 2nd edition (Oxford University Press, Oxford, 2017).
  * (42) T. Soddemann, B. Dünweg, and K. Kremer Phys. Rev. E 68, 046702 (2003).
  * (43) A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972).
  * (44) Y. Rosenfeld and P. Tarazona, Mol. Phys. 95, 141 (1998).
  * (45) T. S. Ingebrigtsen, A. A. Veldhorst, T. B. Schroder, and J. C. Dyre, J. Chem. Phys. 139, 171101 (2013).
  * (46) J. Horbach, W. Kob, K. Binder, and C. A. Angell, Phys. Rev. E 54, R5897 (1996).
  * (47) J. Horbach, W. Kob, and K. Binder, Eur. Phys. J. B 19, 531 (2001).
  * (48) J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rd ed. (Academic Press, 2006).
  * (49) A. L. Thorneywork, D. G. A. L. Aarts, J. Horbach, and R. P. A. Dullens, Soft Matter 12, 4129 (2016).
  * (50) M. Fuchs, W. Götze, and M. R. Mayr, Phys. Rev. E 58, 3384 (1998).
  * (51) W. H. Herschel and R. Bulkley, Kolloid Z. 39, 291 (1926),
  * (52) G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 94, 042605 (2016).
  * (53) H. Risken, The Fokker-Planck Equation (Springer, Berlin, 1989).
  * (54) G. P. Shrivastav and G. Kahl, Soft Matter 17, 8536 (2021).
  * (55) P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Phys. Rev. Lett. 116, 085501 (2016).
  * (56) I. Procaccia, C. Rainone, and M. Singh, Phys. Rev. E 96, 032907 (2017).
  * (57) G. Parisi, I. Procaccia, C. Rainone, and M. Singh, Proc. Natl. Acad. Sci. USA 114, 5577 (2017).
  * (58) M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
  * (59) A. Lemaitre and C. Maloney, J. Stat. Phys. 123, 415 (2006).
  * (60) S. Ganguly, S. Sengupta, P. Sollich, and M. Rao, Phys. Rev. E 87, 042801 (2013).
  * (61) S. Ganguly, S. Sengupta, and P. Sollich, Soft Matter 11, 4517 (2015).
  * (62) M. Baggioli, I. Kriuchevskyi, T. W. Sirk, and A. Zaccone, Phys. Rev. Lett. 127, 015501 (2021).
  * (63) S. Ganguly, P. Nath, J. Horbach, P. Sollich, and S. Karmakar, and S. Sengupta, J. Chem. Phys. 146, 124501 (2017).
  * (64) P. Nath, S. Ganguly, J. Horbach, P. Sollich, S. Karmakar, and S. Sengupta, Proc. Natl. Acad. Sci. USA 115, E4322 (2018).
  * (65) V. S. Reddy, P. Nath, J. Horbach, P. Sollich, and S. Sengupta, Phys. Rev. Lett. 124, 025503 (2020).


2.5cm2.5cm2.5cm2.5cm

# Dark energy constraints from quasar observations

B. Czerny  M.L. Martínez-Aldama  Center for Theoretical Physics, Polish
Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  G. Wojtkowska
Warsaw University Observatory, Al. Ujazdowskie 4, 00-478 Warszawa, Poland  M.
Zajaček  Center for Theoretical Physics, Polish Academy of Sciences, Al.
Lotników 32/46, 02-668 Warsaw, Poland  P. Marziani  INAF, Osservatorio
Astronomico di Padova, Italy  D. Dultzin  Instituto de Astronomía, UNAM,
Mexico  M. H. Naddaf  Center for Theoretical Physics, Polish Academy of
Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  S. Panda  Center for
Theoretical Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668
Warsaw, Poland  R. Prince  Center for Theoretical Physics, Polish Academy of
Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  R. Przyluski  Space
Research Centre, Polish Academy of Sciences, Bartycka 18a, 00-716 Warsaw,
Poland  M. Ralowski  Astronomical Observatory of the Jagiellonian University,
Orla 171, 30–001 Krakow, Poland  M. Śniegowska  Center for Theoretical
Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland

(Key words: cosmology, dark energy, quasars)

Recent measurements of the parameters of the Concordance Cosmology Model
($\Lambda$CDM) done in the low-redshift Universe with Supernovae Ia/Cepheids,
and in the distant Universe done with Cosmic Microwave Background (CMB) imply
different values for the Hubble constant (67.4 $\pm$ 0.5 [km s-1 Mpc-1] from
Planck vs 74.03 $\pm$ 1.42 [km s-1 Mpc-1] Riess et al. 2019). This Hubble
constant tension implies that either the systematic errors are underestimated,
or the $\Lambda$CDM does not represent well the observed expansion of the
Universe. Since quasars - active galactic nuclei - can be observed in the
nearby Universe up to redshift z $\sim$ 7.5, they are suitable to estimate the
cosmological properties in a large redshift range. Our group develops two
methods based on the observations of quasars in the late Universe up to
redshift z$\sim$4.5, with the objective to determine the expansion rate of the
Universe. These methods do not yet provide an independent measurement of the
Hubble constant since they do not have firm absolute calibration but they
allow to test the $\Lambda$CDM model, and so far no departures from this model
were found.

##  1 Introduction

The cosmological parameters can be estimated from different sets of data at
various redshifts, but if the standard $\Lambda$CDM (Lambda-Cold Dark Matter)
model is valid, they can always be represented by the current (zero redshift)
values. The final results from the Planck mission, based on the analysis of
the Cosmic Microwave Background (CMB) do not indicate any tension with the
standard model, and give the value of the $\Omega_{m}=0.315\pm 0.007$ [1] and
the Hubble constant $H_{0}=67.4\pm 0.5$ [km s-1 Mpc-1]. Many of the
measurements done in the local Universe ($z<10$) are in significant
disagreement with these $\Omega_{m}$ or $H_{0}$ values (e.g. [2]) while other
measurements, also local, are still roughly in agreement with the results from
Planck (e.g. the last results from gravitational waves [3]).

Therefore various probes and methods are needed to confirm, or to reject, the
hypothesis that the $\Lambda$CDM model does not well describe the Universe,
and the evolving dark energy is needed instead of the cosmological constant.
Quasars (QSO) are very attractive cosmological probes, since they cover a
broad range of redshifts, from nearby sources (referred to as Active Galactic
Nuclei, AGN) to most distant objects at redshift above 7 [4, 5]. They also do
not show significant evolution with redshift [6].

##  2 Two methods for using quasars in cosmology

We are currently using two methods of turning quasars into standardizable
candles. The first method is based on the radius-luminosity relation and the
second method is based on super-Eddington sources.

###  2.1 Method based on radius-luminosity relation for the BLR

The reverberation mapping technique is based on the long-term monitoring of a
source in order to determine the time response ($\tau_{\rm BLR}$) of the
emission line to the continuum variations [7]. The most important result from
the reverberation mapping studies is the correlation between the continuum
luminosity ($L$) and the distance ($R\mathrm{{}_{\mathrm{BLR}}}$) where the
emission line is emitted in the broad line region (BLR). This relation is
known as the Radius-Luminosity relation (RL) and it is given approximately by
$R\mathrm{{}_{\mathrm{BLR}}}$$\propto L^{0.5}$. The reverberation mapping
studies require extensive use of telescope time to achieve high quality
results, thus only $\sim 120$ sources have been analyzed with this technique
until date. Most of the monitoring are based on the optical H$\beta$ for low-
redshift sources, while for high redshift regimes, due to the Doppler shift,
the monitoring are focused on the UV emission lines such as Mg ii
$\lambda$2800, C iv $\lambda$1549 and C iii] $\lambda$1909.

For many years the RL relation showed a low scatter ($\sigma_{\rm rms}\sim
0.13$ dex) [8], which ensured its use in the determination of the black hole
mass ($M\mathrm{{}_{BH}}$). However, the inclusion of new sources,
particularly those radiating close to the Eddington limit (high accretion
rates), has led to a much larger scatter, clearly related with the accretion
rate [9, 10]. Some corrections based on the accretion rate [11] and
independent parameters such as the Fe ii strength or the amplitude of
variability [12, 13, 14] (in turn correlated with the accretion rate) have
been proposed to correct this effect, allowing to reduce the scatter.

Figure 1: Top panel: Quasar Hubble diagram using the reverberation mapped
sources. Blue circles and red diamonds correspond to the H$\beta$ and Mg ii
$\lambda$2800 sources, respectively. Black line marks the $\Lambda$CDM model
for flat cosmology, $\Omega_{m}=0.297$, $H_{0}=67.5$ [km s-1 Mpc-1]. Bottom
subpanel shows the residuals. Bottom panel: Confidence contours at 68$\%$
(cyan) and 95$\%$ (blue) for $\Omega_{m}$ and $\Omega_{\Lambda}$ for general
$\Lambda$CDM model based on $\chi^{2}$ fitting where the best $\Omega_{m}$ and
$\Omega_{\Lambda}$ are represented by the yellow symbol.

Besides, the Radius-Luminosity relation offers the possibility to determine
the luminosity independently of the redshift [15, 16], and to estimate the
cosmological parameters. However, in [11] the errors for $\Omega_{m}$ and
$\Omega_{\Lambda}$ based on available H$\beta$ were still very large, despite
the corrections for the accretion rate. In the present paper, we make the
following important modifications. First, we combine the previous sample with
Mg ii $\lambda$2800 reverberation-mapped sources. The new sample includes two
sources monitored by us with The Southern African Large Telescope (SALT)
during 6 years [17, 18]. Next, since the errors of the time delay measurement
are highly asymmetric, we use the $\chi^{2}$ statistics to determine the
cosmological parameters, we use the method of [19] instead of a simple
symmetrization of the errors.

We also modified the approach to the R-L relation in case of the H$\beta$
sample which is extremely heterogeneous. We treated the coefficients of this
relation as arbitrary, and minimized the total $\chi^{2}$ fit to a flat
cosmological model, with Hubble constant fixed at 67.5 [km s-1 Mpc-1]. For Mg
ii $\lambda$2800 sample, we used the R-L parametrization given by Eq. 9 in
[18] since this sample was analysed in a more uniform way. While fitting the
flat cosmology model (both samples combined), we applied the sigma-clipping
approach, and we removed the sources which showed a departure by more than 3
sigma from the best-fit. All the removed sources were from the H$\beta$
sample: Mrk 493, J074352.02+271239.5, Mkn 509, MCG+08-11-011, J142103,
J142043, J141123, and J142052. Thus our total sample has now 120 objects. We
then refitted the Hubble diagram. The best-fit returned the best R-L
parametrization of H$\beta$ sample as $\log L_{5100}=1.489\log\tau_{\rm
corr}-2.222$, where $\tau_{\rm corr}$ is the time delay corrected by the
accretion rate effect [11]. For the flat cosmology, we obtained the best-fit
value $\Omega_{m}=0.297^{+0.060}_{-0.054}$ (see Fig. 1, top panel). This value
is fully consistent with the value $0.3153\pm 0.0073$ from Planck [1] for flat
cosmology, and the error in our new result is much smaller than obtained by
[11], although still much larger than from Planck. This illustrates that the
method could be powerful, if the sample is more uniformly analysed from the
very beginning. If we do not assume a flat cosmology, 2-D contour errors are
still large (see Fig. 1, bottom panel), although considerably smaller than in
[11]. The best-fit Planck values are well within the 1$\sigma$ error, so we do
not see any tension with the results based on Cosmic Microwave Background.

###  2.2 Method based on Super-Eddington sources

Quasars radiating close to the Eddington limit are known as xA-QSO or
super–Eddington sources [20, 21]. This QSO population shows peculiar spectral
and photometric properties, which differentiate them from the rest of the QSO
population and make them easy to identify in catalogs like SDSS or the
upcoming Vera Rubin Observatory’s Legacy Survey of Space and Time (LSST). In
the optical range, they are the strongest Fe ii emitters and do not show a
strong contribution of narrow emission lines such as [O iii]
$\lambda\lambda$4959, 5007 [22]. xA-QSO show the strongest outflows in the
high ionization lines mostly observed in the UV emission lines like C iv
$\lambda$1549 or S iv $\lambda$1397 [23], although in the most extreme cases
the strong radiation forces provoke the presence of outflows in low-ionization
lines such as H$\beta$ or Al iii $\lambda$1860\. According to the
photoionization models, their broad line regions show large densities
(nH$=10^{12-13}$ [cm-3]), low-ionization parameters (log $U<-2$) and high
metallicities ($Z\sim 10Z_{\odot}$) [24, 25, 26, 27]. In addition, Super-
Eddington sources also show remarkably low optical variability and time delays
shorter than the predicted by the RL relation [28]. The UV flux ratios Al iii
$\lambda$1860/Si iii] $\lambda$1892$>$0.5 and C iii] $\lambda$1909/Si iii]
$\lambda$1892$<$1.0 have shown a high effectiveness as selection criteria to
identify xA-QSO sources at high-redshift, while at-low redshift xA-QSO
typically show Fe ii/H$\beta$$>$1.0 [20].

In xA-QSO, despite the rise of the accretion rate, the luminosity saturates
toward a limiting value, since the accretion efficiency decreases. Thus the
ratio luminosity–black hole mass ($L$/$M\mathrm{{}_{BH}}$) does not change and
they can be considered as “Eddington standard candles”. A similarity in the
physical conditions (density, ionization parameter, metallicity) of the BLR is
expected because them belonging to the same population, therefore a
generalization of all of them can be considered [29]. Since the low-ionization
lines are less affected by the strong radiation forces, emission lines like
H$\beta$ and Al iii $\lambda$1860 are excellent candidates for virial
estimators. Based on these assumptions, it is possible to determine the
luminosity distances independently of redshift and get an estimation of the
matter and energy content in the Universe.

The previous xA-QSO sample included $\sim 200$ objects at redshift $z<2.7$
[30]. In order to increase the redshift range, we considered the most recent
edition of the The Sloan Digital Sky Survey Reverberation Mapping (SDSS-RM)
catalog [31]. This catalog includes the automatic measurements of the most
important UV emission lines for 549 sources with 44.1$<{\rm log}L_{1700}<$46.9
[ers s-1] at $0.9<z<4.3$, where $\sim 20\%$ show high accretion rates, so they
can be considered as Super-Eddington candidates. In the xA sources the Fe
iii$\lambda 1914$ shows an important contribution, hence a good deblending of
the C iii] $\lambda$1909 and Fe iii $\lambda 1914$ is required. Unfortunately,
the SDSS catalog does not include Fe iii in their multicomponent fittings, so
not all the sources satisfy the selection criteria to identify them as the xA.
So for the first test, we select the sources based on the Eddington ratio
(L/LEdd$>0.2$) estimated from the Al iii $\lambda$1860 based black hole mass.
After removal of some objects with extreme FWHMAlIII values, our final sample
includes 88 objects at $1<z<4.5$.

Figure 2: Quasar Hubble diagram using the Super-Eddington sources. Diamonds
correspond to the measurements from the H$\beta$, while circles belong to UV
Al iii $\lambda$1860 emission line. Purple, gray, yellow, green and blue
symbols correspond to the xA sources from the SEAMBH project [12], [22], [20],
[25] and SDSS-RM [31] samples, respectively. The best-fit line shows the flat
model, with $H_{0}$ from Planck, and best fit $\Omega_{m}=0.290$.

In order to determine the cosmological constant $\Omega_{m}$ and
$\Omega_{\Lambda}$, we combine the previous xA samples such as the super-
Eddington sources from the super-Eddington accreting massive black holes
(SEAMBHs) project with the most recent measurements from [25]. The quasar
Hubble diagram with the super-Eddington sources is shown in Fig. 2. We adopt
the scaling of the virial estimator to be consistent with Planck $H_{0}$
value, and we assume the flat cosmology. In this case we obtain
$\Omega_{m}=0.290^{+0.048}_{-0.043}$, fully consistent with the Planck results
despite the fact that quasars cover the redshift range from nearby sources to
almost 4.5.

##  3 Discussion

Here we presented the most recent results based on the two methods for
applying quasars to constrain the expansion rate of the Universe. Our methods,
as for now, are not based on absolute scaling so they cannot predict the value
of $H_{0}$. In principle, such an absolute scaling can be achieved. For method
(i) it would require an independent measurement of the dust temperature at the
BLR onset, and the development of a 3-D BLR model, which is in progress (see
e.g. [32]). For method (ii), we would need an absolute scaling of the radius-
luminosity relation, also establishing the mean density of the BLR (see e.g.
[20, 29]). At this stage, we fixed the value of the Hubble constant at the
Planck value and tested, whether the redshift dependence of the luminosity
distance is consistent with the standard $\Lambda$CDM model, and whether the
remaining cosmological parameters derived from quasar data are consistent with
Planck values.

So far, within the available accuracy, our values of $\Omega_{m}$ are fully
consistent with the Planck value for the flat Universe despite the fact that
second method extends up to the redshift 4.5. Thus we do not support the claim
of the tension with the standard model based on Supernovae Ia with absolute
scaling in turn predominantly based on Cepheid stars [33]. Our results from
method (i) are consistent with the tension found by [34] since they claim to
see departures only above the redshift 1.5 - 2, and method (i) does not go
this far. As for the method (ii), we have many sources up to redshift 2.5, but
indeed very few above 2.5, and our method of analysis was not yet optimized by
an outlier removal through sigma-clipping. Further studies are clearly needed
for this method, both with the current data and eventually by increasing the
number of high redshift quasars.

## Acknowledgement

The project is partially based on observations made with the SALT under
programs 2012-2-POL-003, 2013-1-POL-RSA-002, 2013-2-POL-RSA-001, 2014-1-POL-
RSA-001, 2014-2-SCI-004, 2015-1-SCI-006, 2015-2-SCI-017, 2016-1-SCI-011,
2016-2-SCI-024, 2017-1-SCI-009, 2017-2-SCI-033, 2018-1-MLT-004 (PI: B.
Czerny). The authors acknowledge the financial support by the National Science
Centre, Poland, grant No. 2017/26/A/ST9/00756 (Maestro 9), and by the Ministry
of Science and Higher Education (MNiSW) grant DIR/WK/2018/12. The Polish
participation in SALT is funded by grant No. MNiSW DIR/WK/2016/07.

## References

  * [1] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, and et al.  Planck 2018 results. VI. Cosmological parameters.  $\rm A\&A$, 641:A6, September 2020, doi:10.1051/0004-6361/201833910.
  * [2] Adam G. Riess.  The expansion of the Universe is faster than expected.  Nature Reviews Physics, 2(1):10–12, December 2019, doi:10.1038/s42254-019-0137-0.
  * [3] The LIGO Scientific Collaboration, the Virgo Collaboration, B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, and et al.  A gravitational-wave measurement of the Hubble constant following the second observing run of Advanced LIGO and Virgo.  arXiv e-prints, page arXiv:1908.06060, August 2019.
  * [4] D. J. Mortlock, S. J. Warren, B. P. Venemans, M. Patel, P. C. Hewett, . G. McMahon, C. Simpson, and et al.  A luminous quasar at a redshift of z = 7.085.  $\rm Nature$, 474(7353):616–619, June 2011, doi:10.1038/nature10159.
  * [5] E. Bañados, B. P. Venemans, C. Mazzucchelli, E. P. Farina, F. Walter, F. Wang, R. Decarli, and et al.  An 800-million-solar-mass black hole in a significantly neutral Universe at a redshift of 7.5.  $\rm Nature$, 553(7689):473–476, January 2018, doi:10.1038/nature25180.
  * [6] M. Onoue, E. Bañados, C. Mazzucchelli, B. P. Venemans, J.-T. Schindler, F. Walter, J. F. Hennawi, and et al.  No Redshift Evolution in the Broad-line-region Metallicity up to z = 7.54: Deep Near-infrared Spectroscopy of ULAS J1342+0928.  $\rm ApJ$, 898(2):105, August 2020, doi:10.3847/1538-4357/aba193.
  * [7] B. M. Peterson, L. Ferrarese, K. M. Gilbert, S. Kaspi, M. A. Malkan, D. Maoz, D. Merritt, and et al.  Central Masses and Broad-Line Region Sizes of Active Galactic Nuclei. II. A Homogeneous Analysis of a Large Reverberation-Mapping Database.  $\rm ApJ$, 613:682–699, October 2004.
  * [8] E. Kilerci Eser, M. Vestergaard, B. M. Peterson, K. D. Denney, and M. C. Bentz.  On the Scatter in the Radius-Luminosity Relationship for Active Galactic Nuclei.  $\rm ApJ$, 801:8, March 2015, doi:10.1088/0004-637X/801/1/8.
  * [9] P. Du, C. Hu, K.-X. Lu, Y.-K. Huang, C. Cheng, J. Qiu, and SEAMBH Collaboration.  Supermassive Black Holes with High Accretion Rates in Active Galactic Nuclei. IV. H$\beta$ Time Lags and Implications for Super-Eddington Accretion.  $\rm ApJ$, 806:22, June 2015, doi:10.1088/0004-637X/806/1/22.
  * [10] P. Du, Z.-X. Zhang, K. Wang, Y.-K. Huang, Y. Zhang, K.-X. Lu, and SEAMBH Collaboration.  Supermassive Black Holes with High Accretion Rates in Active Galactic Nuclei. IX. 10 New Observations of Reverberation Mapping and Shortened H$\beta$ Lags.  $\rm ApJ$, 856:6, March 2018, doi:10.3847/1538-4357/aaae6b.
  * [11] M. L. Martínez-Aldama, B. Czerny, D. Kawka, V. Karas, S. Panda, M. Zajaček, and P. T. Życki.  Can Reverberation-measured Quasars Be Used for Cosmology?  $\rm ApJ$, 883(2):170, October 2019, doi:10.3847/1538-4357/ab3728.
  * [12] P. Du and J.-M. Wang.  The Radius-Luminosity Relationship Depends on Optical Spectra in Active Galactic Nuclei.  $\rm ApJ$, 886(1):42, November 2019, doi:10.3847/1538-4357/ab4908.
  * [13] E. Dalla Bontà, B. M. Peterson, M. C. Bentz, W. N. Brandt, S. Ciroi, G. De Rosa, and et al.  The Sloan Digital Sky Survey Reverberation Mapping Project: Estimating Masses of Black Holes in Quasars with Single-Epoch Spectroscopy.  arXiv e-prints, page arXiv:2007.02963, July 2020.
  * [14] M. L. Martínez-Aldama, M. Zajaček, B. Czerny, and S. Panda.  Scatter Analysis along the Multidimensional Radius-Luminosity Relations for Reverberation-mapped Mg II Sources.  $\rm ApJ$, 903(2):86, November 2020, doi:10.3847/1538-4357/abb6f8.
  * [15] D. Watson, K. D. Denney, M. Vestergaard, and T. M. Davis.  A New Cosmological Distance Measure Using Active Galactic Nuclei.  $\rm ApJL$, 740:L49, October 2011, doi:10.1088/2041-8205/740/2/L49.
  * [16] M. Haas, R. Chini, M. Ramolla, F. Pozo Nuñez, C. Westhues, and R. et al. Watermann.  Photometric AGN reverberation mapping - an efficient tool for BLR sizes, black hole masses, and host-subtracted AGN luminosities.  $\rm A\&A$, 535:A73, November 2011, doi:10.1051/0004-6361/201117325.
  * [17] B. Czerny, J.-M. Wang, P. Du, K. Hryniewicz, V. Karas, Y.-R. Li, S. Panda, M. Sniegowska, C. Wildy, and Y.-F. Yuan.  Interpretation of Departure from the Broad-line Region Scaling in Active Galactic Nuclei.  $\rm ApJ$, 870:84, January 2019, doi:10.3847/1538-4357/aaf396.
  * [18] M. Zajaček, B. Czerny, M. L. Martinez-Aldama, M. Rałowski, A. Olejak, S. Panda, K. Hryniewicz, and et al.  Time-delay Measurement of Mg II Broad-line Response for the Highly Accreting Quasar HE 0413-4031: Implications for the Mg II-based Radius-Luminosity Relation.  $\rm ApJ$, 896(2):146, June 2020, doi:10.3847/1538-4357/ab94ae.
  * [19] R. Barlow.  Asymmetric Statistical Errors.  arXiv e-prints, page physics/0406120, June 2004.
  * [20] P. Marziani and J. W. Sulentic.  Highly accreting quasars: sample definition and possible cosmological implications.  $\rm MNRAS$, 442(2):1211–1229, August 2014, doi:10.1093/mnras/stu951.
  * [21] J.-M. Wang, P. Du, C. Hu, H. Netzer, J.-M. Bai, K.-X. Lu, S. Kaspi, and SEAMBH Collaboration.  Supermassive Black Holes with High Accretion Rates in Active Galactic Nuclei. II. The Most Luminous Standard Candles in the Universe.  $\rm ApJ$, 793:108, October 2014, doi:10.1088/0004-637X/793/2/108.
  * [22] C. A. Negrete, D. Dultzin, P. Marziani, D. Esparza, J. W. Sulentic, A. del Olmo, M. L. Martínez-Aldama, and et al.  Highly accreting quasars: The SDSS low-redshift catalog.  $\rm A\&A$, 620:A118, Dec 2018, doi:10.1051/0004-6361/201833285.
  * [23] J. W. Sulentic, A. del Olmo, P. Marziani, M. A. Martínez-Carballo, M. D’Onofrio, D. Dultzin, J. Perea, and et al.  What does CIV$\lambda$1549 tell us about the physical driver of the Eigenvector quasar sequence?  $\rm A\&A$, 608:A122, December 2017, doi:10.1051/0004-6361/201630309.
  * [24] M. L. Martínez-Aldama, A. del Olmo, P. Marziani, J. W. Sulentic, C. A. Negrete, D. Dultzin, and et al.  Extreme quasars at high redshift.  $\rm A\&A$, 618:A179, November 2018, doi:10.1051/0004-6361/201833541.
  * [25] M. Śniegowska, P. Marziani, B. Czerny, S. Panda, M. L. Martínez-Aldama, and et al.  High metal content of highly accreting quasars.  arXiv e-prints, page arXiv:2009.14177, September 2020.
  * [26] S. Panda, M. L. Martínez-Aldama, M. Marinello, B. Czerny, P. Marziani, and D. Dultzin.  The CaFe Project: Optical Fe II and Near-infrared Ca II Triplet Emission in Active Galaxies. I. Photoionization Modeling.  $\rm ApJ$, 902(1):76, October 2020, doi:10.3847/1538-4357/abb5b8.
  * [27] S. Panda.  The CaFe Project: Optical Fe II and Near-Infrared Ca II triplet emission in active galaxies: (II) synthetic EWs, co-dependence between cloud sizes and metal content.  arXiv e-prints, page arXiv:2004.13113, April 2020.
  * [28] P. Du, K.-X. Lu, Z.-X. Zhang, Y.-K. Huang, K. Wang, C. Hu, J. Qiu, and SEAMBH Collaboration.  Supermassive Black Holes with High Accretion Rates in Active Galactic Nuclei. V. A New Size-Luminosity Scaling Relation for the Broad-line Region.  $\rm ApJ$, 825:126, July 2016, doi:10.3847/0004-637X/825/2/126.
  * [29] S. Panda, P. Marziani, and B. Czerny.  The Quasar Main Sequence Explained by the Combination of Eddington Ratio, Metallicity, and Orientation.  $\rm ApJ$, 882(2):79, Sep 2019, doi:10.3847/1538-4357/ab3292.
  * [30] D. Dultzin, P. Marziani, J. A. de Diego, C. A. Negrete, A. Del Olmo, M. L. Martínez-Aldama, and et al.  Extreme quasars as distance indicators in cosmology.  Frontiers in Astronomy and Space Sciences, 6:80, January 2020, doi:10.3389/fspas.2019.00080.
  * [31] Y. Shen, P. B. Hall, K. Horne, G. Zhu, I. McGreer, T. Simm, and et al.  The Sloan Digital Sky Survey Reverberation Mapping Project: Sample Characterization.  $\rm ApJS$, 241(2):34, April 2019, doi:10.3847/1538-4365/ab074f.
  * [32] M.-H. Naddaf, B. Czerny, and R. Szczerba.  BLR size in Realistic FRADO Model.  Frontiers in Astronomy and Space Sciences, 7:15, April 2020, doi:10.3389/fspas.2020.00015.
  * [33] A. G. Riess, S. Casertano, W. Yuan, L. Macri, B. Bucciarelli, M. G. Lattanzi, J. W. MacKenty, and et al.  Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant.  $\rm ApJ$, 861:126, July 2018, doi:10.3847/1538-4357/aac82e.
  * [34] E. Lusso, G. Risaliti, E. Nardini, G. Bargiacchi, M. Benetti, S. Bisogni, S. Capozziello, and et al.  Quasars as standard candles. III. Validation of a new sample for cosmological studies.  $\rm A\&A$, 642:A150, October 2020, doi:10.1051/0004-6361/202038899.