# The Lifshitz nature of the transition between the gap and gapless states of
a superconductor

Yuriy Yerin  Dipartimento di Fisica e Geologia, Universitá degli Studi di
Perugia, Via Pascoli, 06123 Perugia, Italy  Caterina Petrillo  Dipartimento di
Fisica e Geologia, Universitá degli Studi di Perugia, Via Pascoli, 06123
Perugia, Italy  A. A. Varlamov  CNR-SPIN, c/o Dipartimento DICII, Universitá
“Tor Vergata”, Viale del Politecnico 1, I-00133 Rome, Italy

###### Abstract

It is demonstrated that the known for a long time transition between the gap
and the gapless states in the Abrikosov-Gor’kov theory of a superconductor
with paramagnetic impurities is of the Lifshitz type, i.e. of the
$2\frac{1}{2}$ order phase transition. We reveal the emergence of a cuspidal
edge at the density of states surface $N(\omega,\Delta_{0})$ ($\Delta_{0}$ is
the value of the superconducting order parameter in the absence of magnetic
impurities) and the occurrence of the catastrophe phenomenon at the transition
point. We study the stability of such a transition with respect to the spatial
fluctuations of the magnetic impurities critical concentration $n_{s}$ and
show that the requirement for validity of its mean field description is
unobtrusive: $\nabla\left({\ln{n_{s}}}\right)\ll\xi^{-1}$ (here $\xi$ is the
superconducting coherence length). Finally, we show that, similarly to the
Lifshitz point, the $2\frac{1}{2}$ order phase transition should be
accompanied by the corresponding singularities. For instance, the
superconducting thermoelectric effect has a giant peak exceeding the normal
value of the Seebeck coefficient by the ratio of the Fermi energy and the
superconducting gap. The concept of the experiment for the confirmation of
$2\frac{1}{2}$ order transition nature is proposed. The obtained theoretical
results can be applied for the explanation of recent experiments with
lightwave-driven gapless superconductivity, for the new interpretation of the
disorder induced transition $s_{\pm}$-$s_{++}$ states via gapless state in
multi-band superconductors, for better understanding of the gapless color
superconductivity in quantum chromodynamics, the string theory.

Introduction. - In 1960, two seminal papers were published almost
simultaneously, which gave rise to new directions in the research fields of
superconductivity and fermiology AG1960 ; Lifshitz1960 .

In the first paper, Abrikosov and Gor’kov (AG), extending the Bardeen-Cooper-
Schrieffer (BCS) theory to the case of a superconducting alloy containing
paramagnetic impurities, demonstrated that the original BCS identification of
the phenomenon of superconductivity with the presence of the gap in the
quasiparticle spectrum is too restrictive, and, under some conditions, gapless
superconductivity can exist. According to the AG theory AG1960 ; Ambegaokar ,
the transition between gap and gapless regimes was governed by the
concentration of paramagnetic impurities and the properties of such
superconducting system were studied in the mean-field approximation. Gapless
superconductivity occurs in the very narrow interval of paramagnetic impurity
concentrations $0.912\,n_{c}<n<n_{c}$, where $n_{c}$ is the concentration that
completely suppresses the supercurrent flow. Later, it was recognized that the
gapless regime in a superconductor can be induced by different mechanisms
breaking the time-reversal symmetry: magnetic field Maki1968 , current
Maki1968 , proximity effect Hauser  and the light Yang2019 . However, the
order of this transition, to the best of our knowledge, was never discussed.

In the second of mentioned above papers I.M. Lifshitz Lifshitz1960  introduced
the notion of phase transition of fractional, $2\frac{1}{2}$, order. Also, it
was pointed out that by varying some external parameter (pressure or
concentration of the isovalent impurities) one can change the number of
components of topological connectivity of the Fermi surface (FS), which is
accompanied, according to the Ehrenfest terminology Jaeger1998 , by the
$2\frac{1}{2}$ order phase transition. Further studies of these, named today
as Lifshitz’s, transitions revealed that they are supplemented by
singularities in various properties of the system Varlamov1989 ; Blanter1994 .

The ideas proposed 60 years ago remain still requested in modern studies on
the stability of current-biased superconducting wires (see OVKG20  and
references therein), transformations of the complex heavy fermion Fermi
surfaces due to magnetic field effects Pfau2017PRL ; Pourret2019JPSJ , etc.
Moreover, the concept of a connection between the topological properties of
the different materials exhibiting gapless states and the occurrence of the
exotic Lifshitz transitions was recently discussed in literature based on very
general topological arguments. Examples are given by Dirac and Weyl materials,
and even more exotic systems (see the reviews Volovik1 ; Volovik2 ).

In this article we aim at framing these concepts into a unified description
and show that the known for a long time transition between gap and gapless
superconducting states is the phase transition of the Lifshitz type, i.e. of
the $2\frac{1}{2}$ order. This will be proved by a very simple approach, in
spirit of the fundamental paper Lifshitz1960 , just analyzing the properties
of the free energy in a superconductor containing paramagnetic impurities.

Further, we study the requirements on the homogeneity of the paramagnetic
impurities concentration, which is necessary for the validity of the standard
mean field approximation in the description of the “gap-gapless” transition
used in AG1960 , and prove the stability of the transition with respect to
these fluctuations.

Finally, we argue that such a transition would be accompanied by the
appearance of singularities in several properties, in particular an anomalous
growth of the thermoelectric effect (seeGalperin1973 ; Galperin1974 ;
Zavaritskii1974 ; Ginzburg1978 ) close to the critical concentration, which is
valuable for the experimental verification of the proposed connection.

It is important to note that we confine our study to the case of a s-wave
isotropic superconductor and do not consider unconventional and exotic pairing
symmetries.

Free energy and phase transition. - We start from the expression for the free
energy close to the transition between gapless and gap regimes at $T=0$ (see
Maki1968 ; Skalski1964 ), that is

${F_{s-n}}\\!=\\!-\frac{{N\left(0\right){{\Delta}^{2}}}}{2}\left\\{\begin{gathered}1-\frac{\pi}{2}\zeta+\frac{2}{3}{\zeta^{2}},{\text{
}}\zeta\leqslant 1\hfill\\\
1\\!-\zeta\arcsin{\zeta^{-1}}\\!+\\!{\zeta^{2}}\left({1\\!-\\!\sqrt{1\\!-\\!{\zeta^{-2}}}}\right)\hfill\\\
-\frac{1}{3}{\zeta^{2}}\left({1-{{\left({1-{\zeta^{-2}}}\right)}^{3/2}}}\right),{\text{
}}\zeta>1{\text{ }}\hfill\\\ \end{gathered}\right.$ (1)

$\Delta=\Delta(\tau_{s})$ is the order parameter in the presence of impurities
($\Delta\in\mathbb{R}$) and $N(0)=\frac{mp_{F}}{\pi^{2}\hbar^{3}}$ is the
density of states (DOS) at the Fermi level. The parameter

$\zeta=(\tau_{s}\Delta)^{-1},$ (2)

with $\tau_{s}$ as the electron spin-flip scattering time due to the presence
of paramagnetic impurities, governs the phase transition between the gap and
gapless states. Namely, when $\zeta<1$ the energy gap $\Delta_{g}$ in the
quasiparticle spectrum of a superconductor has a nonzero value, while for
$\zeta\geq 1$ the energy gap remains identically equal to zero and the gapless
state is realized. At the same time the order parameter $\Delta$ is different
from zero and the phenomenon of supercurrent flow occurs. The critical point
$\zeta=1$ separates the gap and the gapless states. We remark that the authors
of AG1960  were the first who pointed out at the importance of making a
distinction between the order parameter $\Delta$ and the energy gap
$\Delta_{g}$ existing in the quasiparticles spectrum.

Figure 1: The second derivative of Eq. (1). The kink is clearly observed at
$\zeta=1$. Purple and cyan colors in the background of the plot illustrate
separation between gap and gapless states respectively. Weak blur near
$\zeta=1$ represents smearing of the transition due to spacial fluctuations of
the magnetic impurities concentration (see the corresponding section of the
paper).

To elucidate what is the order of the phase transition, we first studied the
behavior of the free energy (Eq. 1) and its derivatives over the parameter
$\zeta$ that drives the transition. It turns out that the free energy together
with its first and second derivatives remain continuous function at the
transition point $\zeta=1$. However, the plot of the second derivative
$\partial^{2}F_{s\\!-\\!n}/\partial\zeta^{2}$ unambiguously shows the kink at
$\zeta=1$ (Fig. 1). Moreover, from the expression of the third derivative

$\frac{{{\partial^{3}}{F_{s-n}}}}{{\partial{\zeta^{3}}}}=N\left(0\right){\Delta^{2}}\left\\{\begin{gathered}0,{\text{
}}\zeta\leqslant 1\hfill\\\
\frac{1}{{{\zeta^{4}}\sqrt{{\zeta^{2}}-1}}},{\text{ }}\zeta>1{\text{
}}.\hfill\\\ \end{gathered}\right.$ (3)

one can see the occurrence of the characteristic discontinuity in it with the
square root singularity from the gapless side. I.e., the situation is
completely analogous to the Lifshitz $2\frac{1}{2}$ order phase transitions in
metals. The analogy is also confirmed by DOS dependence on the parameter
$\zeta$ driving the transition. It was shown Maki1968 ; Ambegaokar ;
Skalski1964 that the quasiparticle DOS of the superconductor
$N_{s}\left(\omega\right)$ remains finite at $\omega=0$ and has a typical cusp
for $2\frac{1}{2}$ order phase transition at $\zeta=1$

$N_{s}\left(0\right)=N\left(0\right)\frac{\sqrt{\zeta+1}}{\zeta}\sqrt{\zeta-1}.$
(4)

The appropriate interpretation of the gap-gapless transition can be given
studying the transformation of the surface $N\left(\omega,\Delta_{0}\right)$
in the phase space $\omega$-$\Delta_{0}$, where $\Delta_{0}$ is the value of
the superconducting order parameter in the absence of magnetic impurities. For
this purpose we start with the general expression AG1960 ; Ambegaokar

$N\left(\omega,\zeta\right)=N\left(0\right){\zeta^{-1}}\operatorname{Im}u,$
(5)

where $u$ is given by

$\frac{\omega}{{\Delta}}=u\left({1-\frac{\zeta}{{\sqrt{1-{u^{2}}}}}}\right),$
(6)

and the expression which determines the order parameter $\Delta$ at $T=0$
AG1960 ; Maki1968

$\ln\left({\frac{\Delta}{{{\Delta_{0}}}}}\right)=\left\\{\begin{gathered}-\frac{\pi}{4}\zeta,{\text{
}}\zeta\leqslant 1\hfill\\\
-\operatorname{arccosh}\zeta-\frac{1}{2}\left({\zeta\arcsin{\zeta^{-1}}-\sqrt{1-{\zeta^{-2}}}}\right),{\text{
}}\hfill\\\ \zeta>1.{\text{ }}\hfill\\\ \end{gathered}\right.$ (7)

Figure 2: Evolution of the DOS in the $\omega$-$\Delta_{0}$ space from the gap
state with $\zeta=0.1$ (a), through the collapse of the energy gap $\zeta=1$
(b), to the gapless state with $\zeta=1.75$ (c). Arbitrary units for $\omega$
and $\Delta_{0}$ are used in the legends.

Based on Eqs. (5) and (7) one can track the evolution of the function
$N(\omega,\Delta_{0})$ in the $\omega$-$\Delta_{0}$ phase space via three
stages that are characterized by dissimilar surfaces for $\zeta<1$, $\zeta=1$
and $\zeta>1$ (see Fig. 2). The first stage corresponds to the gap state with
$\zeta<1$ and with the characteristic narrowing hollow between two glued
sheets of the DOS surfaces at $\Delta_{g}=0$ (Fig. 2a). The collapse of the
energy gap when $\zeta=1$, and the subsequent emergence of a topological
feature in the form of the pleat known as a cuspidal edge at $\omega=0$, is
shown in Fig. 2b. This feature allows to speculate about the occurrence of the
catastrophe phenomenon in the $\omega$-$\Delta_{0}$ space over the gap-gapless
phase transition Arnold1 ; Arnold2 ; Saji ; Manfredo . Finally, the last stage
with $\zeta>1$ corresponds to the gapless state with the gradual degradation
of the DOS curved surface to a plane for $\zeta\to\infty$ (Fig. 2c).

In this representation one can freely “travel” over the each surface
$N(\omega,\Delta_{0})$ by changing the variables $\omega$ and $\Delta_{0}$
while keeping $\zeta=\zeta(\Delta_{0},\tau_{s})=const$ and adjusting the value
of $\tau_{s}$ for each $\Delta_{0}$ to satisfy the constancy of the given
value of $\zeta$. One should not be surprised by the manipulating of
$\tau_{s}$ for each $\Delta_{0}$ in order to carry $\zeta=$ fixe at whole
considered surface. In some sense during the study the Lifshitz transition the
experimentalists do the same, for example, investigation of the anomalous
behavior of thermopower in ${\text{L}}{{\text{i}}_{{\text{1 -
x}}}}{\text{M}}{{\text{g}}_{\text{x}}}$ alloy in a dependence of $x$ close to
transition at $x=$ 0.19 (see Ref. Egorov, ).

Therefore, one can conclude that while the Lifshitz transition is governed by
the parameter $z=\mu-\mu_{c}$ ($\mu$ is the chemical potential) Blanter1994 ;
Volovik1 ; Volovik2 ; Volovik3 ; Varlamov1 , the driving parameter of the
transition under consideration is $\zeta-1$.

Smearing of the transition due to spacial fluctuations of the magnetic
impurities concentration. - In the case of the order parameter varying in
space, Eq. (1) for the free energy can be generalized by adding heuristically
the corresponding gradient term, like in the Ginzburg-Landau theory, that is

$\displaystyle F$ $\displaystyle=$
$\displaystyle-\frac{N(0)}{2}\left\\{a\Delta^{2}+\frac{1}{4m}(\nabla\Delta)^{2}\right\\}$
$\displaystyle=$
$\displaystyle-\frac{N(0)}{2}\left\\{a\Delta^{2}+\frac{1}{4m}\left(\frac{d\Delta}{d\zeta}\right)^{2}\left(\frac{d\zeta}{dn_{s}}\right)^{2}(\nabla
n_{s})^{2}\right\\},$

where $a=1-\frac{\pi}{2}\zeta+\frac{2}{3}\zeta^{2}$ and we recall that $n_{s}$
is the concentration of magnetic impurities.

Here we already attributed the variation of the value of order parameter to
the spacial inhomogeneity of the paramagnetic impurities distribution,
elucidating the corresponding gradient in the last term of Eq. (The Lifshitz
nature of the transition between the gap and gapless states of a
superconductor). One can neglect the impurities concentration fluctuations
until the contribution of the “kinetic energy” remains small in comparison to
the superconducting condensation energy.

Correspondingly, comparing the second term in Eq. (The Lifshitz nature of the
transition between the gap and gapless states of a superconductor) with the
first one and using the expression determing the order parameter $\Delta$ at
$T=0$ given by Eq. (7) we find (see [SM, ]) that the fluctuations of the
impurities concentration remain insignificant until

$\frac{\nabla
n_{s}}{n_{s}}\ll\frac{1}{\xi}\left(1-\frac{\pi}{2}\zeta+\frac{2}{3}\zeta^{2}\right)\left(\frac{4}{\pi\zeta}-1\right).$
(9)

Here $\xi$ is the superconducting coherence length. This evaluation is valid
close to the transition point $\zeta=1$ (indeed, the limit of superconductor
without paramagnetic impurities $\zeta\rightarrow 0$ does not make sense in
such consideration), i.e.

$\frac{d\left[\ln n_{s}(r)\right]}{dr}\ll\frac{1}{\xi}.$ (10)

Thermoelectric effect. - It is well known that the Lifshitz transition in
normal metals is accompanied by a giant asymmetric peak in the Seebeck
coefficient Egorov1983 ; Vaks1981 ; Varlamov1985 . Despite the opinion
prevailing in the early period of the study of superconductivity concerning
the vanishing of all conventional thermoelectric properties, today we know
that a wide variety of interesting thermoelectric effects can exist in
superconductors Van1982 . Among them is the quantization of the magnetic flux
passing through the loop consisting of two different superconductors whose
junctions are at different temperatures. As demonstrated in Galperin1974  the
correction to the integer number of flux quanta appears to depend on the
temperature difference and thermoelectric coefficients of the superconductors
in their normal state. Hence, one could expect the giant growth of this effect
when one of the ring legs is close to the gap-gapless transition.

In order to demonstrate this we will calculate the corresponding quasiparticle
contribution to the thermoelectric coefficient following the scheme proposed
by Ambegaokar and Griffin to calculate the corresponding thermal conductivity
(see Ref. Ambegaokar, ). We perform our calculation in the assumption of
validity of the weak enough scattering and applicability of the Born
approximation. The thermoelectric coefficient $\alpha$ relating the
quasiparticle current to the temperature gradient, can be expressed in the
form

$\begin{gathered}\alpha=-\frac{{eN\left(0\right)v_{F}^{2}}}{{{4T^{2}}}}\int\limits_{-\infty}^{+\infty}\frac{\omega
d\omega}{\cosh^{2}\left({\frac{{\beta\omega}}{2}}\right)}\times\hfill\\\
\frac{h\left({\omega,\Delta,\zeta}\right)}{{\operatorname{Im}\left\\{{\Omega\left({\omega,\Delta,\zeta}\right)+\frac{{\text{i}}}{{2{\tau_{tr}}}}+\frac{{\text{i}}}{{{\tau_{s}}}}\left[{1-h\left({\omega,\Delta,\zeta}\right)}\right]}\right\\}}},\hfill\\\
\end{gathered}$ (11)

where $\tau_{tr}$ is the transport collision time AG1960  that enters in the
conductivity of a normal metal, $v_{F}$ is the Fermi velocity and
$\beta=\frac{1}{{{k_{B}}T}}$ is the inverse temperature. The functions
${h\left({\omega,\Delta,\zeta}\right)}$ and
${\Omega\left({\omega,\Delta,\zeta}\right)}$ are given by

$\displaystyle\frac{{\Omega\left({\omega,\Delta,\zeta}\right)}}{\Delta}$
$\displaystyle=$ $\displaystyle\sqrt{{u^{2}}-1}-{\text{i}}\zeta,$ (12)
$\displaystyle h\left({\omega,\Delta,\zeta}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[{1+\frac{{{{\left|u\right|}^{2}}-1}}{{\left|{{u^{2}}-1}\right|}}}\right],$
(13)

where we recall that parameter $u$ is defined by Eq. (6).

Figure 3: Thermoelectric coefficient as a function of $\zeta$ for different
values of $\tau_{s}\Delta_{0}$ for the given dimensionless inverse temperature
$\beta\Delta_{0}=50$ (up) and for different $\beta\Delta_{0}$ for the fixed
$\tau_{s}\Delta_{0}=2.25$ (bottom). For both plots $\tau_{tr}\Delta_{0}=1$.

At temperatures close to zero (large values of $\beta$) the main contribution
to the integral in Eq. (11) comes from the low frequencies domain
$\omega\lesssim\beta^{-1}$. The numerical analysis of the Eq. (11) shows the
dramatic enhancement of the thermoelectric coefficient approaching the
transition from the gap side (see Fig. 3). The absolute value of the peak
linearly decreases with the temperature, still remaining $\epsilon_{F}/\Delta$
times larger than the background value.

In the immediate vicinity of the transition one can expand the parameter $u$
and consequently the functions
$\Omega\left({\omega,\left|\Delta\right|,\zeta}\right)$ and
$h\left({\omega,\left|\Delta\right|,\zeta}\right)$ for small values of
$\omega$. In result one can obtain the asymptotic behavior of the
thermoelectric coefficient close to the phase transition for both the gap and
the gapless states.

From the gap side ($\zeta\rightarrow 1-$) one finds (see Supplemental Material
[SM, ]) that the thermoelectric coefficient takes the form

$\alpha=\frac{{4\sqrt{2}{\pi^{2}}}}{3e}\frac{T}{\Delta}\frac{\sigma_{n}\tau_{s}}{\tau_{s}+2\tau_{tr}}\sqrt{1-{\zeta^{\frac{1}{3}}}},$
(14)

where $\sigma_{n}=\frac{2}{3}N(0)v_{F}^{2}e^{2}\tau_{tr}$. Eq. (14) determines
the magnitude of the Seebeck coefficient in the gap state. Recalling that the
value of Seebeck coefficient in the normal metal is
$S_{n}=\frac{\pi^{2}k_{B}}{3e}\frac{T}{E_{F}}$ one can find that $S_{g}$ is
giant with respect to the latter by the parameter $E_{F}/\Delta$:

${S_{g}}=\frac{\alpha}{{{\sigma_{n}}}}=\frac{{{2^{5/2}}{\tau_{s}}}}{{{\tau_{s}}+2{\tau_{tr}}}}\sqrt{1-{\zeta^{\frac{1}{3}}}}\left({\frac{{{E_{F}}}}{\Delta}}\right){S_{n}}.$
(15)

When performing the same procedure from the gapless side of the transition it
can be disappointing to find $\alpha\equiv 0$. Formally this is related to the
oddness of the integrand function over $\omega$ (see Ref. SM ) in this region.
Yet, the obtained result does not mean that the thermoelectric coefficient
here turns identically zero: in our expansions we did not retain terms of the
order $\omega/E_{F}$, hence the thermoelectric effect from the right of
transition point can be comparable to its normal background.

The results of the numerical calculations of the thermoelectric coefficient
based on Eq. (11) are shown in Fig. 3. For the evaluation of $\alpha$ we used
the dependence of the order parameter modulus as a function of $\zeta$ at zero
temperature given by Eq. (7) and, for large values of $\beta$ (i.e. in the
vicinity of $T=0$), we assumed the temperature variation of $\Delta$ to be
very weak and approximated by Eq. (21) in Ref. SM, .

There are several remarkable hallmarks of the found effect. First of all, the
thermoelectric coefficient has a giant magnitude in the gap region. Second,
the peak is asymmetric, and, third, the peak is shifted from the transition
point into the gap domain ($\zeta<1$) when the temperature increases ($\beta$
decreases). All these features are characteristic also for the Seebeck signal
behaviour close to the $2\frac{1}{2}$ phase transitions Varlamov1989  and can
be considered as the smoking gun for the experimental verification of the
proposed phenomenon.

We should note that a similar strong enhancement of the thermoelectric
coefficient in the presence of magnetic impurities was theoretically predicted
in Ref. Digor, ; Kalenkov, . However, the authors of Refs. Digor, ; Kalenkov,
did not relate the revealed giant thermoeffect to the manifestation of the
$2\frac{1}{2}$ phase transition. They specified that this phenomenon is caused
by violation of the symmetry between electron-like and hole-like excitations
due to formation of the subgap Andreev bound states in the vicinity of
magnetic impurities Kalenkov . Although different assumptions have been used
in Refs. Digor,  and Kalenkov,  for the calculation of the thermoelectric
coefficient the effect of its enhancement remain the same on the qualitative
level.

From the experimental point of view the detection of such a phase transition
can be performed by means of placing in magnetic field a ring, one half of
which is a gap superconductor with the concentration of magnetic impurities
close to the transition value and the other half is an arbitrary
superconductor. In this case, provided superconducting contacts are kept at
different but low temperatures, anomaly strong thermoelectric current is
induced inside the ring and the measured magnetic flux should considerably
deviate from the integer values of the magnetic flux quantum $\Phi_{0}$
Zavaritskii1974 . An alternative method for detecting the discussed transition
can be a jump in the derivative of the specific heat capacity with respect to
the impurity concentration.

The effect of impurity scattering of the DOS dependence on energy and the
anisotropy degree has been investigated in Ref. Hohenberg,  for the case of
superconductors with the anisotropic gap. It was shown that presence of a
relatively small concentration of impurities leads to the isotropisation of
the DOS and decrease of its smearing over energy. With the further increase of
the concentration of impurities the region of smearing shrinks to zero. This
means that the superconductor becomes effectively isotropic. A full analysis
of this problem is outside of the scope of the present paper and therefore
left for future studies.

Conclusions. - We have demonstrated that the known for a long time transition
between the gap and the gapless states of a superconducting alloy with
paramagnetic impurities is the phase transition of the $2\frac{1}{2}$ order.
We have shown that the mean-field approximation used in the Abrikosov-Gor’kov
theory AG1960  is very stable: fluctuations of the impurities concentration
remain irrelevant in the logarithmic scale. Finally, such a phase transition
can be detected by the giant (by the parameter $E_{F}/\Delta$) thermoelectric
effect possessing the characteristic features which would clearly distinguish
it from others. We have proposed experiments for detection of such an effect
and the subsequent confirmation of the $2\frac{1}{2}$ phase transition nature.

Our theoretical results may help to take a fresh look at recent experiments
with light-wave-driven gapless superconductivity Yang2019 , for the new
interpretation of the theoretically predicted disorder induced transition
$s_{\pm}$-$s_{++}$ states via gapless phase in multi-band superconductors
Barzykin ; Efremov and can be useful for the understanding gapless color
superconductivity in quantum chromodynamics and the string theory Alford . In
the case of a dirty multi-band superconductor with increasing of the
nonmagnetic impurities concentration, one of the gaps is seen to close,
leading to a finite residual DOS, followed by a reopening of the gap. Such a
behavior allows to speculate about the Lifshitz origin of $s_{\pm}$-$s_{++}$
transition. For a color superconductor it was shown that, at zero temperature
and small values of the strange quark mass, the ground state of neutral quark
matter corresponds to the so-called color-flavor-locked phase. At some
critical value of the strange quark mass, there is a transition to the gapless
color-flavor-locked phase, where the energy gap in the quasiparticle spectrum
is not mandatory Alford ; Ruester . As in the case of multi-band
superconductivity one can again speculate about the emergence of the Lifshitz
nature of the transition in the phase diagram of the neutral quark matter.

Acknowledgements. - A.V. is grateful to Yu. Galperin, A. I. Buzdin, Aviad
Frydman, S. Frolov, S. Bergeret, and O. Dobrovolskiy for valuable and fruitful
discussions. C.P. and Y.Y. acknowledge support by the CarESS project. A.V.
acknowledges the financial support under the STSM COST Action CA16218
Nanoscale Coherent Hybrid Devices For Superconducting Quantum Technologies.

## References

  * (1) A.A. Abrikosov, L.P.Gor’kov,Contribution to the theory of superconducting alloys with paramagnetic impurities, Sov. Phys. JETP 12, 1243 (1961).
  * (2) I.M. Lifshitz, Anomalies of electron characteristics of a metal in the high pressure region, Zh. Eksp. teor. Fiz. 38, 1569 (1960) [Soviet Physics JETP 11, 1130].
  * (3) Vinay Ambegaokar and Allan Griffin, Theory of the Thermal Conductivity of Superconducting Alloys with Paramagnetic Impurities, Phys. Rev. 137, A1151 (1965) doi:10.1103/PhysRev.137.A1151.
  * (4) K. Maki, Gapless Superconductivity, Chapter 18, ”Superconductivity” Edited by Parks, 1968.
  * (5) J. J. Hauser, Gapless Superconductivity Induced by the Proximity Effect, Phys. Rev. 164, 558 (1967) doi:10.1103/PhysRev.164.558.
  * (6) X. Yang, C. Vaswani, C. Sundahl, M. Mootz, L. Luo, J. H. Kang, I. E. Perakis, C. B. Eom and J. Wang, Lightwave-driven gapless superconductivity and forbidden quantum beats by terahertz symmetry breaking, Nature Photonics 13, 707 (2018) doi:10.1038/s41566-019-0470-y.
  * (7) G. Jaeger, The Ehrenfest Classification of Phase Transitions: Introduction and Evolution, Arch. Hist. Exact Sci. 53, 51 (1998) doi:10.1007/s004070050021.
  * (8) A.A. Varlamov, V. Egorov, and A. Pantsulaya, Kinetic properties of metals near electronic topological transitions $(2\frac{1}{2}$-order transitions), Adv. Phys. 38, 469 (1989) doi:10.1080/00018738900101132.
  * (9) Y.M. Blanter, M.I. Kaganov, A.V. Pantsulaya, and A.A. Varlamov, The theory of electronic topological transitions, Phys. Rep. 245, 159 (1994) doi:10.1016/0370-1573(94)90103-1.
  * (10) Yu.N. Ovchinnikov, A.A. Varlamov, G.J. Kimmel, and A. Glatz, Instabilities of the normal state in current-biased narrow superconducting strips, Phys. Rev. B 101, 014511 (2020) doi:10.1103/PhysRevB.101.014511.
  * (11) H. Pfau, R. Daou, S. Friedemann, S. Karbassi, S. Ghannadzadeh, R. Küchler, S. Hamann, A. Steppke, D. Sun, M. König, A.P. Mackenzie, K. Kliemt, C. Krellner, and M. Brando, Cascade of Magnetic-Field-Induced Lifshitz Transitions in the Ferromagnetic Kondo Lattice Material ${\mathrm{YbNi}}_{4}{\mathrm{P}}_{2}$, Phys. Rev. Lett. 119, 126402 (2017) doi:10.1103/PhysRevLett.119.126402.
  * (12) A. Pourret, S.G. Sharapov, T.D. Matsuda, G. Knebel, G. Zwicknagl, and A.A. Varlamov, Transport Spectroscopy of the Field Induced Cascade of Lifshitz Transitions in YbRh2Si2, J. Phys. Soc. Jpn. 88, 104702 (2019) doi:10.7566/JPSJ.88.104702.
  * (13) G.E. Volovik, Exotic Lifshitz transitions in topological materials, Phys.Usp. 61, 89 (2018) doi:10.3367/UFNe.2017.01.038218.
  * (14) G.E. Volovik, Topological Lifshitz transitions, Low Temp. Phys. 43, 47 (2017) doi:10.1063/1.4974185.
  * (15) Yu.M. Galperin, V.L. Gurevich, V.I. Kozub, Acoustoelectric and Thermoelectric Effects in Superconductors, Sov. Phys. JETP Lett. 17, 476 (1973).
  * (16) Yu.M. Galperin, V.L. Gurevich, V.I. Kozub, Nonlinear acoustic effects in superconductors, Sov. Phys. JETP 65, 1045 (1974).
  * (17) N.V. Zavaritskii, Observation of superconducting current excited in a superconductor by heat flow, Pis’ma Zh. Eksp. Teor. Fiz. 19, 205 (1974) [JETP Lett. 19, 126 (1974)].
  * (18) V.L. Ginzburg , G.F. Zharkov, A.A. Sobyanin, Thermoelectric phenomena in superconductors and thermomechanical circulation effect in a superfluid liquid, Pis’ma Zh. Eksp. Teor. Fiz. 20, 223 (1974) [JETP Lett. 20, 97 (1974)].
  * (19) S.Skalski, O. Betbeder-Mattibet, and P.R. Weiss, Properties of Superconducting Alloys Containing Paramagnetic Impurities Phys. Rev. 136, A1500, (1964) doi:10.1103/PhysRev.136.A1500.
  * (20) Kuang Zhang and G.E. Volovik,Lifshitz transitions via the type-II Dirac and type-II Weyl points, Pis’ma ZhETF 105, 504 (2017) [JETP Lett. 105, 519 (2017)] doi:10.1134/S0021364017080094.
  * (21) A. A. Varlamov, Y. M. Galperin, S. G. Sharapov, Yuriy Yerin, Concise guide for electronic topological transitions, Low Temperature Physics 47, 672 (2021) doi:10.1063/10.0005556.
  * (22) V.I. Arnold, Alexander Varchenko, S.M. Gusein-Zade, Singularities of Differentiable Maps. Volume I: The Classification of Critical Points Caustics and Wave Fronts, (Monographs in mathematics ; Vol. 82), Birkhauser Boston Basel Stuttgart (1985) doi:10.1007/978-1-4612-5154-5.
  * (23) V.I. Arnold, Basic Concepts. In: Arnol’d V.I. (eds) Dynamical Systems V. Encyclopaedia of Mathematical Sciences, vol 5., Springer-Verlag Berlin Heidelberg (1994) doi:10.1007/978-3-642-57884-7_5.
  * (24) Kentaro Saji, Masaaki Umehara, Kotaro Yamada, The geometry of fronts, Annals of Mathematics 169, 491 (2009) doi:10.4007/annals.2009.169.491.
  * (25) Manfredo Perdigao do Carmo, Differential geometry of curves and surfaces. Upper Saddle River, N.J.: Prentice-Hall (1976).
  * (26) V. S. Egorov and A. N. Fedorov, Thermopower of lithium-magnesium alloys at the 2 1/2-order transition , Sov. Phys. JETP 58, 959 (1983).
  * (27) See the Supplemental Material for the detailed derivation of Eqs. (9), (10) and (14) and the low-temperature behavior of the order parameter.
  * (28) V. S. Egorov and A.N. Fedorov, Thermopower of lithium-magnesium alloys at the 2 1/2 -order transition, Zh. Eksp. Teor Phys. 85, 1647 (1983) [Sov. Phys. JETP 58, 959 (1983)].
  * (29) V. G. Vaks, A.V. Trefilov, and S.V. Fomichev, Singularities of resistivity and thermoelectric power of metals in phase transitions of order 2$\frac{1}{2}$, Th. Eksp. Teor. Fiz. 80, 1613 (1981) [Sov.Phys. JETP 53, 830 (1981)].
  * (30) A. A. Varlamov and A.V. Pantsulaya, Anomalous kinetic properties of metals near the Lifshitz topological transition, Sov. Phys. JETP 62, 1263 (1985).
  * (31) D. J. Van Harlingen, Thermoelectric effects in the superconducting state, Physica B+C 109–110, 1710 (1982) doi:10.1016/0378-4363(82)90195-4.
  * (32) D. F. Digor and L. Z. Kon, Thermoelectric effect in superconducting alloys, Theoretical and Mathematical Physics 97, 1137 (1993) doi:10.1007/BF01014806.
  * (33) Mikhail S. Kalenkov, Andrei D. Zaikin, and Leonid S. Kuzmin, Theory of a Large Thermoelectric Effect in Superconductors Doped with Magnetic Impurities, Phys. Rev. Lett. 109, 147004 (2012) doi:10.1103/PhysRevLett.109.147004.
  * (34) P. Hohenberg, Anisotropic Superconductors with Nonmagnetic Impurities, Sov. Phys. JETP 18, 834 (1964).
  * (35) Victor Barzykin, Magnetic-field-induced gapless state in multiband superconductors, Phys. Rev. B 79, 134517 (2009) doi:10.1103/PhysRevB.79.134517.
  * (36) D. V. Efremov, M. M. Korshunov, O. V. Dolgov, A. A. Golubov, P. J. Hirschfeld, Disorder-induced transition between s± and s++ states in two-band superconductors, Phys. Rev. B 84, 180512(R) (2011) doi:10.1103/PhysRevB.84.180512.
  * (37) Mark Alford, Jürgen Berges, and Krishna Rajagopal, Gapless Color Superconductivity, Phys. Rev. Lett. 84, 598 (2000) doi:10.1103/PhysRevLett.84.598
  * (38) Stefan B. Ruester, Verena Werth, Michael Buballa, Igor A. Shovkovy, Dirk H. Rischke, The phase diagram of neutral quark matter: Self-consistent treatment of quark masses, Phys. Rev. D 72, 034004 (2005) doi:10.1103/PhysRevD.72.034004.

##  I Supplemental material for “Topological phase transition between the gap
and the gapless superconductors”

##  II Spacial fluctuations of the magnetic impurities concentration

Let us perform its evaluation, for simplicity, from the “gap side” of the
phase transition. The first derivative in the “kinetic energy” term of Eq. (9)
can be easily obtained by direct differentiation of Eq. (10) in the main
paper:

$\left(\frac{d\Delta}{d\zeta}\right)=-\frac{\pi}{4}\Delta.$ (16)

What concerns the derivative $d\zeta/dn_{s}$ its calculation is more delicate
since $\zeta=(\tau_{s}\Delta)^{-1}$. The scattering lifetime $\tau_{s}$ is
determined by the integral over the solid angle $\Omega$

$\frac{1}{{{\tau_{s}}}}=\left[N(0)\frac{S(S+1)}{(2S+1)^{2}}\int{\left|{{f_{+}}-{f_{-}}}\right|^{2}d\Omega}\right]{n_{s}}=A{n_{s}},$
(17)

where $n_{s}$ is the concentration of the magnetic impurities and $f_{+}$ and
$f_{-}$ are the scattering amplitudes of an electron with a total angular
momentum $S+1/2$ and $S-1/2$.

Based on Eq. (17) and the fact that the order parameter $\Delta$ depends on
$n_{s}$ we obtain

$\frac{{d\zeta}}{{d{n_{s}}}}=\frac{A}{{\Delta}}-\frac{A}{{{{\Delta}^{2}}}}\frac{{d\Delta}}{{d{n_{s}}}}=\frac{A}{{\Delta}}-\frac{{{An_{s}}}}{{{{\Delta}^{2}}}}\frac{{d\Delta}}{{d\zeta}}\frac{{d\zeta}}{{d{n_{s}}}}.$
(18)

Taking into account Eq. (16) one finds

$\frac{d\zeta}{dn_{s}}=\frac{1}{n_{s}}\frac{\zeta}{1-\frac{\pi\zeta}{4}}.$
(19)

Relating the Cooper pair mass to the coherence length as $\xi^{2}=1/(4ma)$ and
returning to Eq. (9) in the main paper one finds that the kinetic energy term
in the gap domain ($\zeta<1$) is expressed as

$\frac{(\nabla\Delta)^{2}}{4m}=\frac{\pi^{2}\xi^{2}\Delta^{2}}{16}\frac{\zeta^{2}}{(1-\frac{\pi\zeta}{4})^{2}}\left(\frac{\nabla
n_{s}}{n_{s}}\right)^{2}.$ (20)

##  III Low temperature behavior of the order parameter

When the temperature is slightly above $T=0$ the temperature dependence of
order parameter $\Delta$ is given by expressions

$\Delta\left(T\right)=\left\\{\begin{gathered}\Delta\left(0\right)-\sqrt{\frac{{2\pi}}{3}}\frac{{\Delta_{g}^{1/6}}}{{{{\left({\zeta\Delta\left(0\right)}\right)}^{2/3}}}}{T^{3/2}}{\left({1-\frac{\pi}{4}\zeta}\right)^{-1}}{e^{-\beta{\Delta_{g}}}},{\text{
}}\zeta<1\hfill\\\
\Delta\left(0\right)-\frac{{{2^{2/3}}}}{{\sqrt{3}}}\frac{{{T^{5/3}}}}{{{{\left({\Delta\left(0\right)}\right)}^{2/3}}}}{\left({1-\frac{\pi}{4}}\right)^{-1}}\Gamma\left({\frac{2}{3}}\right){\rm
Z}\left({\frac{5}{3}}\right)\left({1-{2^{-2/3}}}\right),{\text{
}}\zeta=1{\text{ }}\hfill\\\
\Delta\left(0\right)-\frac{{{\pi^{2}}}}{6}{\left({1-{\zeta^{-2}}}\right)^{1/2}}{\left({1-\frac{1}{2}\sqrt{1-{\zeta^{-2}}}-\frac{\zeta}{2}\arcsin{\zeta^{-1}}}\right)^{-1}}\frac{{{T^{2}}}}{{\Delta\left(0\right)}},{\text{
}}\zeta>1,{\text{ }}\hfill\\\ \end{gathered}\right.$ (21)

where we recall $\Delta_{g}$ is the energy gap, $\Delta\left(0\right)$ is the
value of the parameter at $T=0$, $\Gamma(x)$ is the gamma function and ${\rm
Z}(x)$ is the Riemann zeta function defined by means the capital letter ${\rm
Z}$ to avoid a confusion with the driving parameter $\zeta$.

##  IV Asymptotic expressions for the thermoelectric coefficient

In the vicinity of the zero temperature or for the large values of $\beta$ the
contribution to the Eq. (13) in the main paper gives the low order
frequencies. Such a restriction allows to obtain several useful asymptotics
for the thermoelectric coefficients from the gap and the gapless side of the
phase transition. The starting point is the approximated expression for the
parameter $u$ in the case of the small $\omega$.

###  IV.1 Gap state

For the gap regime, where $\zeta<1$ we have

$\frac{{\omega-\Delta_{g}}}{{\Delta}}=-\frac{3}{2}{\zeta^{-\frac{2}{3}}}{\left({1-{\zeta^{\frac{2}{3}}}}\right)^{\frac{1}{2}}}{\left({u-{u_{0}}}\right)^{2}},$
(22)

where

$\frac{\Delta_{g}}{{\Delta}}={\left({1-{\zeta^{\frac{2}{3}}}}\right)^{\frac{3}{2}}},$
(23)

and

${u_{0}}={\left({1-{\zeta^{\frac{2}{3}}}}\right)^{\frac{1}{2}}}.$ (24)

Substitution of Eqs. (23)-(24) into Eq. (22) yields an equation for $u$ with
the solution for $\Delta_{g}=0$, i.e. when $\zeta\to 1$

$\begin{gathered}u=\sqrt{1-{\zeta^{\frac{2}{3}}}}\pm\frac{1}{3}\frac{{\sqrt{6{\zeta^{\frac{2}{3}}}\left[{{{\left({1-{\zeta^{\frac{2}{3}}}}\right)}^{\frac{3}{2}}}-\frac{\omega}{{\Delta}}}\right]}}}{{{{\left({1-{\zeta^{\frac{2}{3}}}}\right)}^{\frac{1}{4}}}}}\hfill\\\
\approx\sqrt{2\left({1-{\zeta^{\frac{1}{3}}}}\right)}\pm\frac{1}{3}{\text{i}}\frac{{\sqrt{6\left[{\frac{\omega}{{\Delta}}-2\sqrt{2}{{\left({1-{\zeta^{\frac{1}{3}}}}\right)}^{\frac{3}{2}}}}\right]}}}{{{2^{\frac{1}{4}}}{{\left({1-{\zeta^{\frac{1}{3}}}}\right)}^{\frac{1}{4}}}}}\hfill\\\
\end{gathered}$ (25)

Based on Eqs. (25) for the parameter $u$ one can write the expression for
functions ${\Omega\left({\omega,\Delta,\zeta}\right)}$ and
${h\left({\omega,\Delta,\zeta}\right)}$ that are entered in Eq. (13) in the
main paper for the thermoelectric coefficient in the main text. Introducing a
new parameter $z={\sqrt{1-{\zeta^{\frac{1}{3}}}}}$ near the the phase
transition we have

$\frac{{\Omega\left({\omega,\Delta,\zeta}\right)}}{\Delta}=\sqrt{2{z^{2}}+\frac{{\sqrt{2}w}}{{3z}}-1+{\text{i}}\frac{2}{3}{2^{\frac{1}{4}}}\sqrt{6zw}}-{\text{i}},$
(26)

and

$h\left({\omega,\Delta,\zeta}\right)=\frac{1}{2}\left[{1+\frac{{2{z^{2}}+\frac{{\sqrt{2}w}}{{3z}}-1}}{{\sqrt{{{\left({2{z^{2}}+\frac{{\sqrt{2}w}}{{3z}}-1}\right)}^{2}}+\frac{{8\sqrt{2}}}{3}zw}}}}\right],$
(27)

where $w={\frac{\omega}{{\Delta}}-2\sqrt{2}{z^{3}}}$ and for the extraction of
the square root of a complex number in Eq. (26) the well-known formula is
applied

$\sqrt{a+{\text{i}}b}=\pm\sqrt{\frac{{\sqrt{{a^{2}}+{b^{2}}}+a}}{2}}\pm{\text{i}}\operatorname{sgn}b\sqrt{\frac{{\sqrt{{a^{2}}+{b^{2}}}-a}}{2}}.$
(28)

Using Eq. (26) and (27) one can expand in series for small $\omega$ the part
of the integrand in Eq. (13) in the main paper

$\begin{gathered}\frac{{h\left({\omega,\Delta,\zeta}\right)}}{{\operatorname{Im}\left\\{{\Omega\left({\omega,\Delta,\zeta}\right)+\frac{{\text{i}}}{{2{\tau_{tr}}}}+\frac{{\text{i}}}{{{\tau_{s}}}}\left({1-h\left({\omega,\Delta,\zeta}\right)}\right)}\right\\}}}\approx\hfill\\\
{\Upsilon_{0}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)+{\Upsilon_{1}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)\frac{\omega}{{\Delta}}\hfill\\\
\end{gathered}$ (29)

where $\Upsilon_{0}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)$ and
$\Upsilon_{1}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)$ are some function
that we do not present explicitly due to their very cumbersome expressions.
However, one can also perform the expansion in series of this function for
$z=0$ or $(\zeta=1)$ to simplify further analytical calculations

$\Upsilon_{1}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)\approx\frac{{4\sqrt{2}}}{3}\frac{{{\tau_{s}}{\tau_{tr}}}}{{{\tau_{s}}+2{\tau_{tr}}}}z.$
(30)

Therefore, combining Eqs. (29) and (30) finally we obtain asymptotic
expression for the thermoelectric coefficient close to the phase transition
from the gap side

$\begin{gathered}\alpha=\frac{{8\sqrt{2}}}{3}\frac{{eN\left(0\right)Tv_{F}^{2}}}{{\Delta}}\frac{{{\tau_{s}}{\tau_{tr}}}}{{{\tau_{s}}+2{\tau_{tr}}}}z\int\limits_{-\infty}^{+\infty}{\frac{{{\omega^{2}}d\omega}}{{{{\cosh}^{2}}\left({\frac{\omega}{{2T}}}\right)}}}=\hfill\\\
\qquad\frac{{4\sqrt{2}{\pi^{2}}}}{9}\frac{{eN\left(0\right)Tv_{F}^{2}}}{{\Delta}}\frac{{{\tau_{s}}{\tau_{tr}}}}{{{\tau_{s}}+2{\tau_{tr}}}}z=\hfill\\\
\qquad\frac{{4\sqrt{2}{\pi^{2}}}}{9}\frac{{eN\left(0\right)Tv_{F}^{2}}}{{\Delta}}\frac{{{\tau_{s}}{\tau_{tr}}}}{{{\tau_{s}}+2{\tau_{tr}}}}\sqrt{1-{\zeta^{\frac{1}{3}}}}\hfill\\\
\end{gathered}$ (31)

###  IV.2 Gapless state

In the case of the gapless regime, where $\zeta>1$ the expansion of $u$ is
given by

$u={\text{i}}\sqrt{{\zeta^{2}}-1}+{\zeta^{2}}{\left({{\zeta^{2}}-1}\right)^{-1}}\frac{\omega}{{\Delta}}+...$
(32)

This allows to obtain in the same way the expression for functions
${\Omega\left({\omega,\Delta,\zeta}\right)}$ and
${h\left({\omega,\Delta,\zeta}\right)}$ that are entered in Eq. (13) in the
main paper for the thermoelectric coefficient in the main text. As in the
previous case introducing a new parameter $z={\sqrt{{\zeta-1}}}$ near the the
phase transition one can write

$\frac{{\Omega\left({\omega,\Delta,\zeta}\right)}}{\Delta}=\sqrt{2z+\frac{1}{4}\frac{1}{{{z^{2}}}}\frac{{{\omega^{2}}}}{{{{\Delta}^{2}}}}-1+{\text{i}}\frac{{\sqrt{2}}}{{\sqrt{z}}}\frac{\omega}{{\Delta}}}-{\text{i}},$
(33)

and

$h\left({\omega,\Delta,\zeta}\right)=\frac{1}{2}\left[{1+\frac{{2z+\frac{1}{{4{z^{2}}}}\frac{{{\omega^{2}}}}{{{{\Delta}^{2}}}}-1}}{{\sqrt{\sqrt{{{\left({2z+\frac{1}{4}\frac{1}{{{z^{2}}}}\frac{{{\omega^{2}}}}{{{{\Delta}^{2}}}}-1}\right)}^{2}}+\frac{2}{z}\frac{{{\omega^{2}}}}{{{{\Delta}^{2}}}}}}}}}\right],$
(34)

Based on Eq. (33) and (34) we expand in series for small $\omega$ the part of
the integrand in Eq. (13) in the main paper

$\begin{gathered}\frac{{h\left({\omega,\Delta,\zeta}\right)}}{{\operatorname{Im}\left\\{{\Omega\left({\omega,\Delta,\zeta}\right)+\frac{{\text{i}}}{{2{\tau_{tr}}}}+\frac{{\text{i}}}{{{\tau_{s}}}}\left({1-h\left({\omega,\left|\psi\right|,\zeta}\right)}\right)}\right\\}}}\approx\hfill\\\
{\Theta_{0}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)+{\Theta_{1}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)\frac{{{\omega^{2}}}}{{{{\Delta}^{2}}}}.\hfill\\\
\end{gathered}$ (35)

Due to long expressions for functions
${\Theta_{0}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)$ and
${\Theta_{1}}\left({\zeta,{\tau_{tr}},{\tau_{s}}}\right)$ we do not provide
them in an explicit form. Nevertheless, since the expansion in series given by
Eq. (35) contains only the even degree of $\omega$ it is easy to understand
that the integrand in Eq. (11) in the main paper is the odd function of
$\omega$ and, hence, the integral is equal to zero.


Arnaldo, Cuesta-Infante, Colmenar, Risco-Martín, AyalaBoosting the 3D Thermal-
Aware Floorplaning Problem through a M-W PMOEA

C/. Capitán, 39, 28300, Aranjuez (Madrid), Spain. E-mail:
<EMAIL_ADDRESS>
Ignacio Arnaldo is supported by Spanish Government Avanza Competitividad
I+D+I: TSI-020100-2010-962 project. The work has also been supported by
Spanish Government grants TIN 2008-00508 and MEC CONSOLIDER CSD00C-07-20811

# Boosting the 3D Thermal-Aware Floorplanning Problem through a Master-Worker
Parallel MOEA

Ignacio Arnaldo  1  Alfredo Cuesta-Infante  2  J. Manuel Colmenar  2  José L.
Risco-Martín  1  José L. Ayala  1  11affiliationmark:  DACYA, Complutense
University of Madrid22affiliationmark:  CES Felipe II, Aranjuez, Complutense
University of Madrid

###### Abstract

The increasing transistor scale integration poses, among others, the thermal-
aware floorplanning problem; consisting of how to place the hardware
components in order to reduce overheating by dissipation. Due to the huge
amount of feasible floorplans, most of the solutions found in the literature
include an evolutionary algorithm for, either partially or completely,
carrying out the task of floorplanning. Evolutionary algorithms usually have a
bottleneck in the fitness evaluation. In the problem of thermal-aware
floorplanning, the layout evaluation by the thermal model takes 99.5% of the
computational time for the best floorplanning algorithm proposed so far. The
contribution of this paper is to present a parallelization of this evaluation
phase in a master$-$worker model to achieve a dramatic speed-up of the
thermal-aware floorplanning process. Exhaustive experimentation was done over
three dimensional integrated circuits, with 48 and 128 cores, outperforming
previous published works.

###### keywords:

Thermal-aware floorplanning, optimization, multi-objective evolutionary
algorithm, parallelization

##  1 Introduction

Consumers continuously demand faster applications, smaller devices and
recently also ubiquitous computing. So far, developments in materials and
technology have allowed processor manufacturers to provide chips that attained
the expected serial performance. As we approach to the limits of
miniaturization, these demands become harder to accomplish. In order to remain
competitive, industry has moved to parallel architectures such as integrating
more cores in a die, data-parallel execution units, additional register sets
for hardware threads, bigger caches and more independent memory controllers to
increase memory bandwidth. For instance, multi-core general purpose computers
are being shipped for years and data-centers implement heterogeneous many-core
systems. Multi-processor systems-on-chip (MPSoCs) are nowadays also considered
as many-core systems. Up to now, the top core integration silicon CPU chip is
proposed by Intel Labs with an experimental Single-chip Cloud Computer (SCC),
a research microprocessor containing 48 Intel Architecture cores [1]. Also,
novel 3D multi-processor chips have been recently presented [2] as the
alternative to provide the required area of integration and to reduce the
communication delay among the large number of cores.

While the fabrication techniques have driven the integration of an increased
number of transistors to provide the required throughput, these improvements
have posed major problems regarding the operating temperature, directly
related to the power density [3]. As temperature increases, the carrier
mobility degrades, the leakage power consumption increases, gradient
temperatures appear on the surface creating electro-migrations and the
lifetime of the chip decreases exponentially, all in all reducing dramatically
the reliability of the system [4]. In addition, the specific placement of the
functional units on the chip surface (floorplan) also affects to the
temperature distribution because of the diffusive nature of heat [5]. Besides,
in the 3D configuration, the power density increases with the number of
layers. This effect is even more negative due to the problematic cooling of
inner layers of the 3D stack. The impact of power density in the
microprocessor is augmented due to the dielectric insulation layers inserted
between active layers. The reason is that the thermal conductivity of the
formers is very low compared to silicon and metal.

Increasing the chip area to reduce the power density has two shortcomings: it
is costly and requires to solve all the geometric constraints. Static external
coolers reduce the temperature of the chip surface by a constant factor but do
not reduce the temperature gradient across the chip. Instead, thermal-aware
floorplanning algorithms attempt to place functional units in order to achieve
a satisfactory temperature distribution; thus they tackle with both heat
dissipation and component placement at a time.

Floorplanning proposals are frequently formulated as combinatorial
optimization problems that can be smoothly fit to genetic algorithms (GA).

Broadly speaking, a GA performs a heuristic search throughout the solution
space inspired on Darwin’s principle of Natural Selection: The basic features
are:
(i) each iteration tests a small number of solutions (known as population)
compared to the cardinality of the solution space,
(ii) each solution (referred to as individual) is represented in a way
suitable both for evaluation and for producing the subset of the next
generation,
(iii) the next generation is obtained applying genetic operators such as
crossover, mutation and selection, and
(iv) there is a fitness function that evaluates the individuals.

Early floorplanning solutions tackled with GA proposed representations such as
Polish notation [6], combined bucket array [7] and O-trees [8] that are not
satisfactory in the thermal-aware problem because they were engineered to
minimize area. On the contrary, the hottest elements should be spread and
placed as far as possible for reducing the heating produced by closer hot
units. Thus, in 2D, works like [9] decreased the peak temperature using
genetic algorithms, and [10] using simulated annealing; whereas on 3D stacked
systems linear programming and simulated annealing combinations may be found
[11].

These works solve a single-objective optimization that takes into account only
the impact that temperature has on reliability. Therefore, they cannot provide
a thermally optimal solution with a minimum impact on the area of the chip and
the delay due to wiring. A more comprehensive approach in 3D systems is
presented in [12] and [13]. Their proposal consists of a Multi-Objective
Evolutionary Algorithm (MOEA) that tackles with the thermal-aware problem
(optimal placement of blocks) and also with the performance of the system
(minimum wire length delay) satisfying the topological constraints. On the
other hand, they show a critical bottleneck in the evaluation phase due to the
complexity of the computation, and can lead to very long execution times when
complex 3D many-cores architectures are considered.

Table 1: Execution time of the different methods involved in the algorithm
proposed in [12].

Method/Operator | Time (ms)
---|---
Evaluation | 31264762
Selection | 257
Reduction | 251
Mutation | 57.6
Crossover | 1

Our contribution in this paper is the parallelization of the thermal-aware
floorplanner proposed in [12] and [13], with the aim of reducing the
optimization execution time.

Evolutionary algorithms (EA) are intrinsically parallel but it is in the
fitness evaluation where more speed-up can be gained. Table 1 shows the
complete execution time of the sequential version of the algorithm proposed in
[12] until a solution is reached. The evaluation and reduction phases
correspond to methods devoted to compute the fitness and ranking of the
individuals, the rest implement the genetic operators of selection, crossover
and mutation. In this case, fitness is computed using a thermal model that
takes $83.1\%$ of the execution time. Due to the simple representation chosen
for the candidate solutions, a decodification phase is required before the
fitness computation. Adding the decodification phase time and the feasibility
verification time, the evaluation of individuals takes $99.5\%$ of the total
execution time. It is then clear that this task is by far the most time
consuming which justifies the necessity and effort of parallelization.

An EA is usually parallelized at two different levels: definition of
population or fitness evaluation [14]. At the former, the population is split
in a number of non-overlapping subpopulations that evolve independently but
with a probability of interaction. The two most popular models are Islands and
Grid models [14]. In Islands, some individuals are allowed to migrate with a
given frequency. There is a rich variety of Islands topologies, being the most
frequent rings, $n$-dimensional meshes and stars. When migrating, the worst
$k$ individuals in the destination are replaced by the new-comers, which are
the best $k$ individuals in their original island. In grids, each individual
is placed in a cell of a one- or two-dimensional grid. Genetic operations take
place in a small neighborhood of a given individual and their implementation
is straightforward on clusters.

The fitness evaluation parallelization is a much simpler and intuitive
approach. All the genetic operations are performed sequentially over the whole
population but, once a new generation is obtained, the individuals are
evaluated in parallel.

Regarding parallel MOEAs; they were early analyzed in [15]. Shortly
afterwards, the master-slave paradigm was employed in [16].

Since our baseline sequential algorithm presents a high computational load in
the evaluation of individuals, this paper proposes a master-worker algorithm
to parallelize that phase. Our approach was tested with a set of master-worker
configurations, ranging from 1 to 9 workers, as well as the sequential
algorithm, in two experimental multi-core platforms with 48 and 128 3D stacked
core processors each. Speedup review, validation of the proposal, study of
convergence and thermal analysis are presented in this paper. Results suggest
that a new representation could improve future algorithms.

The rest of the paper is organized as follows. The parallel Multi-Objective
Evolutionary Algorithm proposal is presented in Section 2. Experimental
results are shown and discussed in Section 3. Finally, conclusions and future
work are detailed in Section 4.

##  2 Parallel Multi-Objective Evolutionary Algorithm

In this section we present a parallel thermal-aware floorplanner capable of
optimizing many-core heterogeneous platforms under a master-worker schema.

###  2.1 Details of the MOEA

In [12], a Multi-Objective Evolutionary Algorithm (MOEA) is proposed for the
floorplanning of 3D stacked multi-processor single-chips. This kind of chips
consists of a number of layers of fixed area where the functional units
(processors, memories and interconnection blocks) must be placed. Our approach
considers this scenario performing a thermal-aware placement of the different
components while the wiring delay is also minimized. Moreover, it overcomes
[12] with a parallel implementation of the floorplanner that avoids the
constraints imposed to the optimization time.

Block placement problem.  All the components that model the functional units
of the many-core system must be placed in the 3D stack, which imposes the
physical boundaries of maximum length $L$, width $W$ and height $H$. Each
component is represented by a block $B_{i}$ with length $l_{i}$, width
$w_{i}$, height $h_{i}$ and it is denoted by its left-bottom-back corner, with
coordinates $(x_{i},y_{i},z_{i})$, taken the left-bottom-back corner of the
chip as origin of coordinates; where $0\leq x_{i}\leq L-l_{i}$, $0\leq
y_{i}\leq W-w_{i}$, $0\leq z_{i}\leq H-h_{i}$. These blocks cannot overlap. A
schematic representation is shown in Figure 1.

Figure 1: Block representation.

Blocks are placed sequentially. Since each component incorporates its
coordinates, this method leads to a floorplan whose components are not
necessarilly adjacent; unlike the state of the art works [6], [7] or [8]. This
represents a great advantage because cores, which are the most likely to
increase the temperature, can be placed explicitly separated, thus reducing
their impact in the overall temperature of the 3D chip. In order to select the
best placement coordinate $r_{i}=(x_{i},y_{i},z_{i})$ for block $B_{i}$, given
those already placed $B_{j},j<i$, a dominance relationship is established.
Therefore, a set of objective functions that evaluate the fitness, as well as
a suitable representation and appropriate genetic operators, must be derived
for the MOEA approach. The solution is obtained using a Non-dominated Sorting
Genetic Algorithm (NSGA-II) implementation [17].

Fitness.  There are three objective functions. The first objective $J_{1}$ is
the number of topological constraints violated by $B_{i}$ with respect to the
already placed $B_{j}$. The second objective is the wire length, approximated
by the Manhattan distance between blocks with coordinates $r_{i}$ and $r_{j}$,
$J_{2}=|x_{i}-x_{j}|+|y_{i}-y_{j}|+|z_{i}-z_{j}|.$ Finally the thermal impact
is measured through the power consumed by the unitary cells of the chip. A
thermal model that considers the power density of such cells and their
neighbors is used as an approximation to the steady state of the more accurate
thermal model that includes non-linear and differential equations. We evaluate
the thermal response of a given individual with the following model:

$J_{3}=\sum_{i<j\in 1..n}(dp_{i}*dp_{j})/(d_{ij})$ (1)

where $dp$ is the density power of the block considered and $d_{ij}$ is the
euclidean distance between blocks. This model has been shown to be accurate
enough and close to the non-linear simulation [18].

Representation.  Each individual is a possible configuration of the system,
and each configuration is represented through a chromosome that stores the
order in which blocks are being placed. The chromosome is an array whose
components contain an identifier of the functional unit that is going to be
placed. For instance, if 3 cores $(C_{1},C_{2},C_{3})$ and 3 memories
$(L_{1},L_{2},L_{3})$ are to be placed, the search space will have cardinality
$6!$. A possible chromosome would be

$[C_{1},L_{2},C_{3},L_{3},C_{2},L_{1}]$

meaning that $C_{1}$ will be fist placed, then $L_{2}$ and so on.

This decodification of the chromosome requires considering the size and
boundaries of the previously placed blocks. Therefore, the more the number of
components, the more the decodification execution time. As we will describe in
the future work section of the paper, a different representation encoding the
location of the components could help reducing the evaluation phase.

Operators.

Figure 2: Cycle crossover and Swap mutation for permutations of six elements.

Selection was carried out by tournament, taking two random individuals and
selecting the best among them. Individuals are mated in order to produce
offspring. Crossover must take into account that all the components must
appear once and only once in the chromosome. The so called _cycle crossover_
assures that the resulting chromosomes are just permutations of the parents.
Mutation consists of swapping the content of two positions inside the
chromosome or in the rotation of a randomly chosen component. Both cycle
crossover and swap mutation are depicted in Figure 2.

###  2.2 Details of Parallelization

In this paper, we propose a parallel implementation of the MOEA described in
the previous section using the master-worker model. As we have previously
shown, the evaluation phase of the algorithm takes over 99% of the execution
time. This is due to both the fact that the thermal response of all the
individuals of the population has to be evaluated in every generation of the
process, and the decodification of each individual, previous to its
evaluation.

The master-worker model satisfies our needs because, even though the fitness
is based on a simplified thermal model, the computational cost of this
evaluation increases quadratically with the number of components. Therefore,
it is interesting to exploit the fact that evolutionary algorithms are
intrinsically parallel and carry out the evaluation of the population in a
concurrent manner. Figure 3 depicts the approach used in this work. The master
distributes the population among $n$ workers, splitting the computational load
in $n$ ways so it does not carry out any evaluation. Once workers have
finished their task, they send the outcome together with the received
population subset to the master. Although the algorithm stops and waits for
all workers to finish, it is clearly much faster than the sequential execution
as long as each subset was large enough for compensating communication times.

We propose a multi-threaded implementation where only the master executes the
main thread of the algorithm. Since only workers execute the evaluation of
different subsets of the population, it is expected to obtain a speed-up
similar to the number of cores in the processor that executes the algorithm.

Figure 3: Master-worker configuration.

##  3 Experimental results

The experimental work will analyze the speedup obtained with the parallel
version of the MOEA while making clear that the quality of the solutions
remains the same. We will also analyze the thermal optimization achieved by
the floorplanner.

In order to evaluate our floorplanner, we study two heterogeneous 3D
architectures where every stacked layer is based on the Niagara platform. They
differ, from each other, in the number of cores. The first architecture is
composed of 48 processor cores: 32 SPARC and 12 Power6 cores. Adding 72
memories and 6 crossbar for inter-processor communication they sum up a total
of 126 components. In the second architecture, 128 cores are included: 96
SPARC plus 32 Power6. In addition, 192 memories and 16 crossbars are
considered, therefore 336 components need to be placed in this scenario. The
floorplanner will place the processors, the local memories and the crossbars
in 4 and 9 layers respectively. Both architectures represent the current and
the nearly future state-of-the-art in 3D many-core integration.

###  3.1 Speedup Analysis

In the first analysis, we studied the speedup obtained with the parallel
version of our floorplanner. We aim to find the optimal number of workers
leading to the maximum speedup. To this end, we perform a parametric sweep of
the number of workers, from 1 to 9, both in the 48 and 128 core scenarios. In
order to obtain the execution time of these optimizations, we run five times
each one of the worker configurations, obtaining the average execution time
and speedup for both scenarios. The experiments were carried out in a
dedicated Intel Core-i5 machine, a 4-core processor, running at 2.80GHz.

We set a fixed population size of 100 individuals and 250 generations
evolution as the MOEA parameters for the 48 core scenario optimization. Table
2 shows the average execution times and corresponding speedups for these runs,
with a number of workers ranging from 1 to 9.

Table 2: Average execution times and speedups obtained in the 48 cores
scenario.

# workers | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---
time (s) | 24171 | 12398 | 8904 | 6918 | 6380 | 6421 | 6665 | 6502 | 6386
speedup | 1 | 1.95 | 2.72 | 3.49 | 3.79 | 3.76 | 3.63 | 3.72 | 3.79

Figure 4 shows the obtained speedups in the 48 core scenario. It shows see
that the speedup increases almost linearly until the number of workers reaches
the number of cores of the processor we used for these optimizations (4-core
processor).

In our master-worker scheme, the execution time of each worker depends on the
particular evaluation time of the set of individuals that was assigned to the
worker. Then, if the worker receives a set of individuals that need more time
to be evaluated, the worker will slow down. On the contrary, if the
individuals need less time to be evaluated, the worker will speed up.

Therefore, in configurations from 2 to 4 workers, the optimization follows
this behavior: the population is divided into as many sets as number of
workers, then the set of individuals is sent to a different worker and the
evaluation begins. Once workers finish their evaluation task, they wait until
the slowest worker ends up, because the master synchronizes all the workers
until the next generation. As a result, the slowest worker establishes the
execution time of the evaluation of each generation, and the processor cores
that run faster workers will be idle waiting for synchronization.

Each worker runs in a different thread, and the load assigned to one thread
cannot be divided into different processor cores. As a consequence, if the
number of workers is higher than 4, the operating system scheduler will
distribute the execution of the worker threads among the 4 cores. Then, the
cores will swap between the threads, therefore advancing on the execution of
each one. As can be seen in Figure 4, results for 5 workers and above present
an asymptotic trend on speedup because the usage of resources is maximized.
The particular case of the 5 workers configuration obtains the maximum speedup
because maximizes the resources occupation with the lowest number of threads.

These results confirm that the parallelization of the evaluation phase in the
master-worker scheme contributes to the best speedup gains. In addition, there
is no remarkable penalty due to the parallelization, because the speedup
values above 4 cores tend to be similar.

Figure 4: Average speedup values obtained in the 48 cores scenario.

In order to strengthen this hypothesis, the same tests were run for the 128
cores scenario, where the evaluation time for each individual is much longer.
Here, the number of generations of the MOEA has to be scaled up because the
number of components to be placed has been increased. Therefore we consider a
number of generations equal to the total number of components, which is 336:
128 cores, 192 memories and 16 crossbars. The population size remains 100
individuals. Table 3 shows the average execution time and speedup for this
scenario from 1 to 9 workers configuration. Figure 5 displays the speedup
trend for these data.

Table 3: Average execution times and speedups obtained in the 128 cores
scenario.

# workers | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---
time (s) | 634070 | 432160 | 264318 | 203749 | 198977 | 248527 | 235313 | 250021 | 221485
speedup | 1 | 1.47 | 2.14 | 3.11 | 3.19 | 2.55 | 2.69 | 2.54 | 2.86

Figure 5: Average speedup values obtained in the 128 cores scenario.

As was shown, the 128 cores scenario presents the same behavior as the 48
cores one. The resources of the CPU are maximized from the 5 workers
configuration on, and higher numbers of workers obtain similar speedup values.
However, the performance improvement is lower than in the 48 cores
configuration. This behavior occurs because the individual evaluation time is
much higher in this 128 cores scenario, and the execution time of the threads
does not differ so much. In the 48 cores case, the processor slots available
due to the different execution time between threads allow the evaluation of
more individuals than in the 128 cores configuration. As a consequence, the
workers queue, waiting for processor cores, advance more in their execution,
obtaining higher speedups. On the contrary, the workers of the 128 cores
scenario are not able to exploit the processor free slots to evaluate as many
individuals as in the 48 cores case, therefore obtaining lower speedup values.

###  3.2 Validation of solutions

Conceptually it is clear that the parallelization of the fitness evaluation
should lead to the same results than the sequential version of the algorithm.
Nevertheless, we have included in this section a brief consideration regarding
the validation of the parallelization that assures such a baseline. Thus, in
order to show that the quality of the solutions proposed by the parallel
version of the floorplanner remains the same than in the sequential version,
we compare the front of non-dominated solutions obtained with the sequential
version of the algorithm with the front obtained with the parallel version
using 4 and 5 workers in the 48 cores scenario. Figure 6 shows the fronts of
non-dominated solutions returned by the floorplanner in these cases.

Figure 6: Non-dominated fronts of solutions returned by the floorplanner in
the 48 cores scenario.

In Figure 7, we compare the front of non-dominated solutions proposed by the
floorplanner working with 4 and 5 nodes for the 128 cores platform.

Figure 7: Non-dominated fronts of solutions returned by the floorplanner in
the 128 cores scenario.

For every run of the algorithm, the returned non-dominated front is different.
In fact each execution explores a different region of the solution space. We
can see that the fronts cross each other in at least one point. Therefore,
none of the returned fronts dominates the others.

Since EA are intrinsically heuristic, two executions will not produce exactly
the same results. Hence, in order to prove that our proposal is valid it is
necessary to define a measure that analyzes the outputs (solution sets) both
from sequential and parallel executions. Such a measure is usually referred to
as _Indicator_ $I$. In this work, the _Hypervolume_ indicator, proposed by
Zitzler and Thiele [19], has been used. The hypervolume $I(A)$ measures the
total amount of the objective space that has been ‘covered’ by the solution
set $A$; returning the hypervolume of that portion of the objective space that
is weakly dominated by $A$. To this end, the objective space must be bounded.
Otherwise a reference point that must be at least weakly dominated by all
solutions in $A$ is used. Finally higher values of $I$ correspond to higher
quality of the measured set.

Figure 8: Hypervolumes for 48 cores (left) and 128 cores(right) after 30
optimization runs. Both hypervolumes measured in the sequential and in the
parallel execution, the latter with 4 and 5 workers. The central line is the
median, the edges of the box are the 25th and 75th percentiles, the whiskers
extend to the most extreme data points not considered outliers, which are
plotted individually (+ mark).

The comparison has been carried out between the sequential execution, the
4-workers and the 5-workers versions of the parallel implementation. This
choice was motivated because these configurations had obtained the highest
non-saturated speed-ups. Results are shown in Figure 8. All were obtained
after running 30 optimizations, each one with 250 generations in the 48 cores
and 366 generation in the 128 cores. As expected, the three boxplots inside
each picture show a similar outcome; with 25th and 75th percentiles almost
identical within the 48 and 128 core plots.

###  3.3 Convergence of the MOEA

The main benefit of parallelization is the considerable reduction in the
fitness evaluation time and, consequently, in the the whole procedure of
finding a good floorplan. In addition, we can take advantage of such a speed
up for carrying out tests in order to detect possible weak points in the
algorithm that, using the sequential version, would take months to complete.

As the population in an EA evolves, it is desirable to keep their diversity.
Otherwise the exploration of the solution space will be guided towards a
region, avoiding others which might be more promising. The analysis of the
convergence is a straightforward method for verifying that diversity is
maintained. At the same time, it reveals whether the EA is well engineered or
there is room for improvement. A slow convergence with good results is usually
due to a poor representation or not appropriate genetic operators.

Although this work tackles a three objective problem, convergence is studied
only for $J_{2}$ (wire length) and $J_{3}$ (thermal response) because
$J_{1}>0$ means that none floorplan satisfies the constraints.

Thus, values $J_{2}$ and $J_{3}$ of feasible solutions are extracted in arrays
$W_{r,g}$ and $T_{r,g}$ respectively, one for each optimization run $r=1\ldots
30$ and each generation $g=1\ldots 250$ for 48 cores and $g=1\ldots 366$ for
128. Then six matrices are constructed in the following way:

$\begin{array}[]{ccc}{{\textbf{W}}_{{\textbf{min}}}}(r,g)=\min\left\\{{{W_{r,g}}}\right\\},&{{\textbf{W}}_{{\textbf{mean}}}}(r,g)=\text{mean}\left\\{{{W_{r,g}}},\right\\}&{{\textbf{W}}_{{\textbf{max}}}}(r,g)=\max\left\\{{{W_{r,g}}}\right\\},\\\
{{\textbf{T}}_{{\textbf{min}}}}(r,g)=\min\left\\{{{T_{r,g}}}\right\\},&{{\textbf{T}}_{{\textbf{mean}}}}(r,g)=\text{mean}\left\\{{{T_{r,g}}},\right\\}&{{\textbf{T}}_{{\textbf{max}}}}(r,g)=\max\left\\{{{T_{r,g}}}\right\\}.\\\
\end{array}$

This procedure is done for both 48 and 128 core configurations. Finally, all
six matrices are scaled between the minimum and maximum wire length and
thermal response respectively, and plotted as shown in Figure 9.

Figure 9: Convergence evolution of the minimum, mean and maximum values for
objectives $J_{2}$ (wire length) and $J_{3}$ (thermal response), considering
only feasible individuals; for 48 cores (above) and 128 cores (below).

The left-most pictures, corresponding to the minimum values in each generation
and optimization run of $J_{2}$ and $J_{3}$, show a decreasing behavior
eventually reaching the global minimum in almost all optimization runs for 48
cores and, at much slower rate, for 128 cores. The middle pictures show
convergence of the mean values of $J_{2}$ and $J_{3}$. Their trend is
decreasing, starting at $1/2$ of the upper bound down to $1/4$ in the best
cases and $1/3$ in the worst. Finally, the convergence of maximum values of
$J_{2}$ and $J_{3}$ is shown in the right-most pictures. Values decrease down
to $1/3$ of the upper bound in most of the optimization runs. To the light of
these results it is clear that generations are better fitted as evolution
advances. Moreover, since the mean of the last generation is quite close to
the maximum in both objectives, less fitted individuals still have a
considerable probability of being chosen, attesting that diversity is
maintained. On the other hand, the slow convergence of the mean, especially in
the larger problem of 128 cores, indicates that the representation could be
better engineered. Genetic operators were discarded as reason because the
offspring is always valid. Thus, no reconstruction is needed, which might lead
to repeat certain schemas.

###  3.4 Thermal analysis

Finally, we analyze the thermal optimization obtained with our algorithm. The
floorplanner works with a fixed die size and aims to minimize both the total
wire length and the maximum temperature of the chip. As we want to perform a
thermal optimization of the described platforms, we need to provide the power
consumption values and the areas of the different elements of the
architectures as inputs to the thermal-aware floorplanner. In [20] we find
that the power consumption of the sparc is 4W at 1.4GHz. In the case of the
Power6, we find that 2.6W is the estimated power dissipation [21]. We consider
the following areas: $3.24mm^{2}$ and $1.5mm^{2}$ for the sparc and Power6
respectively (see [20] and [21]). The power consumption values and areas of
the memories are found with the CACTI software [22].

48-cores configuration.  We compare an optimized configuration of the 48-cores
heterogeneous platform to the 48-cores homogeneous platform represented in
Figure 10. In this baseline configuration, an original architecture composed
of 12 cores is replicated in all the layers. As a consequence, the SPARC cores
(SPC) are placed above the others producing hotspots. On the other hand,
Figure 11 shows the thermal maps of the different layers of a non-dominated
solution returned by the thermal-aware floorplanner. This figure shows an
optimized placement of the SPARC cores (SPC), Power6 cores (P6), memories (L2)
and crossbars (Cross) achieved by the floorplanner. In this configuration the
hottest elements (SPARC cores) are generally placed in the borders of the chip
and in the outer layers, separated as much as possible. In fact, the
floorplanner avoids placing cores above the others as vertical heat spread is
also taken into account. The crossbars are placed in intermediate layers to
minimize the wire length.

Figure 10: Thermal map of the 4 layers of the baseline configuration of the 48
cores platform.

Figure 11: Thermal map of the 4 layers of a non-dominated solution of the 48
cores platform.

The metrics considered for the thermal analysis of these two platforms are the
maximum and mean temperature of the chip and the maximum thermal gradient. In
Table 4 we present the thermal response of these two different configurations.
These results show that our floorplanner proposes thermally optimized
configurations. The peak temperature of 411.82K found in the original
configuration is reduced to 345.30K while the mean temperature is reduced in
12.54K. We can see that the maximum thermal gradient of the optimized
configuration is reduced from 109.75K to 31.81K. Therefore, not only the
temperature of the chip is reduced but it is also more evenly distributed. On
the other hand, the wire length of the optimized configuration is a 2.11%
greater than the original which translates into a small performance penalty.

Table 4: Thermal response of the 48 cores configurations.

| Wire L. | $T_{MAX}$ | $T_{MEAN}$ | $Grad_{MAX}$
---|---|---|---|---
48 baseline | 6733 | 411.82 K | 344.29 K | 109.75
48 opt. | 6875 | 345.30 K | 331.75 K | 31.81

128-cores configuration.  For this larger configuration, we analyze one of the
optimal floorplans obtained with our parallel implementation. Figure 12 shows
the thermal map of the chosen solution. As for the 48-cores platform, we can
see that the SPARC cores tend to be placed in the outer layers and in the
borders of the chip. The memories and the crossbars are placed in the inner
layers. This way both the chip temperature and the wire length are minimized.
Nevertheless hotspots appear in this configuration. Table 5 shows the thermal
response of an optimized configuration of the 128 cores platform. The hotspot
visible in the first layer of the chip corresponds to the peak temperature of
the chip reaching 396.84K. The mean temperature is 362.50K while the maximum
thermal gradient is 75.80K. Further research and simulations with cooling
techniques are required to study the feasibility of these architectures.

Figure 12: Thermal map of the 9 layers of a non-dominated solution of the 128
cores platform. Table 5: Thermal response of the 128 cores configuration.

| Wire L. | $T_{MAX}$ | $T_{MEAN}$ | $Grad_{MAX}$
---|---|---|---|---
128 opt. | 31587 | 396.84 K | 362.50 K | 75.80

##  4 Conclusions and Future Work

Current and short term future state-of-the-art in 3D many-core integration
requires thermal-aware floorplanning techniques able to reduce peak and mean
temperatures. However, current techniques that take into account thermal
issues spend the most of the execution time dealing with decoding and
evaluation of solutions.

This work has proposed a parallel implementation of a thermal-aware Multi-
objective Evolutionary Algorithm for 3D floorplanning using a master-worker
scheme. This model has provided optimized configurations for systems composed
of 48 and 128 heterogeneous processor cores.

We have shown that the highest speedup values are obtained when the number of
workers is closer to the number of cores of the processor that runs the
algorithm. In our experiments, run on a 4-core processor, we have obtained
maximum speedup values of 3.79 and 3.19 respectively for the 48 and 128 core
test configurations by selecting 5 workers in the optimization algorithm.
Furthermore, the parallelization presented in this work has made possible the
study of the convergence of the floorplanner. The performed analysis shows
that convergence evolves successfully in our experiments.

As a future work, we aim to overcome the drawbacks of the floorplanner
presented in this work. First, we plan to replace the current approximated
thermal model with a validated thermal simulator. In fact, the implemented
thermal model is motivated by its low computational cost but it might not be
obtaining accurately the thermal behavior of the different individuals. A
recent work [23] has shown that the integration of a thermal simulator in the
floorplanning process leads to a simultaneous reduction of the peak
temperature and the wire length. Furthermore, the model proposed in the
referred work is claimed to be, not only more accurate, but also faster than
the approximated thermal model.

Another major improvement could be achieved with the use of a more suitable
representation of the solutions. In fact, most of the floorplanning proposals
are based on representations that require time consuming heuristics to decode
the solutions. For instance, in our work, the decoding of the solutions
together with the evaluation remains the bottleneck of the optimization
process. A new representation allowing a direct mapping of the individuals
into configurations of the architecture would alleviate the computational cost
of the algorithm as the decoding step would be avoided. Furthermore, it would
eliminate heuristics that might limit the exploration space and cause
premature convergence problems. Thus, such a representation would be more
suitable for fixed-outline floorplanning problems.

To propose a tool consistent with the state of the art of 3D chip design,
thermal-aware floorplanners have to be implemented in accordance with current
thermal simulators that split the IC into thermal cells (as done by 3D-ICE
[24]). This way, the coding of the solutions has to respect the grid-like
structure used by this kind of simulators. Therefore, the thermal error due to
the different cell sizes used in the optimization and validation processes is
eliminated. Thus, a better thermal optimization can be achieved. The use of a
grid-like representation together with the removal of placement heuristics
makes the process well adapted for execution in massively parallel
architectures such as GPUs.

A deeper study of new representations is our current and future work, and the
preliminary results are very promising.

## References

  * [1] Intel. http://techresearch.intel.com/ProjectDetails.aspx?Id=1 2012.
  * [2] Loh GH, Xie Y. 3D stacked microprocessor: Are we there yet? _Micro, IEEE_ may-june 2010; 30(3):60 –64, 10.1109/MM.2010.45.
  * [3] Borkar S. Design challenges of technology scaling. _Micro, IEEE_ 1999; 19(4), 10.1109/40.782564.
  * [4] Srinivasan J, Adve S, Bose P, Rivers J. The impact of technology scaling on lifetime reliability. _Dependable Systems and Networks, 2004 International Conference on_ , 2004; 177 – 186, 10.1109/DSN.2004.1311888.
  * [5] Sankaranarayanan K, _et al._. A case for thermal-aware floorplanning at the microarchitectural level. _Journal of Instruction-Level Parallelism_ 2005; 7(1):8–16.
  * [6] Berntsson J, Tang M. A slicing structure representation for the multi-layer floorplan layout problem. _EvoWorkshops_ , 2004; 188–197.
  * [7] Cong J, Wei J, Zhang Y. A thermal-driven floorplanning algorithm for 3D ICs. _Proceedings of the 2004 IEEE/ACM International conference on Computer-aided design_ , ICCAD ’04, IEEE Computer Society: Washington, DC, USA, 2004; 306–313, http://dx.doi.org/10.1109/ICCAD.2004.1382591.
  * [8] Tang M, _et al._. A memetic algorithm for VLSI floorplanning. _IEEE Transactions on Systems, Man, and Cybernetics, Part B_ 2007; 37(1):62–69.
  * [9] Hung WL, Xie Y, Vijaykrishnan N, Addo-Quaye C, Theocharides T, Irwin MJ. Thermal-aware floorplanning using genetic algorithms. _Proceedings of the 6th International Symposium on Quality of Electronic Design_ , ISQED ’05, IEEE Computer Society: Washington, DC, USA, 2005; 634–639, http://dx.doi.org/10.1109/ISQED.2005.122.
  * [10] Han Y, Koren I. Simulated annealing based temperature aware floorplanning. _J. Low Power Electronics_ 2007; 3(2):141–155, http://dx.doi.org/10.1166/jolpe.2007.128.
  * [11] Healy M, _et al._. Multiobjective microarchitectural floorplanning for 2D and 3D ICs. _CADICS, IEEE Transactions on_ 2007; 26(1):38–52, 10.1109/TCAD.2006.883925.
  * [12] Cuesta Gómez D, Risco Martín JL, Ayala JL, Hidalgo JI. A combination of evolutionary algorithm and mathematical programming for the 3D thermal-aware floorplanning problem. _Proceedings of the 13th annual conference on Genetic and evolutionary computation_ , GECCO ’11, ACM: New York, NY, USA, 2011; 1731–1738, http://doi.acm.org/10.1145/2001576.2001809.
  * [13] Arnaldo I, Risco-Martín JL, Ayala JL, Hidalgo JI. Power profiling-guided floorplanner for thermal optimization in 3D multiprocessor architectures. _Proceedings of the 21st international conference on Integrated circuit and system design: power and timing modeling, optimization, and simulation_ , PATMOS’11, Springer-Verlag: Berlin, Heidelberg, 2011; 11–21.
  * [14] Cantú-Paz E. A survey of parallel genetic algorithms. _Calculateurs Paralleles Reseaux et Systems Repartis_ 1998; 10:3289–3293.
  * [15] Van Veldhuizen D, Zydallis J, Lamont G. Considerations in engineering parallel multiobjective evolutionary algorithms. _Evolutionary Computation, IEEE Transactions on_ april 2003; 7(2):144 – 173, 10.1109/TEVC.2003.810751.
  * [16] Lopez Jaimes A, Coello Coello C. MRMOGA: parallel evolutionary multiobjective optimization using multiple resolutions. _Evolutionary Computation, 2005\. The 2005 IEEE Congress on_ , vol. 3, 2005; 2294 – 2301 Vol. 3, 10.1109/CEC.2005.1554980.
  * [17] Deb K, _et al._. A fast and elitist multiobjective genetic algorithm: NSGA-II. _IEEE Transactions on Evolutionary Computation_ 2002; 6(2):182–197.
  * [18] Paci G, _et al._. Exploring temperature-aware design in low-power mpsocs. _International journal of embedded systems_ 2007; 3(1):43–51.
  * [19] Zitzler E, Thiele L. Multiobjective Evolutionary Algorithms: A Comparative Case Study and the Strength Pareto Approach. _IEEE Transactions on Evolutionary Computation_ 1999; 3(4):257–271.
  * [20] OpenSPARC. http://www.opensparc.net/pubs/preszo/07/n2isscc.pdf 2007.
  * [21] IBM. http://www.ibm.com/systems/support/tools/estimator/energy 2012.
  * [22] HPlabs. www.hpl.hp.com/research/cacti/ 2012.
  * [23] Arnaldo I, Vincenzi A, Ayala JL, Risco JL, Hidalgo JI, Ruggiero M, Atienza D. Fast and Scalable Temperature-driven Floorplan Design in 3D MPSoCs. _13th Latin American Test Workshop (LATW)_ , 2012.
  * [24] Sridhar A, Vincenzi A, Ruggiero M, Brunschwiler T, Atienza D. 3D-ICE: Fast compact transient thermal modeling for 3D ICs with inter-tier liquid cooling. _Computer-Aided Design (ICCAD), 2010 IEEE/ACM International Conference on_ , 2010; 463 –470, 10.1109/ICCAD.2010.5653749.