# Weight update in ferroelectric memristors with identical and non-identical
pulses

Suzanne Lancaster<EMAIL_ADDRESS>[  Maximilien Remillieux  [
Moritz Engl  [  Viktor Havel  [  Cláudia Silva  [  Xuetao Wang  [  Thomas
Mikolajick  [  Stefan Slesazeck  [

###### Abstract

Ferroelectric tunnel junctions (FTJs) are a class of memristor which promise
low-power, scalable, field-driven analog operation. In order to harness their
full potential, operation with identical pulses is targeted. In this paper,
several weight update schemes for FTJs are investigated, using either non-
identical or identical pulses, and with time delays between the pulses ranging
from 1 $\mu$s to 10 s. Experimentally, a method for achieving non-linear
weight update with identical pulses at long programming delays is demonstrated
by limiting the switching current via a series resistor. Simulations show that
this concept can be expanded to achieve weight update in a 1T1C cell by
limiting the switching current through a transistor operating in sub-threshold
or saturation mode. This leads to a maximum linearity in the weight update of
87%. The scaling behaviour of this scheme as devices are reduced to the sub-
micron range is investigated, and the linearity slightly reduces in scaled
devices to a maximum of 81%. The origin of the remaining non-linearity is
discussed, as well as how this can be overcome.

NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a, 01187 Dresden, Germany INL] École
Supérieure de Chimie Physique Électronique de Lyon, 3 Rue Victor Grignard,
69100 Villeurbanne, France [NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a, 01187
Dresden, Germany NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a, 01187 Dresden,
Germany NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a, 01187 Dresden, Germany
NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a, 01187 Dresden, Germany NaMLab]
NaMLab gGmbH, Nöthnitzer Str. 64a, 01187 Dresden, Germany NaMLab] NaMLab
gGmbH, Nöthnitzer Str. 64a, 01187 Dresden, Germany [IHM]Institute of
Semiconductors and Microsystems, Technische Universität Dresden, Nöthnitzer
Str 64., 01187 Dresden, Germany NaMLab] NaMLab gGmbH, Nöthnitzer Str. 64a,
01187 Dresden, Germany

##  1 Introduction

With the discovery of ferroelectricity in HfO2 1, the possibility of
integrating ferroelectric devices with CMOS circuitry became much more
promising. Ferroelectric devices based on HfO2 demonstrate fast switching
speeds 2, 3, power efficiency 4 and dense integration suitability 5. As such,
they show great promise for applications such as neuromorphic 6, 7 or edge
computing 8, and logic-in-memory 9, 10. For these use cases, it is
particularly interesting to consider multilevel memristive device operation,
which enables synaptic operation 11 and increases bit density 12. In
particular, linear weight update schemes are desirable in order to improve
accuracy in neuromorphic processors 13. One of the main hindrances here is the
application of switching pulses on-chip: in order to set a device into states
of different levels, non-identical pulses are generally required 14.

In the most ideal case, a ferroelectric device should update its weight when
subjected to identical switching pulses. This would facilitate synaptic
behaviour 15 and simplify programming schemes in crossbar configurations 16.
Due to the field-driven switching of ferroelectrics, weight update is
classically unobtainable in this manner. A certain field will be applied
across the ferroelectric on a given pulse, and all domains with a coercive
field less than the applied field at the pulse frequency will switch, so that
if backswitching is neglected, there will be no switchable domains when
applying the next pulse. Polarization switching in hafnia ferroelectrics is
understood to be governed by nucleation-limited switching 17, i.e. domain
nucleation dominates the switching and is limited by the field required to
nucleate oppositely-polarized domains. For a longer switching time, the
required voltage for domain nucleation is lower; conversely, for a higher
voltage, domain nucleation will happen within a shorter pulse width. This
leads to the well-known time-voltage tradeoff in ferroelectrics 18 and
indicates that increasing the voltage or pulse width is necessary for the
nucleation of additional domains in multilevel programming schemes. In
practice, diverse observations in the literature suggest that weight update
with identical pulses can be achieved through extrinsic, rather than intrinsic
effects 19, 20, 21.

In Ferroelectric Field Effect Transistors (FeFETs), an accumulative weight
update effect has been observed with identical pulses 22, while later
experiments on write disturb indicated that this is not intrinsic to
ferroelectric hafnia but may be stack-dependent 23. Similarly in Ferroelectric
Tunnel Junctions (FTJs), it was shown that weight update can be achieved with
identical pulses at a time delay between switching pulses of 1$\mu$s 19, but
weight update greatly diminishes with increasing time delay between identical
pulses 20. The stronger weight update at short time delays is proposed to
originate from extrinsic effects such as interfacial charge trapping 21,
residual fields when switching with fast pulses 20, or ionic motion in metal-
oxide electrodes. While it is possible to design circuitry to set multilevel
states with non-identical pulses 24, this adds additional constraints in terms
of chip area, energy, and programming algorithms 25. In that respect,
multilevel operation using identical pulses is a much more favourable
approach.

In this article, we will present different operating schemes for ferroelectric
devices, using an FTJ as a model device. FTJs can be considered as memristors
since the polarization modulates the device resistance and thus the current
density at an applied read voltage 6, 26, 27. First, we will revisit switching
with non-identical pulses as well as published approaches to switching with
identical pulses and discuss their limitations, i.e. the necessity of strict
timing conditions. Then, we will demonstrate a novel approach which allows
weight update in ferroelectric devices using identical pulses on arbitrary
time scales, when a resistor of a sufficiently high resistance is connected in
series during the switching pulse. In this case, weight update is achieved
with identical pulses, although it is non-linear, and the full memory window
(MW) is limited by the RC constant of the system. A final scheme is proposed
and demonstrated via simulations, where a weight update with identical pulses
is performed while applying a current compliance to the switching circuit,
which should allow a linear weight update spanning the maximum MW. These
schemes greatly simplify the on-chip weight update of ferroelectric devices.

##  2 Weight update with non-identical pulses

###  2.1 Characterization

Figure 1: (a) Current-voltage and (b) polarization-voltage curves (orange
lines) for the 110 nm diameter FTJ devices used in these experiments, measured
at 20 kHz. Purple diamonds represent the current density measured with pulses
of 50 ms, 1.3 V, for pulses of 10 $\mu$s at the voltages shown on the x-axis.
The devices were woken-up at 100 kHz, 4.5 V, 1e4 cycles.

In these experiments, bilayer FTJs were used, with either 11 nm ferroelectric
Hf0.5Zr${0.5}$O2 (HZO) and a 1 nm tunnel barrier (FTJ-1nm) to minimize
switching and read voltages 28 or 12 nm HZO and a 2 nm barrier 27 (FTJ-2nm).
The full description of sample processing is given in the methods section.
Figure 1a shows the I-V at 20 kHz (after wake-up at 4.5 V, 1e4 cycles, 100 kHz
square pulses) curves measured during switching on an FTJ-1nm device. A higher
background current is measured on the biased electrode due to an additional
capacitive charging of the electrical setup, which needs to be considered when
modelling the switching behaviour. Figure 1b shows the P-V loop extracted from
(a) using a PUND technique (orange line, extracted from the I-V of the
grounded electrode). The purple diamonds represent the FTJ current density
measured with rectangular read pulses of 1.3 V, 50 ms after programming at the
indicated voltages with 10 $\mu$s programming pulses. It is clear that the FTJ
read-out current follows the trend of the P-V loop, which indicates that the
read-out current and polarization follow a linear relationship, as shown
previously for FTJ-2nm devices 19.

###  2.2 Varying pulse width and time

The understanding of how a ferroelectric reacts to increasing or decreasing
voltages or pulse widths is paramount to harnessing its synaptic function 27.
A first approximation for appropriate partial switching pulses can be found
from switching kinetics measurements, where the polarization switched under
pulses of a given width and amplitude can be found. Then, one parameter (pulse
width/pulse amplitude) can be fixed and the switching behaviour can be
explored when varying the free parameter. Examples are plotted here for fixing
the pulse width and increasing the voltage (fig 2a&b) or fixing the voltage
and increasing the pulse width (2c&d).

Figure 2: (a), (b) Current density vs. applied voltage for pulsed weight
update on an FTJ device using positive or negative pulses of increasing
voltage magnitude and 10 $\mu$s pulse width. (c), (d) Current density vs.
pulse time for pulsed weight update of an FTJ device using positive or
negative pulses of increasing time width and fixed voltage amplitudes - 2.5 V
(squares) or 3 V (stars). The FTJ state is read in between each applied pulse,
and the pulse train is shown inset.

In addition, fig 2c&d compare the switching behaviour at two fixed voltage
amplitudes - a larger fixed voltage can be applied in order to reduce the
pulse time needed to switch the ferroelectric, shifting the curve to smaller
times. This voltage-time trade-off is inherent to hafnia ferroelectrics 18
since they follow a nucleation-limited switching behaviour 29 as has also been
shown for FTJ devices 19. The specific voltage-time relationship depends
sensitively on material process parameters 30 and is therefore highly device-
specific. FTJ devices in particular show an asymmetry in the weight update in
positive and negative polarities, related to an internal field introduced by
the dielectric layer 31, which modifies the applied bias in each switching
direction and can be observed as an asymmetry in the switching I-V curves in
positive and negative polarities (fig. 1a). Despite a smaller negative
coercive voltage, the voltage/time window for resetting the device is larger
than for setting, due to a broader switching distribution in the negative
polarity. If the switching amplitude is not large enough, it may be difficult
to fully reset the device, even at long pulse times (see fig2d).

###  2.3 Generating an automata scheme

Figure 3: (a) FTJ current densities measured for programming voltages with an
interval of 0.1 V, for positive and negative programming pulses. Red (blue)
dashed lines indicate the current density reached with positive (negative)
programming pulses, and purple lines indicate four states with a <5% overlap.
(b) A programming scheme extracted from (a) for programming four states with
asymmetric positive and negative voltages. (c) The corresponding device
current densities measured when applying this programming scheme.

Based off of figure 2, it is possible to produce pairs of pulses in positive
and negative polarities which reach the same polarization state. This can be
described in an automata model which defines several states and the required
pulse parameters to move between them. This can form the basis of a finite
state machine and has applications in logic-in-memory circuits, where the low
readout currents of the FTJ device can be exploited for low-power operation 9.

In order to generate an automata model based on a ferroelectric device, the
experiments shown in 2a&b were first repeated with a voltage step of 0.1 V.
The results are plotted in 3a. The dashed horizontal lines represent the
current densities achieved for voltage pulses with positive voltages (red) and
negative voltages (blue). To define the automata, four equidistant states were
targeted. States where positive and negative voltages give an overlapping
current density (within $\pm 5\%$) were identified, and the most equidistant
of those were chosen as the targeted states, highlighted with purple
horizontal lines.

Figure 3b shows the programming scheme used to write the four states during
programming; the corresponding current densities are plotted in figure 3c. The
programming scheme could alternatively be written as a state machine. This
defines four deterministic states with a separation of 1 pA/$\mu$$m^{2}$ and
the positive/negative voltages required to reach them. Here it should be noted
that due to the ferroelectric switching mechanism, for intermediate states it
is still necessary to know the current state of the device in order to know
whether to apply a positive or negative switching pulse. During real device
operation, this can be circumvented by first fully switching the device in one
direction, and then partially switching back, in order to reach intermediate
states without an initial read operation.

##  3 Weight update with identical pulses

###  3.1 Time-correlated switching pulses

Figure 4: (a)Pulse trains for measuring time-correlated switching with n
switching pulses. (b) Switched polarization measured for n pulses of 10 $\mu$s
at different voltages and for a delay time between pulses of 1 $\mu$s (solid
lines) or 10 ms (dashed lines), on FTJ-2nm.

The term ´time-correlated switching pulses’ was coined by Siannas et al. 20 to
describe the effect of a cumulative weight update behaviour when the time
between switching pulses is kept very small. Such measurements have been
demonstrated previously on our FTJ-2nm devices 19 and are revisited here. A
schematic for this type of switching is shown in figure 4. The switching pulse
width tpls is fixed at 2 $\mu$s, and the switching pulse amplitude Vswitch is
varied for different measurement series. N = [1, 100] pulses are applied with
either a 1 $\mu$s or 10 ms delay (tdel) between pulses (where 10 ms delays are
chosen to represent switching on biological timescales of 100 Hz 32). The
polarization is measured with an ND (Negative-Down) triangular pulse set; this
value is then compared to the switching measured on a full switching pulse at
the start of each measurement. The results are plotted in figure 4b, where
solid (dashed) lines represent the polarization switched with a tdel of 1
$\mu$s (10 ms).

To explain the effect of time-correlated pulses, it has previously been
posited that the effective field in the capacitor does not drop to zero in
very short delay times due to the capacitor discharging time 20. It should be
noted that ionic effects may also play a role for HZO on oxide electrodes,
given the lower time constant of ionic motion 33. However, the results shown
here clearly reveal a continued switching effect for time delays far beyond
the discharging time, which can be estimated from the device RC constant
($\sim$ 1$e^{-6}$ s, based on the measured device capacitance 34). On longer
timescales, the additional weight update can be attributed to a fluid imprint
effect 35 whereby injected carriers modify the probability of domain
nucleation on timescales lower than the charge redistribution time 21. If
charges injected at the interface are particularly stable, this may eventually
lead to the accumulative switching behaviour which has been observed 22.
However, such weight update will generally be strongly time-dependent and will
lead to negligible weight update at time intervals between pulses of seconds
or longer.

In summary, in the context of pulse-delay we can distinguish between three
switching regimes and corresponding models. In the case of very short pulse
delays, NLS-like switching kinetics which treat the pulse train as a
continuous pulse can appropriately describe the weight update, where the
integral of accumulated pulse duration (pulse width x number of pulses) is
decisive for weight update, independently of the individual pulse width. That
is, the capacitors remain in a ‘excited’ state during the (< $\mu$s) pulse
delay time.

In contrast, for very long pulse delays (> s) this ‘exited’ state decays
before the next pulse arrives (see Figure 5a). Under such a condition pulses
can be treated as independent and one can assume that all domains with a
coercive field below the actual effective switching field (depending on pulse
width according the NLS mode and non ideal properties of the stack e.g. like
charge trappingl) will be flipped. Thus, allowing the use of a Preisach model
(time independent model) to simulate the behaviour of the device Any
additional similar pulse would not cause additional weight update.

Finally, there is a third regime at the transition between the two extremes.
The separation between this three regimes in terms of pulse delay and its
dependence on pulse width and pulse amplitude depends strongly not only on the
ferroelectric material properties, but the full device stack (including
interlayers and electrodes) and might strongly differ for the different device
concepts (i.e. FeFET, FeCAP or FTJ). Moreover, endurance cycling causing wake-
up and fatigue, or imprint during storage 36 might further change the devices
operation regime. Therefore, so far no satisfactory linear weight update using
identical voltage pulses was found.

###  3.2 Current-limited switching with resistor in series

In order to achieve weight update with identical pulses and arbitrary time
delays, current limiting schemes were applied. When a voltage pulse, V(t), is
applied to a ferroelectric, the current, I(t) which flows has three
components:

$I(t)=A\frac{dP(t)}{dt}+C_{FE}\frac{dV_{FE}(t)}{dt}+I_{leak}(t)$ (1)

where the first term corresponds to the switching current, proportional to the
area of the capacitor, A and the change in polarization $\frac{dP(t)}{dt}$;
the second term is the dielectric displacement current and is proportional to
the ferroelectric capacitance, CFE and the voltage drop over the ferroelectric
$\frac{dV_{FE}(t)}{dt}$; and Ileak is the leakage through the stack.

During polarization switching, a current proportional to $\frac{dP(t)}{dt}$
will flow. The principle behind this type of weight update using identical
pulses is to resitrict the polarization switching charge QP for a single
pulse, by controlling both the polarization current flow IP through the device
during switching as well as the respective switching time tP. In the ideal
case the polarization charge update yields then QP = IP * tP.

In the simplest case, the current control can be realized by a resistor, e.g.
with a resistor much larger than the parasitic resistance of the setup and
electrodes, i.e. $R_{ser}\gg R_{par}$ being inserted in series with the
ferroelectric capacitor. Then, the effective voltage $V_{FTJ}$ that is applied
to the FTJ is limited by the externally applied voltage $V_{app}$ and the
resistance $R_{ser}$, according to:

$V_{FTJ}=V_{app}-R_{ser}*I_{fe}(t),$ (2)

where $I_{FTJ}(t)$, is the capacitor current according to equation 1 with the
capacitance $C_{FE}$ replaced by $C_{FTJ}$. In that way, the field across the
ferroelectric will no longer be determined by the externally applied field
(minus a field drop across the dielectric), but rather will be self-limited by
the polarization current flow through the stack and the external resistance
37. In the equation the $V_{FTJ}$ is considered to be the voltage drop over
the whole device, i.e. the voltage drop over the ferroelectric will be lower
38. In this case the Preisach Model is extracted for the bilayer system and
thus incorporating the influence of the additional dielectric layer e.g. a
shift in the coercive field to higher values and a combined capacitance of the
ferroelectric and dielectric $C_{FTJ}$ 7.

Figure 5: (a) The weight update of an FTJ device when applying identical
pulses with an arbitrary time delay between pulses. The switching parameters
(shown inset) are 10 $\mu$s width and 2.4 V amplitude. (b) A schematic circuit
diagram showing the insertion of a resistor in series with the FTJ. The
resistor is inserted between the voltage source and the device. (c) The weight
update of an FTJ device when applying identical pulses with an arbitrary time
delay between pulses, for the circuit shown in (b). The pulse parameters
(shown inset) are 10 $\mu$s width and 4.5 V amplitude. The weight update is
plotted for a series resistance of 15 k$\Omega$ (squares) and 12 k$\Omega$
(diamonds). Measurements were performed on FTJ-1nm. (d) Simulated FTJ weight
update behaviour using the FTJ VerilogA model calibrated on FTJ-1nm data,
using the RC circuit sketched in (b).

When applying an ideal square switching pulse to this ’RC’ circuit, the
external voltage will first drop over the resistor and the capacitor needs to
charge to reach the respective coercive voltage $V_{C}$ before polarization
switching can occur 39. In addition, especially in the case of an FTJ a
background leakage current $I_{leak}$ will always flow through the capacitor.
Thus, switching will dominate on the pulse plateau, when
$\frac{dV_{app}}{dt}=0$ and with the series resistance chosen such that:

$\frac{V_{app}-V_{max}}{I_{DE}(V_{max})+I_{leak}(V_{max})+I_{switch,max}}<R_{ser}\lesssim\frac{V_{app}-V_{min}}{I_{DE}(V_{min})+I_{leak}(V_{min})}.$
(3)

Here, $V_{min}$ and $V_{max}$ denote the lower and the upper boundary of the
coercive voltage distribution, respectively. Thus, within these limits,
$R_{ser}$ will control the rate at which the polarization switches. Note, a
higher value of $R_{ser}$ will greater restrict the ferroelectric switching,
which could be compensated on the other hand by applying a higher voltage
$V_{app}$ from outside. If $R_{ser}$ is too large, the capacitor cannot charge
and polarization switching will not occur. It should be noted that such weight
update schemes have been successfully implemented in resistive random-access
memory (RRAM) systems 40.

The above discussed concept was experimentally verified. For comparison,
Figure 5a shows the behaviour of an FTJ-1nm device when identical pulses are
applied without connecting an external resistor. The pulse parameters are 2.4
V, 10 $\mu$s. The FTJ state is read out after each of 5 pulses. It can be seen
that a large weight update ($\sim 2pA/\mu m^{2}$, or 50$\%$ of the MW) occurs
on the first pulse. Note that in contrast to the experiments shown in the
previous section, the time delay between pulses is on the order of seconds,
when no accumulative effects due to a residual field in the capacitor should
be present 20. Thus, under this condition the FTJ-1nm is operated with fully
independent switching pulses where no cumulative weight update occurs. All
domains with a coercive voltage below the applied voltage will switch, and
subsequent similar pulses barely change the state, with a strongly diminishing
effect from the second pulse. This behaviour can be described by the Preisach
model 41.

In contrast, in figure 5c, full switching pulses (4.5 V, 10 $\mu$s) are
applied to an FTJ-1nm device, with resistors of different values connected in
series before the biased electrode, according to the scheme shown in figure
5b. Clearly, a more gradual weight update can be observed even thought the
switching conditions are still in a Preisach-like or non-cumulative regime.
This results from the fact, that in this current-limiting approach the
effective voltage over the ferroelectric layer increases with subsequent
pulses. Obviously, for a broader switching peak the range of the effectively
applied voltage for switching domains with lowest ($V_{min}$) towards highest
($V_{max}$) coercive voltage increases. The extension of this voltage range
$V_{max}-V_{min}$ affects the required voltage drop over $R_{ser}$ and thus,
while approaching $V_{max}$ the switching current limit lowers consecutively.
Furthermore note, that an asymmetric weight update occurs in the devices. In
negative polarity, the weight update is more abrupt but slowly saturating, due
to the lower Vc and broader switching distribution in negative polarities.
Then, as the resistor value is increased, the weight update becomes more
linear but the total polarization which can be reached decreases. This can be
attributed again to an RC effect, whereby the voltage drop over the
ferroelectric is limited by the large resistance and the capacitance of the
setup (including device and parasitics), and does not reach a large enough
value to switch all domains within the time width of the pulse.

In the example of the FTJ used here, the read current is very small compared
to the switching current ($I_{read}<I_{switch,max}$), so the series resistance
does not interfere with the read operation for any value of Rser.

###  3.3 Strategy towards linear weight update

In order to achieve efficient, accurate training in a neural network, a
synaptic device should linearly change its conductance state during its weight
update 15. As seen in the previous example, the series resistor method, while
successful in achieving gradual weight update with identical pulses, may lead
to either a non-linear behaviour or a significantly reduced MW. Furthermore,
resistor values need to be specifically chosen depending on the device size,
device capacitance and desired weight update properties. Finally, a general
issue of this series-resistor based current limit is the fact that the
controlled switching current decreases for domains that exhibit a larger
coercive voltage, due to the fact that the current limit is voltage-dependent.
In order to improve this non-linearity, a voltage-independent current limit
would be desirable.

It has been demonstrated via current-limited PFM measurements on BiFeO3 that
programming with an increasing current compliance allows multilevel switching
42. In that work, the compliance was continuously increased in order to
control the number of switching domains, and the scheme has recently been
demonstrated on ferroelectric HfAlOx FeFETs 43.

Figure 6: (a) 1T1C circuit modelled for controlling the current during
switching. Switching pulses of 4.5 V, 10 $\mu$s are applied on the plate line
(PL) while a voltage applied on the word line (WL) limits the switching
current. (b) Polarization and (c) read current density of an FTJ device
simulated for 8 pulses at different WL voltages, for a 9500 $\mu$m2 device.
(d) Polarization and (e) read current density of an FTJ device simulated for 8
pulses at different WL voltages, for a 0.25 $\mu$m2 device.

Similarly, to increase the linearity of the weight update, we propose an
improved solution, whereby the current reaching the ferroelectric capacitor is
limited by the drain current of a transistor that is operated either in sub-
threshold or saturation mode. In this case, the capacitor is connected between
drain and the pulse voltage source (e.g. in a typical 1T1C FeRAM
configuration), and the drain current is independent of the drain voltage. In
contrast to the method demonstrated above, no changing current needs to be
applied to the ferroelectric device: rather, the word line (WL) voltage of the
transistor is kept constant, significantly simplifying the application of this
scheme.

In this case, experiments on large capacitors have not yet been enabled due to
the difficulty of applying a current compliance in pulsed operation using
commercial probe stations. Instead, the principle is demonstrated via Cadence
simulations of simple bitcells. The ferroelectric device is modelled using a
Preisach approach 9, 7.

In order to prove this concept, we prove via circuit simulations that
multilevel behaviour can be achieved with identical pulses by applying a
current compliance on-chip through an access transistor. In order to calibrate
our simulation deck in a first step the experimental results utilizing a
series resistor as discussed in the previous section have been repeated in
simulation. The FTJ model was calibrated to exactly replicate the PV-hystereis
curve of the FTJ-1nm device under test (see Figure S1). Then the same circuit
as shown in Figure 3b was implemented. Thereby, an additional parasitic cable
capacitance of 500 pF was connected towards ground in order to replicate the
cable capacitance from our experimental setup. 5d reveals a good agreement
between measurement and simulation. In this simulation identical pulses with
an amplitude of 4.5 V and a width of 10 $\mu$s have been applied. In between
each switching pulse a smaller read voltage pulse was utilized to monitor the
FTJ read current. The small differences in the weight update behavior between
experiment and simulation can be attributed to small deviations in the exact
PV-curve between simulation and experiment, and the occurrence of potential
imprint effects that so far are not included in our model 21, 36.

In order to demonstrate improved weight update linearity, the same calibrated
FTJ model was then implemented into a 1T1C circuit (Figure 6a). Similar as in
the case above, identical switching pulses have then been applied. This
simulation was performed for varying gate voltages, thus resulting in
different drain current limits. The resulting weight update behavior is
depicted in Figures 6b-c for a large device as measured above, and in Figures
6d-e for a scaled device with nm side lengths. It can be clearly seen that the
linearity of weight update can be considerably improved while retaining a
large MW. Moreover, utilizing the transistor as current limit allows a
flexible tuning of the slope of the weight update. Figures 6c and 6e report
linearity factors (LF) which describe the linearity of the weight update
behaviour, calculated by:

$LF=\frac{J_{7}-J_{0}}{J_{3}-J_{0}}/2$ (4)

where Jn represents the current density on the nth switching pulse and J0
represents the initial device current density i.e. in the fully Off state. The
parameter LF should be 1 for a linearly switching device and can be expressed
as a percentage. LF reaches 0.87 (87%) for the larger simulated device ($\mu$m
side lengths), and 0.81 (81%) for the smaller simulated device (nm side
lengths), compared to 65% (79%) measured for the switching with a 12 k$\Omega$
(12 k$\Omega$) resistor in series, respectively. Furthermore, the larger
device shows a LF of 0.83 (83%) for switching covering the whole MW, whereas
the highest LF for switching with the series resistance method covers a more
limited MW.

The remaining non-linearity can be attributed to the large distribution of the
switching voltages of the individual domains within our device. Even though
the current limit is now independent on the voltage drop at the transistors
drain terminal, the required charging of the dielectric background capacitance
of the FTJ consumes an increasing amount of charge with increasing effective
switching voltage. Hence, in order to obtain perfect linearization of the
weight update, a very narrow ferroelectric switching peak is desirable. These
results show that this is even more critical for devices scaled down to the nm
regime.

##  4 Conclusions

In conclusion, we have presented here several methods for achieving weight
update in devices based on ferroelectric capacitors, with an FTJ device as a
model example. It should be stressed that since the weight update occurs in
the ferroelectric film, these methods will be applicable (with possible
circuit modifications) in other device types, including FeCAP, FeFET and
potential other device types.

To achieve weight update with non-identical pulses, a method was presented to
find an automata (or state machine) model with 4 states, which could also be
used for a higher number of states. Here, deterministic current levels could
be defined in an FTJ device, and the required pulses to transition between
them. These require complex circuitry to apply on-chip, where time and/or
voltage of the applied pulses should be dynamically modified 24.

Besides that, two methods for achieving weight update with identical pulses
and arbitrary time delays were presented. The first, demonstrated
experimentally, involves adding a resistor in series with the ferroelectric
device. Depending on the resistor value, a linearity of up to 79% can be
achieved, but owing to an RC effect, higher resistor values limit the memory
window. Nonetheless, several distinct levels could be reached. This operation
is highly different from the case of identical pulses and short fixed time
delays between pulses, where if the delay time between pulses is increased,
the weight update diminishes 20, 21.The time delay between pulses again is on
the order of seconds, which makes the pulse sequence highly applicable for a
number of neuromorphic computing elements 15 including integrate-and-fire
neurons 7 or as a synaptic element.

The final method, demonstrated via Cadence simulations, involves performing
the weight update of a capacitor under an external current compliance. Here,
the switching is limited by the drain current of a transistor that is operated
either in sub-threshold or saturation mode. We have demonstrated that this
works for different device sizes and can achieve linearity up to 87% covering
most of the memory window. With engineering of the device switching
distribution, this linearity is expected to improve further.

##  5 Experimental section

###  5.1 Sample preparation

FTJ capacitors were fabricated by sputter deposition, atomic layer deposition
(ALD), metal evaporation and dry etching (inductively coupled plasma, ICP).
Sputtered metallic electrodes were employed consisting of W (30 nm)/TiN (5
nm)/TiAl (7.5 nm)/TiN (5 nm), forming TiAlN, on the BE and TiAlN on the TE.
Al2O3 was deposited via ALD (precursor: TMA, oxidant: H2O, T = 260°C) followed
by 10.5 nm HZO (precursors: HyALD, ZyALD, oxidant: O3, T = 280°C). The samples
were annealed at 450°C in Ar atmosphere for 30 mins. Capacitors were defined
via the metal evaporation of Ti/Pt (10/25 nm) which was used as a hard mask
for rate-controlled ICP etching of the TE in BCl3 chemistry. Chlorine
chemistry is needed to etch the TiAlN electrodes but the rate must be measured
since this chemistry can also etch HZO 44.

###  5.2 Electrical measurements

All electrical analyses were performed on a Cascade Microtech probe station
using a Keithley 4200 Parameter Analyser with Keithley 4225 RPMs which allow
pulsed measurements down to a pulse width of 20 ns. P-V loops were extracted
using a PUND sequence at 20 kHz. For multilevel experiments, pulse rise and
fall times were chosen as 500 ns, which is one order of magnitude larger than
the resolution of the measurement equipment and should not be affected by any
internal parasitics. Read-out pulses had a width of 50 ms, which was chosen to
further avoid influence of RC effects which delay the displacement current in
large capacitors, obscuring the leakage current in the read-out.

###  5.3 Switching simulations

A VerilogA file modelling the ferroelectric behaviour with a Preisach approach
was constructed (available upon reasonable request) and used in Cadence
simulations. The FTJ parameters were fitted from experimental data and
verified by comparison of simulated and measured P-V loops (Figure S1). The
modelled bitcells are shown schematically in Figure 5b and 6a. These were
constructed of the FTJ device model as well as a standard voltage source
(analogLib), resistor (Primlib) and transistor (ne5).

Acknowledgements

This work was supported by the DFG Priority Program ´Memristec’ through the
project ReLoFeMris (441909639), and by the European Commission through the
project FIXIT (No. GA: 101135398). T.M. was financially supported out of the
Saxonian State budget approved by the delegates of the Saxon State Parliament.
Views and opinions expressed are those of the authors only and do not
necessarily reflect those of the European Union.

## References

  * Böscke et al. 2011 Böscke, T.; Müller, J.; Bräuhaus, D.; Schröder, U.; Böttger, U. Ferroelectricity in hafnium oxide thin films. _Applied Physics Letters_ 2011, _99_
  * Francois et al. 2022 Francois, T.; Coignus, J.; Makosiej, A.; Giraud, B.; Carabasse, C.; Barbot, J.; Martin, S.; Castellani, N.; Magis, T.; Grampeix, H., et al. High-performance operation and solder reflow compatibility in BEOL-integrated 16-kb HfO 2: Si-based 1T-1C FeRAM arrays. _IEEE Transactions on Electron Devices_ 2022, _69_ , 2108–2114
  * Dahan et al. 2023 Dahan, M. M.; Mulaosmanovic, H.; Levit, O.; Duünkel, S.; Beyer, S.; Yalon, E. Sub-Nanosecond Switching of Si: HfO2 Ferroelectric Field-Effect Transistor. _Nano Letters_ 2023, _23_ , 1395–1400
  * Covi et al. 2021 Covi, E.; Duong, Q. T.; Lancaster, S.; Havel, V.; Coignus, J.; Barbot, J.; Richter, O.; Klein, P.; Chicca, E.; Grenouillet, L., et al. Ferroelectric tunneling junctions for edge computing. 2021 IEEE International Symposium on Circuits and Systems (ISCAS). 2021; pp 1–5
  * De et al. 2022 De, S.; Müller, F.; Thunder, S.; Abdulazhanov, S.; Laleni, N.; Lederer, M.; Ali, T.; Raffel, Y.; Dünkel, S.; Mojumder, S., et al. 28 nm HKMG-based current limited FeFET crossbar-array for inference application. _IEEE Transactions on Electron Devices_ 2022, _69_ , 7194–7198
  * Ryu et al. 2019 Ryu, H.; Wu, H.; Rao, F.; Zhu, W. Ferroelectric tunneling junctions based on aluminum oxide/zirconium-doped hafnium oxide for neuromorphic computing. _Scientific reports_ 2019, _9_ , 20383
  * Gibertini et al. 2022 Gibertini, P.; Fehlings, L.; Lancaster, S.; Duong, Q.; Mikolajick, T.; Dubourdieu, C.; Slesazeck, S.; Covi, E.; Deshpande, V. A Ferroelectric Tunnel Junction-based Integrate-and-Fire Neuron. 2022 29th IEEE International Conference on Electronics, Circuits and Systems (ICECS). 2022; pp 1–4
  * O’Connor et al. 2018 O’Connor, I.; Cantan, M.; Marchand, C.; Vilquin, B.; Slesazeck, S.; Breyer, E. T.; Mulaosmanovic, H.; Mikolajick, T.; Giraud, B.; Noël, J.-P., et al. Prospects for energy-efficient edge computing with integrated HfO 2-based ferroelectric devices. 2018 IFIP/IEEE International Conference on Very Large Scale Integration (VLSI-SoC). 2018; pp 180–183
  * Reuben et al. 2023 Reuben, J.; Fey, D.; Lancaster, S.; Slesazeck, S. A Low-Power Ternary Adder Using Ferroelectric Tunnel Junctions. _Electronics_ 2023, _12_ , 1163
  * Yang 2020 Yang, R. In-memory computing with ferroelectrics. _Nature Electronics_ 2020, _3_ , 237–238
  * Bégon-Lours et al. 2022 Bégon-Lours, L.; Halter, M.; Puglisi, F. M.; Benatti, L.; Falcone, D. F.; Popoff, Y.; Dávila Pineda, D.; Sousa, M.; Offrein, B. J. Scaled, Ferroelectric Memristive Synapse for Back-End-of-Line Integration with Neuromorphic Hardware. _Advanced Electronic Materials_ 2022, _8_ , 2101395
  * Goh et al. 2021 Goh, Y.; Hwang, J.; Kim, M.; Jung, M.; Lim, S.; Jung, S.-O.; Jeon, S. High performance and self-rectifying hafnia-based ferroelectric tunnel junction for neuromorphic computing and TCAM applications. 2021 IEEE International Electron Devices Meeting (IEDM). 2021; pp 17–2
  * Moon et al. 2019 Moon, K.; Lim, S.; Park, J.; Sung, C.; Oh, S.; Woo, J.; Lee, J.; Hwang, H. RRAM-based synapse devices for neuromorphic systems. _Faraday discussions_ 2019, _213_ , 421–451
  * Oh et al. 2017 Oh, S.; Kim, T.; Kwak, M.; Song, J.; Woo, J.; Jeon, S.; Yoo, I. K.; Hwang, H. HfZrO x-based ferroelectric synapse device with 32 levels of conductance states for neuromorphic applications. _IEEE Electron Device Letters_ 2017, _38_ , 732–735
  * Ielmini and Ambrogio 2019 Ielmini, D.; Ambrogio, S. Emerging neuromorphic devices. _Nanotechnology_ 2019, _31_ , 092001
  * Gokmen and Vlasov 2016 Gokmen, T.; Vlasov, Y. Acceleration of deep neural network training with resistive cross-point devices: Design considerations. _Frontiers in neuroscience_ 2016, 333
  * Tagantsev et al. 2002 Tagantsev, A. K.; Stolichnov, I.; Setter, N.; Cross, J. S.; Tsukada, M. Non-Kolmogorov-Avrami switching kinetics in ferroelectric thin films. _Physical Review B_ 2002, _66_ , 214109
  * Mulaosmanovic et al. 2020 Mulaosmanovic, H.; Müller, F.; Lederer, M.; Ali, T.; Hoffmann, R.; Seidel, K.; Zhou, H.; Ocker, J.; Mueller, S.; Dünkel, S., et al. Interplay between switching and retention in HfO 2-based ferroelectric FETs. _IEEE Transactions on Electron Devices_ 2020, _67_ , 3466–3471
  * Lancaster et al. 2022 Lancaster, S.; Mikolajick, T.; Slesazeck, S. A multi-pulse wakeup scheme for on-chip operation of devices based on ferroelectric doped HfO2 thin films. _Applied Physics Letters_ 2022, _120_
  * 20 Siannas, N.; Zacharaki, C.; Tsipas, P.; Kim, D. J.; Hamouda, W.; Istrate, C.; Pintilie, L.; Schmidbauer, M.; Dubourdieu, C.; Dimoulas, A. Electronic Synapses Enabled by an Epitaxial SrTiO3-$\delta$/Hf0. 5Zr0. 5O2 Ferroelectric Field-Effect Memristor Integrated on Silicon. _Advanced Functional Materials_ 2311767
  * Vecchi et al. 2024 Vecchi, S.; Puglisi, F. M.; Appelt, P.; Guido, R.; Wang, X.; Slesazeck, S.; Mikolajick, T.; Lancaster, S. Evaluation of Imprint and Multi-Level Dynamics in Ferroelectric Capacitors. _Advanced Electronic Materials_ 2024, 2400204
  * Mulaosmanovic et al. 2020 Mulaosmanovic, H.; Dünkel, S.; Trentzsch, M.; Beyer, S.; Breyer, E. T.; Mikolajick, T.; Slesazeck, S. Investigation of accumulative switching in ferroelectric FETs: Enabling universal modeling of the switching behavior. _IEEE Transactions on Electron Devices_ 2020, _67_ , 5804–5809
  * Hoffmann et al. 2022 Hoffmann, M.; Tan, A. J.; Shanker, N.; Liao, Y.-H.; Wang, L.-C.; Bae, J.-H.; Hu, C.; Salahuddin, S. Write Disturb-Free Ferroelectric FETs With Non-Accumulative Switching Dynamics. _IEEE Electron Device Letters_ 2022, _43_ , 2097–2100
  * Narayanan et al. 2022 Narayanan, S.; Covi, E.; Havel, V.; Frenkel, C.; Lancaster, S.; Duong, Q.; Slesazeck, S.; Mikolajick, T.; Payvand, M.; Indiveri, G. A 120dB Programmable-Range On-Chip Pulse Generator for Characterizing Ferroelectric Devices. 2022 IEEE International Symposium on Circuits and Systems (ISCAS). 2022; pp 717–721
  * Zhou et al. 2020 Zhou, H.; Ocker, J.; Padovani, A.; Pesic, M.; Trentzsch, M.; Dünkel, S.; Mulaosmanovic, H.; Slesazeck, S.; Larcher, L.; Beyer, S., et al. Application and benefits of target programming algorithms for ferroelectric HfO 2 transistors. 2020 IEEE International Electron Devices Meeting (IEDM). 2020; pp 18–6
  * Chanthbouala et al. 2012 Chanthbouala, A.; Crassous, A.; Garcia, V.; Bouzehouane, K.; Fusil, S.; Moya, X.; Allibe, J.; Dlubak, B.; Grollier, J.; Xavier, S., et al. Solid-state memories based on ferroelectric tunnel junctions. _Nature nanotechnology_ 2012, _7_ , 101–104
  * Max et al. 2020 Max, B.; Hoffmann, M.; Mulaosmanovic, H.; Slesazeck, S.; Mikolajick, T. Hafnia-based double-layer ferroelectric tunnel junctions as artificial synapses for neuromorphic computing. _ACS Applied Electronic Materials_ 2020, _2_ , 4023–4033
  * Lancaster et al. 2023 Lancaster, S.; Segatto, M.; Silva, C.; Max, B.; Mikolajick, T.; Esseni, D.; Driussi, F.; Slesazeck, S. Reducing the tunneling barrier thickness of bilayer ferroelectric tunnel junctions with metallic electrodes. 2023 Device Research Conference (DRC). 2023; pp 1–2
  * Zhou et al. 2013 Zhou, D.; Xu, J.; Li, Q.; Guan, Y.; Cao, F.; Dong, X.; Müller, J.; Schenk, T.; Schröder, U. Wake-up effects in Si-doped hafnium oxide ferroelectric thin films. _Applied Physics Letters_ 2013, _103_
  * Hsain et al. 2022 Hsain, H. A.; Lee, Y.; Lancaster, S.; Materano, M.; Alcala, R.; Xu, B.; Mikolajick, T.; Schroeder, U.; Parsons, G. N.; Jones, J. L. Role of Oxygen Source on Buried Interfaces in Atomic-Layer-Deposited Ferroelectric Hafnia–Zirconia Thin Films. _ACS Applied Materials & Interfaces_ 2022, _14_ , 42232–42244
  * Lancaster et al. 2022 Lancaster, S.; Lomenzo, P. D.; Engl, M.; Xu, B.; Mikolajick, T.; Schroeder, U.; Slesazeck, S. Investigating charge trapping in ferroelectric thin films through transient measurements. _Frontiers in Nanotechnology_ 2022, _4_ , 939822
  * Jackson et al. 2013 Jackson, B. L.; Rajendran, B.; Corrado, G. S.; Breitwisch, M.; Burr, G. W.; Cheek, R.; Gopalakrishnan, K.; Raoux, S.; Rettner, C. T.; Padilla, A., et al. Nanoscale electronic synapses using phase change devices. _ACM Journal on Emerging Technologies in Computing Systems (JETC)_ 2013, _9_ , 1–20
  * Halter et al. 2023 Halter, M.; Bégon-Lours, L.; Sousa, M.; Popoff, Y.; Drechsler, U.; Bragaglia, V.; Offrein, B. J. A multi-timescale synaptic weight based on ferroelectric hafnium zirconium oxide. _Communications Materials_ 2023, _4_ , 14
  * Benatti et al. 2023 Benatti, L.; Vecchi, S.; Puglisi, F. M. Linking the intrinsic electrical response of ferroelectric devices to material properties by means of impedance spectroscopy. _IEEE Transactions on Device and Materials Reliability_ 2023,
  * Buragohain et al. 2019 Buragohain, P.; Erickson, A.; Kariuki, P.; Mittmann, T.; Richter, C.; Lomenzo, P. D.; Lu, H.; Schenk, T.; Mikolajick, T.; Schroeder, U., et al. Fluid imprint and inertial switching in ferroelectric La: HfO2 capacitors. _ACS applied materials & interfaces_ 2019, _11_ , 35115–35121
  * Barbot et al. 2023 Barbot, J.; Fontanini, R.; Segatto, M.; Coignus, J.; Triozon, F.; Carabasse, C.; Bedjaoui, M.; Andrieu, F.; Esseni, D.; Grenouillet, L. Dynamics of polarization loss and imprint in bilayer ferroelectric tunnel junctions. _Journal of Applied Physics_ 2023, _134_
  * Massarotto et al. 2024 Massarotto, M.; Driussi, F.; Bucovaz, M.; Affanni, A.; Lancaster, S.; Slesazeck, S.; Mikolajick, T.; Esseni, D. Analysis and Compensation of the Series Resistance Effects on the Characteristics of Ferroelectric Capacitors. 2024 IEEE 36th International Conference on Microelectronic Test Structures (ICMTS). 2024; pp 1–5
  * Park et al. 2021 Park, H. W.; Hyun, S. D.; Lee, I. S.; Lee, S. H.; Lee, Y. B.; Oh, M.; Kim, B. Y.; Ryoo, S. G.; Hwang, C. S. Polarizing and depolarizing charge injection through a thin dielectric layer in a ferroelectric–dielectric bilayer. _Nanoscale_ 2021, _13_ , 2556–2572
  * Si et al. 2019 Si, M.; Lyu, X.; Shrestha, P. R.; Sun, X.; Wang, H.; Cheung, K. P.; Ye, P. D. Ultrafast measurements of polarization switching dynamics on ferroelectric and anti-ferroelectric hafnium zirconium oxide. _Applied Physics Letters_ 2019, _115_
  * Moon et al. 2017 Moon, K.; Kwak, M.; Park, J.; Lee, D.; Hwang, H. Improved conductance linearity and conductance ratio of 1T2R synapse device for neuromorphic systems. _IEEE Electron Device Letters_ 2017, _38_ , 1023–1026
  * Bartic et al. 2001 Bartic, A. T.; Wouters, D. J.; Maes, H. E.; Rickes, J. T.; Waser, R. M. Preisach model for the simulation of ferroelectric capacitors. _Journal of Applied Physics_ 2001, _89_ , 3420–3425
  * Lee et al. 2012 Lee, D.; Yang, S. M.; Kim, T. H.; Jeon, B. C.; Kim, Y. S.; Yoon, J.-G.; Lee, H. N.; Baek, S. H.; Eom, C. B.; Noh, T. W. Multilevel data storage memory using deterministic polarization control. _Advanced Materials_ 2012, _24_ , 402–406
  * Kim and Lee 2024 Kim, I.-J.; Lee, J.-S. Unlocking large memory windows and 16-level data per cell memory operations in hafnia-based ferroelectric transistors. _Science Advances_ 2024, _10_ , eadn1345
  * Mauersberger et al. 2021 Mauersberger, T.; Trommer, J.; Sharma, S.; Knaut, M.; Pohl, D.; Rellinghaus, B.; Mikolajick, T.; Heinzig, A. Single-step reactive ion etching process for device integration of hafnium-zirconium-oxide (HZO)/titanium nitride (TiN) stacks. _Semiconductor Science and Technology_ 2021, _36_ , 095025


# Building a testable shear viscosity across the QCD phase diagram

Emma McLaughlin  Department of Physics, Columbia University, 538 West 120th
Street, New York, NY 10027, USA  Jacob Rose  Department of Physics, Dr. Karl-
Remeis Observatory - Astronomical Institute, Bamberg D-96049, Germany  Travis
Dore  Illinois Center for Advanced Studies of the Universe, Department of
Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Paolo Parotto  University of Wuppertal, Department of Physics, Wuppertal
D-42097, Germany  Claudia Ratti  Department of Physics, University of Houston,
Houston, TX 77204, USA  Jacquelyn Noronha-Hostler  Illinois Center for
Advanced Studies of the Universe, Department of Physics, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

Current experiments at the Relativistic Heavy Ion Collider (RHIC) are probing
finite baryon densities where the shear viscosity to enthalpy ratio $\eta T/w$
of the Quark Gluon Plasma remains unknown. We use the Hadron Resonance Gas
(HRG) model with the most up-to-date hadron list to calculate $\eta T/w$ at
low temperatures and at finite baryon densities $\rho_{B}$. We then match
$\eta T/w$ to a QCD-based shear viscosity calculation within the deconfined
phase to create a table across $\left\\{T,\mu_{B}\right\\}$ for different
cross-over and critical point scenarios at a specified location. We find that
these new $\eta T/w(T,\mu_{B})$ values would require initial conditions at
significantly larger $\rho_{B}$, compared to ideal hydrodynamic trajectories,
in order to reach the same freeze-out point.

##  I Introduction

Relativistic heavy-ion collisions at the Large Hadron Collider (LHC) and
Relativistic Heavy-Ion Collider (RHIC) have successfully reproduced the phase
transition from the Quark Gluon Plasma into a hadron resonance gas in the
laboratory. Since 2006 it has been understood that this phase transition is a
smooth cross-over Aoki  _et al._ (2006); Borsanyi  _et al._ (2010a, b);
Bazavov  _et al._ (2014); Borsanyi  _et al._ (2020); Ratti (2018); Guenther
(2020). For this reason, the definition of the pseudo-critical temperature
$T_{pc}$ is not unique, but rather depends on the observable considered. If
Quantum Chromodynamics (QCD) exhibits a critical point at large baryon
densities, one expects the different pseudo-critical temperatures to converge
to the critical temperature $T_{c}$ Endrodi  _et al._ (2011); Rougemont  _et
al._ (2017); Grefa  _et al._ (2021); Bellwied  _et al._ (2019); Noronha-
Hostler  _et al._ (2016a) at the critical point. However, it also possible
that there is no critical point and only a cross-over is seen at larger baryon
densities.

One crucial achievement of heavy-ion collisions has been the discovery that
the Quark Gluon Plasma is a nearly perfect fluid with a shear viscosity to
entropy density ratio $\eta/s$ which at $\mu_{B}=0$ is close to the KSS limit
Kovtun  _et al._ (2005) (that we now know is no longer a strict bound Kats and
Petrov (2009); Brigante  _et al._ (2008a, b); Buchel  _et al._ (2009);
Brigante  _et al._ (2008b); Critelli  _et al._ (2014); Finazzo  _et al._
(2016)). Due to the extremely short lifetime of the produced Quark Gluon
Plasma, relativistic viscous hydrodynamic calculations have been vital to
understand the dynamics of this cross-over phase transition. It is common
practice in the field to use $\eta/s$ as a free parameter, and extract its
value from experimental data by comparing flow harmonics to theoretical
predictions from relativistic hydrodynamics Romatschke and Romatschke (2007);
Bozek (2012); Heinz and Snellings (2013); Luzum and Petersen (2014); Niemi
_et al._ (2016); Noronha-Hostler  _et al._ (2016b); McDonald  _et al._ (2017);
Bernhard  _et al._ (2019); Alba  _et al._ (2018a). However, since there are
many other free parameters that influence the extraction of $\eta/s$, this is
not a trivial task.

Ideally, $\eta/s$ would be calculated directly from lattice QCD but this turns
out to be an ill-posed problem, as it involves inversion methods over a
discrete set of lattice points for the correlator of the energy-momentum
tensor. Besides, the signal for this specific correlator is dominated by the
high-$\omega$ part of the spectral function Aarts and Martinez Resco (2002);
Borsanyi  _et al._ (2018), which makes the inversion even harder. For these
reasons, one must turn to alternative approaches such as the hadron resonance
gas (HRG) model Noronha-Hostler  _et al._ (2009); Pal (2010); Khvorostukhin
_et al._ (2010); Tawfik and Wahba (2010); Alba  _et al._ (2015); Ratti  _et
al._ (2011); Tiwari  _et al._ (2012); Noronha-Hostler  _et al._ (2012); Kadam
and Mishra (2014, 2015, 2016); Kadam and Pawar (2019); Mohapatra  _et al._
(2019), transport theory Wesp  _et al._ (2011); Ozvenchuk  _et al._ (2013);
Rose  _et al._ (2018); Rais  _et al._ (2020), holography Kovtun  _et al._
(2005), Color String Percolation Model Sahu and Sahoo (2020), linear sigma
model Heffernan  _et al._ (2020) or QCD-motivated alternatives Haas  _et al._
(2014); Christiansen  _et al._ (2015); Dubla  _et al._ (2018); Ghiglieri  _et
al._ (2018), to name a few. These calculations are often only performed at
$\mu_{B}=0$. Hence, they would be appropriate primarily for LHC and top RHIC
energies, but not for large baryon densities.

While calculations currently exist from some of these models at finite baryon
densities, there are difficulties in systematically studying $\eta
T/w(T,\mu_{B})$111Note that at finite $\mu_{B}$ the enthalpy $w=\varepsilon+p$
is used rather than entropy for different types of phase transitions (cross-
over vs. a critical point within the dynamical H or B university classes)
since some of these models are only applicable within the hadronic phase (e.g.
HRG model Denicol  _et al._ (2013); Kadam and Mishra (2014), UrQMD Demir and
Bass (2009), and SMASH Rose  _et al._ (2018) calculations), or they are
beholden to dynamic university class of the model itself (e.g. holography in
Universality class B Rougemont  _et al._ (2017); Critelli  _et al._ (2017),
which is also limited to a constant $\eta/s(T,\mu_{B})=0.08$ for current
studies at finite $\mu_{B}$ 222Alternatives to $\eta/s(T,\mu_{B}=0)=0.08$
exist Kats and Petrov (2009); Brigante  _et al._ (2008a, b); Buchel  _et al._
(2009) at vanishing baryon densities where an action is formulated that allows
for derivatives up to the 4th order. However, these calculations have not yet
been incorporated into non-conformal AdS. Since non-conformality is vital for
understanding the QCD phase transition, the current framework cannot provide
non-trivial information about $\eta/s(T,\mu_{B})$ at a critical point.).
Furthermore, at larger baryon densities repulsive interactions have been shown
to be more relevant Vovchenko  _et al._ (2017a); Huovinen and Petreczky
(2018), which should be taken into account in the hadronic phase.

In this paper we establish a framework for combining $\eta T/w(T,\mu_{B})$ in
the confined phase – from an interacting hadron resonance gas (HRG) – and in
the deconfined phase – with a QCD-based assumption – in the presence of a
cross-over phase transition. The key element is that we use the finite
$\mu_{B}$ behavior from the HRG model to construct $\eta T/w(T,\mu_{B})$
across the phase diagram relevant to heavy-ion collisions. We note that the
framework is generic enough to be easily updated if more realistic QCD-based
calculations become available at finite $\mu_{B}$. We use an interacting
hadron resonance gas with the most up-to-date particle resonance list from the
Particle Data Group (PDG16+) Alba  _et al._ (2017), which was shown to be a
reasonable fit compared to lattice QCD data, can describe net-Kaon
fluctuations Bellwied  _et al._ (2019), off-diagonal susceptibilities Bellwied
_et al._ (2020), thermal fits Alba  _et al._ (2020), and works well within
relativistic viscous hydrodynamic calculations Alba  _et al._ (2018a). We base
the treatment of the deconfined phase on the parametrized version of
$\eta/s(T,\mu_{B}=0)$ from Ref. Dubla  _et al._ (2018), which is adjusted to
match the hadron resonance gas model at the phase transition assuming a
minimum value $\eta/s\sim 0.08$. In this work we construct four different
profiles for $\eta T/w(T,\mu_{B})$: we first consider the cases of a smooth or
sharp crossover, then the case with a critical point, first at
$\left\\{T,\mu_{B}\right\\}=\left\\{143,350\right\\}$ MeV – to match to the
publicly available BEST collaboration EoS Parotto  _et al._ (2020), then at
$\left\\{T,\mu_{B}\right\\}=\left\\{89,724\right\\}$ MeV – matching the
prediction from holography in Ref. Critelli  _et al._ (2017)333We remind the
reader that, while the result in Ref. Critelli  _et al._ (2017) is a
prediction on the location of the critical point, the BEST collaboration EoS
allows the user to pick its location, and the one mentioned here is the one
used in Ref. Parotto  _et al._ (2020) for illustration purposes..

##  II Hadron Resonance Gas Model

With the HRG model one can calculate the pressure, energy density and total
particle density of species $i$ assuming that the hadrons are point like
particles:

$\displaystyle\frac{p(T,\mu_{B},\mu_{S},\mu_{Q})}{T^{4}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i}(-1)^{B_{i}+1}\frac{g_{i}}{2\pi^{2}}\int\limits_{0}^{\infty}p^{2}\ln{[1+(-1)^{B_{i}+1}e^{(-\frac{\sqrt{p^{2}+m_{i}^{2}}}{T}+\tilde{\mu}_{i})}]}dp$
(1) $\displaystyle\frac{\epsilon(T,\mu_{B},\mu_{S},\mu_{Q})}{T^{4}}$
$\displaystyle=$
$\displaystyle\sum\limits_{i}\frac{g_{i}}{2\pi^{2}}\int\limits_{0}^{\infty}\frac{p^{2}\sqrt{p^{2}+m_{i}^{2}}}{(-1)^{B_{i}+1}+e^{(-\frac{\sqrt{p^{2}+m_{i}^{2}}}{T}+\tilde{\mu}_{i})}}dp$
(2) $\displaystyle\frac{n_{i}(T,\mu_{B},\mu_{S},\mu_{Q})}{T^{3}}$
$\displaystyle=$
$\displaystyle\frac{g_{i}}{2\pi^{2}}\int\limits_{0}^{\infty}p^{2}\left[\exp\left(\frac{\sqrt{p^{2}+m_{i}^{2}}}{T}-\tilde{\mu}_{i}\right)+(-1)^{B_{i}-1}\right]^{-1}dp$
(3)

where

$\tilde{\mu}_{i}\equiv B_{i}\mu_{B}/T+S_{i}\mu_{S}/T+Q_{i}\mu_{Q}/T$ (4)

and $g_{i}$ is the degeneracy of each hadron, $m_{i}$ is the mass, and
$B_{i}$, $S_{i}$, and $Q_{i}$ are the baryon number, strangeness and electric
charge carried by each hadron. Additionally, $\mu_{B}$, $\mu_{S}$, and
$\mu_{Q}$ are the corresponding chemical potentials for each conserved charge.

The other thermodynamic quantities follow

$\displaystyle s(T,\mu_{B},\mu_{S},\mu_{Q})$ $\displaystyle=$
$\displaystyle\frac{\partial p(T,\mu_{B},\mu_{S},\mu_{Q})}{\partial T}$
$\displaystyle\frac{\rho_{i}(T,\mu_{B},\mu_{S},\mu_{Q})}{T^{3}}$
$\displaystyle=$ $\displaystyle\frac{\partial
p(T,\mu_{B},\mu_{S},\mu_{Q})/T^{4}}{\partial(\mu_{i}/T)}$ $\displaystyle
p+\epsilon$ $\displaystyle=$ $\displaystyle sT+\sum_{i=B,S,Q}\mu_{i}\rho_{i}$
(5)

where $s$ is the entropy, $\rho_{i}$ where $i=B,S,Q$ is the net density, and
$\epsilon$ is the energy density. Note that for the rest of the paper we will
abbreviate $\mu_{B},\mu_{S},\mu_{Q}\Rightarrow\mu_{i}$.

In this paper we first compare two different lists of hadrons from the
Particle Data Group, one from 2005 (PDG05) and another developed in Ref. Alba
_et al._ (2017) from 2016 that includes all *-**** states (PDG16+).

###  II.1 Excluded Volume

One method for taking into account repulsive interactions is the excluded
volume approach Rischke  _et al._ (1991) wherein each hadron is delegated a
volume $v$, and the excluded volume pressure $p_{v}$ can then be written as

$\frac{p_{v}(T,\mu_{i})}{T}=n(T,\mu_{i})\exp(-vp_{v}(T,\mu_{i})/T),$ (6)

which can be solved analytically using the Lambert W function

$p_{v}(T,\mu_{i})=\frac{T}{v}W(v\,n(T,\mu_{i}))\,,$ (7)

where $n(T,\mu_{i})=\sum_{i}n_{i}(T,\mu_{i})$ is defined in Eq. (3). The
remaining thermodynamic quantities can be obtained using the thermodynamic
relationships from Eq. (5) and are denoted with a subscript $v$ for excluded
volume. Here we are assuming that all particles have the same volume. In fact,
one could relax that assumption as in Ref. Albright  _et al._ (2014) or even
consider a multicomponent van der Waals Vovchenko  _et al._ (2017b), but we
leave this for future work.

We obtain the effective hard-core volume from $v=4\cdot 4\pi r^{3}/3$, where
$r$ is the effective core radius. The point-like thermodynamic properties in
Eqs. (1-5) can be then restored with $r\rightarrow 0$.

Figure 1:  Comparison of the PDDG05 and PDG16+ list with various excluded
volumes to $s/T^{3}$ results from lattice QCD Borsanyi  _et al._ (2010a). The
band for the PDG16+ denotes the minimum and maximum $r$ that fits within
constraints from Lattice QCD.  Figure 2:  Pressure comparison to lattice QCD
results from Borsanyi  _et al._ (2010a) with an ideal vs. excluded volume HRG
model.

In the excluded volume calculations, there is always the issue of constraining
the volume size, since it is an unconstrained free parameter. In an attempt to
constrain the hard core radius, we vary $r$ such that is fits within lattice
QCD error-bars up to $T\sim 150$ MeV using the entropy from lattice QCD
results shown in Fig. 1 for the entropy and in Fig. 2 for the pressure. We
find that, for the PDG05, there are so few resonances that the largest
possible excluded volume that we can use corresponds to a radius $r=0.1$ fm,
whereas for the PDG16+ we have significantly more resonances and we can
accommodate an excluded volume with a radius up to $r=0.25$ fm. Beyond this
value, we can no longer reproduce the lattice QCD results within error bars. A
quick comment on this result compared to previous papers Gorenstein  _et al._
(2007); Noronha-Hostler  _et al._ (2012) is that here we use the most up-to-
date lattice QCD results, whereas in previous work Gorenstein  _et al._ (2007)
the lattice QCD results were not yet in the continuum limit and were,
therefore, much lower in the confined phase (see e.g., Fig 18 from Borsanyi
_et al._ (2010a)) yielding a larger extracted excluded volume. With current
continuum extrapolated lattice results and the most up-to-date particle list,
we do not find a need for a large excluded volume. Additionally, in Ref.
Noronha-Hostler  _et al._ (2012) Hagedorn resonances were added beyond the
PDG, which required a larger excluded volume. This result still holds: if more
resonances were eventually measured, a larger excluded would become necessary.
While in this work we only explore a single excluded volume, other types of
interactions and/or volumes that vary with the hadron mass or flavor content
Alba  _et al._ (2018b) would lead to different conclusions and may be relevant
at large $\mu_{B}$. We leave this extension for future work.

###  II.2 Shear Viscosity

Figure 3:  Ratio $\eta/s$ as a function of the temperature. The excluded
volume result for the PDG05 list and $r=0.1$ fm (black line) is compared to
the results obtained with the PDG16+ list with $r=0.1$ fm (blue solid line)
and $r=0.25$ fm (blue dashed line).

At vanishing baryon density, an excluded volume description has been used in
Refs. Gorenstein  _et al._ (2007); Noronha-Hostler  _et al._ (2012) to
estimate $\eta(T)$ within the hadronic phase. We extend here this formalism to
finite $\mu_{B}$:

$\eta^{HRG}=\frac{5}{64\sqrt{8}}\frac{1}{r^{2}}\frac{1}{n}\sum\limits_{i}n_{i}\frac{\int\limits_{0}^{\infty}p^{3}e^{-\sqrt{p^{2}+m_{i}^{2}}/T+\tilde{\mu}_{i}}dp}{\int\limits_{0}^{\infty}p^{2}e^{-\sqrt{p^{2}+m_{i}^{2}}/T+\tilde{\mu}_{i}}dp}\,\,,$
(8)

assuming – as mentioned earlier – a single radius for all species. Note that,
in this case, the following relationship between number densities holds:

$\frac{n_{i,v}}{n_{v}}=\frac{n_{i}}{n}\,\,,$ (9)

i.e., one can exchange the excluded volume number density ratio and the ideal
gas number density ratio. If this assumption is relaxed, $n_{i,v}$ and $n_{v}$
must be used in Eq. (8).

After constraining the radius $r$, it is possible to calculate $\eta/s$ as
shown in Fig. 3. Using the PDG05 list, we find that $\eta/s(T)$ is extremely
large and is not remotely able to reach the $1/4\pi$ estimate from the KSS
bound Kovtun  _et al._ (2005). Because of the larger number of states present
in the PDG16+ list, a significant decrease in $\eta/s(T)$ occurs, which brings
it down to $\eta/s(T)\sim 0.8$ at T=165 MeV and $\eta/s(T)\sim 0.34$ at T=199
MeV.

Our next step is to connect to the deconfined phase $\eta^{QGP}/s$ from Refs.
Christiansen  _et al._ (2015); Dubla  _et al._ (2018). In order to avoid
confusion, we first define the following temperatures:

  * •

$T_{min}(\mu_{B})$: Pseudo-critical temperature, i.e., the temperature where
the minimum of $\eta T/w$ is reached for a fixed $\mu_{B}$;

  * •

$T_{sw}$: Switching temperature where the connection between $\eta^{HRG}/s$
and $\eta^{QGP}/s$ occurs. This might correspond to a lower temperature
$T_{sw}\leq T_{min}$;

  * •

$T_{ch}(\mu_{B})$: Chiral transition line temperature, which varies with
$\mu_{B}$. At vanishing chemical potential one has
$T_{ch,0}=T_{ch}(\mu_{B}=0)$;

  * •

$T_{YM}$: a parameter in the definition of $\eta^{QGP}/s$ from Refs.
Christiansen  _et al._ (2015); Dubla  _et al._ (2018);

  * •

$T_{c}$: critical temperature, in the case where there is a critical point
(the critical point has a corresponding $\mu_{B,c}$ as well).

Since the value $\eta/s(T)\sim 0.8$ is still larger than what one would expect
from comparisons of viscous hydrodynamics to experimental data Bernhard  _et
al._ (2016), we propose the following solution:

  * •

Assume a higher pseudo-critical temperature $T_{min}$ for shear viscosity.
This is not an unreasonable assumption in light of Refs. Christiansen  _et
al._ (2015); Rougemont  _et al._ (2017); Dubla  _et al._ (2018); Everett  _et
al._ (2020), wherein the minimum $\eta/s(T)$ is found around $T_{min}\sim 200$
MeV. This leads to a lower $\eta/s(T)$ in the hadronic phase of $\sim 0.34$;

  * •

Renormalize the minimum to $\eta/s(T_{min})=0.08$;

  * •

Utilize the parametrized $\eta/s(T)$ from Ref. Dubla  _et al._ (2018) above
$T_{sw}$, which we also normalize to $\eta/s(T_{min})=0.08$. Then, in the
deconfined phase we have:

$\quad\eta^{QGP}/s=f\left[\frac{a}{\alpha^{\gamma}_{s,HQ}(cT/T_{YM})}+\frac{b}{(T/T_{YM})^{\delta}}\right]\,\,,$
(10)

where the function $\alpha_{s,HQ}(z)$ reads:

$\alpha_{s,HQ}(z)=\frac{1}{\beta_{0}}\frac{z^{2}-1}{z^{2}logz^{2}}\,\,,$ (11)

the constants $a=0.2$, $b=0.15$, $c=0.79$, and the exponents $\delta=5.1$,
$\gamma=1.6$, $\beta_{0}=11-2n_{f}/3$ (where $n_{f}$ is the number of
flavors), and $f$ is an overall normalization constant that we need to ensure
the correct normalization at $T=T_{min}$ and $\mu_{B}=0$.

This allows us to build a $\eta/s(T)$ at $\mu_{B}=0$ that has a temperature
dependence at least motivated by the HRG model (even if the overall magnitude
is not).

We must also implement a finite $\mu_{B}$ dependence, which is somewhat non-
trivial because it depends on the existence and location of a critical point.
Subtle differences appear when only a cross-over or a critical point are
present in the phase diagram. We also must state quite clearly, that we do
_not_ incorporate any critical scaling of the shear viscosity within this
framework. However, since the shear viscosity scales as $\xi^{(4-d)/19}$ where
$d$ is the number of spatial dimensions (in comparison the bulk viscosity
scales as $\xi^{3}$ Monnai  _et al._ (2017)), we argue that criticality would
not considerably affect the shear viscosity. What in turn will affect the
shear viscosity is the exact profile of the line of pseudo-critical
temperatures $T_{min}(\mu_{B})$, i.e. the location of the minimum of $\eta
T/w$ across the phase diagram. For this reason, in the following we consider
very different choices depending on whether only a cross-over or a critical
point are present.

###  II.3 Cross-over

The implementation of a cross-over phase transition is straightforward. We
test two different approaches: a sharp cross-over vs. a smooth cross-over. In
both cases we connect to $\eta^{QGP}/s$.

One of the main purposes of this work is to construct a temperature-dependent
shear viscosity over enthalpy ratio $\eta T/w$ (recall $\eta T/w=\eta
T/(\varepsilon+p)$, which simplifies to $\eta/s$ at $\mu_{B}=0$). In fact, we
would like to construct a shear viscosity that can vary across the entire
4-dimensional QCD phase diagram relevant to heavy-ion collisions, i.e.
spanning over $\left\\{T,\mu_{B},\mu_{S},\mu_{Q}\right\\}$. Unfortunately,
such a model does not exist in the QGP phase beyond Ref. Rougemont  _et al._
(2017), although the authors of Ref. Christiansen  _et al._ (2015) may
eventually incorporate a dependence on the baryon chemical potential, which
would be quite useful. For this reason, in this paper we let the behavior of
the HRG model drive the $\mu_{B}$-dependence of the shear viscosity as follows
(a flow chart of this procedure is shown in Fig. 4). We reserve for future
work the extension to the full $\left\\{\mu_{B},\mu_{S},\mu_{Q}\right\\}$,
which we can anticipate will further influence hydrodynamic trajectories
Monnai  _et al._ (2021).

Figure 4:  Flow chart of our $\eta T/w\left\\{T,\mu_{B}\right\\}$ procedure
for a cross-over.

First, we calculate the excluded volume shear viscosity across
$\left\\{T,\mu_{B}\right\\}$. The enthalpy is then calculated once again in
the HRG model with the mass-dependent volume corrections. With the simple
assumption $\mu_{S}=\mu_{Q}=0$ one finds that, at large $\mu_{B}$, there is a
suppression in $\eta T/w$, which has already been demonstrated in a number of
previous papers Denicol  _et al._ (2013); Kadam and Mishra (2014).

As previously discussed, we do not have a good theory to describe the shear
viscosity of the QGP phase across $\left\\{T,\mu_{B}\right\\}$. Thus, we
simply take the same parameterization form as in Eq. (10), and vary the
overall magnitude as well as the switching temperature with $\mu_{B}$.
Everything remains the same as in Eq. (10), except for $T_{YM}$ and $f$. In
order to determine $T_{YM}(\mu_{B})$ and $f(\mu_{B})$ we need to determine the
temperature at which the matching of the $\eta T/w$ in the confined and
deconfined phases occurs. We use a Taylor expansion for the switching
temperature such that:

$T_{sw}=T_{sw,0}\left(1-\kappa_{2}\Big{(}\frac{\mu_{B}}{T_{ch,0}}\Big{)}^{2}-\kappa_{4}\Big{(}\frac{\mu_{B}}{T_{ch,0}}\Big{)}^{4}\right).$
(12)

Note that this has the same format as what is used in the lattice QCD
community for the chiral phase transition. The values we employ for the chiral
phase transition are Bellwied  _et al._ (2015); Borsanyi  _et al._ (2020):
$T_{sw,0}=T_{ch,0}=158$, $\kappa_{2}=0.0149$, and $\kappa_{4}=0.00032$.
However, there is no reason to assume that the minimum of $\eta T/w$ should
occur exactly at the chiral transition temperature $T_{ch}$; we then study two
different scenarios with a cross-over phase transition. First, we define a
“smooth" cross-over in which we identify $T_{sw}$ exactly with the chiral
transition line, which in turn leads to a flat region at $T_{sw}<T<T_{min}$.
Second, we contrast this with a “sharp" cross-over where we set
$T_{sw}=T_{min}$, which leads to a clear kink in $\eta T/w$ at $T_{min}$.

Transition Parameter Values
---
| Transition | $T_{sw,0}$ | $T_{ch,0}$ | $\kappa_{2}$ | $\kappa_{4}$
CO smooth | Crossover | 158 | 158 | 0.0149 | 0.00032
CO sharp | Crossover | 199 | 158 | 0.0149 | 0.00032
CP(143,350) | Crossover | 199 | 199 | 0.08 | 0.005
First Order | 158 | 158 | 0.0149 | 0.001
CP(89,724) | Crossover | 199 | 199 | 0.038 | 0.0005
First Order | 158 | 158 | 0.0149 | 0.00032
Table 1: The parameters that define the pseudo-critical temperature where the
minimum of $\eta T/w$ occurs, based on Eq. (12.) Figure 5: Temperature
dependence of $\eta T/w(T,\mu_{B})$ for all 4 phase diagrams along slices of
constant $\mu_{B}$

The final results for $\eta T/w(T,\mu_{B})$ are shown in Fig. 5. We can
clearly see that, at vanishing chemical potential, the $\eta T/w$ has a longer
flat region in the case of a “smooth” crossover. This small difference at
$\mu_{B}=0$ has large effects at finite $\mu_{B}$ as well. In fact, the smooth
crossover leads to a much larger $\eta T/w$ across all $\mu_{B}$ values,
because the matching to the HRG occurs at lower temperatures, hence leading to
a strong $\mu_{B}$-dependence. In contrast, matching the HRG model results at
higher temperatures (as in the case of a sharp cross-over) leads to a weaker
$\mu_{B}$-dependence, which in turn significantly decreases $\eta T/w$ at
large $\mu_{B}$. This can be seen even more clearly in Fig. 6, where the
minimum value of $\eta T/w$ is approximately twice as large for the smooth
cross-over than for the sharp one at around $\mu_{B}\sim 600$ MeV. This
indicates that the exact temperature for the switch from the deconfined to a
confined phase has a significant effect on $\eta T/w$. In other words, the
transition line for the minimum of $\eta T/w$, specifically how that varies
with $\mu_{B}$, plays a large role in the overall magnitude of $\eta T/w$.

Figure 6: Minimum $\eta T/w(T,\mu_{B})$ as a function of $\mu_{B}$

###  II.4 Critical point

In the case of a critical point, we need two separate $T_{pc}$ descriptions:
one for the cross-over regime and one for the first-order phase transition
line. At low $\mu_{B}$, in the cross-over phase, we use $T_{0}=199$ MeV from
Eq. (12), but adjust the $\kappa_{2}$ and $\kappa_{4}$ parameters to ensure
that the minimum of $\eta T/w$ hits the critical point. Once the first-order
phase transition is reached, we readjust the parameters in Eq. (12) such that
the first-order phase transition falls onto the chiral transition line
described by $T_{0}=158$ MeV, $\kappa_{2}=0.0149$, and $\kappa_{4}=0.00032$.

The results for our critical point scenarios are shown in Fig. 5. The closer
the critical point lies to $\mu_{B}=0$, the more steeply the minimum of $\eta
T/w$ must drop to reach the critical point. Since the low-T region of
$\eta^{HRG}T/w$ has a stronger dependence on $\mu_{B}$, the magnitude of $\eta
T/w$ must increase dramatically. We conclude that if a critical point is
present at low $\mu_{B}$, the shear viscosity (at least with our setup) must
be quite large. In contrast, if the critical point is far from the $\mu_{B}=0$
axis, such as in Ref. Critelli  _et al._ (2017), then the minimum of $\eta
T/w$ can decrease more gradually across $\mu_{B}$, which in turn leads to a
smaller overall $\eta T/w$ at finite $\mu_{B}$. It is suggestive to compare
the minimum of $\eta T/w$ across $\mu_{B}$ for our two setups with a critical
point, as shown in Fig. 6. We observe that the critical point at low $\mu_{B}$
has a significantly larger minimum $\eta T/w$. The kink in Fig. 6 arises at
the critical point because the transition line changes to the chiral
transition one once the first-order regime is reached.

Finally, in Fig. 7 we show density plots of $\eta T/w$ to illustrate its
different profiles. In the two cross-over setups, one can see clear
differences in the trajectories of the minimum of $\eta T/w$ across the phase
diagram. In the two cases with a critical point (shown as a red dot), one can
clearly see the difference in the transitions lines between the cross-over and
the first order phase transition line.

|
---|---
|
Figure 7: Profile of $\eta T/w$ for our four different scenarios. The cross-
over line is shown in red-dashed, the critical point is a red dot, and the
solid black line is the first-order phase transition line. Values of $\eta
T/w>5$ are shown in white.

##  III Fluctuating Hydrodynamics at $\sqrt{s_{NN}}=19.6$ GeV

Because of event-by-event fluctuations from the initial conditions, large
variations in the local temperature in a quickly expanding and cooling QGP
droplet, and entropy creation due to viscosity, it can be somewhat misleading
to argue that a specific beam energy only probes a single isentropic (where
the ratio of entropy over baryon number $s/\rho_{B}=\text{const}$) trajectory
Dore  _et al._ (2020). Rather the expanding, cooling fireball probes a range
of temperatures over time as shown in Ref. Shen and Schenke (2019).
Furthermore, it has been argued that applying different cuts in rapidity could
provide a method for fine tuning $\mu_{B}$ at a fixed beam energy Brewer  _et
al._ (2018). Thus, in order to determine the $\eta T/w(T,\mu_{B})$ ranges
explored at a set beam energy, we first determine the path across
$\left\\{T,\mu_{B}\right\\}$ taken by AuAu collision with energy
$\sqrt{s_{NN}}=19.6$ GeV.

First, we compare the commonly used isentropes calculated from lattice QCD to
the $T-\mu_{B}$ evolution of a single hydrodynamic event from Ref. Shen and
Schenke (2019) at $\sqrt{s_{NN}}=19.6$ GeV. From hydrodynamics it is possible
to calculate the evolution of the mean temperature $\langle T\rangle$ and
chemical potential $\langle\mu_{B}\rangle$, which is shown in Fig. 8.
Realistically, entropy is produced because of the finite viscosity of the
Quark Gluon Plasma (both bulk and shear viscosity contribute to entropy
production), so one does not expect the relationship
$s/\rho_{B}=\text{const.}$ to actually hold. Close to the phase transition a
minimum of $\eta T/w$ is present, so it is not surprising that exactly at this
point the lattice QCD isentrope most closely matches the hydrodynamic
calculation. It appears that the lattice QCD isentropes can provide reasonable
estimates at the phase transition, but not at higher temperatures.

Figure 8: The temperature and chemical potential dependence of a single
hydrodynamic event from Ref. Shen and Schenke (2019) at $\sqrt{s_{NN}}=19.6$
GeV for the mean temperature and $\mu_{B}$ compared to one standard deviation
from the mean at mid-rapidity $|\eta_{s}|<0.5$. The isentrope at the same beam
energy taken from lattice QCD data is also shown in blue.

Additionally, in Ref. Brewer  _et al._ (2018) it was pointed out that forward
rapidity can be used to fine tune $\mu_{B}$ in order to refine the search for
the critical point. This does indeed appear to be the case, as shown in Fig.
9, so we expect to also see changes in the $\eta T/w(T,\mu_{B})$ ranges
scanned at forward rapidity. In Fig. 9, two different rapidity cuts are shown:
$|\eta_{s}|<0.5$ and $0.5<|\eta_{s}|<1$ (note that rapidity is typically
denoted $\eta_{s}$, so we add the subscript to indicate spacetime rapidity and
avoid confusion with viscosity).

Figure 9: The temperature and chemical potential dependence of a single
hydrodynamic event from Ref. Shen and Schenke (2019) at $\sqrt{s_{NN}}=19.6$
GeV for the mean temperature and $\mu_{B}$, compared to one standard deviation
from the mean at mid-rapidity $|\eta_{s}|<0.5$ (black) and $0.5<\eta_{s}<1$
(red).

Finally, in Fig. 10 we compare the values of $\eta T/w$ probed across $\pm 1$
standard deviations from the mean $T$ and $\mu_{B}$ at the rapidity cuts
$|\eta_{s}|<0.5$ and $0.5<|\eta_{s}|<1$. Here we only compare our setups with
a crossover and that with the critical point from Ref. Parotto  _et al._
(2020). We find that, because the average $\mu_{B}$ is still relatively low
(approximately $\mu_{B}\sim 200-300$ MeV), the location of the minimum $\eta
T/w$ is relatively compatible between the different cases. Because of the
longer flat region, the smooth cross-over actually sees a smaller $\eta T/w$
over a longer period of time. In contrast, in the case with a critical point
we observe a larger $\eta T/w$. We find that, at forward rapidity, the Quark
Gluon Plasma sees a slightly higher $\eta T/w$ at around $\tau=1-4$ fm. Thus
we expect that, while forward rapidity does reach large $\mu_{B}$, we also
anticipate that it would be farther from equilibrium.

Figure 10: Approximation of $\eta T/w$ probed within relativistic viscous
hydrodynamics vs. time (using the $\left(T,\mu_{B}\right)$ from a single
hydrodynamic background from Ref. Shen and Schenke (2019))

.

##  IV Influence on $T,\mu_{B}$ trajectories

In this section we give a qualitative idea of the effect of using different
$\eta T/w$ profiles on the trajectories taken through the QCD phase diagram by
the expanding hydrodynamic system. Here, we only explore the effect of a
crossover transition. The hydrodynamic setup is identical to Ref. Dore  _et
al._ (2020), other than a replacement of the EoS with the one developed in
Ref. Noronha-Hostler  _et al._ (2019) with $\mu_{S}=\mu_{Q}=0$. This is done
in order to compare the effect specifically with a crossover transition. The
EoS used here is directly obtained from lattice QCD through a Taylor expansion
of the pressure in powers of $\mu_{B}/T$. The bulk viscosity has the form

$\frac{\zeta T}{w}=36\times\frac{1/3-c_{s}^{2}}{8\pi}\,\,,$ (13)

so that all of its functional dependence on thermodynamic variables enters
through the speed of sound. The hydrodynamic equations exhibit Bjorken
symmetry and are readily solved as a set of coupled ODE’s. We compare the
sharp vs. smooth cross-over form of $\eta T/w$ to either a constant $\eta
T/w=0.08$ or an ideal hydrodynamic expansion. For an ideal fluid, entropy is
conserved so one can describe that expansion through isentropic trajectories.

Figure 11: Comparison of $\left\\{T,\mu_{B}\right\\}$ trajectories with
various initial conditions for $s/\rho_{B}$ with either $\eta T/w=0$ (ideal
hydrodynamics), $\eta T/w=0.08$, or $\eta T/w$ following our “sharp” or
“smooth” cross-over transitions.

The effect on the trajectories can be seen in Fig. 11, where they all begin at
a temperature of $T_{i}=500$ MeV, and a chemical potential of
$\mu_{i}=[200,400,800]$ MeV. The green lines are isentropes produced by
running ideal hydrodynamics given the same initial conditions. The results are
consistent with what is shown in Fig. 5. At lower T and $\mu_{B}$, the smooth
and sharp shear viscosities behave very similarly, and one would expect the
associated hydrodynamic trajectories to look similar if initialized in this
region. Hydrodynamic runs that begin at large $\mu_{B}$, on the other hand,
are more strongly influenced by the actual profile of $\eta T/w$, which leads
to different trajectoriesDore  _et al._ (2020).

We note that the effect of $\eta T/w$ is mainly to shift the trajectories
towards $\mu_{B}=0$, such that a very small initial $s/\rho_{B}$ is necessary
in order to see dramatic effects around the phase transition. For example, we
reach here a similar freeze-out point with the initial condition
$s/\rho_{B}=13$ for both crossover $\eta T/w$ profiles, as the hydrodynamic
simulation shown in Fig. 8 where the isentropic trajectory corresponded to
$s/\rho_{B}=48.4$. This is not surprising, because ideal hydrodynamics has
maximum entropy, whereas a viscous system produces entropy over time, so that
entropy production drives the system from a smaller to a larger $s/\rho_{B}$.

We also compare to a constant $\eta T/w=0.08$ to show the overall effect of
having a dependence on $T$ and $\mu_{B}$. It can be seen in Fig. 11 that, the
higher the initial chemical potentials, the more the associated trajectories
deviate from each other. It becomes apparent that, at higher $\mu_{B}$, one
becomes more sensitive to the $T$ and $\mu_{B}$ dependence of $\eta T/w$ as
well. Nonetheless, any viscosity (even constant $\eta T/w=0.08$) causes
substantial deviations from the isentropic lines. We conclude that heavy-ion
collisions at low beam energies are more affected by the $T$ and $\mu_{B}$
dependence of $\eta T/w$ and we anticipate large deviations from isentropic
trajectories, especially if the $\eta T/w$ profile deviates substantially from
a constant value. There are some indications of this from Bayesian analyses
Auvinen  _et al._ (2018), which imply that extreme caution should be exercised
about any conclusions drawn from work based on ideal hydrodynamics at large
baryon densities (see Ref. Dexheimer  _et al._ (2020) for further discussion).

##  V Conclusions

We used a hadron resonance gas picture with a state-of-the-art list of
resonances to study the QCD shear viscosity at large baryon densities. We
developed a phenomenological approach to produce profiles of $\eta
T/w\left(T,\mu_{B}\right)$ on the QCD phase diagram, which can be used in
relativistic viscous hydrodynamic codes simulating collisions at the RHIC Beam
Energy Scan. We considered four variations: a “smooth” or “sharp” crossover,
and a critical point at $\\{T_{c},\mu_{B}^{c}\\}=\\{143,350\\}$ MeV or
$\\{T_{c},\mu_{B}^{c}\\}=\\{89,725\\}$ MeV. We have demonstrated that the
trajectories followed by hydrodynamic evolution are sensitive to different
behavior in the shear viscosity. It is clear that a shear viscosity with a
functional dependence on thermodynamic variables may probe different areas of
the phase diagram than a constant one, especially if the differences are large
at early times. We note that this framework does not incorporate critical
scaling in $\eta T/w\left(T,\mu_{B}\right)$, since this is expected to be
extremely small Monnai  _et al._ (2017), but does ensure that its inflection
point (i.e. minimum of $\eta T/w\left(T,\mu_{B}\right)$) does converge to the
critical point.

Using this framework we found that the minimum $\eta
T/w\left(T,\mu_{B}\right)$ increases with $\mu_{B}$ and that this effect is
strongly influenced by the exact location of the transition line one chooses
for the minimum of $\eta T/w\left(T,\mu_{B}\right)$. For instance, a critical
point at low $\mu_{B}$ would dramatically increase the magnitude of $\eta T/w$
compared to a critical point at larger $\mu_{B}$. We showed that the
trajectories from heavy-ion collisions that pass through the QCD phase diagram
would have significant deviations from isentropic lines, the larger the
chemical potential is. The effect of viscosity is mainly to push the
trajectories towards lower $\mu_{B}$ at freeze-out. This implies that a
significantly smaller value of $s/\rho_{B}$ is needed as an initial condition,
in order to reach the same point of freeze-out compared to an ideal
hydrodynamics system. We only consider a cross-over phase transition for this
particular part of the study but the presence of a critical point would only
enhance this effect due to the critical scaling of bulk viscosity, as shown
already in Ref. Dore  _et al._ (2020). We also expect early times ($\tau\sim
1-4$ fm/c) to be the regime probing the minimum of $\eta T/w$.

##  VI Acknowledgements

This material is based upon work supported by the National Science Foundation
under grant no. PHY-1654219 and by the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics, within the framework of the Beam Energy
Scan Theory (BEST) Topical Collaboration. We also acknowledge the support from
the Center of Advanced Computing and Data Systems at the University of
Houston. J.N.H. acknowledges support from the US-DOE Nuclear Science Grant No.
DE-SC0020633 and from the Illinois Campus Cluster, a computing resource that
is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with
the National Center for Supercomputing Applications (NCSA), and which is
supported by funds from the University of Illinois at Urbana-Champaign. E.M.
was supported by the National Science Foundation via grant PHY-1560077. P.P.
acknowledges support by the DFG grant SFB/TR55.

## References

  * Aoki  _et al._ (2006) Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, Phys. Lett. B 643, 46 (2006), arXiv:hep-lat/0609068 .
  * Borsanyi  _et al._ (2010a) S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, JHEP 11, 077 (2010a), arXiv:1007.2580 [hep-lat] .
  * Borsanyi  _et al._ (2010b) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo (Wuppertal-Budapest), JHEP 09, 073 (2010b), arXiv:1005.3508 [hep-lat] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Borsanyi  _et al._ (2020) S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 125, 052001 (2020), arXiv:2002.02821 [hep-lat] .
  * Ratti (2018) C. Ratti, Rept. Prog. Phys. 81, 084301 (2018), arXiv:1804.07810 [hep-lat] .
  * Guenther (2020) J. N. Guenther, (2020), arXiv:2010.15503 [hep-lat] .
  * Endrodi  _et al._ (2011) G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, JHEP 04, 001 (2011), arXiv:1102.1356 [hep-lat] .
  * Rougemont  _et al._ (2017) R. Rougemont, R. Critelli, J. Noronha-Hostler, J. Noronha, and C. Ratti, Phys. Rev. D 96, 014032 (2017), arXiv:1704.05558 [hep-ph] .
  * Grefa  _et al._ (2021) J. Grefa, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, (2021), arXiv:2102.12042 [nucl-th] .
  * Bellwied  _et al._ (2019) R. Bellwied, J. Noronha-Hostler, P. Parotto, I. Portillo Vazquez, C. Ratti, and J. M. Stafford, Phys. Rev. C 99, 034912 (2019), arXiv:1805.00088 [hep-ph] .
  * Noronha-Hostler  _et al._ (2016a) J. Noronha-Hostler, R. Bellwied, J. Gunther, P. Parotto, A. Pasztor, I. Portillo Vazquez, and C. Ratti, (2016a), arXiv:1607.02527 [hep-ph] .
  * Kovtun  _et al._ (2005) P. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005), arXiv:hep-th/0405231 .
  * Kats and Petrov (2009) Y. Kats and P. Petrov, JHEP 01, 044 (2009), arXiv:0712.0743 [hep-th] .
  * Brigante  _et al._ (2008a) M. Brigante, H. Liu, R. C. Myers, S. Shenker, and S. Yaida, Phys. Rev. D 77, 126006 (2008a), arXiv:0712.0805 [hep-th] .
  * Brigante  _et al._ (2008b) M. Brigante, H. Liu, R. C. Myers, S. Shenker, and S. Yaida, Phys. Rev. Lett. 100, 191601 (2008b), arXiv:0802.3318 [hep-th] .
  * Buchel  _et al._ (2009) A. Buchel, R. C. Myers, and A. Sinha, JHEP 03, 084 (2009), arXiv:0812.2521 [hep-th] .
  * Critelli  _et al._ (2014) R. Critelli, S. I. Finazzo, M. Zaniboni, and J. Noronha, Phys. Rev. D 90, 066006 (2014), arXiv:1406.6019 [hep-th] .
  * Finazzo  _et al._ (2016) S. I. Finazzo, R. Critelli, R. Rougemont, and J. Noronha, Phys. Rev. D 94, 054020 (2016), [Erratum: Phys.Rev.D 96, 019903 (2017)], arXiv:1605.06061 [hep-ph] .
  * Romatschke and Romatschke (2007) P. Romatschke and U. Romatschke, Phys. Rev. Lett. 99, 172301 (2007), arXiv:0706.1522 [nucl-th] .
  * Bozek (2012) P. Bozek, Phys. Rev. C 85, 034901 (2012), arXiv:1110.6742 [nucl-th] .
  * Heinz and Snellings (2013) U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Luzum and Petersen (2014) M. Luzum and H. Petersen, J. Phys. G 41, 063102 (2014), arXiv:1312.5503 [nucl-th] .
  * Niemi  _et al._ (2016) H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph] .
  * Noronha-Hostler  _et al._ (2016b) J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 93, 034912 (2016b), arXiv:1511.06289 [nucl-th] .
  * McDonald  _et al._ (2017) S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Phys. Rev. C 95, 064913 (2017), arXiv:1609.02958 [hep-ph] .
  * Bernhard  _et al._ (2019) J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nature Phys. 15, 1113 (2019).
  * Alba  _et al._ (2018a) P. Alba, V. Mantovani Sarti, J. Noronha, J. Noronha-Hostler, P. Parotto, I. Portillo Vazquez, and C. Ratti, Phys. Rev. C 98, 034909 (2018a), arXiv:1711.05207 [nucl-th] .
  * Aarts and Martinez Resco (2002) G. Aarts and J. M. Martinez Resco, JHEP 04, 053 (2002), arXiv:hep-ph/0203177 .
  * Borsanyi  _et al._ (2018) S. Borsanyi, Z. Fodor, M. Giordano, S. D. Katz, A. Pasztor, C. Ratti, A. Schaefer, K. K. Szabo, and C. Toth, Balint, Phys. Rev. D 98, 014512 (2018), arXiv:1802.07718 [hep-lat] .
  * Noronha-Hostler  _et al._ (2009) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. Lett. 103, 172302 (2009), arXiv:0811.1571 [nucl-th] .
  * Pal (2010) S. Pal, Phys. Lett. B 684, 211 (2010), arXiv:1001.1585 [nucl-th] .
  * Khvorostukhin  _et al._ (2010) A. S. Khvorostukhin, V. D. Toneev, and D. N. Voskresensky, Nucl. Phys. A 845, 106 (2010), arXiv:1003.3531 [nucl-th] .
  * Tawfik and Wahba (2010) A. Tawfik and M. Wahba, Annalen Phys. 522, 849 (2010), arXiv:1005.3946 [hep-ph] .
  * Alba  _et al._ (2015) P. Alba, R. Bellwied, M. Bluhm, V. Mantovani Sarti, M. Nahrgang, and C. Ratti, Phys. Rev. C 92, 064910 (2015), arXiv:1504.03262 [hep-ph] .
  * Ratti  _et al._ (2011) C. Ratti, S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo (Wuppertal-Budapest), Nucl. Phys. A 855, 253 (2011), arXiv:1012.5215 [hep-lat] .
  * Tiwari  _et al._ (2012) S. K. Tiwari, P. K. Srivastava, and C. P. Singh, Phys. Rev. C 85, 014908 (2012), arXiv:1111.2406 [hep-ph] .
  * Noronha-Hostler  _et al._ (2012) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. C 86, 024913 (2012), arXiv:1206.5138 [nucl-th] .
  * Kadam and Mishra (2014) G. P. Kadam and H. Mishra, Nucl. Phys. A 934, 133 (2014), arXiv:1408.6329 [hep-ph] .
  * Kadam and Mishra (2015) G. P. Kadam and H. Mishra, Phys. Rev. C 92, 035203 (2015), arXiv:1506.04613 [hep-ph] .
  * Kadam and Mishra (2016) G. P. Kadam and H. Mishra, Phys. Rev. C 93, 025205 (2016), arXiv:1509.06998 [hep-ph] .
  * Kadam and Pawar (2019) G. Kadam and S. Pawar, Adv. High Energy Phys. 2019, 6795041 (2019), arXiv:1802.01942 [hep-ph] .
  * Mohapatra  _et al._ (2019) R. K. Mohapatra, H. Mishra, S. Dash, and B. K. Nandi, (2019), arXiv:1901.07238 [hep-ph] .
  * Wesp  _et al._ (2011) C. Wesp, A. El, F. Reining, Z. Xu, I. Bouras, and C. Greiner, Phys. Rev. C 84, 054911 (2011), arXiv:1106.4306 [hep-ph] .
  * Ozvenchuk  _et al._ (2013) V. Ozvenchuk, O. Linnyk, M. I. Gorenstein, E. L. Bratkovskaya, and W. Cassing, Phys. Rev. C 87, 064903 (2013), arXiv:1212.5393 [hep-ph] .
  * Rose  _et al._ (2018) J. B. Rose, J. M. Torres-Rincon, A. Schäfer, D. R. Oliinychenko, and H. Petersen, Phys. Rev. C 97, 055204 (2018), arXiv:1709.03826 [nucl-th] .
  * Rais  _et al._ (2020) J. Rais, K. Gallmeister, and C. Greiner, Phys. Rev. D 102, 036009 (2020), arXiv:1909.04522 [hep-ph] .
  * Sahu and Sahoo (2020) D. Sahu and R. Sahoo, (2020), arXiv:2006.04185 [hep-ph] .
  * Heffernan  _et al._ (2020) M. Heffernan, S. Jeon, and C. Gale, Phys. Rev. C 102, 034906 (2020), arXiv:2005.12793 [hep-ph] .
  * Haas  _et al._ (2014) M. Haas, L. Fister, and J. M. Pawlowski, Phys. Rev. D 90, 091501 (2014), arXiv:1308.4960 [hep-ph] .
  * Christiansen  _et al._ (2015) N. Christiansen, M. Haas, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. Lett. 115, 112002 (2015), arXiv:1411.7986 [hep-ph] .
  * Dubla  _et al._ (2018) A. Dubla, S. Masciocchi, J. M. Pawlowski, B. Schenke, C. Shen, and J. Stachel, Nucl. Phys. A 979, 251 (2018), arXiv:1805.02985 [nucl-th] .
  * Ghiglieri  _et al._ (2018) J. Ghiglieri, G. D. Moore, and D. Teaney, JHEP 03, 179 (2018), arXiv:1802.09535 [hep-ph] .
  * Denicol  _et al._ (2013) G. S. Denicol, C. Gale, S. Jeon, and J. Noronha, Phys. Rev. C 88, 064901 (2013), arXiv:1308.1923 [nucl-th] .
  * Demir and Bass (2009) N. Demir and S. A. Bass, Phys. Rev. Lett. 102, 172302 (2009), arXiv:0812.2422 [nucl-th] .
  * Critelli  _et al._ (2017) R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 96, 096026 (2017), arXiv:1706.00455 [nucl-th] .
  * Vovchenko  _et al._ (2017a) V. Vovchenko, M. I. Gorenstein, and H. Stoecker, Phys. Rev. Lett. 118, 182301 (2017a), arXiv:1609.03975 [hep-ph] .
  * Huovinen and Petreczky (2018) P. Huovinen and P. Petreczky, Phys. Lett. B 777, 125 (2018), arXiv:1708.00879 [hep-ph] .
  * Alba  _et al._ (2017) P. Alba _et al._ , Phys. Rev. D 96, 034517 (2017), arXiv:1702.01113 [hep-lat] .
  * Bellwied  _et al._ (2020) R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther, J. Noronha-Hostler, P. Parotto, A. Pasztor, C. Ratti, and J. M. Stafford, Phys. Rev. D 101, 034506 (2020), arXiv:1910.14592 [hep-lat] .
  * Alba  _et al._ (2020) P. Alba, V. M. Sarti, J. Noronha-Hostler, P. Parotto, I. Portillo-Vazquez, C. Ratti, and J. M. Stafford, Phys. Rev. C 101, 054905 (2020), arXiv:2002.12395 [hep-ph] .
  * Parotto  _et al._ (2020) P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
  * Rischke  _et al._ (1991) D. H. Rischke, M. I. Gorenstein, H. Stoecker, and W. Greiner, Z. Phys. C 51, 485 (1991).
  * Albright  _et al._ (2014) M. Albright, J. Kapusta, and C. Young, Phys. Rev. C 90, 024915 (2014), arXiv:1404.7540 [nucl-th] .
  * Vovchenko  _et al._ (2017b) V. Vovchenko, A. Motornenko, P. Alba, M. I. Gorenstein, L. M. Satarov, and H. Stoecker, Phys. Rev. C 96, 045202 (2017b), arXiv:1707.09215 [nucl-th] .
  * Gorenstein  _et al._ (2007) M. I. Gorenstein, M. Hauer, and O. N. Moroz,  , 214 (2007), arXiv:0708.0137 [nucl-th] .
  * Alba  _et al._ (2018b) P. Alba, V. Vovchenko, M. I. Gorenstein, and H. Stoecker, Nucl. Phys. A 974, 22 (2018b), arXiv:1606.06542 [hep-ph] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Everett  _et al._ (2020) D. Everett _et al._ (JETSCAPE), (2020), arXiv:2011.01430 [hep-ph] .
  * Monnai  _et al._ (2017) A. Monnai, S. Mukherjee, and Y. Yin, Phys. Rev. C 95, 034902 (2017), arXiv:1606.00771 [nucl-th] .
  * Monnai  _et al._ (2021) A. Monnai, B. Schenke, and C. Shen, (2021), arXiv:2101.11591 [nucl-th] .
  * Bellwied  _et al._ (2015) R. Bellwied, S. Borsanyi, Z. Fodor, J. Günther, S. D. Katz, C. Ratti, and K. K. Szabo, Phys. Lett. B 751, 559 (2015), arXiv:1507.07510 [hep-lat] .
  * Dore  _et al._ (2020) T. Dore, J. Noronha-Hostler, and E. McLaughlin, Phys. Rev. D 102, 074017 (2020), arXiv:2007.15083 [nucl-th] .
  * Shen and Schenke (2019) C. Shen and B. Schenke, Nucl. Phys. A 982, 411 (2019), arXiv:1807.05141 [nucl-th] .
  * Brewer  _et al._ (2018) J. Brewer, S. Mukherjee, K. Rajagopal, and Y. Yin, Phys. Rev. C 98, 061901 (2018), arXiv:1804.10215 [hep-ph] .
  * Noronha-Hostler  _et al._ (2019) J. Noronha-Hostler, P. Parotto, C. Ratti, and J. M. Stafford, Phys. Rev. C 100, 064910 (2019), arXiv:1902.06723 [hep-ph] .
  * Auvinen  _et al._ (2018) J. Auvinen, J. E. Bernhard, S. A. Bass, and I. Karpenko, Phys. Rev. C 97, 044905 (2018), arXiv:1706.03666 [hep-ph] .
  * Dexheimer  _et al._ (2020) V. Dexheimer, J. Noronha, J. Noronha-Hostler, C. Ratti, and N. Yunes, (2020), arXiv:2010.08834 [nucl-th] .
  * Bazavov  _et al._ (2019) A. Bazavov _et al._ (HotQCD), Phys. Lett. B 795, 15 (2019), arXiv:1812.08235 [hep-lat] .