# Superradiant Superconductivity

G. Baskaran  The Institute of Mathematical Sciences, C.I.T. Campus, Chennai
600 113, India
Perimeter Institute for Theoretical Physics, Waterloo, Ontario, Canada N2L 2Y5

###### Abstract

We suggest possibility of Dicke superradiance in superconductors. The
necessary 2-level atoms are identified with Anderson pseudo spins in k-space,
seeing a k-dependent self consistent mean field. A way to couple these 2-level
bose atoms to a macroscopically excited coherent boson mode and create a novel
nonequilibrium superradiant superconductivity (SRSC) is suggested. Our
coherence transfer mechanism offers a hope to realize transient
superconductivity, even at room temperatures, in the pseudo gap phase of
certain underdoped cuprates. Recent experiments are briefly discussed in the
light of our theory. Quantum entanglement, QCP and superfluorescence
properties follow.

Introduction Superconductivity is a remarkable macroscopic manifestation of
quantum mechanics. A rich physics and phenomenology, including Meissner and
Josephson effects are parts of superconductivity BCS . Dicke’s Superradiance
Dicke1954  is another macroscopic manifestation, exhibited by a collection of
2-level atoms interacting with a single boson mode. The coupled system can
develop quantum coherence, enhanced emission properties and complex dynamics.
Certain phenomena in NMR, ESR, optics and cold atoms are related to
superradiance.

In the present work we suggest a way to combine superconductivity and
superradiance, We call the resultant non equilibrium state as superradiant
superconductivity (SRSC). In our proposal a macroscopically occupied long
wavelength single boson mode interacts with a collection of independent
2-level atoms located in k-space and creates a Dicke superradiant situation,
under certain conditions. In this state certain deformation of Cooper pair
wave function is entangled with a coherent external bosonic mode.

Interaction of coherent electromagnetic radiation and ultrasound with
superconductor is a well studied subjectMicrowaveExpt ; Eliashberg ;
OwenScalapino ; KumarSinha ; McIntoshLindesay . Our proposal of SRSC may have
relevance to some known results. An exciting recent development is
experimental observation of transient superconductivity well above Tc, induced
by certain femtosecond laser pulses, in the pseudo gap phase of cuprates
liscLBCO ; liscYBCO .

In a pioneering theoretical work Eliashberg Eliashberg  in 1970 showed that
microwave induced quasi particle redistribution self consistently enhances gap
values and Jc. Works by Scalapino, Owen and Chang OwenScalapino , also focused
on quasi particle redistribution. In a later theory in 1994, McIntosh and
Lindesey McIntoshLindesay  showed that stimulated emission and reabsorption of
photon by correlated electron pairs play a fundamental role in
superconductivity enhancement. This key insight is one of the triggers for our
proposal. Interestingly, in 1968, there was a theoretical suggestion
KumarSinha  for photon induced room temperature superconductivity.

In what follows, we start with an ideal BCS superconductor and show how Dicke
superradiance emerges, when the wavelength of the macroscopically occupied
external single boson mode $\lambda\geq L$, the sample size $L$. Then we
discuss how our mechanism could gnerate transient superconductivity abouve Tc
and discuss recent experiments liscLBCO  and (see note liscYBCO  in the light
of our mechanism.

In our work we make the tacit assumption that there are suitable relaxation
processes involving quasiparticles and phonons that drains energy to the heat
bath efficiently to avoid heating. At the same time some energy gets pumped to
the electronic sub system to help reach a new non equilibrium coherent state
for a short time scale. It is the nature of non equilibrium coherent state
that we are after. To achieve this we assume that the coherent state of the
single boson mode is long lived and does not radiate away its energy. It
exchanges its quanta with the electron subsystem only and gets quantum
entangled. Ours is an equilibrium statistical mechanics approximation tailored
to get a glimpse of a remarkable non equilibrium situation.

Model. To develop our theory we follow Anderson’s pseudo spin formulation of
BCS theory PWApseudoSpin . It helps us to view BCS mean field eigen states as
a k-space lattice containing 2-level bose atoms and free fermions. Consider
time reversed single particle states $({\bf k}\uparrow,{\bf-k}\downarrow)$,
with empty state written as $|0\rangle_{\bf k}$. To generate complete Fock
space, we need only 4 states in each $({\bf k}\uparrow,{\bf-k}\downarrow)$ :
i) $|0\rangle_{\bf k}$, ii)
$c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow}|0\rangle_{\bf k}$, iii)
$c^{\dagger}_{k\uparrow}|0\rangle_{\bf k}$ and iv)
$c^{\dagger}_{-k\downarrow}|0\rangle_{\bf k}$. BCS interaction mixes only the
0 and 2-fermion states. Resulting ground and excited paired states are two
orthogonal states:
$|g\rangle_{k}\equiv(u_{k}+v_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow})|0\rangle_{k}$
and
$|e\rangle_{k}\equiv(u_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow}-v_{k})|0\rangle_{k}$.
We call these 2-level bosonic states as Anderson atom or A-atom. A-atom
carries zero total momentum. Single fermion states
$c^{\dagger}_{k\uparrow}|0\rangle_{\bf k}$ and
$c^{\dagger}_{-k\downarrow}|0\rangle_{\bf k}$, in $({\bf
k}\uparrow,{\bf-k}\downarrow)$ remain unaffected by BCS interaction.

An A-atom close to fermi suface is special (see note note2 ). It is a coherent
superposition of 0 and 2-electron states. Consequent non zero value of the
product $u_{k}v_{k}$, around the fermi surface quantifies superconductivity.

BCS mean field Hamiltonian has the familiar form:

$H_{mf}=\sum\varepsilon_{k}\alpha^{\dagger}_{k\sigma}\alpha_{k\sigma},$ (1)

where, Bogoliubov quasi particle operators $\alpha^{\dagger}_{k\sigma}\equiv
u_{k}c^{\dagger}_{k\sigma}+\sigma v_{k}c_{-k-\sigma}$ and
$\alpha_{k\sigma}\equiv u_{k}^{*}c_{k\sigma}+\sigma
v_{k}^{*}c^{\dagger}_{-k-\sigma}$. The quasi particle energy
$\varepsilon_{k}\equiv\sqrt{(\frac{\hbar
k^{2}}{2m}-\mu)^{2}+\triangle_{k}^{2}}$.

Complete set of BCS mean field eigen states can be written as product over all
states, $({\bf k}\uparrow,{\bf-k}\downarrow)$, each containing either an
A-atom in the ground or excited state or a single upspin or down spin fermion
state. Bogoliubov quasi particle operators have very simple action on the BCS
eigen states. BCS vacuum,
$|BCS\rangle=\prod_{k}(u_{k}+v_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow})|0\rangle$
is annihilated by the annihilation operator, $\alpha_{q\sigma}|BCS\rangle=0$.

Bogoliubov creation operator, while acting on the BCS ground state, removes an
A-atom and replaces it by a fermion:
$\alpha^{\dagger}_{q\uparrow}|BCS\rangle=c^{\dagger}_{q\uparrow}\prod_{k\neq
q}(u_{k}+v_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow})|0\rangle$ and
$\alpha^{\dagger}_{-q\downarrow}|BCS\rangle=c^{\dagger}_{-q\downarrow}\prod_{k\neq
q}(u_{k}+v_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow})|0\rangle$.

What is the operator that excites an A-atom ? Pair of Bogoliubov quasi
particle operators
$\alpha^{\dagger}_{q\uparrow}\alpha^{\dagger}_{-q\downarrow}$, with total
momentum zero and total spin projection zero, acting within $({\bf
q}\uparrow,{\bf-q}\downarrow)$ excites an A-atom:
$\alpha^{\dagger}_{q\uparrow}\alpha^{\dagger}_{-q\downarrow}|BCS\rangle=(u_{q}c^{\dagger}_{q\uparrow}c^{\dagger}_{-q\downarrow}-v_{q})\prod_{k\neq
q}(u_{k}+v_{k}c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow})|0\rangle$.

The 2-level (bosonic) A-atom subspace can be studied using pseudo spin (Pauli)
operators. Pseudo spin operators (see note note3 ) are defined as,
$\sigma^{z}_{k}\equiv(1-\alpha^{\dagger}_{k\uparrow}\alpha_{k\uparrow}-\alpha^{\dagger}_{-k\downarrow}\alpha_{-k\downarrow}),\leavevmode\nobreak\
\sigma^{+}_{k}\equiv\alpha^{\dagger}_{k\uparrow}\alpha^{\dagger}_{-k\downarrow}$
and $\sigma^{-}_{k}\equiv\alpha_{-k\downarrow}\alpha_{k\uparrow}$.

The BCS mean field Hamiltonian (equation 1) in the boson subspace takes a
suggestive form :

$H_{mf}=-\sum\varepsilon_{k}\sigma_{k}^{z},$ (2)

It describes a collection of non-interacting pseudo spins in the prsence of a
k-dependent magnetic field of magnitude $\varepsilon_{k}$. Energy level
separation of a 2-level A-atom is 2$\varepsilon_{k}$. Notice that long range
interaction in k-space in the BCS Hamiltonian leads to free spins in the mean
field description, but in the presence of a self consistent mean field of
magnitude $\varepsilon_{k}$ in k-space. In our pseudo spin basis BCS ground
state is a fully aligned ferromagnet, while in Andersons basis pseudo spins
twist to form Bloch wall across the fermi surface in k-space (see note note2
).

Now we consider a simple way to couple A-atoms selectively to a single
external boson mode, with creation and annihilation operators
($b^{\dagger},b$). Interaction of electrons with this mode, in the long wave
length (zero momentum transfer) limit, $\lambda>>L$, where $L$ is the size of
the sample, has a simple form:

$H_{int}=\frac{1}{{\sqrt{N}}}\sum
B_{k}(c^{\dagger}_{k\sigma}c_{k\sigma}+H.c.)(b+b^{\dagger})$ (3)

Here Bk is a momentum dependent coupling constant and N $\sim$ number of
electrons in the interaction region. In terms of Bogoliubov quasiparticle
operators,

$\displaystyle H_{int}=\frac{1}{{\sqrt{N}}}\sum
B_{k}(u_{k}^{2}-B_{-k}v_{k}^{2})\alpha^{\dagger}_{k\sigma}\alpha_{k\sigma}(b+b^{\dagger})+\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ (4) $\displaystyle+$
$\displaystyle\frac{1}{{\sqrt{N}}}\sum(B_{k}+B_{-k})u_{k}v_{k}(\alpha^{\dagger}_{k\uparrow}\alpha^{\dagger}_{-k\downarrow}+H.c.)(b+b^{\dagger})$

We ignore non resonant terms using rotating wave approximation. Further
quasiparticle number operators can be also be taken care of using Hartree type
of approximations. We are left with the important pair annihilation and
creation terms:

$H_{int}\approx\frac{1}{{\sqrt{N}}}\sum
B_{k}u_{k}v_{k}(\alpha^{\dagger}_{k\uparrow}\alpha^{\dagger}_{-k\downarrow}b+\alpha_{-k\downarrow}\alpha_{k\uparrow}b^{\dagger})$
(5)

Interms of pseudo spin operators it takes the form
$H_{int}\approx\frac{1}{{\sqrt{N}}}\sum
B_{k}u_{k}v_{k}(\sigma^{+}_{k}b+\sigma^{-}_{k}b^{\dagger})$. Thus the final
form of the Hamiltonian of the superconductor interacting with a single boson
mode is:

$H=\hbar\omega_{0}(b^{\dagger}b+\frac{1}{2})-\sum\varepsilon_{k}\sigma_{k}^{z}+\frac{1}{{\sqrt{N}}}\sum\lambda_{k}(\sigma_{k}^{+}b+\sigma_{k}^{-}b^{\dagger})$
(6)

where $\lambda_{k}\equiv(B_{k}+B_{-k})u_{k}v_{k}$. Equation 5 is a generalized
Dicke Hamiltonian Dicke1954 , where 2-level atoms in k-space have a
k-dependent energy level separation, The sum, N${}_{t}\equiv$ N∗ \+ Nboson of
number of excited N∗ atoms and number of photons Nboson, commutes with the
Hamiltonian (equation 6)2$\varepsilon_{k}$.

Finding a Dicke like Hamiltonian is a key result of our paper, from which
several consequences follow.

Notice that A-atom-boson mode coupling
$\lambda_{k}\equiv(B_{k}+B_{-k})u_{k}v_{k}$ is appreciable only in regions
where the product $u_{k}v_{k}$ is appreciable. That is, possibility of
superrandiace is intimately connected with pairing phenomenon. The matrix
element Bk = - B-kfor electron-photon coupling. And Bk = + B-k for electron-
acoustic phonon coupling BCS . Thus in simple geometries, $\lambda_{k}=0$ for
electron-electromagnetic radiation coupling.

Our restriction to bosonic subspace and our effective Hamiltonian is a good
low temperature approximation because i) kBT << $\Delta_{0}$, the minimum
superconducting gap and density of thermal fermionic quasi particles is small
and ii) when $\lambda>>L$, the boson mode excites only the A-atoms.

More importantly, we have ignored back reaction, i.e., self consistent
modification of uk, vk or gap function $\Delta_{k}$, arising from interaction
with the boson mode. We will see later that selfconsistent modification
reinforces superradiant superconductivity.

To illustrate superradiance, consider a simple Dicke Hamiltonian, with
identical two level atoms in resonance with the boson mode,
$H_{D}=\hbar\omega_{0}(b^{\dagger}b+\frac{1}{2})-\frac{\hbar\omega_{0}}{2}\sum_{i}\sigma^{z}_{i}+\frac{g}{\sqrt{N}}\sum_{i}(b^{\dagger}\sigma^{-}_{i}+b\sigma^{+}_{i})$.
For every value of Nt there is an unique ground state, a nodeless in phase
superposition of degenerate states with real positive coefficients. The ground
state is a superradiant state capable of undergoing a spontaneous emission
with an emission strength that scales as N${}^{2}_{t}$.

For our purpose consider a superconductor at T = 0 in the presence of a
macroscopically occupied single boson mode $|N_{b}\rangle$ and allow the
coupled system to evolve in time. When $\hbar\omega_{0}$ start increasing
towards $\triangle_{0}$, minimum of the two quasi particle gap, a set of
k-points which are near resonance with energy of a boson quanta actively
participate in superradiance and modify the ground state wave function. Net
density of these active A-atoms depend on quasi particle density of states and
the coupling constant $\lambda_{k}$.

Dicke Hamiltonian, equation 6, admits Bethe Ansatz solutionBetheAnsatz  for
k-independent $\lambda_{k}=\lambda_{0}$. Using the approximation,
$\lambda_{k}\approx\lambda_{0}$ for our set of near resonant A-atoms, our
ground state wave function has Bethe Ansatz form:

$\displaystyle|SRSC\rangle$ $\displaystyle\sim$
$\displaystyle(b^{\dagger}+\sum_{k}w_{k}\leavevmode\nobreak\
\alpha^{\dagger}_{k\uparrow}\alpha^{\dagger}_{-k\downarrow})^{N_{b}}|BCS\rangle\otimes|0_{b}\rangle$
(7) $\displaystyle\equiv$
$\displaystyle(b^{\dagger}+\sum_{k}w_{k}\leavevmode\nobreak\
\sigma^{+}_{k})^{N_{b}}|BCS\rangle\otimes|0_{b}\rangle$

Here $|0_{b}\rangle$ is the vacuum of the single boson mode. Superradiance
mixes (hybridizes or entangles) two nearly degenerate neutral modes. One is
the single mode external Bose oscillator. Second is a coherent sum of zero
momentum Bogoliubov pair excitations,
$\sum_{k}w_{k}\alpha^{\dagger}_{k\uparrow}\alpha^{\dagger}_{-k\downarrow}$; or
equivalently an Anderson pseudo spin wave packet mode in k-space. It is easy
to show that the second boson mode is a dynamic deformation of the Cooper pair
wave function (in the relative coordinate of the two electrons, characterized
by wk). The center of mass degree of freedom of the Cooper pairs, and hence
the phase of superconducting order parameter is not directly influenced by
superradiance phenomenon.

Superradiance effect in an s-wave superconductor is maximum, when the boson
frequency $\hbar\omega_{0}$ passes through minimum gap 2$\Delta_{o}$, where
quasi particle density of states has a maximum. If a superconductor supports
excited Cooper bound states below $2\Delta_{0}$, depending on the symmetry of
the excited states, there will be enhanced superradiance around these bound
state energies.

It follows from our work that one should be able to see i) a well known
quantum phase transition HeppLieb , as a function $\omega,\lambda$ ii)
enhanced quantum entanglement TobiasEntanglement  around the transition point
and iii) superfluorescence Superfluorescence .

Application to Pseudogap Phase of Cuprates Having theoretically suggested
possibility of Dicke superradiance in a BCS superconductor, we will address
recent experimental observation liscLBCO ; liscYBCO  of femtosecond laser
induced transient superconductivity in the pseudogap normal state of some
cuprates. In the two experiments two different Cu-O bond stretching modes are
resonantly excited by an 80 meV ( $\sim$ 20 THz) femtosecond laser. In view of
resonance, laser pumps its energy and coherence to the infrared phonon mode.
Electronic subsystem receives its energy and coherence from the infrared mode.
We have a phonon-photon polariton Hamiltonian:

$H=\hbar\omega_{0}a^{\dagger}a+\hbar\omega_{0}b^{\dagger}b+g(a^{\dagger}b+H.c.)$
(8)

Here $(a^{\dagger},a)$ and $(b^{\dagger},b)$ are the photon and phonon
operators respectively. As wavelengths of 20 THz infrared radiation and the
optic modes is $\sim 150$ microns, we will approximate the wavelengths by size
of the sample. The phonon optical mode coupling ‘g’ is of the order of 10 meV.
This coupling will lead to interesting Rabi oscillation between two modes,
after the femtosecond photon pulse impinges on the superconducting crystal.

It is safe to assume that Cu-O stretching lattice modes in both experiments
modulate i) site energy and ii) the hopping matrix element ‘t’ of the tight
binding electronic Hamiltonian for cuprates. To leading order in the normal
coordinate displacement u of this mode we have $t=t_{0}+\frac{\partial
t}{\partial u}|_{0}u\equiv
t_{0}+\alpha_{t}(b^{\dagger}+b){\rm\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and}$ to is the value of hopping
integral in the absence of resonant excitation of Cu-O stretching infrared
mode.

As far as the pseudogap normal state of cuprates is concerned, there are
experimental evidences Ong ; RamanPseudoGap ; STM  and theoretical support
BZAEmeryKivelson  that this metallic state has substantial pairing amplitude
and a strong phase fluctuations. It is well described as a 2D vortex liquid
above a Kosterlitz-Thouless transition point. This is borne out by Nernst
effect Ong , Raman effect RamanPseudoGap , among other experiments. In what
follows we propose an effective Hamilotonian that is an expression of the fact
that pseudogap phase gap supports local superconductivity. We assumes presence
of equal density of positive and negative vortices that are quasi static and
spatially random. That is, the thermal vortices are slowly moving compared to
time scale of interest to us.

Our effective Hamiltonian for pseudogap normal state has the form:

$H_{\rm normal}=\sum\varepsilon_{m}\alpha^{\dagger}_{\sigma}\alpha_{m\sigma}$
(9)

The index m denotes eigen modes of Bogolibov quasiparticle operator
$\alpha^{\dagger}_{m\sigma}\equiv u_{m}c^{\dagger}_{m\sigma}+\sigma
v_{m}c_{-m-\sigma}$ . In view of presence of disordered vortices in the
background, single particle eigen modes are not Bloch states; some localized
and rest extended. Our conclusions hold good even for the d-wave symmetry
situation in cuprates.

In the absence of external magnetic field we have pairs of degenerate single
particle eigen states (m$\uparrow$, - m $\downarrow$), connected by time
reversal symmetry. As in the BCS case, a pair subspace (m$\uparrow$, - m
$\downarrow$) is occupied by A-atom in its ground or excited state or an
unpaired fermion. By using same arguments as in the BCS case, bosonic
excitation sector in the normal state pseudo gap phase, coupled to the single
phonon mode gives us a Dicke type pseudo spin Hamiltonian:

$H=\hbar\omega_{0}(b^{\dagger}b+\frac{1}{2})+\sum\varepsilon_{m}\sigma^{z}_{m}+\frac{1}{{\sqrt{N}}}\sum\lambda_{i}(\sigma^{+}_{m}b+\sigma^{-}_{m}b^{\dagger})$
(10)

Here the operator
$\alpha^{\dagger}_{m\uparrow}\alpha^{\dagger}_{{-m}\downarrow}\equiv\sigma^{+}_{m}$
excites an A-atom.

In terms of A-atom and ferminic quasi particle there is a key difference
between the BCS supercodnuctor and the cuprate superconductors above Tc. In a
standard BCS superconductor, the pair subspace $({\bf
k}\uparrow,{\bf-k}\downarrow)$, is dominated by fermionic quasi particles and
nearly vanishing density of A-atom. Whereas, in the pseudogap, which exists
over a wide temperature range above Tc, the pair subspace (m$\uparrow$, -m
$\downarrow$) is dominated by ground and excited A-atoms and nearly vanishing
density of fermions. This makes pseudogap phase special and susceptible for
transient superconductivity.

To understand how superradiance induces transient superconductivity in the
pseudo gap phase, we have to go beyond our model Hamiltonian (equation 10) and
consider selfconsistent modification of um and vm’s. We offer a feed back
mechanism. Qualitatively it is as follows. A subspace (m$\uparrow$, - m
$\downarrow$) contains A-atoms with high probability, in ground or excited
states; fermions with low probability. A fraction of excited A-atoms are in
resonance with the macroscopically occupied phonon mode. In view of
macroscopic occupancy, the boson mode stimulates the near resonant excited
A-atom to emit a boson and reach its ground state. In the process we create an
excess population of ground state A-atoms. Increase in density of ground state
A-atoms means increased superconducting correlation (increase in magnitude of
ukvk); consequently an increase in superradiance interaction. Thus there is a
positive feedback, which could establishes a transient long range
superconducting order. As pseudo gap phase extends to room temperatures in
some of the underdoped cupraets, our mechamism offers a possibility to observe
room temperature transient superconductivity. This is one more incentive for
authors of reference liscYBCO  to confirm their exciting observations.

To establish superconductivity in the normal state of a Kosterlitz Thouless
superconductor, what we need is only a spatial reorganization of random
thermal vortices into either i) a fluid of bound vortex-antivortex pairs as in
Kosterlitz-Thouless phase or ii) an ordered lattice of positive and negative
vortices (see note note5 ). The increased pairing correlation from
superradiance increases the core energy of the thermal vortices and a
corresponding increase of vortex pair binding energy. Resulting increase in
population of paired vortices help create transient superconductivity.

In addition to superconductors, it will be interesting look for superradiant
superfluidity in pairing dominatated fermion systems: superfluid He3, cold
atoms, heavy nucleii and nuclear matter.

Acknowledgement I thank - N Kumar, K P Sinha, R K Shankar and R Nityananda for
early discussions on photoinduced superconductivity; P W Anderson and N P Ong
for an encouraging discussion; N P Ong for bringing to my attention reference
RamanPseudoGap ; B. Keimer for an encouraging information liscYBCO ; DAE,
India for a Raja Ramanna Fellowship. This research was supported by Perimeter
Institute for Theoretical Physics.

## References

  * (1) J. Bardeen, J R Schreiffer and L. Cooper, Phys. Rev., 108, 1175 (1957); Introduction to Superconductivity, M. Tinkham (Dover, NY 2004)
  * (2) R. H. Dicke, Phys. Rev., 93, 99 (1954); Super-radiance, M.G. Benedict et al., (IOP Publishing, Bristol 1996)
  * (3) P. W. Anderson and A. H. Dayem, Phys. Rev. Lett., 13, 195 (1964); A. F. G. Wyatt et al., Phys. Rev. Lett., 16, 1166 (1966); A. H. Dayem, J. J. Wiegand, Phys. Rev., 155, 419 (1967); R. Escudero and H.J.T. Smith, Phys. Rev. B31, 2725 (1985); S. I. Vedeneev, D. K. Maude, and J. M. Byrne, Phys. Rev. B78, 052509 (2008)
  * (4) G. M. Eliashberg, JETP Letters, 11, 114 (1970); B. I. Ivlev and G. M. Eliashberg, JETP Letters 13, 333 (1971)
  * (5) C.S. Owen and D.J. Scalapino, Phys. Rev Lett. 25, 1559 (1972)
  * (6) N Kumar and K P Sinha, Phys. Rev., 1̱74, 482 (1968)
  * (7) D. R. McIntosh and J. Lindesay, Phys. Rev., B50, 15852 (1994)
  * (8) D. Fausti et al., Science, 331, 6014 (2011)
  * (9) S. Kaiser et al., arXiv:1205.466 v2. According to version 3, in view of a calibration error, part of the claim needs to be verified; experiments are being repeated. However, signal for transient superconductivitty, seen as an appearance of c-axis plasma edge, remains robust (B. Keimer, private communication)
  * (10) P.W. Anderson, Phys. Rev.,112, 1900 (1958)
  * (11) In circuit QED, a collective degree of freedom of a Josephson junction is called Josephson atom (see for example, M. Devoret, S. Girvin and R. Schoelkopf, Ann. Phys. (Leipzig), 16, 767 .(2007)). A-atom is different - it fills k-space and is a bulk property of the superconductor.
  * (12) Our pseudo spin ${\vec{\sigma}}_{\bf k}$ is related to Anderson’s pseudo spin, $\tau^{z}_{k}\equiv(1-c^{\dagger}_{k\uparrow}c_{k\uparrow}-c^{\dagger}_{-k\downarrow}c_{-k\downarrow}),\leavevmode\nobreak\ \tau^{+}_{k}\equiv c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow}$ and $\tau^{-}_{k}\equiv c_{-k\downarrow}c_{k\uparrow}$ by a ${\bf k}$-dependent rotation of quantization direction to $(\theta,\phi_{k})$, where $u_{k}\equiv\cos\frac{\phi_{k}}{2}$ and $v_{k}\equiv\sin\frac{\phi_{k}}{2}e^{i\theta}$.
  * (13) M. Gaudin, J. Phys. (Paris), 37, 1087 (1976); A. Kundu, J. Phys. A: Math. Gen., 37, L281 (2004); J. Dukelsky et al., Phys. Rev. Lett., 93, 050403 (2004). Depending on wk, pseudo spin wave packet mode may have a overlap with the Higgs amplitude mode of the superconducting order parameter.
  * (14) N. Lambert, C. Emary, and T. Brandes, Phys. Rev. Lett., 92, 073602 (2004)
  * (15) K. Hepp and E. Lieb, Annals of Physics, 76, 360 (1973); Y. K. Wang and F. T. Hioe, Phys. Rev., A7, 831 (1973)
  * (16) R. Bonifacio, L. A. Lugiato, Phys. Rev., A11 1507 (1975)
  * (17) Z.A. Xu et al., Nature 406, 486 (2000); Y. Wang, Lu Li and N.P. Ong, Phys. Rev. B73, 024510 (2006)
  * (18) I. Iguchi, T. Yamaguchi and A. Sugimoto NATURE, 412, 420 (2001); C.V. Parker et al., Nature, 468, 677 (2010) vanishing phase coherence J. Corson et al., Nature 398, 221-223 (18 March 1999)
  * (19) A. Dubroka et al., Phys. Rev. Lett. 106, 047006 (2011)
  * (20) G. Baskaran, Z. Zou and P. W. Anderson, Sol. St. Commn., 63, 973 (1987); V. Emery and S. Kivelson, Nature, 374, 434 (1995)
  * (21) I thank P.W. Anderson for suggesting this possibility.


11institutetext: 1 Center for Imaging Science, 2 Biomedical Engineering,
Rochester Institute of Technology, Rochester, NY USA
11email<EMAIL_ADDRESS>
3 Bioengineering Graduate Program, 4 Electrical Engineering and Computer
Science, 5 Information and Telecommunication Technology Center, University of
Kansas, Lawrence, KS, USA

# CNN-based Cardiac Motion Extraction to Generate Deformable Geometric Left
Ventricle Myocardial Models from Cine MRI

Roshan Reddy Upendra${}^{\textrm{{\char 0\relax}}}$ 11 Brian Jamison Wentz
3355 Richard Simon 22 Suzanne M. Shontz Cristian A. Linte 3344551122

###### Abstract

Patient-specific left ventricle (LV) myocardial models have the potential to
be used in a variety of clinical scenarios for improved diagnosis and
treatment plans. Cine cardiac magnetic resonance (MR) imaging provides high
resolution images to reconstruct patient-specific geometric models of the LV
myocardium. With the advent of deep learning, accurate segmentation of cardiac
chambers from cine cardiac MR images and unsupervised learning for image
registration for cardiac motion estimation on a large number of image datasets
is attainable. Here, we propose a deep leaning-based framework for the
development of patient-specific geometric models of LV myocardium from cine
cardiac MR images, using the Automated Cardiac Diagnosis Challenge (ACDC)
dataset. We use the deformation field estimated from the VoxelMorph-based
convolutional neural network (CNN) to propagate the isosurface mesh and volume
mesh of the end-diastole (ED) frame to the subsequent frames of the cardiac
cycle. We assess the CNN-based propagated models against segmented models at
each cardiac phase, as well as models propagated using another traditional
nonrigid image registration technique.

###### Keywords:

Patient-specific Modeling Deep Learning Image Registration Cine Cardiac MRI

##  1 Introduction

To reduce the morbidity and mortality associated with cardiovascular diseases
(CVDs) [3], and to improve their treatment, it is crucial to detect and
predict the progression of the diseases at an early stage. In a clinical set-
up, population-based metrics, including measurements of cardiac wall motion,
ventricular volumes, cardiac chamber flow patterns, etc., derived from cardiac
imaging are used for diagnosis, prognosis and therapy planning.

In recent years, image-based computational models have been increasingly used
to study ventricular mechanics associated with various cardiac conditions. A
comprehensive review of patient-specific cardiovascular modeling and its
applications is described in [18]. Cardiovascular patient-specific modeling
includes a geometric representation of some or all cardiac chambers of the
patient’s anatomy and is derived from different imaging modalities [8].

The construction of patient-specific geometric models entails several steps:
clinical imaging, segmentation and geometry reconstruction, and spatial
discretization (i.e., mesh generation) [13]. For example, Bello et al. [2]
presented a deep learning based framework for human survival prediction for
patients diagnosed with pulmonary hypertension using cine cardiac MR images.
Here, the authors employ a 4D spatio-temporal B-spline image registration
method to estimate the deformation field at each voxel and at each timeframe.
The estimated deformation field was used to propagate the ED surface mesh of
the right ventricle (RV), reconstructed from the segmentation map, to the rest
of the timeframes of a particular subject. Cardiac MRI is a current gold
standard to assess global (ventricle volume and ejection fraction) and
regional (kinematics and contractility) function of the heart under various
diseases. In particular, cardiac MRI enables the generation of high quality
myocardial models, which can, in turn, be used to identify reduced function.

In this work, we propose a deep learning-based pipeline to develop patient-
specific geometric models of the LV myocardium from cine cardiac MR images
(Fig. 1). These models may be used to conduct various simulations, such as
assessing myocardial viability. In our previous work [19], we introduced a
preliminary, proof of concept, CNN-based 4D deformable registration method for
cardiac motion estimation from cine cardiac MR images, using the ACDC dataset
[4]. Here, we demonstrate the use of the CNN-based 4D deformable registration
technique to build dynamic patient-specific LV myocardial models across
subjects with different pathologies, namely normal, dilated cardiomyopathy
(DCM), hypertrophic cardiomyopathy (HCM) and subjects with prior myocardial
infarctions (MINF). Following segmentation of the ED cardiac frame, we
generate both isosurface and volume LV meshes, which we then propagate through
the cardiac cycle using the CNN-based registration fields. In addition, we
demonstrate the generation of dynamic LV volume meshes depicting the heart at
various cardiac phases by warping a patient-specific ED volume mesh based on
the registration-based propagated surface meshes. Lastly, we compare these
meshes to those obtained by directly propagating the ED volume mesh using the
CNN-based deformation fields.

##  2 Methodology

###  2.1 Cardiac MRI Data

We use the 2017 ACDC dataset that was acquired from real clinical exams. The
dataset is composed of cine cardiac MR images from $150$ subjects, divided
into five equally-distributed subgroups: normal, MINF, DCM, HCM and abnormal
RV. The MR image acquisitions were obtained using two different MR scanners of
$1.5$ T and $3.0$ T magnetic strength. These series of short axis slices cover
the LV from base to apex such that one image is captured every $5$ mm to $10$
mm with a spatial resolution of $1.37$ mm2/pixel to $1.68$ mm2/pixel.

###  2.2 Image Preprocessing

We first correct for the inherent slice misalignments that occur during the
cine cardiac MR image acquisition. We train a modified version of the U-Net
model [14] to segment the cardiac chambers, namely LV blood-pool, LV
myocardium and RV blood-pool, from 2D cardiac MR images. We identify the LV
blood-pool center, i.e., the centroid of the predicted segmentation mask and
stack the 2D cardiac MR slices collinearly to obtain slice misalignment
corrected 3D images [19, 7].

Figure 1: Overview of the proposed CNN-based workflow to generate patient-
specific LV myocardial geometric model.

###  2.3 Deformable Image Registration

####  2.3.1 CNN-based Image Registration.

We leverage our 4D deformable registration method described in [19] which
employs the VoxelMorph [1] framework to determine the optical flow
representation between the slice misalignment corrected 3D images. The CNN is
trained using the following loss function:

${L}={L}_{\text{similarity}}+\lambda{L}_{\text{smooth}},$ (1)

where ${L}_{\text{similarity}}$ is the mean squared error (MSE) between the
target frame and the warped frame, ${L}_{\text{smooth}}$ is the smoothing loss
function to spatially smooth the registration field, and $\lambda$ is the
regularization parameter, which is set to $10^{-3}$ in our experiments.
Inspired by Zhu et al. [20], we use the Laplacian operator in the smoothing
loss function as it considers both global and local properties of the
objective function $y=x^{2}$ instead of the traditional gradient operator
which considers only the local properties of the function $y=x^{2}$. A
detailed comparison of both these smoothing loss functions with respect to
cardiac motion estimation from cine MR images is found in [19].

The 4D cine cardiac MRI datasets are composed of $28$ to $40$ 3D image frames
that cover the complete cardiac cycle. For this discussion, we shall refer to
the 3D images as $I_{ED}$, $I_{ED+1}$,…,$I_{ED+N_{T}-1}$ where $I_{ED}$ is the
end-diastole image frame, and $N_{T}$ is the total number of 3D images. We
employ the fixed reference frame registration method, wherein the task is to
find an optical flow representation between the image pairs
$\\{(I_{ED},I_{ED+t})\\}_{t=1,2,3,...,N_{T}-1}$.

During training, we use $110$ of the total $150$ MR image dataset for
training, $10$ for validation and the remaining $30$ for testing. The CNN for
cardiac motion estimation is trained using an Adam optimizer with a learning
rate of $10^{-4}$, halved at every $10^{th}$ epoch for $50$ epochs. Both, the
U-Net model used for slice misalignment correction and VoxelMorph network
trained to estimate cardiac motion were trained on NVIDIA RTX 2080 Ti GPU.

####  2.3.2 Conventional Image Registration.

We compare the performance of the VoxelMorph framework with that of the
B-spline free form deformation (FFD) nonrigid image registration algorithm
[15]. This iterative intensity-based image registration method was implemented
using SimpleElastix [12, 9], which enables a variety of image-registration
algorithms in different programming languages. The FFD algorithm was set to
use the adaptive stochastic gradient descent method as the optimizer, MSE as
the similarity measure and binding energy as the regularization function. The
FFD-based image registration was optimized in $500$ iterations, while sampling
$2048$ random points per iteration, on an Intel(R) Core(TM) i9-9900K CPU.

###  2.4 Mesh Generation and Propagation

We use the manual segmentation map of the ED frame to generate isosurface
meshes. The slice thickness of each MRI image slice is $5$ mm to $10$ mm,
however, in order to obtain good quality meshes, the segmentation maps were
resampled to a slice thickness of $1$ mm. We use the Lewiner marching cubes
[11] algorithm to generate the meshes from the resampled segmentation maps of
the ED frames, and then simplification techniques, such as vertex
simplification and edge collapse, were performed using MeshLab $2020.07$ [5].
The simplification techniques are repeated multiple times to reduce the number
of vertices until the mesh has been fully decimated while preserving the
anatomical integrity and aspect ratio of the isosurface meshes.

Volume meshes of the initial surface meshes at the end-diastolic phases for
four patients with various heart conditions were generated based on the
decimated patient-specific surface meshes using Tetgen 1.6 [17]. In
particular, a constrained Delaunay mesh generation algorithm was used to
generate tetrahedral meshes based on the triangulated surface meshes. Steiner
points were added within the boundary of the surface mesh so that the
tetrahedra maintained a radius-edge ratio of $1.01$ and a maximum volume of
$9$ mm3 as needed for generation of valid meshes [17].

Mesh quality assessment was performed on the ED volume meshes utilizing the
scaled Jacobian metric, which ranges between $-1$ to $+1$, where $+1$
indicates an ideal equilateral tetrahedron, while negative and zero scaled
Jacobian values indicate inverted and degenerate tetrahedral elements,
respectively. Tetrahedra with a scaled Jacobian greater than or equal to $0.2$
are considered acceptable [10]. The ED volume mesh has a minimum scaled
Jacobian value of $0.078$, which demonstrates a valid, non-tangled mesh.
However, the end-systole (ES) phase mesh contains some lower quality elements
indicated by lower minimum scaled Jacobian values.

To demonstrate the VoxelMorph-based motion extraction and propagation to build
patient-specific LV myocardial models, we generate two sets of volume meshes
at each cardiac frame for each patient in each pathology group (Fig. 2).

Figure 2: Pipeline to generate dynamic volume meshes (at cardiac frames (ED +
k)) by direct CNN-based propagation, as well as volume mesh warping based on
dynamic boundary meshes.

The first set is produced by propagating the volume meshes at the ED frame to
all the subsequent frames of the cardiac cycle using the deformation field
estimated by the VoxelMorph-based registration method. For the second set, the
ED volume mesh generated with Tetgen was used to generate the volume meshes
corresponding to the other cardiac phases. We employed the log barrier-based
mesh warping (LBWARP) method [16] to deform the ED volume mesh onto the target
surface mesh for the new cardiac phase (Fig. 3). The method computes new
positions for the interior vertices in the ED volume mesh, while maintaining
the mesh topology and point-to-point correspondence [16].

Figure 3: LV volume meshes at three cardiac phases (a) end-diastole; (b) end-
systole; and (c) mid-diastole generated using LBWARP.

Briefly, LBWARP first calculates a set of local weights for each interior
vertex in the initial (ED) volume mesh based on the relative inverse distances
from each of its neighbors, which specify the representation of each interior
vertex in terms of its neighbors. Next, the vertices in the ED surface mesh
are mapped onto the new surface boundary. Finally, the interior vertices in
the ED volume mesh are then repositioned to reflect the updated positions of
the boundary nodes, while maintaining edge connectivity and point-to-point
correspondence, and ultimately yielding the volume meshes that correspond to
each new cardiac phase.

##  3 Results and Discussion

To evaluate the registration performance, the LV isosurface (generated from
the ED image segmentation map) is propagated to all the subsequent cardiac
frames using the deformation field estimated by FFD and VoxelMorph. We then
compare these isosurfaces to those directly generated by segmenting all
cardiac image frames using a modified U-Net model [14] (Section 2.2), which we
refer to as the “silver standard”.

Table 1 summarizes the performance of the FFD and VoxelMorph registration by
assessing the Dice score and mean absolute distance (MAD) between the
propagated and directly segmented (i.e., “silver standard”) isosurfaces.

Fig. 4 illustrates the distance between the three sets of isosurfces
(segmented, CNN-propagated and FFD-propagated) for one patient from each
pathology. The MAD between the surfaces is less than 2 mm at all frames, with
the CNN-propagated isosurfaces being closest to the “silver standard”
segmented surfaces.

Table 1: Mean Dice score (%) and mean absolute distance (MAD) (mm) between FFD and segmentation (FFD-SEG), CNN and segmentation (CNN-SEG), and FFD and CNN (FFD-CNN) results. Statistically significant differences were evaluated using the t-test (* for p $<$ 0.1 and ** for p $<$ 0.05). | Normal | MINF | DCM | HCM
---|---|---|---|---
| Dice | MAD | Dice | MAD | Dice | MAD | Dice | MAD
FFD-Segmentation | 74.80 | 1.53 | 77.69 | 1.09 | 80.41 | 0.91 | 77.39 | 1.97
CNN-Segmentation | 80.41** | 1.15 | 81.21* | 0.87 | 83.39* | 0.91 | 82.46* | 1.09
FFD-CNN | 77.81 | 1.13 | 82.12 | 0.75 | 81.67 | 0.97 | 77.34 | 1.77

Figure 4: MAD between FFD- and CNN-propagated, and segmented (i.e., “silver
standard”) isosurfaces at all cardiac frames for all patient pathologies.

Figure 5: Mean node-to-node distance at each cardiac frame between the CNN-
propagated and LBWARP-generated volume meshes (left); mean (std-dev) node
distance across all frames for each patient pathology (right).

As mentioned in Section 2.4 and shown in Fig. 2, we generate two sets of
volume meshes at each frame of the cardiac cycle. Fig. 5 shows the mean node
distance between the two sets of volume meshes across all cardiac frames for
one subject in each of the four pathologies. Fig. 5 shows the mean node
distance between the two sets of volume meshes at each frame of the cardiac
cycle for the four subjects. It can be observed that the two sets of volume
meshes are in close agreement with each other, exhibiting a mesh-to-mesh
distance within 2 mm.

We also briefly investigated the effect of using initial-to-final frame vs.
adjacent frame-to-frame registration to extract the cardiac motion throughout
the cycle. Although the sequential registration method estimates smaller
deformation between two consecutive, adjacent image frames compared to the
larger deformations estimated by the initial-to-final frame registration,
their concatenation across several frames accumulates considerable
registration errors. As such, when using these concatenated registration-
predicted deformation fields to propagate the ED isosurfaces and volume meshes
to the subsequent cardiac phases, the Dice score and MAD between the
propagated and segmented geometries rapidly deteriorate, along with the
quality of the propagated surface and volume meshes.

Moreover, although the proposed VoxelMorph-based cardiac motion extraction
method can capture the frame-to-frame motion with sufficient accuracy, as
shown in this work, our ongoing and future efforts are focused on further
improving the algorithm by imposing diffeomorphic deformations [6]. This
improvement will help maintain a high quality of the meshes and prevent mesh
tangling and element degeneration, especially for the systolic phases.

##  4 Conclusion

In this work, we show that the proposed deep learning framework can be used to
build LV myocardial geometric models. The proposed framework is not limited to
any pathology and can be extended to LV and RV blood-pool geometry.

## Acknowledgments

This work was supported by grants from the National Science Foundation (Award
No. OAC 1808530, OAC 1808553 & CCF 1717894) and the National Institutes of
Health (Award No. R35GM128877).

## References

  * [1] Balakrishnan, G., Zhao, A., Sabuncu, M.R., Guttag, J., Dalca, A.V.: Voxelmorph: A learning framework for deformable medical image registration. IEEE Transactions on Medical Imaging 38(8), 1788–1800 (2019)
  * [2] Bello, G.A., Dawes, T.J., Duan, J., Biffi, C., de Marvao, A., Howard, L.S., Gibbs, J.S.R., Wilkins, M.R., Cook, S.A., Rueckert, D., et al.: Deep-learning cardiac motion analysis for human survival prediction. Nature Machine Intelligence 1(2), 95–104 (2019)
  * [3] Benjamin, E.J., Blaha, M.J., Chiuve, S.E., Cushman, M., Das, S.R., Deo, R., Floyd, J., Fornage, M., Gillespie, C., Isasi, C., et al.: Heart disease and stroke statistics-2017 update: a report from the american heart association. Circulation 135(10), e146–e603 (2017)
  * [4] Bernard, O., Lalande, A., Zotti, C., Cervenansky, F., Yang, X., Heng, P.A., Cetin, I., Lekadir, K., Camara, O., Ballester, M.A.G., et al.: Deep learning techniques for automatic MRI cardiac multi-structures segmentation and diagnosis: is the problem solved? IEEE Transactions on Medical Imaging 37(11), 2514–2525 (2018)
  * [5] Cignoni, P., Callieri, M., Corsini, M., Dellepiane, M., Ganovelli, F., Ranzuglia, G.: Meshlab: an open-source mesh processing tool. In: Eurographics Italian chapter conference. vol. 2008, pp. 129–136. Salerno, Italy (2008)
  * [6] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.: Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces. Medical image analysis 57, 226–236 (2019)
  * [7] Dangi, S., Linte, C.A., Yaniv, Z.: Cine cardiac MRI slice misalignment correction towards full 3D left ventricle segmentation. In: Medical Imaging 2018: Image-Guided Procedures, Robotic Interventions, and Modeling. vol. 10576, p. 1057607. International Society for Optics and Photonics (2018)
  * [8] Gray, R.A., Pathmanathan, P.: Patient-specific cardiovascular computational modeling: diversity of personalization and challenges. Journal of cardiovascular translational research 11(2), 80–88 (2018)
  * [9] Klein, S., Staring, M., Murphy, K., Viergever, M.A., Pluim, J.P.: Elastix: a toolbox for intensity-based medical image registration. IEEE transactions on medical imaging 29(1), 196–205 (2009)
  * [10] Knupp, P.M.: Algebraic mesh quality metrics for unstructured initial meshes. Finite Elements in Analysis and Design 39(3), 217–241 (2003)
  * [11] Lewiner, T., Lopes, H., Vieira, A.W., Tavares, G.: Efficient implementation of marching cubes’ cases with topological guarantees. Journal of graphics tools 8(2), 1–15 (2003)
  * [12] Marstal, K., Berendsen, F., Staring, M., Klein, S.: Simpleelastix: A user-friendly, multi-lingual library for medical image registration. In: Proceedings of the IEEE conference on computer vision and pattern recognition workshops. pp. 134–142 (2016)
  * [13] Morris, P.D., Narracott, A., von Tengg-Kobligk, H., Soto, D.A.S., Hsiao, S., Lungu, A., Evans, P., Bressloff, N.W., Lawford, P.V., Hose, D.R., et al.: Computational fluid dynamics modelling in cardiovascular medicine. Heart 102(1), 18–28 (2016)
  * [14] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 234–241. Springer (2015)
  * [15] Rueckert, D., Sonoda, L.I., Hayes, C., Hill, D.L., Leach, M.O., Hawkes, D.J.: Nonrigid registration using free-form deformations: application to breast mr images. IEEE transactions on medical imaging 18(8), 712–721 (1999)
  * [16] Shontz, S.M., Vavasis, S.A.: A mesh warping algorithm based on weighted laplacian smoothing. In: IMR. pp. 147–158 (2003)
  * [17] Si, H.: Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software (TOMS) 41(2), 1–36 (2015)
  * [18] Smith, N., de Vecchi, A., McCormick, M., Nordsletten, D., Camara, O., Frangi, A.F., Delingette, H., Sermesant, M., Relan, J., Ayache, N., et al.: euheart: personalized and integrated cardiac care using patient-specific cardiovascular modelling. Interface focus 1(3), 349–364 (2011)
  * [19] Upendra, R.R., Wentz, B.J., Shontz, S.M., Linte, C.A.: A convolutional neural network-based deformable image registration method for cardiac motion estimation from cine cardiac mr images. In: 2020 Computing in Cardiology. pp. 1–4. IEEE (2020)
  * [20] Zhu, Y., Zhou Sr, Z., Liao Sr, G., Yuan, K.: New loss functions for medical image registration based on Voxelmorph. In: Medical Imaging 2020: Image Processing. vol. 11313, p. 113132E. International Society for Optics and Photonics (2020)