# Joule heating effects in high transparency Josephson junctions

Matti Tomi  QTF Centre of Excellence, Department of Applied Physics, Aalto
University, P.O. Box 15100, FI-00076 Aalto, Finland  Mikhail R. Samatov  HSE
University, 101000 Moscow, Russia  Andrey S. Vasenko  HSE University, 101000
Moscow, Russia  I.E. Tamm Department of Theoretical Physics, P.N. Lebedev
Physical Institute, Russian Academy of Sciences, 119991 Moscow, Russia  Antti
Laitinen  QTF Centre of Excellence, Department of Applied Physics, Aalto
University, P.O. Box 15100, FI-00076 Aalto, Finland  Pertti Hakonen  QTF
Centre of Excellence, Department of Applied Physics, Aalto University, P.O.
Box 15100, FI-00076 Aalto, Finland  Low Temperature Laboratory, Department of
Applied Physics, Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland
Dmitry S. Golubev  QTF Centre of Excellence, Department of Applied Physics,
Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland

###### Abstract

We study, both theoretically and experimentally, the features on the current-
voltage characteristic of a highly transparent Josephson junction caused by
transition of the superconducting leads to the normal state. These features
appear due to the suppression of the Andreev excess current. We show that by
tracing the dependence of the voltage, at which the transition occurs, on the
bath temperature and by analyzing the suppression of the excess current by the
bias voltage one can recover the temperature dependence of the heat flow out
of the junction. We verify theory predictions by fabricating two highly
transparent superconductor-graphene-superconductor (SGS) Josephson junctions
with suspended and non-suspended graphene as a non-superconducting section
between Al leads. Applying the above mentioned technique we show that the
cooling power of the suspended junction depends on the bath temperature as
$\propto T_{\rm bath}^{3.1}$ close to the superconducting critical
temperature.

###### pacs:

##  I Introduction

Highly transparent Josephson junctions possess unique properties
distinguishing them from the usual low-transparency junctions, in which two
superconducting leads are separated by a tunnel barrier. For example,
transparent junctions have a non-sinusoidal current-phase relation [1, 2] and
host Andreev levels, which can be probed by microwave spectroscopy [3] and can
be potentially used for quantum information processing [4]. These junctions
also exhibit a characteristic pattern of multiple Andreev reflection (MAR)
peaks in the differential conductance at sub-gap bias voltages $eV<2\Delta$,
where $\Delta$ is the superconducting gap [5, 6, 7, 8, 9]. Two superconducting
leads of a transparent Josephson junction are usually connected by a short
section made of a non-superconducting material, which forms good electric
contacts to them. A variety of materials have been used for this purpose:
carbon nanotubes [10, 11, 12], graphene [13, 14], InAs nanowires [15, 16, 3],
2d [17] and 3d [18, 19] topological insulators, etc.

Here we consider yet another characteristic feature of highly transparent
Josephson junctions – the excess current. More specifically, we investigate,
both theoretically and experimentally, the suppression of the excess current
by Joule heating. It is known that the I-V curve of a transparent junction at
high bias voltage approaches the asymptotic form $I=V/R+I_{\rm exc}$, where
$R$ is the junction resistance in the normal state and $I_{\rm exc}$ is the
excess current. This current is proportional to the value of the gap $\Delta$
in the leads and originates from Andreev reflection [20]. Joule heating of the
superconducting leads suppresses the gap $\Delta$. At some bias point the
temperature of the leads reaches the critical temperature of the
superconducting phase transition $T_{C}$, and both the gap and the excess
current vanish. This bias point is evidenced by a dip in the differential
conductance of the junction. Such dips are often observed in Josephson
junctions, see e.g. Refs. [21, 22, 9]. They have been systematically
investigated by Choi et al [23] in highly transparent superconductor -
graphene - superconductor (SGS) Josephson junctions on silicon oxide
substrate, who have made an important observation that the Joule heating power
$P_{C}$, at which the leads become normal, does not depend on the gate voltage
controlling the junction resistance. Similar dips in $dI/dV$ are also observed
in normal metal - superconductor (NS) junctions [24, 25]. Positions of the
dips in the differential conductance provide information about temperature
dependence of the cooling power, i.e. the heat flow out of the junction, which
is important for bolometry applications. For example, by tracing dip positions
at different bath temperatures $T_{\rm bath}$, Choi et al have shown that in
an SGS junction with Al/Ti leads on the Si/SiO2 substrate $P_{C}$
approximately scales as $P_{C}\propto T_{C}^{4}-T_{\rm bath}^{4}$.

In spite of their frequent observation, detailed theoretical description of
the dips in $dI/dV$ described above is still lacking. Historically, several
models have been proposed in literature to explain their origin: multiple
reflections of quasiparticles between the bulk of the superconducting lead and
the contact area [21], suppression of superconductivity in the leads by high
current density [22] and by magnetic field induced by the current [24, 25]. At
present, it is understood that Joule heating is the main cause of the above-
gap features in junctions with the resistance exceeding that of the leads. The
effect of Joule heating on the positions of the dips has been discussed, for
example, in Refs. [24, 23]. In a broader context, it has been shown that Joule
heating suppresses the gap in the leads of a the superconductor - normal metal
- superconductor (SNS) junction [26], induces hysteresis in the current -
voltage characteristics of SNS junctions [27] and of superconducting nanowires
[28], etc.

Here we develop a theoretical model of Joule heating in high transparency
Josephson junctions and derive a simple analytical expression for the
differential conductance in the vicinity of a high bias dip. We demonstrate
that the critical power $P_{C}$ predominantly depends on the properties of the
leads and the substrate rather than on the junction resistance. We also test
the theory predictions experimentally by studying highly transparent SGS
Josephson junctions with aluminum leads. We show, in particular, that for a
suspended SGS junction the cooling power scales as $P\propto
T_{0}^{3.1}-T_{\rm bath}^{3.1}$, where $T_{0}$ is the temperature of the
superconducting leads at the contacts with graphene. Graphene is known to be a
very promising material for bolometry [29, 30, 31, 32, 33, 34, 35, 36].
Recently, bolometers based on SGS junctions similar to ours and having the
noise equivalent power at the level $10^{-20}-10^{-18}$ W/$\sqrt{\rm Hz}$ in
the microwave frequency range have been demonstrated [37, 38]. We argue that
the technique based on the suppression of the excess current by Joule heating
may provide an additional tool for calibration of such bolometers because it
allows one to obtain the dependence of the cooling power on temperature for a
given device in a certain temperature range.

The paper is organized as follows: in Sec. II we present the theoretical
model, in Sec. III we discuss the experiment and in Sec. IV we summarize our
results.

##  II Model

In this section we present the theoretical model. We assume that the bias
voltage $V$ applied to the junction is positive and sufficiently large,
$eV>2\Delta$. In this limit the I-V curve of a transparent Josephson junction
can be well approximated as

$I(V)=\frac{V}{R}+\alpha\frac{\Delta_{L}(T_{L})+\Delta_{R}(T_{R})}{2eR}.$ (1)

Here the second term in the right hand side is the excess current and
$\Delta_{L}(T_{L})$ and $\Delta_{R}(T_{R})$ are the superconducting gaps in
the left and the right leads. They depend on the temperatures of the leads at
the contacts with the central section, $T_{L}$ and $T_{R}$. For a short
Josephson junction with the distance between the leads, $L_{\rm jct}$, much
shorter than the superconducting coherence length $\xi$, i.e. for $L_{\rm
jct}\ll\xi$, the pre-factor $\alpha$ is expressed as [39]

$\alpha=\frac{\sum_{j}\frac{\tau_{j}^{2}}{1-\tau_{j}}\bigg{(}1-\frac{\tau_{j}^{2}}{2\sqrt{1-\tau_{j}}(2-\tau_{j})}\ln{\frac{1+\sqrt{1-\tau_{j}}}{1-\sqrt{1-\tau_{j}}}}\bigg{)}}{\sum_{j}\tau_{j}}.$
(2)

Here $\tau_{j}$ is the transmission probability of $j$-th conducting channel
and the sums run over all channels. The parameter $\alpha$ varies between 0
and $8/3$ and approaches its maximum value in a highly transparent junction
with channel transmissions close to 1. In a tunnel junction with very small
transmissions $\tau_{j}\ll 1$ one finds $\alpha\ll 1$ and the excess current
in Eq. (1) vanishes. Two other important cases - a normal metal diffusive wire
(SNS junction) and a clean short and wide strip of graphene - are
characterized by Dorokhov’s distribution of transmission probabilities [42,
43], which results in $\alpha={\pi^{2}}/{4}-1\approx 1.4674$. Finally, in a
junction having a long diffusive section with $L_{\rm jct}>\xi$ the pre-factor
$\alpha$ decreases with the length as [44]

$\displaystyle\alpha=2.47\,{\xi}/{L_{\rm jct}}.$ (3)

We first consider a simple case of two identical superconducting leads, which
are symmetrically coupled to the central normal section. Accordingly, we put
$T_{L}=T_{R}=T_{0}$ and $\Delta_{L}(T_{L})=\Delta_{R}(T_{R})=\Delta(T_{0})$.
In order to find the dependence of the temperature $T_{0}$ on the bias voltage
$V$ we need to solve the heat balance equation

$\displaystyle P(T_{0},T_{\rm bath})=IV,$ (4)

where we have introduced the bath temperature $T_{\rm bath}$ and the cooling
power of the junction $P(T_{0},T_{\rm bath})$, which is, in general, unknown.
The Joule heating power is partially carried away to the superconducting
leads, while the rest of it goes to the substrate phonons and possibly to some
other cooling channels, see Fig. 1(a). Here we assume that the junction
resistance is much larger than the resistances of the leads and, hence, the
Josephson critical current of the junction is much lower than the depairing
current in the leads. In this case, at certain bias voltage $V_{C}$ the
contact temperature becomes equal to the critical temperature of the
superconductor, $T_{0}=T_{C}$. At this point the differential conductance has
a dip. Tracing the position of this dip at different bath temperatures, one
can restore the cooling power $P(T_{C},T_{\rm bath})=I(V_{C})V_{C}$ as a
function of $T_{\rm bath}$.

Figure 1: (a) High transparency Josephson junction with normal-conductor
section connecting the two identical superconducting leads, fully residing in
contact with the substrate. Part of the Joule heat goes to the substrate
($P_{\rm substrate}$) and the rest of it – to the left ($P_{SL}$) and to the
right ($P_{SR}$) superconducting leads. $T_{L}$ and $T_{R}$ are the
temperatures of the superconducting leads at the contacts with the central
section. (b) Josephson junction with a suspended central part, in which the
Joule heat escapes only through the superconducting leads.

We now find the shape of the dip close to the critical voltage value $V_{C}$.
In the vicinity of this bias point one can linearize Eq. (4) writing it in the
form

$\displaystyle(T_{0}-T_{C})G_{C}^{\rm th}=IV-I(V_{C})V_{C}.$ (5)

Here $G_{C}^{\rm th}=\partial P(T_{0},T_{\rm bath})/\partial
T_{0}|_{T_{0}=T_{C}}$ is the thermal conductance of the leads at the critical
temperature. According to the theory by Bardeen, Cooper and Schrieffer (BCS)
[45], in the vicinity of the critical temperature the superconducting gap
behaves as

$\displaystyle\Delta(T_{0})=k_{B}T_{C}\sqrt{\frac{8\pi^{2}}{7\zeta(3)}\frac{T_{C}-T_{0}}{T_{C}}},$
(6)

where $\zeta(x)$ denotes Riemann’s zeta function. Combining this expression
with the expression for the I-V curve (1), we transform Eq. (5) to the form

$\displaystyle(T_{0}-T_{C})G_{C}^{\rm
th}=\frac{V^{2}-V_{C}^{2}}{R}+\frac{k_{B}T_{C}V}{eR}\sqrt{\frac{8\pi^{2}\alpha^{2}}{7\zeta(3)}\frac{T_{C}-T_{0}}{T_{C}}}.$

Solving this equation for $T_{0}$, substituting the result in the expression
for the current (1) and taking the derivative of it over $V$, we obtain a
simple analytical expression for the differential conductance

$\displaystyle\frac{dI}{dV}=\frac{1}{R}-\frac{4\theta(V_{C}-V)}{bR}\left(\frac{(b-4)V}{\sqrt{bV_{C}^{2}-(b-4)V^{2}}}+2\right).$
(7)

Here $\theta(x)$ is the Heaviside step function and we have introduced the
dimensionless parameter

$\displaystyle b=\frac{14\zeta(3)}{\pi^{2}\alpha^{2}}\frac{e^{2}RG_{C}^{\rm
th}}{k_{B}^{2}T_{C}}.$ (8)

The approximation (7) is formally valid if $b\gg 1$ and
$|V_{C}-V|\lesssim\alpha k_{B}T_{C}/e$. However, Eq. (7) can also be used
outside this voltage interval because the non-linear correction to $dI/dV$
almost vanishes there. Eq. (7) predicts strong dip in the differential
conductance with the negative minimum value, $-1/R$, which is achieved at
$V=V_{C}$. In the experiment, the dip is smeared by the inhomogeneity of the
hot spot in the contact area between the leads and the central section, by the
inverse proximity effect suppressing the gap close to the contact, etc.

Figure 2: (a) Dependence of the contact temperature $T_{0}$ on the bias
voltage. (b) Solid line is the I-V curve of the junction, red dashed line
shows Ohm’s law $I=V/R$, and blue line – Ohm’s law with added excess current,
$I=V/R+\alpha\Delta(0)/eR$. The superconducting branch at $V=0$ is not shown.
(c) Solid line shows the differential conductance obtianed by solving Eq. (12)
numerically, and red dashed line represents the approximate expression (7).
The following parameters have been chosen for this simulation: $R=100$
$\Omega$, $R_{L_{S}}=2$ $\Omega$, $T_{C}=1.2$ K, $\Delta(0)=182$ $\mu$eV,
$\alpha=0.5$, $T_{\rm bath}=10$ mK.

In reality the two superconducting leads of the junction are never fully
identical. They may have different critical temperatures, different geometry,
and different couplings to the central section. As a result, a single dip in
$dI/dV$ splits into two dips occurring at different bias voltages $V_{CL}$ and
$V_{CR}$, at which the left and the right leads switch to the normal state. In
this case Eq. (4) cannot be solved analytically. However, in the limit $b\gg
1$ one can rather accurately approximate the solution by the sum of two
independent dips,

$\displaystyle\frac{dI}{dV}$ $\displaystyle=$
$\displaystyle\frac{1}{R}-\frac{2\theta(V_{CL}-|V|)}{bR}\left(\frac{(b-1)|V|}{\sqrt{bV_{CL}^{2}-(b-1)V^{2}}}+1\right)$
(9)
$\displaystyle-\,\frac{2\theta(V_{CR}-|V|)}{bR}\left(\frac{(b-1)|V|}{\sqrt{bV_{CR}^{2}-(b-1)V^{2}}}+1\right).$

For a symmetric case $V_{CL}=V_{CR}$ and for $b\gg 1$ Eq. (9) coincides with
the exact Eq. (7) everywhere except a narrow interval close to the minimum of
the dip. Eq. (9) is valid for both positive and negative bias and it assumes
that the two dips are sufficiently close so that one can use the same
parameter $b$ (8) for both of them.

###  II.1 Suspended symmetric junction

In this subsection we consider a suspended junction in a symmetric
configuration, in which the only cooling mechanism is the removal of heat
through the superconducting leads, see Fig. 1b. This model can be easily
solved numerically, and it allows us to test the analytical approximation (7).
By symmetry, the power $IV/2$ is dissipated in each lead. We assume that the
right lead has the shape of a one-dimensional wire with the cross-sectional
area $S$, although the final expression for the cooling power of the junction
(13) given below, does not depend on geometry of the lead thanks to the
Wiedemann–Franz law relating heat and electric currents. Then, the temperature
profile in the right lead is determined by the heat diffusion equation

$S\kappa_{S}(T)\frac{dT}{dx}=-\frac{IV}{2},$ (10)

where the themal conductivity of the lead is [46],

$\kappa_{S}(T)=\frac{4k_{B}^{2}\sigma}{e^{2}}T\int^{\infty}_{\frac{\Delta(T)}{2k_{B}T}}\frac{x^{2}}{\cosh^{2}x}dx.$
(11)

Here we have introduced the conductivity of the superconducting material in
the normal state $\sigma$. We assume that the contact with the non-
superconducting central part is located at $x=0$ so that $T(0)=T_{0}$, and at
$x=L$ the superconducting wire contacts the bulk lead kept at the bath
temperature, $T(L)=T_{\rm bath}$. Integrating Eq. (10) from $0$ to $L$, we
obtain the equation for the contact temperature $T_{0}$,

$\displaystyle P_{S}(T_{0},T_{\rm bath})=I(V)V,$ (12)

where

$\displaystyle P_{S}=\frac{8k_{B}^{2}}{e^{2}R_{\rm leads}}\int_{T_{\rm
bath}}^{T_{0}}dT\,T\int_{\frac{\Delta(T)}{2k_{B}T}}^{\infty}dx\frac{x^{2}}{\cosh^{2}x}$
(13)

is the power carried away by quasiparticles through both superconducting
leads. In Eq. (12) we introduced the total resistance of the two leads in the
normal state $R_{\rm leads}=2L/\sigma S$. In this model, the heat conductance
at critical temperature is determined by Wiedemann–Franz law, $G_{C}^{\rm
th}=2\pi^{2}k_{B}^{2}T_{C}/3e^{2}R_{\rm leads}$, and the dimensionless
parameter (8) acquires the form

$\displaystyle b=\frac{28\zeta(3)}{3\alpha^{2}}\frac{R}{R_{\rm leads}}.$ (14)

Next, from Eq. (12) we find the value of the critical voltage $V_{C}$ at low
bath temperatures $T_{\rm bath}\ll T_{C}$,

$\displaystyle V_{C}\approx\sqrt{RP_{S}(T_{C},0)}=0.829\sqrt{\frac{R}{R_{\rm
leads}}}\frac{\Delta(0)}{e},$ (15)

where we have numerically evaluated the cooling power $P_{S}(T_{C},0)$
assuming the BCS dependence $\Delta(T)$. Assuming that the excess current is
fully suppressed for $V>V_{C}$, one finds the critical power as
$P_{C}=V_{C}^{2}/R$, which gives

$\displaystyle P_{C}=P_{S}(T_{C},0)=0.687\left({\Delta^{2}(0)}/{e^{2}R_{\rm
leads}}\right).$ (16)

Interestingly, Eqs. (15,16) agree with many experimetal observations reported
in the literature. Indeed, experiments in finite magnetic field [22] have
shown that the critical voltage is proportional to the superconducting gap,
$V_{C}\propto\Delta(0)$, and the experiments with tip-shaped superconducting
electrodes [24, 25] demonstrated that it is proportional to the square root of
the interfacial resistance $R_{IF}$ ($V_{C}\propto\sqrt{R_{IF}}$). Eq. (16)
predicts that the critical power does not depend on the junction resistance in
agreement with the experiment [23]. We also note that at sufficiently low
temperatures Eq. (15) may also be applicable to junctions lying on the
substrate if the overheated leads have much higher heat conductance than the
substrate.

In Fig. 2 we present numerically exact results obtained from Eq. (12). We have
chosen $\alpha=0.5$ and approximated the I-V curve as follows:
$I(V)=(1+\alpha/2)V/R$ for $eV<2\Delta$, and $I(V)=(V+\alpha\Delta/e)/R$ for
$eV>2\Delta$. We have ignored MAR features at sub-gap voltages because we
focus at voltages above the double gap, $eV>2\Delta$. Voltage dependence of
the contact temperature $T_{0}(V)$ and of the current $I(V)$ are shown in
Figs. 2(a) and 2(b) respectively. In Fig. 2(c) we compare the exact
differential conductance with the approximate expression (7) and observe
perfect agreement between them for all voltages $eV>2\Delta$.

Finally, we provide approximations for the quasiparticle cooling power (13) at
low and at high temperatures. Eq. (13) can be re-written in the form

$\displaystyle P_{S}=F(T_{0})-F(T_{\rm bath}),$ (17)

where at low temperatures $k_{B}T\ll\Delta(0)$ the function $F(T)$ is
exponentially small,

$\displaystyle F(T)=\frac{4\Delta(0)k_{B}T}{e^{2}R_{\rm
leads}}e^{-{\Delta(0)}/{k_{B}T}},$ (18)

while at temperatures slightly below $T_{C}$ and for $T>T_{C}$ it takes the
form

$\displaystyle F(T)=\frac{\pi^{2}k_{B}^{2}T^{2}}{3e^{2}R_{\rm leads}}$
$\displaystyle-\,\frac{2}{15}\left(\frac{8\pi^{2}}{7\zeta(3)}\right)^{3/2}\frac{k_{B}^{2}T_{C}^{2}}{e^{2}R_{\rm
leads}}{\rm Re}\,\left[\left(1-\frac{T}{T_{C}}\right)^{{5}/{2}}\right].$ (19)

The first term in the right hand side of Eq. (19) corresponds to the
Wiedemann–Franz heat current through the normal leads, and the second term
describes a small correction to it caused by superconductivity onset at
temperatures below $T_{C}$.

##  III Experiment

Our back-gated devices were fabricated using exfoliated graphene on a SiO2/Si
wafer. The substrate was highly boron-doped Si (p++/B, $\rho<10$ m$\Omega$cm)
which remains conducting down to mK temperatures. A layer of dry oxide with a
thickness of $250-270$ nm was grown thermally on the substrate at a
temperature of 1000∘C. Before depositing graphene, alignment markers were
patterned on the chips using optical lithography and evaporation of Ti and Au.

For suspended samples, we made graphene exfoliation on LOR resist [47]. After
making the electrical contacts using electron beam lithography, the area to be
suspended was irradiated using a high dose of electrons, developed in ethyl
lactate, and washed in hexane. The electrical leads were made of
aluminum/titanium sandwich structures. Titanium was evaporated at ultra high
vacuum conditions which resulted in high transparency contacts to graphene.
The presence of Ti lowered the superconducting transition temperature of the
leads: in our second sample structure with a trilayer sandwich structure
10/50/10 nm of Ti/Al/Ti, the critical temperature was suppressed to half of
the bulk value. The normal state resistance of the Al leads is $<1$
$\Omega/\mu$m, which means that the Joule heating by the measurement current
will mostly take place in the SGS junction. For details of the structure and
dimensions of the two samples see Table I.

Figure 3: Samples structure with the suspended (etched) part visible in the center between the vertical borders of the LOR resist seen as the darker background color. The Al/Ti measurement leads are seen as horizontal lighter grey stripes, the connection of which to the measurement system is displayed in the upper right corner. The green arrow depicts the electrical current flow in our experiments, while the red arrows illustrate the heat flows to the metallic leads as assumed in our analysis. | type | $L$ | $W$ | leads |  $d$ (nm) | $W_{\mathrm{lead}}$ | $N$ |  $T_{C}$ (K)
---|---|---|---|---|---|---|---|---
NS | ML | 0.30 | 6.0 | Ti/Al/Ti | 10/50/10 | 0.6 | 2 | 0.58
S | BL | 0.50 | 1.0 | Ti/Al | 10/70 | 1.2 | 4 | 0.77
Table 1: Parameters of the two studied samples: NS - nonsuspended, S -
suspended; the types ML and BL refer to monolayer and bilayer graphene,
respectively. $L$ and $W$ denote the sample length and width in micrometers,
respectively. Column $d$ indicates the thickness of the metal layers
evaporated for the contacts, while $W_{\mathrm{lead}}$ denotes the lead width
in $\mu$m and $N$ is the effective number of leads for hat transport. Critical
temperature $T_{C}$ represents an average value of the superconducting
transition temperatures observed on different leads.

The suspended sample (S in Table I) with four measurement leads is illustrated
in Fig. 3. The width of measurement leads is 1.2 $\mu$m, i.e. the contact
areas are more than two times larger than the size of the actual suspended
part of the graphene. The bright looking sections of the horizontal leads in
Fig. 3 are also suspended. Thus, the Joule heating generated in graphene has
to propagate along suspended metallic leads over a distance of at least two
microns before any cooling by the substrate through the LOR resist may take
place. The arrows in the figure illustrate the flow of electrical current as
well as the main directions of Joule heat escape from the sample. The critical
temperature was slightly different for the four leads and the obtained $T_{C}$
values varied between $T_{C}=760-780$ mK (expected BCS gap $\Delta=115\dots
118$ $\mu$eV).

Figure 4: (a) Measured $IV$ curve of the sample NS at gate voltage $V_{g}=30$
V (electron concentration $n\simeq 2\times 10^{12}$ cm-2) with the normal
state resistance $R=128$ $\Omega$ (indicated by the red line). Parameters of
the sample: excess current above the double gap voltage $I_{\rm
exc}^{(1)}=455$ nA, excess current at high bias $I_{\rm exc}^{(2)}=216$ nA,
$\Delta=87.5$ $\mu$eV, the parameter $\alpha=1.5$, critical current
$I_{C}=260$ nA, re-trapping current $I_{r}=175$ nA. (b) $dI/dV$ at several
values of $V_{g}$. The inset shows the critical power $P_{C}$ (red circles)
and the normal state resistance $R$ (black line) as functions of $V_{g}$.
Figure 5: Heating effects in a suspended bilayer SGS sample with the following
parameters: normal state resistance $R=2.2$ k$\Omega$, maximum excess current
$I_{\rm exc}=26.5$ nA, superconducting gap $\Delta\approx 90$ $\mu$eV,
$\alpha=0.646$. (a) Differential resistance as a function of power at
different bath temperatures (successive curves have been shifted vertically by
20 $\Omega$. For clarity, we define the power as $P=IV\times{\rm sign}V$, so
that it has different signs for positive and negative bias voltages. (b)
Critical powers $P_{C}$, i.e. positions of the resistance peaks in frame 5(a),
as a function of bath temperature. Continuous lines are fits using the heat
flow power law in Eq. (21) with the exponent $\gamma$ as a fitting parameter.
All these data were measured at $V_{g}=-20$ V which corresponds to a hole
carrier density of $n=2\times 10^{11}$ cm-2. Inset: Simulated results for the
superconducting-to-normal transition due to Joule heating, where the
resistance peaks arise from a reduction of excess current. Solid lines have
been smoothed by a moving average, whereas the dashed blue line displays the
unsmoothed result for $T_{\mathrm{bath}}=100$ mK.

Measurements were carried out using two dilution refrigerators: a BlueFors
prototype cryostat BF-SD125 and a commercial BlueFors BF-LD250. Most of the
measurements were performed on the BF-LD250 cryostat with the lowest
achievable base temperature of approx. 10 mK. Temperatures were measured with
a RuO2 resistor that had been calibrated against a Coulomb blockade
thermometer [48]. An electrically shielded room with filtered power lines was
employed to prevent electromagnetic interference and low-frequency electrical
noise from aﬀecting the weak signals before ampliﬁcation. To limit the amount
of thermal noise entering from room temperature to the sample, all dc lines
were equipped with three-stage $RC$ low-pass filters at the base temperature
of the cryostat, with a cut-off frequency of $f_{co}\simeq 1$ kHz.

The electrical characteristics were obtained using regular low-noise
measurement methods and lock-in techniques. We employed current bias through
$1-100$ M$\Omega$ resistor, selected on the basis of the required current
range (dependent on gate voltage $V_{g}$). The voltage across the sample was
measured using an LI-75 low-noise preamplifier powered from lead batteries.
Differential resistance at small audio frequency (around 70 Hz) was measured
simultaneously with the DC characteristics using a Stanford SR830 lock-in
amplifier.

In Fig. 4(a) we show the I-V curve of the non-suspended SGS junction taken at
low bath temperature $T_{\rm bath}=90$ mK. This junction has rather large
excess current exceeding the Josephson critical current. For this reason, its
suppression at high bias voltage is clearly visible. The experimental I-V
curve resembles the theoretical one plotted in Fig. 2(b), although in the
experiment only one of the leads completely switches to the normal state in
presented range of bias voltages. This junction exhibits large value of the
pre-factor in the excess current, $\alpha=1.5$, which points to its high
transparency. Assuming that the junction is long and diffusive, we can use Eq.
(3) and estimate the ratio $L_{\rm jct}/\xi\approx 1.65$, which gives the
coherence length $\xi\approx 180$ nm. In Fig. 4(b) we plot the differential
resistance of the same sample at several values of the gate voltage. It is
evident that both the critical voltage $V_{C}$ and the junction resistance $R$
(see the inset) significantly vary with $V_{g}$. However, as the inset shows,
the critical power $P_{C}$ only weakly depends on $V_{g}$ and varies from 0.77
nW to 0.96 nW for the presented five values of the gate voltage. This
obervation agrees with that of Ref. [23], and it is also consistent with Eqs.
(15,16). We believe that the remaining weak dependence of $P_{C}$ on $V_{g}$
results from the incomplete suppression of the $V_{g}$-dependent excess
current above $V_{C}$ in this sample.

Figure 6: (a) Excess current versus bias voltage in a suspended SGS junction
measured at the base temperature ($T=10$ mK). (b) Temperature of the leads in
the vicinity of the suspended SGS junction, $T_{0}(V)$. (c) Cooling power of
the leads versus the temperature at the graphene/Al contacts (given in frame
(b)) for the positive (red line) and for the negative (blue line) bias. Dashed
lines are the power law fits assuming $T_{\rm bath}\ll T_{C}$. (d)
Differential conductances of the suspended SGS junction at negative bias
measured at $T_{\rm bath}=10$mK (superconducting leads, black line) and at
$T_{\rm bath}=1$K (normal leads, blue line). Red dashed line shows the fit
with Eq. (9) and the arrow points to the MAR peak corresponding to the voltage
$\Delta/e$.

The superconducting-to-normal transition due to Joule heating in a suspended
bilayer SGS sample is depicted in Fig. 5a for a temperature range of $10\dots
800$ mK ($T_{C}\simeq 770$ mK). Each resistance peak is clearly split into two
parts, presumably because of junction asymmetry. There is also a small
asymmetry between positive and negative bias voltages, giving rise to a total
of four peaks.

Heat flow in mesoscopic devices can typically be characterized by a power law
of the form

$P=\Sigma(T_{L,R}^{\gamma}-T_{\rm bath}^{\gamma})$ (20)

where $\Sigma$ is a constant and $\gamma$ is a characteristic exponent which
depends on the dissipation mechanism. The transition to the normal state
occurs at the critical temperature $T_{L,R}=T_{C}$ and the corresponding
critical powers can be expressed as

$P_{C}=\frac{P_{C,0}}{T_{C}^{\gamma}}(T_{C}^{\gamma}-T_{\rm bath}^{\gamma}),$
(21)

where $P_{C,0}$ is the critical power at $T_{\rm bath}=0$. The observed peak
positions from Fig. 5(a) are plotted in Fig. 5(b), along with the fits using
Eq. (21) and the exponent $\gamma$ as a fit parameter. We obtain
$\gamma\approx 3.1$ for all four peaks. The four resistance traces in Fig.
5(b) converge to two diﬀerent critical temperatures, $\sim 760$ mK and $\sim
780$ mK. This supports the idea that the two leads have slightly different
critical temperatures.

Next, we analyze the dependence of the differential conductance of the
suspended SGS junction on bias voltage. First, we find bias dependent excess
current by subtracting the IV-curves of the junction measured at the base
temperature of $T=10$ mK, when the leads are superconducting, and at 1K, when
the leads become normal, i.e. we define $I_{\rm exc}(V)=I_{\rm 10mK}(V)-I_{\rm
1K}(V)$. This current is plotted in Fig. 6(a). Two red dashes lines show the
maximum and the minimum value of $I_{\rm exc}(V)$, they are equal to $\pm
26.5$ nA and should correspond to the non-suppressed theoretical excess
current $\pm I_{\rm exc}$ defined in Eq. (1). Analyzing multiple Andreev
reflection pattern in $dI/dV$ for this sample (see Fig. 6(d)), we estimate the
superconducting gaps $\Delta_{L}\approx\Delta_{R}=90$ $\mu$eV and find the
junction resistance $R=2195$ $\Omega$. With these values, we obtain the excess
current suppression parameter from Eq. (1) as $\alpha=0.646$. From Eq. (3) we
estimate the ratio $L_{\rm jct}/\xi\approx 3.8$ and $\xi\approx 130$ nm. The
latter value is similar to the one obtained for the non-suspended junction,
which supports the validity of the diffusive model for graphene.

In Fig. 6(b) we plot the effective temperature of the leads $T_{0}(V)$, which
is obtained from the bias dependent excess current plotted in Fig. 6(a) by
numerically solving the equation $|I_{\rm exc}(V)|/I_{\rm
exc}^{\max}=\Delta(T_{0})/\Delta_{0}$ and assuming BCS gap-temperature
dependence, $\Delta(T_{0})\approx\Delta_{0}\tanh(1.737\sqrt{T_{C}/T_{0}-1})$.
We used $T_{C}=770$ mK, which is the average value between the critical
temperatures of the two leads given above. The accuracy of this procedure is
limited by the noise in $I_{\rm exc}(V)$, which stems from the numerical
subtraction of the two experimental I-V curves. For this reason, the splitting
between the two critical voltages $V_{CL}$ and $V_{CR}$, which is clear in
Fig. 5(a), is not visible in Fig. 6(a). Therefore, here we ignore relatively
weak asymetry between the leads and put $T_{L}=T_{R}=T_{0}$. Furthermore, this
type of thermometry only makes sense in the limited range of bias voltages
where the excess current decreases with the bias, but the critical temperature
of the leads is not yet reached. Therefore, only these two voltage intervals
at positive and at negative bias are shown in Fig. 6(b). The dependence
$T_{0}(V)$, obtained in this way, resembles the model prediction shown in Fig.
2(a).

In Fig. 6(c) we plot the Joule heating power $P=IV$, which is equal to the
cooling power of the junction, versus the temperature of the leads $T_{0}(V)$
presented in Fig. 6(b). The red line corresponds to the positive bias and the
blue one – to the negative bias. The dashed lines present the predictions of
Eq. (20) at very low bath temperature, $T_{\rm bath}\to 0$, and with the same
values of the exponent $\gamma=3.1$ and $\Sigma=P_{C,0}/T_{C}^{\gamma}$ as in
the theory curves shown in Fig. 5(b). Namely, we have used $P_{C,0}=4.7$ nW
for positive and $P_{C,0}=4.4$ nW for negative bias, and these values match
those extracted from the positions of the peaks in $dV/dI$ with larger power,
see Fig. 5(b). The agreement between the dashed lines and the experimental
curves is good at relatively high temperatures close to $T_{C}$. Thus, we have
confirmed the validity of Eq. (20) in a broader range of temperatures by
varying the temperature of the leads $T_{0}(V)$ with the bias voltage and
keeping the bath temperature temperature low, $T_{\rm bath}\sim 10$ mK. At
temperatures $T_{0}(V)\lesssim 0.6$ K the power law curves deviate from the
experimetal data. However, there the accuracy of the numerical conversion of
the excess current to temperature becomes very poor and we cannot draw any
conclusion from this discrepancy.

The differential conductances of the suspended SGS junction in the
superconducting ($T_{\rm bath}=10$ mK, black line) and in the normal ($T_{\rm
bath}=1$ K, blue line) states are presented in Fig. 6(d). We show only
negative bias voltages, where the dips associated with the transition of the
leads to the normal state are more prononced. In the same figure we also show
the fit based on Eq. (9) with the red dashed line. The theoretical curve has
been plotted with the paramters earlier extracted from the fits shown in Figs.
5(b) and 6(a), namely, we have used $\alpha=0.646$, $R=2195$ $\Omega$ and
$P_{C,0}=4.4$ nW. With these parameters, we obtain the thermal conductance at
critical temperature $G_{C}^{\rm th}=\gamma P_{C,0}/T_{C}=17.7$ nW/K. To
obtain this value, we have taken the derivative of the cooling power (20) at
$T_{\rm bath}=T_{C}$ and, as before, have set $\gamma=3.1$ and $T_{C}=0.77$ K.
The parameter $b$, defined in Eq. (8), takes the value $b=2.776\times 10^{4}$.
The two critical voltages were adjustable parameters, they were chosen as
$V_{CL}=-2.82$ mV and $V_{CR}=-3.2$ mV. In addition, we also introduced the
non-linear correction to the differential conductance of the junction in the
normal state, which originated from the energy dependent density of states in
graphene and was needed to better reproduce the voltage dependence of the blue
line in Fig. 6(d). Accordingly, in Eq. (9) we have replaced $1/R\to(1+\beta
V)/R$ with $\beta=0.015$ mV-1. We find that Eq. (9), modified in this way,
rather well reproduces the shape of the experimental $dI/dV$-curve everywhere
except at low bias voltages, where the MAR peaks become visible. As expected,
the dips caused by transistion of the leads to the normal state are smeared in
the experimental curve. We also note that Fig. 6(d) is a different
representation of the inset in Fig. 5(b), and that it resembles the
theoretical model plot of Fig. 2(c).

We conclude this section by noticing that the model of Sec. II, which relates
the dips in the differential conductance of a highly transparent Josephson
junction to the suppression of the excess current by Joule heating, well
describes our experimental data.

##  IV Discussion and conclusions

In the previous section we have demonstrated how one can obtain the dependence
of the cooling power of the Josephson junction on the bath temperature and on
the junction temperature by fitting its differential conductance at high bias
voltages above $2\Delta/e$. This simple technique is applicable to any highly
transparent Josephson junction and it might be useful for high-flux
calibration of bolometers containing such junctions.

Our experimental results on graphene deserve further specific discussion. We
find that the critical Joule heating $P_{C}$ becomes independent of the bath
temperature at $T_{\rm bath}\ll T_{C}$, see Fig. 5(b). Such saturation of the
heat flow at low temperature is rather common in systems with weak coupling
between electrons and lattice phonons. We find similar behavior in the
theoretical electronic heat transport model of Sec. II.1. Indeed, for
$T_{0}=T_{C}$ and at $T_{\rm bath}\ll T_{C}$, where the energy gap of the
superconductor does not depend on temperature, the heat power (13) saturates
at the limiting value (16).

According to Figs. 4(a) and 5(b), the saturated low-temperature critical power
amounts to $P_{C}^{\mathrm{NS}}=1$ nW and $P_{C}^{\mathrm{S}}=4$ nW in our
non-suspended and suspended sample, respectively. If we scale these values by
the number of leads (see Table I) and the cross sectional area of the sandwich
conductors, $d\times W_{\mathrm{leads}}$, we obtain the values
$P_{C}^{\mathrm{NS}}/(d\times W_{\mathrm{leads}})=8.67$ kW/m2 and
$P_{C}^{\mathrm{S}}/(d\times W_{\mathrm{leads}})=13.7$ kW/m2. The similarity
of these two critical power densities supports the conclusion that thermal
transport along the leads is the main mechanism of heat relaxation in our
devices.

The critical voltages amount to $V_{C}^{\mathrm{NS}}=0.34$ mV and
$V_{C}^{\mathrm{S}}=3.0$ mV for samples NS and S, respectively. We have
applied Eq. (15) at low $T$ to investigate whether these values scale
according to the model. Since we do not know exactly the lead resistance, we
compare the measured ratio of critical voltages to the theory by assuming that
the lead resistance scales inversely proportional to the cross sectional area
of the leads [49]. With the parameters from Table I, Eq. (15) yields for the
critical voltage ratio $V_{C}^{\mathrm{S}}\large{/}V_{C}^{\mathrm{NS}}=8.5$,
while from the experimental data in Figs. 4 and 5 we obtain
$V_{C}^{\mathrm{S}}\large{/}V_{C}^{\mathrm{NS}}=8.8$. Thus, our theoretical
model accounts for the variation of the critical voltage within a range
covering one order of magnitude.

The observed exponent in the temperature dependence of the power for the
suspended junction, $\gamma=3.1$, is close to the one found by Borzenets et al
for SGS junctions with Pb leads and graphene deposited on Si/SiO2 substrate
[50]. Lead has a large superconducting gap, which should effectively cut off
electronic thermal conductance below 1 K and leave only electron-phonon
coupling for energy relaxation [52]. For SGS junctions with Al/Ti leads on the
Si/SiO2 substrate the exponents $\gamma\approx 4$ [23] and $\gamma=3.1$ [51],
which agrees with our findings, have been previously reported. Power law
thermal transport with the exponent $\gamma=4.5$ has also been observed in
superconducting Nb at low temperatures $k_{B}T\ll\Delta$ [53].

The origin of the exponent $\gamma\approx 3$ observed in our suspended device
is not fully clear. Let us briefly discuss several possible explanations for
this observation. In our experiment the ratio $k_{B}T/\Delta$ extends to large
values, and, therefore, energy relaxation should be dominated by
quasiparticles carrying heat through the leads. For temperatures close to
$T_{C}$ the quasiparticles behave almost as normal state electrons obeying
Wiedemann-Franz law, which leads to $\gamma=2$, see Eq. (19). The increase of
the exponent $\gamma$ from 2 to 3 may be qualitatively explained by several
different effects or by their combination. First, the quasiparticle heat
current (13) is lower and more sensitive to temperature than the normal state
one in the whole interval $T<T_{C}$, including the temperatures close to
$T_{C}$ as Eq. (19) shows. Second, the parameter $R_{\rm leads}$, appearing in
the cooling power (13), is the resistance of the part of an Al lead in which
thermalization between the quasiparticles and the substrate phonons occurs.
Since the corresponding thermalization length $L_{T}$ is reduced with
increasing temperature, the resistance $R_{\rm leads}\propto L_{T}$ is also
reduced. Hence, the dependence of the power (13) on temperature becomes
stronger than it were in the model with fixed $R_{\rm leads}$. Yet another
possibility for the increase of the exponent is the injection of strongly non-
equilibrium quasiparticles with non-thermal distribution in the Al leads at
high bias voltages. The propagation of such quasiparticles should be described
by the full set of kinetic equations with the inclusion of charge imbalance
instead of just one Eq. (10). Quantitative analysis of the effects mentioned
above goes beyond the scope of this paper. At the moment, the similarity of
the exponents observed in our experiment and in Refs. [50, 51] remains
puzzling.

In conclusion, we have analyzed, both theoretically and experimentally, the
effect of Joule heating on the IV characteristics of highly transparent
Josephson junctions. We have shown that the Joule heating supresses the excess
current and induces dips in the differential conductance at bias voltages
where the two leads switch to the normal state. We have derived simple
analytical expressions for the shape of these dips and for their positions in
a specific case of a suspended junction, which is cooled only through the
leads. We have experimentally tested preditions of the theory by studing SGS
junctions in both suspended and non-suspended configurations. Good agreement
between the experiment and the theory is found. In the suspended junction we
observe power law temperature dependence of the heat current with the exponent
$\gamma=3.1$ instead of theoretically expected exponential dependence (18) at
low temperatures and power law scaling with $\gamma=2$ at higher temperatures
(19). Our work also indicates that excess current can be employed for
thermometry in special cases.

## Acknowledgements

Discussions with Manohar Kumar and Alexander Savin are gratefully
acknowledged. This work was supported by the Academy of Finland projects
314448 (BOLOSE), 310086 (LTnoise) and 312295 (CoE, Quantum Technology Finland)
as well as by ERC (grant no. 670743). This research project utilized the Aalto
University OtaNano/LTL infrastructure which is part of European Microkelvin
Platform (funded by European Union’s Horizon 2020 Research and Innovation
Programme Grant No. 824109). A.L. is grateful to Väisälä foundation of the
Finnish Academy of Science and Letters for scholarship. The work at HSE
University (M.R.S. and A.S.V.) was supported by the framework of the Academic
Fund Program at HSE University in 2021 (grant No 21-04-041).

## References

  * [1] W. Haberkorn, H. Knauer, and J. Richter, Phys. Status Solidi A47, K161 (1978).
  * [2] A. A. Golubov, M. Yu. Kupriyanov, and E. Il’ichev, Rev. Mod. Phys. 76, 411 (2004).
  * [3] M. Hays, G. de Lange, K. Serniak, D.J. van Woerkom, D. Bouman, P. Krogstrup, J. Nygard, A. Geresdi, and M.H. Devoret, Phys. Rev. Lett. 121, 047001 (2018).
  * [4] A. Zazunov, V. S. Shumeiko, E. N. Bratus, J. Lantz, and G.Wendin, Phys. Rev. Lett. 90, 087003 (2003).
  * [5] M. Octavio, M. Tinkham, G.E. Blonder, and T.M.Klapwijk, Phys. Rev. B 27, 6739 (1983).
  * [6] J. C. Cuevas, A. Martin-Rodero, and A. Levy Yeyati, Phys. Rev. B 54, 7366 (1996).
  * [7] A. Bardas and D.V. Averin, Phys. Rev. B 56, R8518 (1997).
  * [8] D.V. Averin, G. Wang, and A.S. Vasenko, Phys. Rev. B 102, 144516 (2020).
  * [9] X. Du, I. Skachko, and E. Y. Andrei, Phys. Rev. B 77, 184507 (2008).
  * [10] W. Liang, M. Bockrath, D. Bozovic, J. H. Hafner, M. Tinkham, and H. Park, Nature 411, 665 (2001).
  * [11] P. Jarillo-Herrero, J.A. van Dam, and L.P. Kouwenhoven, Nature 439, 953 (2006).
  * [12] J-D. Pillet, C. H. L. Quay, P. Morfin, C. Bena, A. Levy Yeyati and P. Joyez, Nat. Phys. 6, 965 (2010).
  * [13] V.E. Calado, S. Goswami, G. Nanda, M. Diez, A.R. Akhmerov, K. Watanabe, T. Taniguchi, T.M. Klapwijk, and L.M.K. Vandersypen, Nat. Nanotechnology 10, 761 (2015).
  * [14] G.-H. Lee, S. Kim, S.-H. Jhi, H.-J. Lee, Nature Comm. 6, 6181 (2015).
  * [15] S. Abay, D. Persson, H. Nilsson, F. Wu, H.Q. Xu, M. Fogelström, V. Shumeiko, and P. Delsing, Phys. Rev. B 89, 214508 (2014).
  * [16] E.M. Spanton, M. Deng, S. Vaitiekenas, P. Krogstrup, J. Nygard, C.M. Marcus, and K.A. Moler, Nat. Phys. 13, 1177 (2017).
  * [17] I. Sochnikov, L. Maier, C.A. Watson, J.R. Kirtley, C. Gould, G. Tkachov, E. M. Hankiewicz, C. Brüne, H. Buhmann, L.W. Molenkamp, and K.A. Moler, Phys. Rev. Lett. 114, 066801 (2015).
  * [18] M. Veldhorst, M. Snelder, M. Hoek, T. Gang, V.K. Guduru, X.L. Wang, U. Zeitler, W.G. van der Wiel, A.A. Golubov, H. Hilgenkamp and A. Brinkman, Nat. Mat. 11, 417 (2012).
  * [19] S. Charpentier, L. Galletti, G. Kunakova, R. Arpaia, Y. Song, R. Baghdadi, S.M. Wang, A. Kalaboukhov, E. Olsson, F. Tafuri, D. Golubev, J. Linder, T. Bauch, and F. Lombardi, Nat. Comm. 8, 2019 (2017).
  * [20] G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Phys. Rev. B 25, 4515 (1982).
  * [21] Chanh Nguyen, Herbert Kroemer, and Evelyn L. Hu, Phys. Rev. Lett. 69, 2847 (1992).
  * [22] Peng Xiong, Gang Xiao and R.B. Lebowitz, Phys. Rev. Lett. 71, 1907 (1993).
  * [23] Jae-Hyun Choi, Hu-Jong Lee and Yong-Joo Doh, J. Korean Phys. Soc. 57, 149 (2010).
  * [24] P. S. Westbrook and A. Javan, Phys. Rev. B 59, 14606 (1999).
  * [25] J. A. Gifford, G. J. Zhao, B. C. Li, J. Zhang, D. R. Kim, and T. Y. Chen, J. Appl. Phys. 120, 163901 (2016).
  * [26] K. Flensberg and J. Bindslev Hansen, Phys. Rev. B 40, 8693 (1989).
  * [27] H. Courtois, M. Meschke, J. T. Peltonen, and J. P. Pekola, Phys. Rev. Lett. 101, 067002 (2008).
  * [28] M. Tinkham, J. U. Free, C. N. Lau, and N. Markovic, Phys. Rev. B 68, 134515 (2003).
  * [29] H. Vora, P. Kumaravadivel, B. Nielsen and X. Du, Appl. Phys. Lett. 100, 153507 (2012).
  * [30] K. C. Fong, and K. Schwab, Phys. Rev. X 2, 031006 (2012).
  * [31] J. Yan, et al., Nat. Nanotechnol. 7, 472–478 (2012).
  * [32] C. McKitterick, D. Prober and B. Karasik, J. Appl. Phys. 113, 044512 (2013).
  * [33] D. K. Efetov, et al., Nat. Nanotechnol. 13, 797–801 (2018).
  * [34] Q. Han, et al., Sci. Rep. 3, 3533 (2013).
  * [35] X. Cai, et al., Nat. Nanotechnol. 9, 814–819 (2014).
  * [36] A. E. El Fatimy, et al., Nat. Nanotechnol. 11, 335–338 (2016).
  * [37] G.-H. Lee, D.K. Efetov, W. Jung, L. Ranzani, E.D. Walsh, T.A. Ohki, T. Taniguchi, K. Watanabe, P. Kim, D. Englund and K.C. Fong, Nature 586, 42 (2020).
  * [38] R. Kokkoniemi, J.-P. Girard, D. Hazra, A. Laitinen, J. Govenius, R. E. Lake, I. Sallinen, V. Vesterinen, M. Partanen, J. Y. Tan, K. W. Chan, K. Y. Tan, P. Hakonen and M. Möttönen, Nature 586, 47 (2020).
  * [39] A.V. Zaitsev, Zh. Eksp. Teor. Fiz. 86, 1742 (1984) [Sov. Phys. JETP 59, 1015 (1984)].
  * [40] M. Titov and C. W. J. Beenakker, Phys. Rev. B 74, 041401(R) (2006).
  * [41] P. Joyez, D. Vion, M. Götz, M. H. Devoret, and D. Esteve, J. Supercond. 12, 757 (1999).
  * [42] O. N. Dorokhov, Solid State Commun. 51, 381 (1984).
  * [43] J. Tworzydlo, B. Trauzettel, M. Titov, A. Rycerz, and C.W.J. Beenakker, Phys. Rev. Lett. 96, 246802 (2006).
  * [44] J. C. Cuevas, J. Hammer, J. Kopu, J. K. Viljas, and M. Eschrig, Phys. Rev. B 73, 184505 (2006).
  * [45] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phy. Rev. 108, 1175 (1957).
  * [46] J. Bardeen, G. Rickayzen, and L. Tewordt, Phys. Rev. 113, 982 (1959).
  * [47] N. Tombros, A. Veligura, J. Junesch, J. J. van den Berg, P. J. Zomer, M. Wojtaszek, I. J. V. Marun, H. T. Jonkman, and B. J. van Wees, J. Appl. Phys. 109, 93702 (2011).
  * [48] J. P. Pekola, K. P. Hirvi, J. P. Kauppinen, and M. A. Paalanen, Phys. Rev. Lett. 73, 2903 (1994).
  * [49] The normal state resistance of the Al layer is much smaller than that of the Ti layer. Hence, we also checked the scaling using the Al layer thickness, but the difference is small, and the change would not modify our conclusions. We use the total thickness because the mechanism of the electronic heat transport in the superconducting state is not fully clear.
  * [50] I. V. Borzenets, U. C. Coskun, H. T. Mebrahtu, Y. V. Bomze, A. I. Smirnov, and G. Finkelstein, Phys. Rev. Lett., 111, 027001 (2013).
  * [51] K. Tsumura, N. Furukawa, H. Ito, E. Watanabe, D. Tsuya, and H. Takayanagi, Appl. Phys. Lett. 108, 033109 (2016).
  * [52] J. Voutilainen, A. Fay, P. Häkkinen, J. K. Viljas, T. T. Heikkilä, and P. J. Hakonen, Phys. Rev. B, 84, 045419 (2011).
  * [53] A. Feshchenko, O-P. Saira, J. Peltonen, J. P. Pekola, Sci Rep 7, 41728 (2017).


# Covariant formulation of gluon self-energy in presence of ellipsoidal
anisotropy

Ritesh Ghosha,b<EMAIL_ADDRESS>Bithika Karmakara,b
<EMAIL_ADDRESS>Arghya Mukherjeea,c<EMAIL_ADDRESS>aSaha
Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata - 700064, India  bHomi
Bhabha National Institute, 2nd floor, BARC Training School Complex,
Anushaktinagar, Mumbai, Maharashtra 400094, India  cSchool of Physical
Sciences, National Institute of Science Education and Research, HBNI, Jatni
752050, India

###### Abstract

In this work, a covariant formulation of the gluon self-energy in presence of
ellipsoidal anisotropy is considered. It is shown that the general structure
of the gluon self-energy can be written in terms of six linearly independent
projection tensors. Similar to the spheroidal anisotropy, mass scales can be
introduced for each of the collective modes considering the static limits.
With a simplified ellipsoidal generalization of the Romatschke-Strickland
form, the angular dependencies of the mass scales are studied. It is observed
that, compared to the spheroidal case, additional unstable mode may appear in
presence of ellipsoidal anisotropy depending upon the choice of the
parameters.

##  I Introduction

Over the last few decades, the plethora of knowledge gained through the
nonperturbative lattice QCD simulations Ding:2015ona ; Ratti:2018ksb ;
Bazavov:2019lgz , perturbative hard thermal loop calculations Braaten:1989mz ;
Andersen:2004fp ; Haque:2014rua ; Su:2015esa , effective hydrodynamical
modeling Gale:2013da ; Jeon:2015dfa ; Jaiswal:2016hex  as well as AdS/CFT
inspired studies Son:2007vk ; Nastase:2007kj ; Peschanski:2008xn ;
Sadeghi:2013zma  have paved the way of achieving a remarkable advancement in
our understanding of the experimental data available from the Relativistic
Heavy Ion Collider (RHIC) facility at Brookhaven National Lab and the Large
Hadron Collider (LHC) facility at European Center for Nuclear Research
(CERN)Stock:2010hoa ; Busza:2018rrf . The concerted efforts from different
heavy ion research communities suggest that the deconfined quark-gluon
plasma(QGP) matter produced in the ultrarelativistic heavy ion collision
experiments is most likely to possess substantial deviation from perfect local
isotropic equilibrium Strickland:2013uga . In fact, soon after the initial
nuclear impact, very large pressure anisotropy is expected in the center of
the fireball with even larger anisotropy prevailing in the cooler regions of
the plasma Strickland:2014pga . A promising way for hydrodynamical modeling of
such highly momentum space anisotropic system is to consider the framework of
anisotropic hydrodynamics Alqahtani:2017jwl ; Alqahtani:2017mhy . On the other
hand, incorporating certain classes of anisotropic momentum distributions in
the framework of hard thermal loop perturbation theory, the nonequilibrium
plasma properties can be extracted by studying the collective modes of the
quasipartons Mrowczynski:2004kv ; Mrowczynski:2000ed . A suitable
parametrization for such distributions can be achieved following Refs.
Romatschke:2003ms ; Romatschke:2004jh ; Schenke:2006fz  where the one particle
isotropic momentum space distribution function is deformed by introducing a
directional dependency. This particular Romatschke-Strickland (RS) form has
been widely used in different phenomenological applications such as photon and
dilepton production from anisotropic QGP Bhattacharya:2008up ;
Bhattacharya:2008mv ; Bhattacharya:2015ada , anisotropic heavy quark potential
Nopoush:2017zbu , bottomonia suppression Krouppa:2016jcl , plasma wakes
Mandal:2012wi ; Mandal:2013jla , nuclear modification factor Mandal:2011xn ,
quasiparticle descriptions of particle production Chandra:2016dwy  and so on.

An important aspect of considering the nonequilibrium momentum distributions
in QGP medium is the occurrence of kinetic instabilities. These are, in simple
terms, the collective modes that possess a positive imaginary part in their
mode frequencies resulting an exponential growth in the chromomagneic and
chromoelectric fields. The existence of such Chromo-Weibel instabilities
Strickland:2006cv  can influence the thermalization and isotropization of the
medium Rebhan:2004ur . A recent review with a pedagogical introduction to the
required theoretical tools for such studies can be found in Ref.
Mrowczynski:2016etf . In general, how many unstable modes are possible depends
on the choice of the parameters. For example, in case of spheroidal anisotropy
with Romatschke-Strickland parametrized form, the number of the possible
unstable modes differs for the prolate and the oblate case which are obtained
considering the negative and positive values for the anisotropy parameter
respectively Romatschke:2003ms . Moreover, there will be a directional
dependence as well. In other words, depending on the angle of propagation with
respect to the anisotropy direction, a stable collective mode may become
unstable and can give rise to additional instabilities.

As already mentioned, the large momentum space anisotropy in early stages can
be efficiently incorporated in the aHydro framework. This is because, unlike
usual viscous hydrodynamic set up, here, the dominant anisotropic
contributions to the distribution functions are captured in the leading order
of hydrodynamic expansion. The second order anisotropic hydrodynamics as
developed in Bazow:2013ifa  can take into account arbitrary transverse
expansion and is consistent with the traditional dissipative hydrodynamics
approach in the limiting case of small anisotropy. However, as argued in Ref.
Tinti:2013vba , the azimuthally symmetric ansatz in such approaches involving
a single anisotropy parameter can further be generalized in a systematic way.
More specifically, in the original anisotropic hydrodynamic setup, the leading
order local rest frame distribution is taken to be of RS form and thus the two
components of the pressure in the transverse plane can be different only after
including the second order corrections. A generalization to include three
different pressure components in the leading order of hydrodynamic expansion
has been developed in Refs. Tinti:2013vba ; Nopoush:2014pfa ;
Alqahtani:2015qja  and recently has further been generalized in Ref.
Nopoush:2019vqc . In the hard loop approach too, the generalized ellipsoidal
distribution has been considered in the parton self-energy studies
Kasmaei:2016apv  and implemented in phenomenological applications of photon
and dilepton production Kasmaei:2019ofu ; Kasmaei:2018oag . Obviously, the
consideration of an ellipsoidal anisotropy increases the number of parameters
compared to the spheroidal case, as, apart from one additional anisotropy
parameter dependence, the collective properties of the partons in this case
also possess azimuthal angular dependency. Recently, the gluonic unstable mode
with an ellipsoidal anisotropic momentum distribution has been reported for
the first time in Ref. Kasmaei:2018yrr . It is observed that the growth rate
of the unstable modes can become several times larger than the spheroidal case
for certain propagation directions. Also the number of possible unstable modes
are direction dependent. A convenient method to describe such directional
dependency is to consider the static limits of the collective modes. In that
case, one introduces mass scales corresponding to each of the collective modes
and studies their angular dependency. The occurrence of a negative value in
the mass scale in fact signifies the presence of an instability in the
corresponding mode which are well studied in the spheroidal case. However, for
the ellipsoidal case, such mass scales could not be defined in the formalism
adopted in Kasmaei:2018yrr  where the general structure of the gluon
polarization has not been considered. An important application of gluon self
energy lies in the determination of the perturbative part of the heavy quark
potential Dumitru:2007hy ; Dumitru:2009fy ; Nopoush:2017zbu  which requires a
covariant formulation of the general structure.

The primary objective of this work is to formulate a general structure for the
gluon polarization in presence of ellipsoidal momentum space anisotropy in a
covariant way and study the occurrence of the unstable modes in a similar
fashion as done in case of spheroidal momentum space anisotropy. A suitable
approach for constructing the tensor basis for gluon self energy is to choose
the maximum possible mutually orthogonal set. The choice becomes useful for
the derivation of the effective propagator where contractions among the basis
tensors are involved. Once the effective gluon propagator is derived, one can
obtain the collective modes from its pole. It should be mentioned here that
the gluon collective modes can also be obtained by solving the characteristic
equation without requiring any consideration of the general structure. In
fact, this is the procedure adopted in Ref.Kasmaei:2018yrr . However, the
advantage of considering the general structure is that, the collective modes
are expressed in terms of coordinate independent form factors which is
essential for introducing the mass scales. In fact, the nontrivial angular
dependence in the hard-thermal loop resummed potential enters through these
mass scales Nopoush:2017zbu .

The article is organized as follows. We first describe the formalism for
constructing the general structure for the gluon polarization function in
presence of the ellipsoidal anisotropy. The derived structure is subsequently
used to obtain the effective propagator and the collective modes of the gluon.
Considering to a simple ellipsoidal generalization of the RS form, the
numerical results for the form factors and the mass scales are presented in
the following section. Finally we conclude with a brief summary of the
presented work. The following convention is used throughout the article:
metric $\eta^{\mu\nu}=\mbox{diag}(1,-1,-1,-1)$ and a general momentum four
vector $Q^{\mu}=(q^{0},\bm{q})$ with $Q^{2}$ representing
$\eta_{\mu\nu}Q^{\mu}Q^{\nu}=Q\cdot Q$ and $q=|{\bm{q}}|$.

##  II Formalism

We consider a thermal QCD medium with a general anisotropic momentum
distribution characterized by two independent four vectors $a^{\mu}$ and
$b^{\mu}$ having unit norms. These two, together with the normalized heat bath
velocity ($u^{\mu}$) and the gluon momentum ($P^{\mu}$), can be used to form a
set of ten independent symmetric tensors so that the symmetric gluon
polarization tensor can be expressed as a linear combination of them. A simple
choice for this purpose may be the set of tensors $P^{\mu}P^{\nu}$,
$u^{\mu}u^{\nu}$, $b^{\mu}b^{\nu}$, $a^{\mu}a^{\nu}$,
$P^{\mu}u^{\nu}+P^{\nu}u^{\mu}$, $P^{\mu}b^{\nu}+P^{\nu}b^{\mu}$,
$P^{\mu}a^{\nu}+P^{\nu}a^{\mu}$, $u^{\mu}b^{\nu}+b^{\mu}u^{\nu}$,
$u^{\mu}a^{\nu}+u^{\nu}a^{\mu}$ and $b^{\mu}a^{\nu}+a^{\mu}b^{\nu}$. Notice
that, along with the four vectors, we have not considered the metric
$\eta^{\mu\nu}$ as usually done in isotropic or even in anisotropic case with
one anisotropy direction. This is because, in case of ellipsoidal anisotropy,
the metric itself no longer remains an independent tensor and can be expressed
as a linear combination of the chosen set. Now, the constraints from the
transversality condition $P^{\mu}\Pi_{\mu\nu}=0$ further reduce the number of
independent basis tensors to six. In the rest frame of the heat bath with
$u^{\mu}=(1,0,0,0)$, one of the anisotropy directions can be taken along z,
say $b^{\mu}=(0,0,0,1)$, whereas the other anisotropy direction can be assumed
to lie in the $xz$ plane without any loss of generality. In the following we
discuss a convenient method to obtain the basis tensors in a systematic way.

Let us first consider the general structure of the gauge boson self-energy in
vacuum that can be written as

$\displaystyle\Pi^{\mu\nu}=\left(\eta^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}\right)\Pi(P^{2})=V^{\mu\nu}\Pi(P^{2}).$
(1)

Using the tensor $V^{\mu\nu}$, we obtain $\tilde{u}^{\mu}=V^{\mu\nu}u_{\nu}$
which can be used to construct the first basis tensor given by

$\displaystyle
A^{\mu\nu}=\frac{\tilde{u}^{\mu}\tilde{u}^{\nu}}{\tilde{u}^{2}}.$ (2)

Now, it is useful to define the $A^{\mu\nu}$ subtracted part of $V^{\mu\nu}$
as $U^{\mu\nu}=V^{\mu\nu}-A^{\mu\nu}$ which can be used to obtain
$\tilde{b}^{\mu}$ defined as $\tilde{b}^{\mu}=U^{\mu\nu}b_{\nu}$ such that it
becomes orthogonal to $\tilde{u}^{\mu}$ by construction. Similar to the
earlier case, we obtain our second basis tensor using $\tilde{b}^{\mu}$ as

$\displaystyle
B^{\mu\nu}=\frac{\tilde{b}^{\mu}\tilde{b}^{\nu}}{\tilde{b}^{2}}.$ (3)

Another symmetric tensor that can be constructed intuitively using
$\tilde{b}^{\mu}$ and $\tilde{u}^{\mu}$ together is given by

$\displaystyle
C^{\mu\nu}=\frac{\tilde{u}^{\mu}\tilde{b}^{\nu}+\tilde{b}^{\mu}\tilde{u}^{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{b}^{2}}}.$
(4)

To obtain the rest of the tensors, once again we go through the similar steps:
at first we define $R^{\mu\nu}=U^{\mu\nu}-B^{\mu\nu}$, which, in this case,
can be viewed as the $A^{\mu\nu}$ and the $B^{\mu\nu}$ subtracted part of our
vacuum basis tensor $V^{\mu\nu}$. Then we obtain the $\tilde{a}^{\mu}$ from
$a^{\mu}$ as $\tilde{a}^{\mu}=R^{\mu\nu}a_{\nu}$. This time, the newly
constructed $\tilde{a}^{\mu}$ becomes orthogonal to $\tilde{b}^{\mu}$ as well
as $\tilde{u}^{\mu}$. Note that, all the four vectors of the set
$\tilde{u}^{\mu}$, $\tilde{b}^{\mu}$ and $\tilde{a}^{\mu}$ are orthogonal to
the gluon four momentum $P^{\mu}$ and hence the basis tensors constructed
using them trivially satisfy the transversality condition. Moreover, because
of the orthogonality among the constructed four vectors, the extraction of the
form factors gets simplified considerably. Now, with the constructed set, the
rest of the independent symmetric tensors can also be obtained intuitively as

$\displaystyle D^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{\tilde{a}^{\mu}\tilde{a}^{\nu}}{\tilde{a}^{2}},$ (5)
$\displaystyle E^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{\tilde{u}^{\mu}\tilde{a}^{\nu}+\tilde{a}^{\mu}\tilde{u}^{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{a}^{2}}},$
(6) $\displaystyle F^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{\tilde{a}^{\mu}\tilde{b}^{\nu}+\tilde{b}^{\mu}\tilde{a}^{\nu}}{\sqrt{\tilde{a}^{2}}\sqrt{\tilde{b}^{2}}}.$
(7)

The general structure of the gauge boson self-energy in presence of an
ellipsoidal anisotropic medium can be expressed as a linear combination of the
six basis tensors as

$\displaystyle\Pi^{\mu\nu}$ $\displaystyle=$ $\displaystyle\alpha
A^{\mu\nu}+\beta B^{\mu\nu}+\gamma C^{\mu\nu}+\delta D^{\mu\nu}+\sigma
E^{\mu\nu}+\lambda F^{\mu\nu}~{}.$ (8)

It is worth mentioning here that though tensors like $V^{\mu\nu}$ are used to
obtain the set, once we declare our constructed set of tensors as independent,
all the other symmetric and transverse tensors not belonging to the set become
expressible as their linear combination. For example, it can be shown that

$\displaystyle V^{\mu\nu}$ $\displaystyle=$ $\displaystyle
A^{\mu\nu}+B^{\mu\nu}+D^{\mu\nu}.$ (9)

Now, we can obtain the effective propagator from the Dyson-Schwinger equation

$\displaystyle i\mathcal{D}^{\mu\nu}$ $\displaystyle=$ $\displaystyle
i\mathcal{D}_{0}^{\mu\nu}+i\mathcal{D}_{0}^{\mu\rho}(i\Pi_{\rho\rho^{\prime}})i\mathcal{D}^{\rho^{\prime}\nu},$
(10)

where the inverse bare propagator $\mathcal{D}_{0}^{\mu\nu}$ without the
explicit color indices is given by

$\displaystyle(\mathcal{D}^{-1}_{0})^{\mu\nu}$ $\displaystyle=$
$\displaystyle-P^{2}\eta^{\mu\nu}-\frac{1-\zeta}{\zeta}P^{\mu}P^{\nu},$ (11)

with $\zeta$ representing the gauge fixing parameter. To obtain the effective
propagator, let us first consider the inner product identities among the basis
tensors. To express the identities in a compact form, here we suppress the
Lorentz indices of the basis tensors. Also the tensors $C$, $E$, and $F$ are
considered as a sum of two parts, for example $C=\overline{C}+\underline{C}$
and in similar fashion for the other two where

$\displaystyle\overline{C}_{\nu}^{\mu}=\frac{\tilde{u}^{\mu}\tilde{b}_{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{b}^{2}}},\hskip
14.22636pt\underline{C}_{\nu}^{\mu}=\frac{\tilde{b}^{\mu}\tilde{u}_{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{b}^{2}}},$
$\displaystyle\overline{E}_{\nu}^{\mu}=\frac{\tilde{u}^{\mu}\tilde{a}_{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{a}^{2}}},\hskip
14.22636pt\underline{E}_{\nu}^{\mu}=\frac{\tilde{a}^{\mu}\tilde{u}_{\nu}}{\sqrt{\tilde{u}^{2}}\sqrt{\tilde{a}^{2}}},$
$\displaystyle\overline{F}_{\nu}^{\mu}=\frac{\tilde{a}^{\mu}\tilde{b}_{\nu}}{\sqrt{\tilde{a}^{2}}\sqrt{\tilde{b}^{2}}},\hskip
14.22636pt\underline{F}_{\nu}^{\mu}=\frac{\tilde{b}^{\mu}\tilde{a}_{\nu}}{\sqrt{\tilde{a}^{2}}\sqrt{\tilde{b}^{2}}}.$

With this notation, the inner product between the basis tensors can be written
in a compact form as

(19)

where the composition rule for each entry of the multiplication table is
defined as $r^{\mu}_{\rho}c^{\rho}_{\nu}$ with $r$ and $c$ being members
representing the row and the column respectively. Moreover, further
contraction of the free indices as $r^{\mu}_{\rho}c^{\rho}_{\mu}$ simplifies
to

(27)

which will be useful to extract the coefficients of the basis tensors, i.e.,
the form factors from the polarization function. Using the properties of the
basis tensors from Eq. (19) one can obtain the effective gluon propagator
given by

$\displaystyle\mathcal{D}^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{\beta\delta-\left(\beta+\delta\right)P^{2}-\lambda^{2}+P^{4}}{\Delta}A^{\mu\nu}-\frac{\delta\alpha-\left(\delta+\alpha\right)P^{2}-\sigma^{2}+P^{4}}{\Delta}B^{\mu\nu}-\frac{\gamma\left(P^{2}-\delta\right)+\sigma\lambda}{\Delta}C^{\mu\nu}$
(28) $\displaystyle-$
$\displaystyle\frac{\alpha\beta-\left(\alpha+\beta\right)P^{2}-\gamma^{2}+P^{4}}{\Delta}D^{\mu\nu}-\frac{\sigma\left(P^{2}-\beta\right)+\lambda\gamma}{\Delta}E^{\mu\nu}-\frac{\lambda\left(P^{2}-\alpha\right)+\gamma\sigma}{\Delta}F^{\mu\nu}-\zeta\frac{P^{\mu}P^{\nu}}{P^{4}}.$

where the common denominator of the basis tensors, $\Delta$ is given by

$\displaystyle\Delta$ $\displaystyle=$ $\displaystyle
P^{6}-(\alpha+\beta+\delta)P^{4}-(\gamma^{2}+\sigma^{2}+\lambda^{2}-\alpha\beta-\beta\delta-\delta\alpha)P^{2}+\alpha\lambda^{2}+\beta\sigma^{2}+\delta\gamma^{2}-\alpha\beta\delta-2\gamma\sigma\lambda.$
(29)

Notice that, in presence of ellipsoidal anisotropy, the denominator of the
effective propagator becomes a cubic equation of $P^{2}$ and can be written as
a product of three factors as

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle(P^{2}-\Omega_{0})(P^{2}-\Omega_{+})(P^{2}-\Omega_{-}),$ (30)

where $\Omega_{0}$ and $\Omega_{\pm}$ can be written in terms of the form
factors as

$\displaystyle\Omega_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\alpha+\beta+\delta)-\frac{1}{3}\frac{\varpi}{\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}}}+\frac{1}{3}\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}},$
(31) $\displaystyle\Omega_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\alpha+\beta+\delta)+\frac{1+i\sqrt{3}}{6}\frac{\varpi}{\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}}}-\frac{1-i\sqrt{3}}{6}\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}},$
(32) $\displaystyle\Omega_{-}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\alpha+\beta+\delta)+\frac{1-i\sqrt{3}}{6}\frac{\varpi}{\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}}}-\frac{1+i\sqrt{3}}{6}\Big{(}\frac{\chi+\sqrt{4\varpi^{3}+\chi^{2}}}{2}\Big{)}^{\frac{1}{3}},$
(33)

where $\varpi$ and $\chi$ are given by

$\displaystyle\varpi$ $\displaystyle=$
$\displaystyle\alpha(\beta-\alpha)+\beta(\delta-\beta)+\delta(\alpha-\delta)-3(\gamma^{2}+\lambda^{2}+\sigma^{2}),$
(34) $\displaystyle\chi$ $\displaystyle=$
$\displaystyle(2\alpha-\beta-\delta)(2\beta-\delta-\alpha)(2\delta-\alpha-\beta)+54\gamma\lambda\sigma$
(35) $\displaystyle-$ $\displaystyle
9\alpha(2\lambda^{2}-\sigma^{2}-\gamma^{2})-9\beta(2\sigma^{2}-\gamma^{2}-\lambda^{2})-9\delta(2\gamma^{2}-\lambda^{2}-\sigma^{2}).$

Once the form factors are extracted from the polarization function, the
desired collective modes of the gluon can be obtained from the pole of the
effective propagator.

##  III Results

Figure 1: The real and imaginary parts of the form factors are plotted in
(a)–(f) as a function of $\omega/p$ for three different sets of anistotropy
tuple $\xi=(\xi_{a},\xi_{b})=(0,0),(0,10),$ and $(5,10)$. In each case, the
imaginary parts are shown with comparatively thicker style. The spheroidal set
is shown at fixed $\theta_{p}=\pi/4$ whereas for the ellipsoidal case,
$(\theta_{p},\phi_{p})=(\pi/4,\pi/3)$ is considered.

In this work we consider the hard-loop gluon polarization tensor in a
nonequilibrium QGP medium given by Mrowczynski:2000ed

$\displaystyle\Pi^{\mu\nu}(p^{0}=\omega,{\bm{p}})$ $\displaystyle=$
$\displaystyle
g_{s}^{2}\int\frac{d^{3}{\bm{k}}}{(2\pi)^{3}}\frac{K^{\mu}}{E_{k}}\frac{\partial
f({\bm{k}})}{\partial
K^{\rho}}\Big{[}\eta^{\rho\nu}-\frac{K^{\rho}P^{\nu}}{P\cdot
K+i0^{+}}\Big{]},$ (36)

where $g_{s}$ is the strong coupling constant and $k^{0}=E_{k}$ represents the
energy of the massless partonic degrees of freedom that modify the gluon
dispersion in presence of anisotropic medium. The effective distribution
function $f({\bm{k}})$ is given by

$\displaystyle
f({\bm{k}})=2N_{c}n_{g}({\bm{k}})+N_{f}\big{[}n_{q}({\bm{k}})+n_{\overline{q}}({\bm{k}})\big{]},$
(37)

where $n_{g}({\bm{k}})$ represents the gluon number density whereas
$n_{q}({\bm{k}})$ and $n_{\overline{q}}({\bm{k}})$ are the quark and antiquark
number densities respectively. $N_{c}$ and $N_{f}$ are respectively the number
of colors and quark flavors. A general method of constructing anisotropic
momentum distribution function is to transform the argument of an isotropic
momentum space distribution essentially by introducing deformations in a
parametrized way. A discussion on the generalization of the RS form to
ellipsoidal anisotropies can be found in Refs.Nopoush:2014pfa ;
Kasmaei:2018yrr . In our case we use the ellipsoidal momentum distribution
parametrized as Kasmaei:2016apv ; Kasmaei:2018yrr

$\displaystyle f_{\mbox{aniso}}({\bm{k}})\equiv
f_{\mbox{iso}}\big{(}\frac{1}{\Lambda}\sqrt{{\bm{k}}^{2}+\xi_{a}({\bm{k}}\cdot{\bm{a}})^{2}+\xi_{b}({\bm{k}}\cdot{\bm{b}})^{2}}\big{)},$
(38)

where $\Lambda$ represents a temperature-like scale which, in the equilibrium
limit, corresponds to the temperature. The gluon polarization tensor with such
anisotropic parton distributions can be written as

$\displaystyle\Pi^{\mu\nu}(\omega,{\bm{p}},\xi)$ $\displaystyle=$
$\displaystyle
m_{D}^{2}\int\frac{d\Omega}{4\pi}v^{\mu}\frac{v^{l}+\xi_{a}({\bm{v}}\cdot{\bm{a}})a^{l}+\xi_{b}({\bm{v}}\cdot{\bm{b}})b^{l}}{(1+\xi_{a}({\bm{v}}\cdot{\bm{a}})^{2}+\xi_{b}({\bm{v}}\cdot{\bm{b}})^{2})^{2}}\Big{[}\eta^{\nu
l}-\frac{v^{\nu}P^{l}}{\omega-{\bm{p}}\cdot{\bm{v}}+i0^{+}}\Big{]},$ (39)

where $\xi$ represents the anisotropy tuple $(\xi_{a},\xi_{b})$ and
$m_{D}^{2}=(N_{c}+N_{f}/2)\frac{g_{s}^{2}\Lambda^{2}}{3}$ corresponds to the
QCD Debye mass scale. In the above expression, the parton four velocity
$v^{\mu}=K^{\mu}/k$ and the components of ${\bm{a}}$, ${\bm{b}}$, ${\bm{v}}$,
and ${\bm{p}}$ in the rest frame of the medium are chosen as

$\displaystyle{\bm{a}}$ $\displaystyle=$ $\displaystyle(1,0,0),$ (40)
$\displaystyle{\bm{b}}$ $\displaystyle=$ $\displaystyle(0,0,1),$ (41)
$\displaystyle{\bm{v}}$ $\displaystyle=$
$\displaystyle(\sin\theta_{k}\cos\phi_{k},\sin\theta_{k}\sin\phi_{k},\cos\theta_{k}),$
(42) $\displaystyle{\bm{p}}$ $\displaystyle=$ $\displaystyle
p(\sin\theta_{p}\cos\phi_{p},\sin\theta_{p}\sin\phi_{p},\cos\theta_{p}),$ (43)

with $d\Omega$ representing the differential solid angle corresponding to the
internal angular coordinates $(\theta_{k},\phi_{k})$. It should be mentioned
here that another common choice of reference frame is the parton specific
coordinate as used in Ref. Kasmaei:2018yrr  where one reorients the axes so
that the polar angle is measured with respect to the parton momentum. However,
as we are interested in the evaluation of the form factors, the results do not
depend on any specific choice of the reference frame.

Now, as discussed in the previous section, the gluon self-energy can be
decomposed in terms of six independent basis tensors. Utilizing the
contraction relations among the basis tensors, the corresponding form factors
can be obtained in terms of the self-energy components given in Eq.(39) as a
two dimensional integral over the solid angle. However, because of the
transversality condition and the symmetry under the exchange of the free
indices, the number of independent components of the self-energy reduces to
six. Thus, any chosen set of six independent components will be sufficient to
determine the form factors. It should be mentioned here that, though in our
work, the form factors are obtained numerically using the conventional
quadrature routines, an alternative way involving the hypergeometric expansion
method is expected to be much more efficient Kasmaei:2018yrr  for this
purpose. The real and imaginary parts of the form factors are shown in Fig. 1,
for three different situations, namely the isotropic case with $\xi=(0,0)$,
the spheroidal case with $\xi=(0,10)$ and the ellipsoidal case with
$\xi=(5,10)$. The spheroidal case is shown for $\theta_{p}=\pi/4$ whereas for
the ellipsoidal anisotropy, $\theta_{p}$ and $\phi_{p}$, are chosen as $\pi/4$
and $\pi/3$, respectively. The imaginary parts of the form factors in all the
three cases exist only in the spacelike region. It should be noted that, when
the isotropic medium is considered, $\beta$ and $\delta$ become degenerate
whereas $\gamma$, $\lambda$, and $\sigma$ become zero. This results in two
distinct dispersive modes of the gluon among which one is degenerate, i.e.,
$\Omega_{0}=\alpha$ and $\Omega_{+}=\Omega_{-}=\beta=\delta$. The analytic
expression of the degenerate form factors is given by

$\displaystyle\beta=\delta=\Pi_{T}=\frac{m_{D}^{2}}{2}\frac{\omega^{2}}{p^{2}}\bigg{[}1-\frac{\omega^{2}-p^{2}}{2\omega
p}\ln\frac{\omega+p}{\omega-p}\bigg{]},$ (44)

whereas the dispersion for the other distinct mode can be obtained from

$\displaystyle\alpha=\Pi_{L}=\frac{m_{D}^{2}}{\tilde{u}^{2}}\bigg{[}1-\frac{\omega}{2p}\ln\frac{\omega+p}{\omega-p}\bigg{]},$
(45)

which are the familiar results of the gluon self-energy in isotropic thermal
medium Bellac:2011kqa . Here
$m_{D}^{2}=(N_{c}+N_{f}/2)\frac{g_{s}^{2}T^{2}}{3}$ and
$\tilde{u}^{2}=-p^{2}/P^{2}$. In the presence of spheroidal anisotropy,
$\sigma$ and $\lambda$ remain zero. However, $\beta$ and $\delta$ are no
longer degenerate. In that case the functions corresponding to the dispersive
modes simplify to

$\displaystyle\Omega_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}\alpha+\beta+\sqrt{(\alpha-\beta)^{2}+4\gamma^{2}}\bigg{)},$
(46) $\displaystyle\Omega_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{(}\alpha+\beta-\sqrt{(\alpha-\beta)^{2}+4\gamma^{2}}\bigg{)},$
(47) $\displaystyle\Omega_{-}$ $\displaystyle=$ $\displaystyle\delta.$ (48)

which may be compared with the modes obtained in Ref. Romatschke:2003ms  (see
for example Eq.(43) therein). Though we find a different combination of the
form factors in the expressions due to the different choice of our basis
tensors, it should be noticed that, in our case too, the arguments inside the
square root appear as a sum of two complete squares, thereby allowing similar
interpretation for the collective modes [see for example Sec. VI of
Romatschke:2003ms ]. When ellipsoidal anisotropy is considered, it is observed
that all the form factors are nonzero and they contribute in the gluon
dispersion. However, It should be noticed that, for the fixed set of
parameters chosen for the figure, the values of $\lambda$ and $\sigma$ are an
order of magnitude smaller than the other form factors.

Now, following Ref. Romatschke:2003ms , a mass scale corresponding to a given
form factor, say for example $\alpha$, can be defined in the static limit as

$\displaystyle m_{\alpha}^{2}=\lim_{\omega\rightarrow 0}\alpha.$ (49)

It is evident from Fig. 1 that the imaginary part of each of the form factors
vanishes at $\omega\rightarrow 0$ limit. Thus, the mass scales as defined
above are real quantities. Now, in a similar way, one can define the mass
scales corresponding to the gluon dispersive modes as

$\displaystyle m_{\Omega_{0,\pm}}^{2}=\lim_{\omega\rightarrow
0}\Omega_{0,\pm}(\omega,p,\theta_{p},\phi_{p}).$ (50)

Figure 2: Squared values of the mass scales corresponding to $\Omega_{0}$ and
$\Omega_{\pm}$ are plotted in (a)–(c) as a function of $\theta_{p}$ for three
different values of $\phi_{p}=\\{\pi/2,\pi/4,\pi/6\\}$ (shown in discrete
style) considering $(\xi_{a},\xi_{b})=(5,10)$. The solid line represents the
spheroidal case with $(\xi_{a},\xi_{b})=(0,10)$.

Their variation with the polar angle $\theta_{p}$ is shown in Fig. 2 for three
fixed values of $\phi_{p}=\\{\pi/2,\pi/4,\pi/6\\}$. For each modes, the
corresponding spheroidal version is also shown for comparison. It is evident
from the figures that the azimuthal symmetry of the spheroidal case is now
broken with the introduction of additional anisotropy direction. Consequently,
a nontrivial $\phi_{p}$ dependence can be observed in the mass scales. As can
be seen from the figure, the mass scale corresponding to $\Omega_{0}$ remains
positive throughout the range of $\theta_{p}$ values as also found in the
spheroidal case. Again, similar to the spheroidal anisotropy, negative values
in the mass scale is observed for $\Omega_{\pm}$ which correspond to
instability Romatschke:2003ms . However, it can be noticed that as the value
of $\phi_{p}$ approaches to $\pi/2$, $m^{2}_{\Omega_{\pm}}$ becomes positive
at smaller values of $\theta_{p}$. In other words, with larger deviation from
azimuthal anisotropy direction, the collective modes can be unstable only in
shorter window of $\theta_{p}$ values as compared to the spheroidal case.
However, it should be noted that, the above observation is made with a
particular set of anisotropy parameter where both $\xi_{a}$ and $\xi_{b}$ are
positive. An interesting situation occurs for negative value of anisotropy
parameter as shown in the left panel of Fig. 3. Here, the angular variation of
$m_{\Omega_{-}}^{2}$ is shown with $(\xi_{a},\xi_{b})=(-0.5,-0.9)$. In this
case, it is observed that, unlike the spheroidal case (which remains stable
throughout the $\theta_{p}$ range), the mass scale corresponding to
$\Omega_{-}$ can be negative, indicating an unstable mode. The corresponding
growth rate $\Gamma_{\Omega_{-}}$ can be obtained by solving the dispersion
with the replacement $\omega\rightarrow i\Gamma_{\Omega_{-}}$, i.e.,  by
solving the equation

$\displaystyle(\omega=i\Gamma_{\Omega_{-}})^{2}-p^{2}-\Omega_{-}(\omega=i\Gamma_{\Omega_{-}},p,\theta_{p},\phi_{p})=0.$
(51)

The corresponding solution of the growth rate is shown in the right panel of
Fig. 3 with $\theta_{p}$ and $\phi_{p}$ both fixed at $5\pi/12$. It should be
mentioned here that in this case the growth rate of the unstable mode has
amplitudes similar to the spheroidal case Romatschke:2003ms  whereas, with
positive anisotropy parameters, a several times larger growth rate can be
observed for certain angular values as reported in Ref. Kasmaei:2018yrr .

Figure 3: In the left panel, the squared value of the mass scale corresponding
to $\Omega_{-}$ is plotted as a function of $\theta_{p}$ for three different
values of $\phi_{p}=\\{\pi/4,\pi/3,5\pi/12\\}$ (shown in discrete style)
considering $(\xi_{a},\xi_{b})=(-0.5,-0.9)$. The spheroidal case (shown in
solid style) with $(\xi_{a},\xi_{b})=(0,-0.9)$ is also plotted for comparison.
In the right panel, the growth rate for the $\Omega_{-}$ mode is plotted as a
function of $p/m_{D}$ with $(\xi_{a},\xi_{b})=(-0.5,-0.9)$ and
$(\theta_{p},\phi_{p})=(5\pi/12,5\pi/12)$.

##  IV Summary

In this article, we have studied the gluon polarization in presence of
ellipsoidal momentum-space anisotropy. The momentum distribution function in
our case is parametrized with two anisotropy parameters $\xi_{a}$ and
$\xi_{b}$, represented together as an anisotropy tuple
$\xi=(\xi_{a},\xi_{b})$. This is a simple generalization of the spheroidal RS
form that has been extensively used in the literature. The general structure
of the gluon polarization tensor in presence of such ellipsoidal anisotropy
has been formulated and subsequently, the gluon effective propagator is
obtained. As shown earlier, from the pole of the effective propagator, the
three collective modes can be obtained in terms of the form factors. The
results obtained using our formulation are in agreement with the previous
study incorporating ellipsoidal anisotropy Kasmaei:2018yrr . It should be
mentioned here that it is not mandatory to consider the general structure of
the polarization function to obtain the collective modes, as, those can also
be obtained solving the characteristic equation directly (see for example
Kasmaei:2018yrr ). However, one of the important advantages of considering the
general structure is that, the collective modes, as in our case, can be
expressed in terms of the form factors which do not depend on the choice of
the frame of reference. As a consequence, it is possible to define mass scales
by taking the static limits of the functions characterizing the collective
modes ($\Omega_{0}$ and $\Omega_{\pm}$) in a similar fashion as done in case
of spheroidal anisotropy. The importance of such definition lies in the fact
that, the existence of instability can be inferred systematically by studying
the angular variations of the mass scales. More specifically, for the given
external angles, the negative value of the squared mass indicates that the
corresponding mode is unstable. In our analysis with $\xi_{a}=-0.5$ and
$\xi_{b}=-0.9$, we have observed that, unlike the spheroidal case, the mode
corresponding to $\Omega_{-}$ becomes unstable. The appearance of such
additional unstable mode in presence of ellipsoidal anisotropy may have
important influences on the isotropization of the QGP medium produced in HIC
experiments. As mentioned earlier, the formulation, as developed in this work,
will be particularly useful in the studies concerning the heavy quark
potential in presence of ellipsoidal anisotropy. The usual procedure to obtain
the hard-thermal loop resummed perturbative part of the heavy quark potential
is to consider the Fourier transform of the 00 component of the effective
propagator [obtained in Eq.(28)] in the static limits. Consequently, the
nontrivial angular dependence enters in the potential through the mass scales.
It should be noted that, in this work, only the retarded part of the gluon
self energy is considered. On the other hand, the imaginary part of the
potential can be obtained from the Feynman effective propagator in the real
time Keldysh formalism Nopoush:2017zbu . Due to the azimuthal angular
dependence of the mass scales, a nontrivial modification in the real as well
as in the imaginary parts of the heavy quark potential is expected in presence
of ellipsoidal anisotropy, which will be an interesting future direction.

##  V ACKNOWLEDGMENTS

R. G. is funded by University Grants Commission (UGC). B. K. is funded by
Department of Atomic Energy (DAE), India via the project TPAES. A. M.
acknowledges Najmul Haque for immense encouragement and support. A. M. would
like to acknowledge Science and Engineering Research Board (SERB) for funding.

## References

  * (1) H. T. Ding, F. Karsch and S. Mukherjee, Int. J. Mod. Phys. E 24, no.10, 1530007 (2015) doi:10.1142/S0218301315300076 [arXiv:1504.05274 [hep-lat]].
  * (2) C. Ratti, Rept. Prog. Phys. 81, no.8, 084301 (2018) doi:10.1088/1361-6633/aabb97 [arXiv:1804.07810 [hep-lat]].
  * (3) A. Bazavov et al. [USQCD], Eur. Phys. J. A 55, no.11, 194 (2019) doi:10.1140/epja/i2019-12922-0 [arXiv:1904.09951 [hep-lat]].
  * (4) E. Braaten and R. D. Pisarski, Nucl. Phys. B 337, 569-634 (1990) doi:10.1016/0550-3213(90)90508-B
  * (5) J. O. Andersen and M. Strickland, Annals Phys. 317, 281-353 (2005) doi:10.1016/j.aop.2004.09.017 [arXiv:hep-ph/0404164 [hep-ph]].
  * (6) N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, JHEP 05, 027 (2014) doi:10.1007/JHEP05(2014)027 [arXiv:1402.6907 [hep-ph]].
  * (7) N. Su, Int. J. Mod. Phys. A 30, 1530025 (2015) doi:10.1142/S0217751X15300252 [arXiv:1502.04589 [hep-ph]].
  * (8) C. Gale, S. Jeon and B. Schenke, Int. J. Mod. Phys. A 28, 1340011 (2013) doi:10.1142/S0217751X13400113 [arXiv:1301.5893 [nucl-th]].
  * (9) S. Jeon and U. Heinz, Int. J. Mod. Phys. E 24, no.10, 1530010 (2015) doi:10.1142/S0218301315300106 [arXiv:1503.03931 [hep-ph]].
  * (10) A. Jaiswal and V. Roy, Adv. High Energy Phys. 2016, 9623034 (2016) doi:10.1155/2016/9623034 [arXiv:1605.08694 [nucl-th]].
  * (11) D. T. Son and A. O. Starinets, Ann. Rev. Nucl. Part. Sci. 57, 95-118 (2007) doi:10.1146/annurev.nucl.57.090506.123120 [arXiv:0704.0240 [hep-th]].
  * (12) H. Nastase, [arXiv:0712.0689 [hep-th]].
  * (13) R. Peschanski, Acta Phys. Polon. B 39, 2479-2510 (2008) [arXiv:0804.3210 [hep-ph]].
  * (14) J. Sadeghi, B. Pourhassan and S. Heshmatian, Adv. High Energy Phys. 2013, 759804 (2013) doi:10.1155/2013/759804
  * (15) Relativistic Heavy Ion Physics, R. Stock (Ed.), Springer, Heidelberg, 2010, ISBN: 978-3-642-01538-0 (Print) 978-3-642-01539-7 (Online)” doi:10.1007/978-3-642-01539-7.
  * (16) W. Busza, K. Rajagopal and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339-376 (2018) doi:10.1146/annurev-nucl-101917-020852 [arXiv:1802.04801 [hep-ph]].
  * (17) M. Strickland, Pramana 84, no.5, 671-684 (2015) doi:10.1007/s12043-015-0972-1 [arXiv:1312.2285 [hep-ph]].
  * (18) M. Strickland, Acta Phys. Polon. B 45, no.12, 2355-2394 (2014) doi:10.5506/APhysPolB.45.2355 [arXiv:1410.5786 [nucl-th]].
  * (19) M. Alqahtani, M. Nopoush, R. Ryblewski and M. Strickland, Phys. Rev. Lett. 119, no.4, 042301 (2017) doi:10.1103/PhysRevLett.119.042301 [arXiv:1703.05808 [nucl-th]].
  * (20) M. Alqahtani, M. Nopoush and M. Strickland, Prog. Part. Nucl. Phys. 101, 204-248 (2018) doi:10.1016/j.ppnp.2018.05.004 [arXiv:1712.03282 [nucl-th]].
  * (21) S. Mrowczynski, A. Rebhan and M. Strickland, Phys. Rev. D 70, 025004 (2004) doi:10.1103/PhysRevD.70.025004 [arXiv:hep-ph/0403256 [hep-ph]].
  * (22) S. Mrowczynski and M. H. Thoma, Phys. Rev. D 62, 036011 (2000) doi:10.1103/PhysRevD.62.036011 [arXiv:hep-ph/0001164 [hep-ph]].
  * (23) P. Romatschke and M. Strickland, Phys. Rev. D 68, 036004 (2003) doi:10.1103/PhysRevD.68.036004 [arXiv:hep-ph/0304092 [hep-ph]].
  * (24) P. Romatschke and M. Strickland, Phys. Rev. D 70, 116006 (2004) doi:10.1103/PhysRevD.70.116006 [arXiv:hep-ph/0406188 [hep-ph]].
  * (25) B. Schenke and M. Strickland, Phys. Rev. D 74, 065004 (2006) doi:10.1103/PhysRevD.74.065004 [arXiv:hep-ph/0606160 [hep-ph]].
  * (26) L. Bhattacharya, R. Ryblewski and M. Strickland, Phys. Rev. D 93, no.6, 065005 (2016) doi:10.1103/PhysRevD.93.065005 [arXiv:1507.06605 [hep-ph]].
  * (27) L. Bhattacharya and P. Roy, Phys. Rev. C 78, 064904 (2008) doi:10.1103/PhysRevC.78.064904 [arXiv:0809.4596 [hep-ph]].
  * (28) L. Bhattacharya and P. Roy, Phys. Rev. C 79, 054910 (2009) doi:10.1103/PhysRevC.79.054910 [arXiv:0812.1478 [hep-ph]].
  * (29) M. Nopoush, Y. Guo and M. Strickland, JHEP 09, 063 (2017) doi:10.1007/JHEP09(2017)063 [arXiv:1706.08091 [hep-ph]].
  * (30) B. Krouppa and M. Strickland, Universe 2, no.3, 16 (2016) doi:10.3390/universe2030016 [arXiv:1605.03561 [hep-ph]].
  * (31) M. Mandal and P. Roy, Phys. Rev. D 86, 114002 (2012) doi:10.1103/PhysRevD.86.114002 [arXiv:1310.4657 [hep-ph]].
  * (32) M. Mandal and P. Roy, Phys. Rev. D 88, no.7, 074013 (2013) doi:10.1103/PhysRevD.88.074013 [arXiv:1310.4660 [hep-ph]].
  * (33) M. Mandal, L. Bhattacharya and P. Roy, Phys. Rev. C 84, 044910 (2011) doi:10.1103/PhysRevC.84.044910 [arXiv:1101.5855 [hep-ph]].
  * (34) V. Chandra and V. Sreekanth, Eur. Phys. J. C 77, no.6, 427 (2017) doi:10.1140/epjc/s10052-017-4992-5 [arXiv:1602.07142 [nucl-th]].
  * (35) M. Strickland, Braz. J. Phys. 37, 762-765 (2007) doi:10.1590/S0103-97332007000500021 [arXiv:hep-ph/0611349 [hep-ph]].
  * (36) A. Rebhan, P. Romatschke and M. Strickland, Phys. Rev. Lett. 94, 102303 (2005) doi:10.1103/PhysRevLett.94.102303 [arXiv:hep-ph/0412016 [hep-ph]].
  * (37) S. Mrowczynski, B. Schenke and M. Strickland, Phys. Rept. 682, 1-97 (2017) doi:10.1016/j.physrep.2017.03.003 [arXiv:1603.08946 [hep-ph]].
  * (38) D. Bazow, U. W. Heinz and M. Strickland, Phys. Rev. C 90, no.5, 054910 (2014) doi:10.1103/PhysRevC.90.054910 [arXiv:1311.6720 [nucl-th]].
  * (39) L. Tinti and W. Florkowski, Phys. Rev. C 89, no.3, 034907 (2014) doi:10.1103/PhysRevC.89.034907 [arXiv:1312.6614 [nucl-th]].
  * (40) M. Nopoush, R. Ryblewski and M. Strickland, Phys. Rev. C 90, no.1, 014908 (2014) doi:10.1103/PhysRevC.90.014908 [arXiv:1405.1355 [hep-ph]].
  * (41) M. Alqahtani, M. Nopoush and M. Strickland, Phys. Rev. C 92, no.5, 054910 (2015) doi:10.1103/PhysRevC.92.054910 [arXiv:1509.02913 [hep-ph]].
  * (42) M. Nopoush and M. Strickland, Phys. Rev. C 100, no.1, 014904 (2019) doi:10.1103/PhysRevC.100.014904 [arXiv:1902.03303 [nucl-th]].
  * (43) B. S. Kasmaei, M. Nopoush and M. Strickland, Phys. Rev. D 94, no.12, 125001 (2016) doi:10.1103/PhysRevD.94.125001 [arXiv:1608.06018 [hep-ph]].
  * (44) B. S. Kasmaei and M. Strickland, Phys. Rev. D 102, no.1, 014037 (2020) doi:10.1103/PhysRevD.102.014037 [arXiv:1911.03370 [hep-ph]].
  * (45) B. S. Kasmaei and M. Strickland, Phys. Rev. D 99, no.3, 034015 (2019) doi:10.1103/PhysRevD.99.034015 [arXiv:1811.07486 [hep-ph]].
  * (46) B. S. Kasmaei and M. Strickland, Phys. Rev. D 97, no.5, 054022 (2018) doi:10.1103/PhysRevD.97.054022 [arXiv:1801.00863 [hep-ph]].
  * (47) A. Dumitru, Y. Guo and M. Strickland, Phys. Lett. B 662, 37-42 (2008) doi:10.1016/j.physletb.2008.02.048 [arXiv:0711.4722 [hep-ph]].
  * (48) A. Dumitru, Y. Guo and M. Strickland, Phys. Rev. D 79, 114003 (2009) doi:10.1103/PhysRevD.79.114003 [arXiv:0903.4703 [hep-ph]].
  * (49) M. L. Bellac, “Thermal Field Theory,” doi:10.1017/CBO9780511721700