# Coherent coupling between vortex bound states and magnetic impurities in 2D
layered superconductors

Sunghun Park  Departamento de Física Teórica de la Materia Condensada,
Instituto Nicolás Cabrera and Condensed Matter Physics Center (IFIMAC),
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  Víctor Barrena
Laboratorio de Bajas Temperaturas y Altos Campos Magnéticos, Departamento de
Física de la Materia Condensada, Instituto Nicolás Cabrera and Condensed
Matter Physics Center (IFIMAC), Unidad Asociada UAM-CSIC, Universidad Autónoma
de Madrid, E-28049 Madrid, Spain  Samuel Mañas-Valero  Instituto de Ciencia
Molecular (ICMol), Universidad de Valencia, Catedrático José Beltrán 2, 46980
Paterna, Spain  José J. Baldoví  Instituto de Ciencia Molecular (ICMol),
Universidad de Valencia, Catedrático José Beltrán 2, 46980 Paterna, Spain  Max
Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee
149, D-22761 Hamburg, Germany  Antón Fente  Laboratorio de Bajas Temperaturas
y Altos Campos Magnéticos, Departamento de Física de la Materia Condensada,
Instituto Nicolás Cabrera and Condensed Matter Physics Center (IFIMAC), Unidad
Asociada UAM-CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain
Edwin Herrera  Laboratorio de Bajas Temperaturas y Altos Campos Magnéticos,
Departamento de Física de la Materia Condensada, Instituto Nicolás Cabrera and
Condensed Matter Physics Center (IFIMAC), Unidad Asociada UAM-CSIC,
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  Federico Mompeán
Instituto de Ciencia de Materiales de Madrid, Consejo Superior de
Investigaciones Científicas (ICMM-CSIC), Sor Juana Inés de la Cruz 3, 28049
Madrid, Spain  Mar García-Hernández  Instituto de Ciencia de Materiales de
Madrid, Consejo Superior de Investigaciones Científicas (ICMM-CSIC), Sor Juana
Inés de la Cruz 3, 28049 Madrid, Spain  Ángel Rubio  Max Planck Institute for
the Structure and Dynamics of Matter, Luruper Chaussee 149, D-22761 Hamburg,
Germany  Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy
Facility (ETSF), Universidad del País Vasco CFM CSIC-UPV/EHU-MPC & DIPC,
Avenida Tolosa 72, 20018 San Sebastián, Spain  Eugenio Coronado  Instituto de
Ciencia Molecular (ICMol), Universidad de Valencia, Catedrático José Beltrán
2, 46980 Paterna, Spain  Isabel Guillamón  Laboratorio de Bajas Temperaturas y
Altos Campos Magnéticos, Departamento de Física de la Materia Condensada,
Instituto Nicolás Cabrera and Condensed Matter Physics Center (IFIMAC), Unidad
Asociada UAM-CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain
Alfredo Levy Yeyati  Departamento de Física Teórica de la Materia Condensada,
Instituto Nicolás Cabrera and Condensed Matter Physics Center (IFIMAC),
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  Hermann Suderow
Laboratorio de Bajas Temperaturas y Altos Campos Magnéticos, Departamento de
Física de la Materia Condensada, Instituto Nicolás Cabrera and Condensed
Matter Physics Center (IFIMAC), Unidad Asociada UAM-CSIC, Universidad Autónoma
de Madrid, E-28049 Madrid, Spain

###### Abstract

Bound states in superconductors are expected to exhibit a spatially resolved
electron-hole asymmetry which is the hallmark of their quantum nature. This
asymmetry manifests as oscillations at the Fermi wavelength, which is usually
tiny and thus washed out by thermal broadening or by scattering at defects.
Here we demonstrate theoretically and confirm experimentally that, when
coupled to magnetic impurities, bound states in a vortex core exhibit an
emergent axial electron-hole asymmetry on a much longer scale, set by the
coherence length. We study vortices in 2H-NbSe2 and in 2H-NbSe1.8S0.2 with
magnetic impurities, characterizing these with detailed Hubbard-corrected
density functional calculations. We find that the induced electron-hole
imbalance depends on the band character of the superconducting material. Our
results open interesting prospects for the study of coupled superconducting
bound states.

Figure 1: Interaction between CdGM and YSR states. a The YSR state is
schematically represented by the red and blue dots with arrows representing
the spin. The value of the superconducting gap as a function of the radius is
shown in green for energies above and below the Fermi level. The CdGM states
are schematically shown by colored circles. The difference in energy between
levels, $\Delta^{2}/E_{F}$, with $\Delta$ the superconducting gap and $E_{F}$
the Fermi energy, is strongly exaggerated. Green arrow on top schematically
represents the vortex core size $\xi_{V}$. The interaction between CdGM and
YSR bound states leads to a spatial shift in the CdGM states, which is
different for electron and hole components of the LDOS inside the
superconducting gap, $N(E>0)$ (b) and $N(E<0)$ (c). The difference is shown in
d and the LDOS is larger along the direction of the YSR impurity, as
corresponding to a hole band. For an electron band, we expect a different
shift in b and c, such that it leads to an inversion of the color scale in d.
The YSR impurity is represented by a black dot and the vortex core center by
the crossing point between black lines in b, c and d.

Bound states appear in superconductors at localized perturbations in the
superconducting order parameter. Both Caroli de Gennes Matricon (CdGM) states
at vortex cores Abrikosov1957 ; Caroli1964  and Yu-Shiba-Rusinov (YSR) states
at magnetic impurities Yu1965 ; Shiba1968 ; Rusinov1969  are examples of this
phenomenon. YSR states provide mixed electron-hole excitations that serve to
create the conditions needed for Majorana states. These are expected for
instance at the ends of magnetic chains of YSR atoms PhysRevB.84.195442 ;
Nadj-Perge2014 ; PhysRevLett.115.197204 ; Kimeaar5251 . On the other hand,
CdGM states have been proposed to isolate and manipulate Majoranas in a
topological superconductor PhysRevLett.98.237002 ; doi:10.1146/annurev-
conmatphys-030212-184337 ; Alicea_2012 ; PhysRevLett.100.096407 ; Sun2017 ;
PhysRevB.91.214514 .

The nature of YSR and CdGM states is, however, quite different. YSR states are
spin polarized and appear at a single or a few subgap energies and exhibit
oscillations at the Fermi wavelength $\lambda_{F}$ that can be resolved with
atomic scale local density of states (LDOS) measurements Yazdani1997 ;
Flatte1997 ; Menard2015 ; HEINRICH20181 . By contrast, CdGM states are spin
degenerate and form a quasi-continuum with a level separation
$\Delta^{2}/E_{F}$ (where $\Delta$ is the superconducting gap and $E_{F}$ is
the Fermi energy), which is usually small compared to $\Delta$. Thus, their
discreteness and their mixed electron hole character only appears at very low
temperatures or for $\Delta\approx E_{F}$ and in absence of scattering, in the
so called quantum limit Hayashi1998 . Otherwise, thermal excitations or
defects produce dephasing resulting in an electron-hole symmetric LDOS pattern
at vortex cores.

Thus, in most cases, CdGM states are electron-hole symmetric and their
features in the LDOS extend to much larger distances than those of YSR states.
Here we ask the question if we can build a hybrid quantum system consisting in
vortices close to magnetic impurities and transfer the quantum property of YSR
states, i.e. their electron-hole asymmetry, into the more extended CdGM states
far from the quantum limit.

As we show below, we indeed theoretically predict and experimentally observe
electron-hole asymmetric features in the LDOS of vortices in presence of
magnetic impurities. As we schematically represent in Fig. 1, a magnetic
impurity close to a vortex core induces a coupling between CdGM states with
$n$ and $n\pm 1$ angular momenta. This coupling produces a slight shift of the
charge density of the positive (negative) energy excitations towards
(outwards) the impurity with respect to their mean position, which remains
even away of the quantum limit. The sign of the coupling changes when the
bands at the Fermi energy have a hole (electron) character. The discrete
nature of CdGM states is thus revealed in the vortex core LDOS: the difference
between the electron and hole LDOS would exhibit an axial asymmetry as
illustrated in Fig. 1d, with a larger LDOS close to the position of the
magnetic impurity. We emphasize that, in contrast to the above mentioned
oscillations at the tiny $\lambda_{F}$ scale, the electron-hole asymmetric
feature in the vortex LDOS occurs at a much larger length scale.

Figure 2: Length scales of CdGM and YSR states. a Tunneling conductance from
the center of a vortex (yellow to blue) in 2H-NbSe1.8S0.2. Dashed lines
represent schematically the splitting of the zero bias conductance peak when
leaving the vortex core, a characteristic feature of CdGM states. In the lower
left inset we show a zero bias conductance map. In the lower right panel we
show as white circles the radially averaged zero bias conductance $G$ as a
function of the radius. b Tunneling conductance as a function of the position
at a magnetic Fe impurity (yellow and ochre curves are close to the impurity
center and separated approximately by $\lambda_{F}\approx$ 0.7 nm) and far
from it (blue) in 2H-NbSe1.8S0.2. The tunneling conductance map at the bias
voltage where we observe a maximum of the conductance in the yellow curve is
shown in the lower left panel. In the lower right panel as white circles we
show the spatial dependence of the conductance along a line from the position
of the impurity outwards. Red lines in the lower right insets follows roughly
the expected dependencies for the spatial variation of the LDOS for CdGM (a)
and YSR (b) states, discussed in previous work Hess1989 ; Hess1990 ;
Guillamon2008c ; Menard2015 . White bars in the lower left insets of a, b are
both 10 nm long.

The transition metal dichalcogenide 2D layered superconductor 2H-NbSe2 is the
first material where the vortex LDOS has been measured and one of the few
where the nature of CdGM states has been extensively studied, both in
experiment and theory Hess1989 ; Hess1990 ; Guillamon2008c . YSR impurities
have been also imaged in detail in this material Menard2015 ; Liebhaber2020 ;
Kezilebieke2018 . Vortex cores in 2H-NbSe2 ($T_{c}=7.2$ K) are highly
anisotropic, with a characteristic sixfold star shape. Previous work imaged
YSR impurities and vortex cores at the same time, but did not identify any
particular connection Menard2015 . On the other hand, when doped with S as in
2H-NbSe1.8S0.2 ($T_{c}=6.6$ K) the vortex core CdGM states are in-plane
isotropic, leading to round shaped, symmetric, vortex cores Fente2016 . Using
these two systems, we can study the interaction between YSR and CdGM states
for in-plane isotropic (2H-NbSe1.8S0.2) and anisotropic (2H-NbSe2) vortices.

To produce the YSR states in 2H-NbSe2 and 2H-NbSe1.8S0.2 we introduce Fe
impurities during sample growth (about 150 ppm), as identified after the
experiment using inductively coupled plasma atomic analysis. In this diluted
regime Fe impurities produce practically no changes in the residual
resistivity or $T_{c}$ both for 2H-NbSe1.8S0.2 or 2H-NbSe2. The amount of Fe
impurities is sufficiently small as to leave the superconducting gap and
vortex structure unaffected, but large enough to be easily detected in the
area occupied by a single vortex. We use a Scanning Tunneling Microscope (STM)
to measure the LDOS as a function of the position at 0.8 K. Samples are
cleaved at or below liquid He, to allow for a clean atomically flat surface
and the tip is prepared in-situ Suderow2011 .

It is of interest to first analyze the different length scales associated to
isolated CdGM and YSR states, particularly in the case of
2H-NbSe${}_{1.8}S_{0.2}$ (Fig. 2). As we show in Fig. 2a, CdGM states provide
a zero bias peak at the center of the vortex that decays with distance at a
scale which is generally larger than the coherence length (of approximately 10
nm in 2H-NbSe2 and 7 nm in 2H-NbSe1.8S0.2 as obtained from Hc2(T)) and is
magnetic field dependentKogan2005 ; Fente2016 . The vortex core size at the
magnetic fields considered here is of $\xi_{V}\approx$ 30 nmFente2016 . The
zero bias peak splits when leaving the vortex core, as shown in previous work
Hess1989 ; Hess1990 ; Guillamon2008c ; Gygi1991 ; Clinton1992 ; Hayashi1998 ;
Fente2016 . The six-fold anisotropy characteristic of vortex cores in 2H-NbSe2
is washed out by the S substitutional disorder in 2H-NbSe${}_{1.8}S_{0.2}$. On
the other hand, a YSR state in 2H-NbSe${}_{1.8}S_{0.2}$ is shown in Fig. 2b.
There is a conductance peak within the superconducting gap which changes from
positive to negative bias voltage values at a scale of order of $\lambda_{F}$
(about 0.7 nm, see lower right inset in Fig. 2b and Ref. Menard2015 ). At the
same time, the height of the peak decreases exponentially with distance (lower
right inset of Fig. 2b). Thus we see that the effect of YSR states is
transposed two orders of magnitude in distance, from $\lambda_{F}$ to
$\xi_{V}$, and leads to the electron-hole asymmetric CdGM states we discuss
below.

Figure 3: YSR states at magnetic impurities. Spin-density isosurface of
2H-NbSe2 a and of 2H-NbSe1.8S0.2 b obtained from DFT+U calculations. Red
stands for spin up and blue for spin-down charge density. The isosurface is
plotted for an imbalance by 0.002 of the spin density. The Fe atom (not shown)
is located on the red spot, substituting a Nb atom. Se atoms are orange and Nb
atoms ochre colors. We represent the two uppermost layers (Se and Nb) from the
top. In b there are S atoms in layers below, shown in the Supplementary
Information. Calculations are made on four layers, using $4\times 4\times 2$
sized slabs (see Supplementary Information for further details). c, d Measured
tunneling conductance maps at a Fe impurity in 2H-NbSe2 c and in
2H-NbSe1.8S0.2 d. White arrows provide the crystalline directions of the
atomic Se lattice.

To better understand the magnetism at the Fe impurities we calculated the
electronic structure and the spin-density isosurface, i.e. the difference
between spin-up and spin-down charge densities (Fig. 3), by means of Hubbard-
corrected density functional theory (DFT+U). We constructed a slab model with
a 4 $\times$ 4 $\times$ 2 supercell that contains one isolated magnetic
impurity at a Nb site. To model 2H-NbSe1.8S0.2, we introduced approximately
10% S atoms (115 Se and 13 S atoms) randomly distributed. Computational
methods are detailed in the Supplemental Information. In 2H-NbSe2 we observe
clearly a strong magnetic moment on the Fe atom (Fig. 3a). The same behavior
is observed at the Fe site in 2H-NbSe1.8S0.2 with however a large
antiferromagnetic coupling with immediately neighboring Se atoms that become
spin-polarized (Fig. 3). Importantly, this coupling breaks the six-fold in-
plane symmetry of the Se lattice. The introduction of S atoms lowers the
symmetry from six to three fold. We show in Fig. 3c,d the measured tunneling
conductance map $G(\mathbb{r})=G(x,y)$ on a YSR state on 2H-NbSe2 (Fig. 3c)
and on 2H-NbSe1.8S0.2 (Fig. 3d). We see that the LDOS at YSR impurities is
modified from a six-fold star shape in 2H-NbSe2 to a predominantly three-fold
star shape in 2H-NbSe1.8S0.2 and ascribe this effect to the symmetry breaking
in the lattice induced by the S distribution in 2H-NbSe1.8S0.2, as suggested
by our calculations (Fig. 3a,b). The three fold symmetry is smeared at the
scale of $\xi_{V}$ leading to the observed round vortex cores shown in Fig.
2a.

Figure 4: Vortex cores coupled with YSR states. a Zero bias conductance map of
a vortex in 2H-NbSe1.8S0.2. b Map showing the difference between the
normalized tunneling conductance at positive and negative bias voltages
$\frac{\delta
G(\mathbb{r},V)}{G_{0}}=\frac{G(\mathbb{r},V)-G(\mathbb{r},-V)}{G_{0}}$ at
$|V|=$0.2 mV and in the same field of view as a. c $\frac{\delta
G(\mathbb{r},V)}{G_{0}}$ obtained from the calculation described in the text.
The same quantities are plotted for 2H-NbSe2 in d, e and f. Color scales are
given by the bars on the right side of panels a and d for the zero bias
conductance, and in between panels b, c and e, f for $\frac{\delta
G(\mathbb{r},V)}{G_{0}}$. Fe impurities are marked by black dots. Vortex
centers are at the crossing point between the black lines. White arrows in a,b
give the directions of the Se lattice.

We show vortices in close proximity to YSR impurities in Fig.4a,d. When we
make the difference between images taken at positive and negative bias
voltages, $\frac{\delta
G(\mathbb{r},V)}{G_{0}}=\frac{G(\mathbb{r},V)-G(\mathbb{r},-V)}{G_{0}}$ (with
$G_{0}$ the averaged tunneling conductance for bias voltages above the gap),
we observe that vortex cores are not axially symmetric (Fig.4b,e). In contrast
(as we show in detail in the Supplemental Information), $\frac{\delta
G(\mathbb{r},V)}{G_{0}}$ is axially symmetric in absence of YSR impurities.

As stated above, we trace the broken axial symmetry to the interplay between
vortex and YSR states. We have calculated the perturbation to a rotationally
symmetric vortex induced by a magnetic impurity located close to the vortex
core. As shown in the Supplemental Information, we start with a two
dimensional superconductor described by a Bogoliubov-de Gennes Hamiltonian. We
find discrete energy levels $E_{n}$ and the shape of electron and hole wave
functions $\psi_{n}^{+}$, $\psi_{n}^{-}$ of CdGM vortex bound states (with $n$
the angular momentum number). Magnetic YSR impurities are characterized as
usual by the exchange coupling $J$ at the impurity sites Huang2020 . This
coupling leads to an effective Hamiltonian in the subspace spanned by the
states $\psi_{n-1},\psi_{n},\psi_{n+1}$, with solution $\tilde{\psi}_{n}^{+}$,
$\tilde{\psi}_{n}^{-}$. Without YSR impurities, the vortex core LDOS obtained
is always axially symmetric, as found previously. There are slight, axially
symmetric electron-hole variations at the Fermi wavelength $\lambda_{F}$
scale, which are smeared out due to dephasing except in the quantum limit. The
vortex core LDOS with YSR impurities, obtained from $\tilde{\psi}_{n}^{+}$,
$\tilde{\psi}_{n}^{-}$, is however axially asymmetric. The asymmetry is due to
the spatial shift in the perturbed CdGM states $\tilde{\psi}_{n}^{+}$,
$\tilde{\psi}_{n}^{-}$ and is induced by the mixing between adjacent CdGM
levels ($n+1$ and $n-1$). This asymmetry is roughly given by $\frac{\delta
G(\mathbb{r},V)}{G_{0}}\propto|\tilde{\psi}_{n}^{+}|^{2}-|\tilde{\psi}_{n}^{-}|^{2}\propto\pm
J^{2}e^{-4r_{p}/\xi_{V}}\cos(\theta-\theta_{p})$, where $\theta$ is the polar
angle with respect to the vortex center, $r_{p}$ and $\theta_{p}$ provide the
length and the angle of the line joining the vortex center and the impurity
position and the $\pm$ sign depends on whether the effective mass is negative
or positive, i.e. whether the bands have a hole or an electron like character.
The magnitude of the perturbation decays exponentially with the distance from
the impurity to the vortex center.

As we show in Fig. 4, the observed LDOS asymmetry can be qualitatively
reproduced using our theory (Fig. 4c,f). For that purpose we introduce an
impurity distribution corresponding to the one in the experiments and add the
contribution of each impurity to the asymmetry. Furthermore, we use an
isotropic gap for 2H-NbSe1.8S0.2 and a six-fold anisotropic gap for 2H-NbSe2.
Detailed parameters of the calculation are provided in the Supplemental
Information. Here we highlight that the exchange coupling $J$ is negative,
corresponding to the antiferromagnetic exchange found in Fig. 3a,b and that we
can use the same value for all the impurities. In practice, due to the already
mentioned distance dependence, the asymmetry $\frac{\delta
G(\mathbb{r},V)}{G_{0}}$ is however dominated by the few impurities which are
closest to the vortex core.

Let us note that vortices in 2H-NbSe2, with their characteristic strong six-
fold star shape, present a rather involved shape of the asymmetry (Fig. 4f).
This suggests that the spatial extension of CdGM states determines the overall
shape of the asymmetry.

In all, we conclude from our combined theoretical and experimental work that
YSR states produce electron-hole asymmetric vortex cores. The YSR states allow
visualizing the discrete nature of CdGM levels and their electron-hole
asymmetry is translated to large scales. Our theory also suggests that a
superconductor with predominantly electron band character should lead to an
opposite shift in the LDOS. For example, vortices have been observed in
$\beta-$Bi2Pd Herrera2015 ; PhysRevLett.120.167001 ; lv2017experimental ,
which has predominantly electron character sakano2015topologically ;
iwaya2017full . YSR states in $\beta-$Bi2Pd have been observed
PhysRevLett.120.167001  but their influence on CdGM states has not yet been
addressed.

Bound states in vortex cores have been considered in the past mostly to
address the influence of pair potential disturbances on vortex pinning
PhysRevB.57.5457 . Recent calculations find vortex core states induced by
magnetic or nonmagnetic impurities which result in small modifications of the
superconducting gap parameter PhysRevB.103.024510 . However in this work we
find that vortex positions and gap parameter do not exhibit visible changes at
energies away from the gap edge in presence of magnetic impurities.

Unconventional d-wave or s$\pm$ superconductors often show pair breaking at
atomic impurities and vortices at the same time Fischer2007 ; Massee2015 ;
Allan2013 ; chen2021friedel . Our results suggest that vortex bound states
might be strongly influenced by such impurities, although a theory for the
interplay between YSR and CdGM states within a non s-wave pairing interaction
is required to analyze this further.

Finally, let us note that one can also envisage experiments using atomic
manipulation or deposition to place impurities at certain positions. The
position of vortices is easily modified by changing the magnetic field. Groups
of geometrically arranged YSR impurities can lead to significant spatially
distortions of the LDOS of vortex states.

###### Acknowledgements.

This work was supported by the EU (ERC-StG-679080, ERC-2015-AdG-694097, ERC
AdG Mol-2D 788222 and H2020-MSCAIF-2016-751047), by EU program Cost (CA16218,
Nanocohybri and MOLSPIN CA15128), by the Spanish State Research Agency
(FIS2017-84330-R, RYC-2014-15093, CEX2018-000805-M, MAT2017-89993-R co-
financed by FEDER and CEX2019-000919-M), by the Deutsche
Forschungsgemeinschaft (DFG) under under Germany’s Excellence Strategy -
Cluster of Excellence Advanced Imaging of Matter (AIM) EXC 2056 - 390715994,
RTG 1995 and GRK 2247, by the Comunidad de Madrid through program NANOMAGCOST-
CM (Program No.S2018/NMT-4321) and by the Generalitat Valenciana (Grupos
Consolidados IT1249-19, Prometeo Programme, iDiFEDER/2018/061 and
CDEIGENT/2019/022). AR acknowledges support by the MPI - New York City Center
for Non-Equilibrium Quantum Phenomena. HS and IG acknowledge SEGAINVEX at UAM
for design and construction of STM cryogenic equipment. SP acknowledges Banco
Santander (ref. 102I0112).

## References

  * (1) Abrikosov, A. A.  On the magnetic properties of superconductors of the second group.  _Zh. Eksp. Teor. Fiz._ 32, 1442–1452 (1957).
  * (2) Caroli, C., deGennes, P. G. & Matricon, J.  Bound fermion states on a vortex line in a type II superconductor.  _Physics Letters_ 9, 307 – 309 (1964).
  * (3) Yu, L.  Bound state in superconductors with paramagnetic impurities.  _Acta Physica Sinica_ 21, 75 (1965).
  * (4) Shiba, H.  Classical spins in superconductors.  _Progress of Theoretical Physics_ 40, 435–451 (1968).
  * (5) Rusinov, A. I.  Superconductivity near a paramagnetic impurity.  _J. Exptl. Theoret. Phys Lett._ 9, 85–87 (1968).
  * (6) Choy, T.-P., Edge, J. M., Akhmerov, A. R. & Beenakker, C. W. J.  Majorana fermions emerging from magnetic nanoparticles on a superconductor without spin-orbit coupling.  _Phys. Rev. B_ 84, 195442 (2011).
  * (7) Nadj-Perge, S. _et al._ Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor.  _Science_ 346, 602–607 (2014).
  * (8) Ruby, M. _et al._ End states and subgap structure in proximity-coupled chains of magnetic adatoms.  _Phys. Rev. Lett._ 115, 197204 (2015).
  * (9) Kim, H. _et al._ Toward tailoring Majorana bound states in artificially constructed magnetic atom chains on elemental superconductors.  _Science Advances_ 4, eaar5251 (2018).
  * (10) Bolech, C. J. & Demler, E.  Observing Majorana bound states in $p$-wave superconductors using noise measurements in tunneling experiments.  _Phys. Rev. Lett._ 98, 237002 (2007).
  * (11) Beenakker, C.  Search for Majorana fermions in superconductors.  _Annual Review of Condensed Matter Physics_ 4, 113–136 (2013).
  * (12) Alicea, J.  New directions in the pursuit of Majorana fermions in solid state systems.  _Reports on Progress in Physics_ 75, 076501 (2012).
  * (13) Fu, L. & Kane, C. L.  Superconducting proximity effect and Majorana fermions at the surface of a topological insulator.  _Phys. Rev. Lett._ 100, 096407 (2008).
  * (14) Sun, H.-H. & Jia, J.-F.  Detection of majorana zero mode in the vortex.  _npj Quantum Materials_ 2, 34 (2017).
  * (15) Björnson, K. & Black-Schaffer, A. M.  Probing vortex majorana fermions and topology in semiconductor/superconductor heterostructures.  _Phys. Rev. B_ 91, 214514 (2015).
  * (16) Yazdani, A., Jones, B. A., Lutz, C. P., Crommie, M. F. & Eigler, D. M.  Probing the local effects of magnetic impurities on superconductivity.  _Science_ 275, 1767–1770 (1997).
  * (17) Flatté, M. E. & Byers, J. M.  Local electronic structure of a single magnetic impurity in a superconductor.  _Phys. Rev. Lett._ 78, 3761–3764 (1997).
  * (18) Menard, G. C. _et al._ Coherent long-range magnetic bound states in a superconductor.  _Nat. Phys._ 11, 1013–1016 (2015).
  * (19) Heinrich, B. W., Pascual, J. I. & Franke, K. J.  Single magnetic adsorbates on s-wave superconductors.  _Progress in Surface Science_ 93, 1 – 19 (2018).
  * (20) Hayashi, N., Isoshima, T., Ichioka, M. & Machida, K.  Low-lying quasiparticle excitations around a vortex core in quantum limit.  _Phys. Rev. Lett._ 80, 2921–2924 (1998).
  * (21) Hess, H. F., Robinson, R. B., Dynes, R. C., Valles, J. M. & Waszczak, J. V.  Scanning-tunneling-microscope observation of the Abrikosov flux lattice and the density of states near and inside a fluxoid.  _Phys. Rev. Lett._ 62, 214–216 (1989).
  * (22) Hess, H. F., Robinson, R. B. & Waszczak, J. V.  Vortex-core structure observed with a scanning tunneling microscope.  _Phys. Rev. Lett._ 64, 2711–2714 (1990).
  * (23) Guillamon, I., Suderow, H., Guinea, F. & Vieira, S.  Intrinsic atomic-scale modulations of the superconducting gap of 2H-NbSe2.  _Phys. Rev. B_ 77, 134505 (2008).
  * (24) Liebhaber, E. _et al._ Yu–Shiba–Rusinov states in the charge-density modulated superconductor NbSe2.  _Nano Letters_ 20, 339–344 (2020).
  * (25) Kezilebieke, S., Dvorak, M., Ojanen, T. & Liljeroth, P.  Coupled Yu–Shiba–Rusinov states in molecular dimers on NbSe2.  _Nano Letters_ 18, 2311–2315 (2018).
  * (26) Fente, A. _et al._ Field dependence of the vortex core size probed by scanning tunneling microscopy.  _Phys. Rev. B_ 94, 014517 (2016).
  * (27) Suderow, H., Guillamón, I. & Vieira, S.  Compact very low temperature scanning tunneling microscope with mechanically driven horizontal linear positioning stage.  _Review of Scientific Instruments_ 82 (2011).
  * (28) Kogan, V. G. & Zhelezina, N. V.  Field dependence of the vortex core size.  _Phys. Rev. B_ 71, 134505 (2005).
  * (29) Gygi, F. & Schlüter, M.  Self-consistent electronic structure of a vortex line in a type-II superconductor.  _Phys. Rev. B_ 43, 7609–7621 (1991).
  * (30) Clinton, W. L.  Approximate solutions for the Bogoliubov de Gennes equations: Superconductor normal metal superconductor junctions and the vortex problem.  _Phys. Rev. B_ 46, 5742–5745 (1992).
  * (31) Huang, H. _et al._ Quantum phase transitions and the role of impurity-substrate hybridization in Yu-Shiba-Rusinov states.  _Communications Physics_ 3, 199 (2020).
  * (32) Herrera, E. _et al._ Magnetic field dependence of the density of states in the multiband superconductor $\beta$-Bi2Pd.  _Phys. Rev. B_ 92, 054507 (2015).
  * (33) Choi, D.-J. _et al._ Influence of magnetic ordering between cr adatoms on the Yu-Shiba-Rusinov states of the $\beta\text{$-$}{\mathrm{Bi}}_{2}\mathrm{Pd}$ superconductor.  _Phys. Rev. Lett._ 120, 167001 (2018).
  * (34) Lv, Y.-F. _et al._ Experimental signature of topological superconductivity and Majorana zero modes on $\beta$-Bi2Pd thin films.  _Science Bulletin_ 62, 852–856 (2017).
  * (35) Sakano, M. _et al._ Topologically protected surface states in a centrosymmetric superconductor $\beta$-PdBi2.  _Nature Communications_ 6, 1–7 (2015).
  * (36) Iwaya, K. _et al._ Full-gap superconductivity in spin-polarised surface states of topological semimetal $\beta$-PdBi2.  _Nature Communications_ 8, 1–7 (2017).
  * (37) Larkin, A. I. & Ovchinnikov, Y. N.  Resistance of layered superclean superconductors at low temperatures.  _Phys. Rev. B_ 57, 5457–5465 (1998).
  * (38) Bespalov, A. A. & Plastovets, V. D.  Large spectral gap and impurity-induced states in a two-dimensional Abrikosov vortex.  _Phys. Rev. B_ 103, 024510 (2021).
  * (39) Fischer, Ø., Kugler, M., Maggio-Aprile, I., Berthod, C. & Renner, C.  Scanning tunneling spectroscopy of high-temperature superconductors.  _Rev. Mod. Phys._ 79, 353–419 (2007).
  * (40) Massee, F. _et al._ Imaging atomic-scale effects of high-energy ion irradiation on superconductivity and vortex pinning in Fe(Se,Te).  _Science Advances_ 1, e1500033 (2015).
  * (41) Allan, M. P. _et al._ Anisotropic impurity states, quasiparticle scattering and nematic transport in underdoped Ca(Fe1-xCox)2As2.  _Nat. Phys._ 9, 220–224 (2013).
  * (42) Chen, X. _et al._ Friedel oscillations of vortex bound states under extreme quantum limit in KCa2Fe4As4F2 (2021).  eprint arxiv:2102.12150.


# Hierarchical Block Low-rank Approximation
of Cavity Radiation

Ivan Baburin<EMAIL_ADDRESS>Jonas Ballani<EMAIL_ADDRESS>John W. Peterson<EMAIL_ADDRESS>David Knezevic
<EMAIL_ADDRESS>
###### Abstract

In this paper we examine the use of low-rank approximations for the handling
of cavity radiation boundary conditions in a transient heat transfer problem.
The finite element discretization that arises from cavity radiation involves
far-field degree of freedom coupling and matrices characterized by a mixed
dense/sparse structure, which pose difficulties for the efficiency and
scalability of solvers. Here we consider a special treatment of the cavity
radiation discretization using a block low-rank approximation combined with
hierarchical matrices. We provide an overview of the methodology and discuss
techniques that can be used to improve efficiency within the framework of
hierarchical matrices, including use of the adaptive cross approximation (ACA)
method. We provide a number of numerical results that demonstrate the accuracy
and efficiency of the approach in practical problems, and demonstrate
significant speed-up and memory reduction compared to the more conventional
“dense matrix” approach.

###### keywords:

Cavity Radiation , Hierarchical Matrices , Radiation Boundary , Sparse
Computing

††journal: CMAME

label=., leftmargin=30pt, rightmargin=10pt

[eth]organization=ETH, city=Zürich, country=Switzerland

[akselos]organization=Akselos SA, city=Lausanne, country=Switzerland

##  1 Introduction

Blackbody radiation is a well-studied phenomenon in classical mechanics.
Though challenging to model mathematically (c.f. the “paradox” known as the
“ultraviolet catastrophe” [1]), the study of blackbody radiation at the
beginning of 20th century gave birth to Planck’s law and the concept of
quantization of energy, and led to further fundamental discoveries that laid
the foundations for the field of quantum theory. In the modern era of
scientific computation, there are still many challenges connected to blackbody
radiation, in particular the computation of radiation-based heat transfer. In
this article we consider a more general setting of “graybody” radiation, in
which the emissivity coefficient may be specified on the cavity surface. This
is relevant in many industrial and scientific settings, such as heat transfer
in reactors [2, 3, 4], or radiative heat transfer on celestial surfaces [5, 6,
7, 8].

The primary challenge of cavity radiation from the computational point of view
is that the interactions are “non-local”: in general every surface facet
interacts with every other surface facet, where the interactions are
quantified via the so-called “reflection matrix,” which in general is dense
and of size $n_{\rm facets}\times n_{\rm facets}$. This non-locality results
in a loss of sparsity in the matrices arising from finite element
discretizations, which in turn leads to very high computational cost when
$n_{\rm facets}$ becomes large (e.g. $>$10000$$). The most straightforward
approach — explicitly allocating and computing the dense matrices associated
with cavity radiation — is prohibitively expensive in terms of both memory
utilization and computation time, especially when applied on large or highly
detailed models.

We propose an alternative method of computing cavity radiation which avoids
dense matrix discretizations by exploiting the physical nature of radiation,
namely its fast (inverse-square) decay relative to the distance between
emitting surfaces. Due to this property, many individual entries in the view
factor matrix $F$ and reflection matrix $C$ (for the precise definitions we
refer to Section 2) underlying the computation of the cavity radiation flux
interact only weakly, and thus, if grouped together, can be approximated
reliably using low-rank blocks. Since the grouping or clustering of facets is
strongly linked to the underlying geometry, it is natural to store the $F$ and
$C$ matrices hierarchically, identifying large, weakly interacting matrix
blocks through a geometric clustering criterion. In the analogous context of
boundary element methods, hierarchical matrices have been successfully applied
to large-scale problems, see [9, 10] for comprehensive introductions.

While we were working on the implementation, a very similar approach was
presented by Potter et al. [8] in the context of thermal irradiance on
planetary surfaces, where the radiosity matrix shares many structural
properties with the view factor matrix $F$ from heat transfer. In their paper,
Potter et al. describe an efficient construction of a hierarchical matrix with
low-rank approximation of individual blocks, and provide an application of
their approach to the computation of irradiance on a lunar crater. The
numerical results demonstrate the potential benefits of the hierarchical
approach when applied to practical problems.

The novelty of our approach with regard to Potter et al. is twofold. First,
due to the nature of thermal radiation, we will perform computations not only
with the hierarchical view factor matrix $F$, but also with the inverse of the
reflection matrix $C$. Due to the nonlinear (temperature-dependent materials,
radiation flux terms) and time-dependent character of the heat equation, we
have to perform many applications of $C^{-1}$, but, because $C$ is itself
_independent_ of the current temperature solution, we are able to employ a
block-LU decomposition of $C$. While there is a non-trivial up-front
computational cost involved in the construction of the block-LU decomposition,
this cost is amortized over the many time steps and nonlinear iterations which
are performed during the course of the simulation. We note that the low-rank
structure of the reflection matrix is preserved in its block-LU decomposition,
and this fact is exploited to allow the action of $C^{-1}$ to be computed
quickly and with low memory requirements.

Second, we use the adaptive cross approximation (ACA) [11] to optimize the
construction of the view factor matrix, since we found that it was
significantly faster than the singular value decomposition when constructing a
low-rank approximation of individual blocks. Due to its strong reliability
guarantees, we have opted to use ACA with full pivoting (running linearly in
the number of entries). Further computational gains can be expected when using
ACA with partial pivoting (running linearly in the number of facets), at the
cost of potentially less reliable low-rank approximations due to the use of
partial information. A comprehensive discussion of different cross
approximation techniques can be found in [12].

Note also that in the current work we do not consider occlusion detection of
surface facets, and hence in Section 6 we conduct numerical experiments in
which occlusion detection is not required. Including occlusion detection is a
natural extension of the implementation, e.g. via a ray-tracing approach [13].
We emphasize that incorporating occlusion detection would not change the
block-low rank approach that we present here in any way, since it corresponds
to “zeroing out” entries in the view factor matrix that correspond to occluded
facets, which is equivalent to modifying only the “input” to the block low-
rank algorithm.

##  2 Problem Setting

###  2.1 Transient heat equation

We consider transient heat transfer in a domain $\Omega\subset\mathbb{R}^{3}$
with radiation boundary $\Gamma\subseteq\partial\Omega$ governed by:

$\displaystyle\rho c_{p}\frac{\partial T}{\partial t}$
$\displaystyle=\nabla\cdot(k\nabla T)+f$ $\displaystyle\quad\in\Omega$ (1)
$\displaystyle k\nabla T\cdot\vec{n}$ $\displaystyle=q(T)$
$\displaystyle\quad\in\Gamma$ (2)

where $T$ is the unknown temperature, $q(T)$ is the radiation heat flux, $f$
is the internal heat source, $k$ is the thermal conductivity, $\rho$ is the
material density, and $c_{p}$ is the specific heat. In addition to the
radiation boundary condition (2) shown above, there may in general be other
types of boundary conditions (Dirichlet, Neumann) present, though their
details are not crucial to our discussion. Furthermore, in order to fully
specify the problem, an initial condition for the temperature which applies on
all of $\Omega$ at time $t=0$ must also be provided.

The variational statement associated with (1) is then: find
$T\in\mathcal{H}^{1}(\Omega)$ such that

$\int_{\Omega}\rho c_{p}\frac{\partial T}{\partial
t}v\,\mathrm{d}{V}+\int_{\Omega}k\nabla T\cdot\nabla
v\,\mathrm{d}{V}-\int_{\Gamma}q(T)v\,\mathrm{d}{S}=\int_{\Omega}fv\,\mathrm{d}{V}$
(3)

holds for all test functions $v$ from the Sobolev space
$\mathcal{H}^{1}(\Omega)$. Equation (3) is in general a nonlinear equation for
$T$, both because the flux $q(T)$ depends on $T$ to the fourth power (see
(4)), but also because the thermal conductivity and specific heat may depend
on the temperature. The nonlinearity of (3) plays a crucial role in its
numerical solution, since we need to iterate (via Newton’s method) within each
time step to find the converged solution at the current time, and each
iteration requires the formation and solution of a large linear system of
equations involving the Jacobian matrix, as we will further elaborate in
Section 3.

###  2.2 Radiation flux

In this section, we briefly recall some of the fundamentals of thermal
radiation which are relevant to the current work, as presented in [14]. The
thermal energy flux (amount of energy emitted per unit time, per unit area)
from a so-called “graybody” thermal radiator is given, according to the
Stefan–Boltzmann law of radiation, by

$q=\epsilon\sigma T^{4}$ (4)

where $\sigma\approx 5.669\times 10^{-9}\;\frac{W}{m^{2}K^{4}}$ is the
Stefan–Boltzmann constant, $0<\epsilon\leq 1$ is the surface emissivity
(dimensionless), and $T$ is the absolute surface temperature of the body,
measured in Kelvin. As $\epsilon\rightarrow 1$, we recover the case of a
“blackbody” radiator. In this formulation, we also assume that the emissivity
is independent of the wavelength of the radiation, that is, the surfaces are
so-called “ideal” graybodies [15]. The SI units of the energy flux $q$ are
thus $\frac{W}{m^{2}}$.

In the case of two radiating graybody infinite flat planes with temperatures
$T_{1}$ and $T_{2}$ and surface emissivities $\epsilon_{1}$ and
$\epsilon_{2}$, respectively, it can be shown [15] that the net radiation
energy flux between them is proportional to the difference between their
temperatures raised to the fourth power:

$q=\frac{\sigma\epsilon_{1}\epsilon_{2}(T_{1}^{4}-T_{2}^{4})}{1-(1-\epsilon_{1})(1-\epsilon_{2})}$
(5)

The case of infinite parallel planes is particularly simple because one can
assume that all the energy emitted by plane 1 is absorbed by plane 2, and
vice-versa, which then admits a closed-form solution based on Kirchoff’s law.

In real applications, one must of course consider finite, non-planar surfaces
with non-uniform temperature fields. The standard approach taken in the
general case is thus to discretize the surfaces comprising the radiation
cavity into approximately planar finite regions or “facets,” and then consider
the exchange of energy between individual pairs of facets separately, the
details of which are discussed in §3.2.

##  3 Discretization details

###  3.1 Finite element discretization

First we discretize the three-dimensional domain $\Omega$ into a mesh
$\mathcal{T}_{h}$ with $n_{\text{nodes}}$ nodes comprising hexahedral or
tetrahedral elements with circumspherical diameter $\mathcal{O}(h)$. Our
approach does not explicitly depend on the geometric type of the elements
used, and should in theory work for all standard element types as well as
“hybrid” meshes consisting of a mixture of geometric types. Then, the faces of
the elements in $\mathcal{T}_{h}$ which lie on the radiation boundary $\Gamma$
can be considered as a lower-dimensional “manifold” mesh $\mathcal{M}_{h}$
consisting of $n_{\rm facets}:=|\mathcal{M}_{h}|$ facets, which will all be
allowed to emit and reflect radiation.

We next introduce the standard continuous Lagrange nodal basis
$\\{\varphi_{N}\\}_{N=1}^{n_{\text{nodes}}}$ associated to $\mathcal{T}_{h}$
and state the semi-discrete Galerkin finite element approximation of (3), find
$T_{h}$ such that:

$\int_{\Omega}\rho c_{p}\frac{\partial T_{h}}{\partial
t}\varphi_{N}\,\mathrm{d}{V}+\int_{\Omega}k\nabla
T_{h}\cdot\nabla\varphi_{N}\,\mathrm{d}{V}-\int_{\Gamma}q(T_{h})\varphi_{N}\,\mathrm{d}{S}=\int_{\Omega}f\varphi_{N}\,\mathrm{d}{V}$
(6)

holds for all $N=1,\ldots,n_{\text{nodes}}$, where

$\displaystyle T_{h}(x):=\sum_{M=1}^{n_{\text{nodes}}}T_{M}\varphi_{M}(x)$ (7)

and $T_{M}$ are the unknown coefficients. The argument $x$ in (7) indicates
the spatial dependence of the basis functions (and lack thereof for the
coefficients). The temperature on node $M$ is consequently given by $T_{M}$,
and the associated “radiation power” on node $M$ is henceforth defined as

$\eta_{M}:=T_{M}^{4}$ (8)

The volume integrals in (6) are computed in the usual manner by looping over
the elements of $\mathcal{T}_{h}$ and using standard Gauss quadrature. For the
radiation surface integral, we follow the approach presented in the Abaqus
manual [16], and treat the radiation flux $q_{i}$ as a constant on each facet
$i\in\mathcal{M}_{h}$, so that this term, which we shall denote by $Q_{N}$, is
approximated as

$\displaystyle
Q_{N}:=\int_{\Gamma}q(T_{h})\varphi_{N}(x)\;\text{d}S\approx\sum_{i=1}^{n_{\text{facets}}}q_{i}\int_{\Gamma_{i}}\varphi_{N}(x)\;\text{d}S$
(9)

In the preceding few equations, we have introduced the notational convention
of using $N$, $M$ to index node-based quantities, and $i$, $j$ to index facet-
based quantities, and this approach will be continued throughout the rest of
this section. To complete the description of the finite element
discretization, we need an expression for $q_{i}$, the radiation flux to the
$i$th facet, in terms of all other facets in the radiation cavity. This is
discussed in §3.2.

###  3.2 Radiation flux discretization

The key ingredient in modeling the interaction between discretized facets $i$
and $j$ of an enclosure is the so-called view factor matrix, whose entries are
given by [15]

$F_{ij}=\int_{A_{i}}\int_{A_{j}}\frac{\cos{\phi_{i}}\cos{\phi_{j}}}{\pi
R^{2}}\,\mathrm{d}{A_{j}}\,\mathrm{d}{A_{i}},$ (10)

where $A_{i}$ is the area of facet $i$, $R$ is the distance between facets $i$
and $j$, and $\phi_{i}$ is the angle between the line segment joining the
facet centroids and the plane-normal vector of facet $i$. We note in
particular that the diagonal of the view factor matrix consists of all zeros,
since (10) is only well-defined when $i\neq j$.

In this work we compute entries of $F$ via numerical quadrature on the facets.
The individual entries $F_{ij}$ are non-polynomial and hence may require
specialized quadrature rules. A typical approach, which we follow in our
implementation, is to apply adaptive Gauss quadrature in which the order of
the quadrature rule is increased for facets that are close together (hence
strongly interacting), and decreased for facets that are further away from one
another. In our testing this leads to accurate and computationally efficient
computation of the view factor matrix. However, we note that more
sophisticated integration techniques can be applied to the $F_{ij}$ if needed,
e.g. a detailed discussion for a more efficient treatment of singular
integrals representing the interaction of neighboring facets may be found in
[10].

A secondary quantity, which depends on the view factor matrix, is the so-
called reflection matrix, which is given by

$C_{ij}=\delta_{ij}-\frac{1-\epsilon_{i}}{A_{i}}F_{ij}\;\Leftrightarrow\;C=I-\Lambda
F,$ (11)

where $I$ denotes the identity matrix and $\Lambda$ is a diagonal scaling
matrix with $\Lambda_{ii}=(1-\epsilon_{i})/{A_{i}}$. With the preceding
definitions in mind, the radiation flux to facet $i$ can be written [16] using
the previously defined quantities as

$q_{i}=\frac{\sigma\epsilon_{i}}{A_{i}}\sum_{j}\epsilon_{j}\sum_{k}F_{ik}C^{-1}_{kj}(\eta_{j}-\eta_{i})$
(12)

where the radiation power $\eta_{i}$ on facet $i$ is treated as a constant
given by the nodal averaging/projection formula

$\displaystyle\eta_{i}:=\sum_{M=1}^{n_{\text{nodes}}}\frac{1}{A_{i}}\int_{\Gamma_{i}}\eta_{M}\varphi_{M}(x)\;\text{d}S$
(13)

where $A_{i}$ is the area of facet $i$, whose domain is $\Gamma_{i}$. We
observe that the sum in (13) nominally goes over all nodes $M$ in the finite
element mesh, but in fact only a small number of nodal basis functions
$\varphi_{M}$ are non-zero on any given facet $\Gamma_{i}$ due to the
localized nature of the finite element basis functions. Therefore, many of the
terms in the sum are zero. Associated with the formula in (13), it is
illustrative to define the “projection” operator

$\displaystyle
P_{iM}:=\frac{1}{A_{i}}\int_{\Gamma_{i}}\varphi_{M}(x)\;\text{d}S$ (14)

so that

$\displaystyle\eta_{i}=\sum_{M=1}^{n_{\text{nodes}}}P_{iM}\eta_{M}$ (15)

The projection operator $P_{iM}$ therefore maps node-based quantities to
facet-based quantities, and we can think of it as a sparse, rectangular
$n_{\text{facets}}\times n_{\text{nodes}}$ matrix for simplicity.

Combining all of the preceding ingredients, we finally get the following
expression for $Q_{N}$, originally defined in (9) (note: all sums go from 1 to
$n_{\text{facets}}$ unless otherwise noted)

$\displaystyle Q_{N}$
$\displaystyle:=\int_{\Gamma}q(T_{h})\varphi_{N}(x)\;\text{d}S$
$\displaystyle\approx\sum_{i}q_{i}\int_{\Gamma_{i}}\varphi_{N}(x)\;\text{d}S$
$\displaystyle=\sum_{i}\frac{\sigma\epsilon_{i}}{A_{i}}\sum_{j}\epsilon_{j}\sum_{k}F_{ik}C^{-1}_{kj}(\eta_{j}-\eta_{i})\int_{\Gamma_{i}}\varphi_{N}(x)\;\text{d}S$
$\displaystyle=\sum_{i}\left(\sum_{j}R_{ij}\eta_{j}-\sum_{k}R_{ik}\eta_{i}\right)\left(\frac{1}{A_{i}}\int_{\Gamma_{i}}\varphi_{N}(x)\;\text{d}S\right)$
$\displaystyle=\sum_{i}\sum_{j}\bar{R}_{ij}\eta_{j}P_{iN}$
$\displaystyle=\sum_{i}\sum_{j}\bar{R}_{ij}\left(\sum_{M=1}^{n_{\text{nodes}}}P_{jM}\eta_{M}\right)P_{iN}$
$\displaystyle=\sum_{i}\sum_{j}\sum_{M=1}^{n_{\text{nodes}}}P_{iN}\bar{R}_{ij}P_{jM}\eta_{M}$
(16)

where we have used (14), (15) and, for simplicity of notation, we have defined
the tensors $R$ and $\bar{R}$ as

$\displaystyle R_{ij}$
$\displaystyle:=\sigma\epsilon_{i}\epsilon_{j}\sum_{k}F_{ik}C^{-1}_{kj}$ (17)
$\displaystyle\bar{R}_{ij}$ $\displaystyle:=R_{ij}-\delta_{ij}\sum_{k}R_{ik}$
(18)

We can write the double facet sum from (16) as

$\displaystyle S_{NM}:=\sum_{i}\sum_{j}P_{iN}\bar{R}_{ij}P_{jM}$ (19)

or, employing matrix notation

$\displaystyle S:=P^{\top}\bar{R}P$ (20)

which gives us the following final expression for $Q_{N}$

$\displaystyle Q_{N}=\sum_{M=1}^{n_{\text{nodes}}}S_{NM}\eta_{M}$ (21)

###  3.3 Time discretization and Newton iteration

In the time domain, we discretize (6) using a standard implicit Euler scheme.
Let $T^{n}_{h}$ be the approximate solution at time level $t_{n}$, and
consider a fixed time step of size $\Delta t$ such that $t_{n+1}=t_{n}+\Delta
t$. The $N$th component of the discrete residual vector associated with (6) is
then given by

$\displaystyle\mathscr{R}_{N}(T^{n+1}_{h})$ $\displaystyle:=\int_{\Omega}\rho
c_{p}\left(\frac{T^{n+1}_{h}-T^{n}_{h}}{\Delta
t}\right)\varphi_{N}\,\mathrm{d}{V}+\int_{\Omega}k\nabla
T^{n+1}_{h}\cdot\nabla\varphi_{N}\,\mathrm{d}{V}$
$\displaystyle-Q^{n+1}_{N}-\int_{\Omega}f\varphi_{N}\,\mathrm{d}{V}$ (22)

The solution of the nonlinear system of equations $\mathscr{R}=0$ implied by
(22) is performed at each time step via Newton–Krylov [17, 18] iterations,
which require the assembly of the Jacobian $J$ associated with $\mathscr{R}$
based on the current temperature iterate.

Note that the Jacobian contribution due to the volume integral terms, which we
shall denote by $J^{\text{sparse}}$, will be sparse due to the usual
properties of the Galerkin finite element method. The flux term $Q_{N}^{n+1}$,
on the other hand, is both nonlinear and non-local, as discussed previously.
The Jacobian contribution for this term, $J^{\text{cav}}$, is thus a dense
matrix given by

$J^{\text{cav}}_{NM}:=\frac{\partial Q_{N}}{\partial
T_{M}}=\sum_{K=1}^{n_{\text{nodes}}}{S}_{NK}\frac{\partial}{\partial
T_{M}}\left(T^{4}_{K}\right)=4{S}_{NM}T^{3}_{M}$ (23)

where we have dropped the $n+1$ superscript for clarity. We also emphasize
that the matrix $S$, as defined in (20), forms the core part of
$J^{\text{cav}}$, and also is temperature independent. The total Jacobian

$J=J^{\text{sparse}}+J^{\text{cav}}$ (24)

therefore has both sparse and dense contributions, and requires specialized
handling, which is described in Section 3.4.

###  3.4 Iterative approach

In case of a large number of facets, the matrix size of the (dense) cavity
part $J^{\text{cav}}$ of the Jacobian becomes prohibitive, making direct
solution methods for the solution of the linear system in each nonlinear
iteration impractical. We hence resort to preconditioned iterative methods
(such as GMRES), requiring the implementation of the action of the Jacobian
and a suitable preconditioner onto an arbitrary vector $x$.

Thanks to the additive decomposition of the total Jacobian (24), we may
consider the two terms in the product

$Jx=J^{\text{sparse}}x+J^{\text{cav}}x$ (25)

in isolation. Since it is inexpensive to compute, we assume that
$J^{\text{sparse}}$ is available in fully assembled (sparse) form, and the
resulting product $J^{\text{sparse}}x$ can therefore be computed efficiently
by standard matrix-vector operations. In addition to being conveniently fully
assembled, we have also observed that $J^{\text{sparse}}$ is suitable for use
as a preconditioner. Our preconditioned iteration hence reads

$x\mapsto(J^{\text{sparse}})^{-1}\left(J^{\text{sparse}}x+J^{\text{cav}}x\right)$
(26)

where $(J^{\text{sparse}})^{-1}$ is computed explicitly via a sparse LU
decomposition. We note that more refined choices of preconditioners are
possible (in particular involving additional terms for the cavity part of the
Jacobian), but we have not investigated these in more detail since, in our
examples, we generally observed good performance (small numbers of GMRES
iterations) with this simple choice.

The approach for efficiently storing and computing matrix-vector products with
$J^{\text{cav}}$ relies on a block low-rank approximation of the view factor
matrix $F$, and a subsequent approximation of the inverse of the reflection
matrix $C$; these are described in detail in Sections 4 and 5. If the material
properties (thermal conductivity, specific heat, etc.) are temperature- and
time-independent, then $J^{\text{sparse}}$ is itself also independent of the
temperature, and thus theoretically only needs to be computed once and reused
for all subsequent iterations. Our present implementation is designed to
handle the general case of temperature-dependent material properties, and thus
does not currently take advantage of this particular optimization. This choice
is justified by the solver performance results reported in Section 6, which
show that, especially in the fine grid limit, the time spent recomputing
$J^{\text{sparse}}$ at each iteration is small compared to the time spent in
the computation of $J^{\text{cav}}$.

##  4 Block low-rank approximation of the view factor matrix

The view factor matrix $F$ defined in (10) represents the radiation
interaction between different facets from the mesh $\mathcal{M}_{h}$. The
particular form of (10) immediately gives rise to the following two
observations (cf. [8]):

  1. (i)

Although we do not consider occlusion detection in the current work, matrix
blocks corresponding to occluded groups of facets are identically zero. This
may result in major computational gains, as the corresponding matrix entries
need not be computed or stored explicitly.

  2. (ii)

The radiation interaction decays with increasing distance $R$. Not only does
this make the entries $F_{ij}$ smaller in absolute value, but it also leads to
a weak interaction between groups of facets that are far apart. This enables
data-sparse approximations of the corresponding matrix blocks of the view
factor matrix.

Our focus in this paper is to explore aspect ii using hierarchical block low-
rank approximations of the view factor matrix. The key elements of this
approximation are:

  1. (a)

Geometric clustering of the facets via a $k$-d tree to hierarchically
subdivide the row and column index set $\mathcal{I}=\\{1,...,n\\}$ of $F$.

  2. (b)

Using the subdivision from a, the block index set
$\mathcal{I}\times\mathcal{I}$ is hierarchically subdivided into smaller
blocks until the corresponding facet groups are sufficiently far apart. To
prevent unreasonably small matrix blocks, the subdivision is also stopped if a
user-defined minimal block size has been reached.

  3. (c)

Matrix blocks corresponding to distant groups of facets are approximated by
low-rank matrices. The rest of the blocks are stored as dense matrices.

In the following, we provide a more detailed description of the techniques
used for a, b and c.

###  4.1 Geometric Clustering

To group facets from the mesh $\mathcal{M}_{h}$ into clusters, we use a $k$-d
tree based on the centroids of the individual facets. We note that an
alternative, yet similar approach would be to consider a clustering approach
based on the bounding boxes of the facets.

With the index set $\mathcal{I}:=\\{1,...,n\\}$ corresponding to the centroids
of all facets as a starting point, the clustering proceeds iteratively by
geometrically splitting the point cloud of centroids via a hyperplane. Once
the split has been defined, the iteration continues recursively for the two
disjoint sets of centroids. The procedure for an arbitrary index set
$\tau\subset\mathcal{I}$ is defined in Algorithm 1.

Algorithm 1 BuildIndexTree($\tau\subset\mathcal{I}$)

if $\\#\tau\leq n_{\text{min}}$ then

$T\leftarrow\\{\tau\\}$ $\triangleright$ leaf index set

else

$\tau=:\tau_{\text{left}}\mathbin{\mathaccent
0{\cdot}\cup}\tau_{\text{right}}$ $\triangleright$ split using a hyperplane

$T_{\text{left}}\leftarrow$ BuildIndexTree$(\tau_{\text{left}})$

$T_{\text{right}}\leftarrow$ BuildIndexTree$(\tau_{\text{right}})$

$T\leftarrow T_{\text{left}}\uplus T_{\text{right}}$ $\triangleright$ left and
right subtrees

end if

The initial call to Algorithm 1 reads

$T_{\mathcal{I}}:=\texttt{BuildIndexTree}(\mathcal{I}),$ (27)

where the obtained _index tree_ $T_{\mathcal{I}}$ represents the hierarchical
split of all centroids. We note that in Algorithm 1 the recursion is only
continued in case the cardinality of the index set $\tau$ is larger than a
user-defined minimal leaf size $n_{\text{min}}\geq 1$ (say $n_{\text{min}}\sim
100$). This is to prevent the creation of arbitrarily small matrix blocks
which would negatively impact the performance of the overall hierarchical
approach.

Since the splitting in Algorithm 1 is unaware of the original ordering of the
facets in $\mathcal{M}_{h}$, the nodes $\tau\in T_{\mathcal{I}}$ generally
correspond to non-contiguous index sets. To associate the obtained splitting
to contiguous row and column index sets of the view factor matrix $F$, it is
hence convenient to permute the rows and columns of $F$ accordingly.

###  4.2 Block Subdivision

Using the hierarchical subdivision of the index set $\mathcal{I}$, we next
introduce a block subdivision of $\mathcal{I}\times\mathcal{I}$ defining the
eventual matrix blocks of the view factor matrix $F$. As before, we proceed by
recursively subdividing the block index set $\mathcal{I}\times\mathcal{I}$
into smaller subblocks until a stopping criterion is fulfilled.

Since we expect that groups of facets which are far apart interact only
weakly, we use a standard geometrically-based _admissibility condition_ (cf.
[9]) as a stopping criterion. Given two index sets
$\sigma,\tau\subset\mathcal{I}$, we denote the corresponding matrix block
$\sigma\times\tau$ as _admissible_ if

$\min(\operatorname{diam}(A_{\sigma}),\operatorname{diam}(A_{\tau}))\leq
c\operatorname{dist}(A_{\sigma},A_{\tau}),$ (28)

where $A_{\sigma},A_{\tau}$ denote the point clouds of centroids of the facets
belonging to index sets $\sigma,\tau$, respectively. The recursive procedure
for the block subdivision is defined in Algorithm 2.

Algorithm 2 BuildBlockTree($\sigma,\tau\in T_{\mathcal{I}}$)

if the block $\sigma\times\tau$ is admissible then

$T\leftarrow\\{\sigma\times\tau\\}$ $\triangleright$ leaf block

else

$T\leftarrow\emptyset$

for $\sigma^{\prime}\in\operatorname{child}(\sigma)$ and
$\tau^{\prime}\in\operatorname{child}(\tau)$ do $\triangleright$ recursion to
children

$T_{\text{sub}}\leftarrow\texttt{BuildBlockTree}(\sigma^{\prime},\tau^{\prime}$)

$T\leftarrow T\uplus T_{\text{sub}}$

end for

end if

The initial call to Algorithm 2 reads

$T_{\mathcal{I}\times\mathcal{I}}:=\texttt{BuildBlockTree}(\mathcal{I},\mathcal{I}),$
(29)

where the obtained _block tree_ $T_{\mathcal{I}\times\mathcal{I}}$ represents
the hierarchical block subdivision of the view factor matrix $F$.

######  Example 1.

Consider a simple equidistant arrangement of the centroids of $n$ facets along
a straight line. For a trivial minimal leaf size of $n_{\text{min}}=1$,
Algorithms 1 and 2 would result in a block structure as depicted in Figure 1.

$n=4$

$n=8$

$n=16$

$n=32$

Figure 1: Visualisation of the block structure of the view factor matrix for a
trivial one-dimensional arrangement of $n$ facets. Gray blocks are represented
as low-rank matrices and blue blocks as dense matrices.

######  Example 2.

Consider two finite parallel flat plates of size $L\times L$, separated by a
distance $L$, and let their surfaces be successively refined so that in total
there are $40\times 40$, $64\times 64$, and $106\times 106$ facets,
respectively, on each surface as shown in Fig. 2. Then, using a leaf size of
128, we obtain the block low-rank structure shown in Fig. 3. Note that the
“on-diagonal” blocks are initially represented as dense matrices on the coarse
grids, but they are eventually revealed to have some internal low-rank
structure as the mesh is refined.

Figure 2: Separated parallel flat plates with 3200 facets on the cavity
surface.

(a) $3200$ facets

(b) $8192$ facets

(c) $22472$ facets

Figure 3: View factor matrix block low-rank approximation for successively
refined parallel flat plates. These block structures were computed numerically
using the approach described in this Section.

###  4.3 Adaptive Cross Approximation

For the low-rank approximation of data sparse blocks one can use the singular
value decomposition (SVD) to compute the low-rank decomposition, but it is
well know that this becomes computationally expensive for large blocks. An
alternative approach that is typically more computationally efficient is to
apply adaptive cross approximation (ACA, cf. [11]). Using ACA we can construct
a $k$-rank approximation of a data sparse matrix $A\in\mathbb{R}^{m,n}$ just
in $\mathcal{O}(mnk)$ time, which is significantly faster than classical SVD.

The main idea of ACA is to iteratively determine sets
$\mathcal{R}=\\{i_{1},\ldots,i_{k}\\}$, $\mathcal{C}=\\{j_{1},\ldots,j_{k}\\}$
of row and column indices, respectively, to construct a low-rank approximation
of admissible matrix blocks $X$ of $F$ of the form

$X\approx
X(:,\mathcal{C})\left(X(\mathcal{R},\mathcal{C})\right)^{-1}X(\mathcal{R},:),$
(30)

see Figure 4 for an illustration. The representation (30) may then be further
postprocessed (or truncated) to obtain low-rank factors
$X\approx\tilde{X}:=UV^{\top}$ in standard form.

136247136247$\approx$ Figure 4: Low rank approximation by Adaptive Cross
Approximation (ACA) with row and column indices $\mathcal{R}=\\{2,4,7\\}$ and
$\mathcal{C}=\\{1,3,6\\}$, respectively

There are multiple approaches to how the pivoting should be performed (some of
them provide accuracy guarantees related to best low-rank approximation), but
exact details go beyond the scope of this paper. A more detailed discussion on
the different techniques can be found in [12]. In Section 6 we investigate
different approximation tolerances $\epsilon_{\text{rel}}$ which can be set
for ACA to satisfy:

$\frac{\|X-\tilde{X}\|}{\|X\|}\leq\epsilon_{\text{rel}}$ (31)

for all blocks $X$ as defined above. We will analyze the resulting accuracy
and demonstrate the computational advantage of ACA for practical cavity
radiation problems with respect to the choice of $\epsilon_{\text{rel}}$ in
the solver.

##  5 Cavity Jacobian assembly

The iterative approach introduced in Section 3.4 requires us to efficiently
construct the cavity part $J^{\text{cav}}$ of the Jacobian from (23) and to
implement its action on an arbitrary vector. We recall from (20) that the
matrix-vector product with $J^{\text{cav}}$ is defined through a (sparse)
projection operator $P$ and a dense coupling matrix $\bar{R}$ which depends on
the view factor matrix $F$ and the inverse of the reflection matrix $C$ via
(17), (18).

We now assume that the view factor matrix $F$ has been constructed by the
techniques introduced in Section 4 and is hence available in a hierarchical
block low-rank format. It is then straightforward to obtain the reflection
matrix $C=I-\Lambda F$ from (11) within the same hierarchical format as well,
since the diagonal matrix $\Lambda$ only introduces a simple row scaling in
each matrix block of $F$, and the identity matrix $I$ can be added to the
diagonal (dense) blocks without modification of the block structure.

The key computational task is hence to efficiently compute an (approximate)
factorization of $C$. For simplicity, we construct this factorization within
the same hierarchical block low-rank format as $F$, though alternative choices
are possible (cf. [9]). The remaining task within our iterative framework is
then to implement matrix-vector multiplication with matrices in hierarchical
block low-rank format and arbitrary vectors.

###  5.1 Hierarchical block LU decomposition of the reflection matrix

The factorization of the reflection matrix $C\approx LU$ proceeds recursively
using the block structure defined through Algorithm 2. The first block
subdivision leads to a block system

$\begin{bmatrix}C_{11}&C_{12}\\\
C_{21}&C_{22}\end{bmatrix}=\begin{bmatrix}L_{11}&0\\\
L_{21}&L_{22}\end{bmatrix}\cdot\begin{bmatrix}U_{11}&U_{12}\\\
0&U_{22}\end{bmatrix}$

which can be solved via the following four steps:

  1. 1.

perform the LU decomposition $L_{11}U_{11}=C_{11}$,

  2. 2.

solve for $U_{12}=L_{11}^{-1}C_{12}$ by forward substitution,

  3. 3.

solve for $L_{21}=C_{21}U_{11}^{-1}$ by backward substitution,

  4. 4.

perform the LU decomposition $L_{22}U_{22}=C_{22}-L_{21}U_{12}$.

The LU decompositions in Steps 1 and 4 are computed recursively in case the
corresponding blocks are further subdivided, whereas dense LU factorizations
are employed in leaf blocks. The forward and backward substitutions in Steps 2
and 3 are again implemented by recursive algorithms involving block matrix-
matrix products and additions, cf. [9] for a detailed outline.

It is important to note that matrix additions and multiplications will
formally lead to a rank increase in the corresponding low-rank blocks and
hence need to be coupled with a suitable rank truncation policy. In our
experiments, we have chosen to fix a relative accuracy at the value
$\epsilon_{\text{rel}}$ (equal to the approximation accuracy used for the view
factor matrix $F$) and to determine the ranks adaptively based on the singular
value decay within the low-rank blocks. As a consequence of the truncation,
the computed LU decomposition $C\approx LU$ is only approximate, but within a
prescribed, controllable accuracy.

###  5.2 Matrix-vector assembly

Once the initial factorization of the reflection matrix has been obtained, it
can be re-used within all iterations of the cavity Jacobian assembly. In
particular, the action $x\mapsto C^{-1}x$ is approximated by $x\mapsto
U^{-1}(L^{-1}x)$ which is implemented as a subsequent forward and backward
substitution. Similarly, the assembly of the action $x\mapsto Fx$ is performed
recursively, following the block structure of $F$. Most remarkably, the
matrix-vector assembly does not involve any further approximation but is
performed “exactly” (within the numerical accuracy of floating point
operations).

##  6 Numerical experiments

In this section, we consider radiative heat transfer between a collection of
spheres arranged in a Fibonacci spiral pattern [19]. The spiral configuration
is chosen because it leads to a distribution of both nearby and distant facets
that provide a good illustration of the capabilities of the block low-rank
approach presented in the preceding sections. (We also note that there is
little occlusion between the various spheres when they are arranged in this
manner, which justifies our omission of occlusion detection for this test
case.) Since we include the entire exterior boundary of each sphere in the
radiation cavity, there will be some facets (which we term “isolated” facets)
that cannot “see” any other facets within the cavity, and thus produce
identically zero rows in the view factor matrix.

###  6.1 Problem setup and representative results

The physical setup of the problem consists of a central, heated sphere
surrounded by several ambient temperature spheres of the same size arranged in
a spiral pattern, as shown in Fig. 5. Referring to the labeling system shown
in Fig. 5(a), sphere 1 is located at the origin, and then spheres 2, 3, 5, 7,
9, 11, and 13 are offset by $(\pm a_{n},\pm a_{n})$, from the previous
sphere’s location, where $a_{n}$ is the $n$th Fibonacci number, and the signs
are chosen so that a counter-clockwise spiral is formed. Finally, the
remaining five spheres are placed at the center of the quarter-circular arcs
joining the previous spheres. The spheres are meshed with four-node
tetrahedral elements, and have 88, 240, 432, 720, 1248, 1992, and 2664
exterior facets at Levels 1–7, respectively, as shown in Fig. 5(b). The total
number of mesh nodes, elements, and cavity facets (including all 13 spheres)
for each mesh level is given in Table 1.

(a) Spiral configuration

(b) Mesh refinement levels

Figure 5: Spheres (a) arranged in Fibonacci spiral configuration and (b)
showing the different mesh refinement levels used.

This geometric configuration provides several different opportunities for low-
rank approximation. First, because the distance between the spheres increases
rapidly as one proceeds around the spiral, there will be cavity facets in
relatively close proximity (near the center of the spiral) as well as far
apart, relative to the mesh spacing (facet sizes) present in the mesh. Second,
because the spheres are arranged in a planar spiral, there are some facets
(around 7% of the total number of facets) on the “top” and “bottom” of each
sphere, as well as along the exterior of the outer arm of the spiral, which
cannot “see” any other facets, i.e. are isolated.

Table 1: Number of facets, nodes, and elements for each mesh level in the Fibonacci spheres problem. We use a first-order Lagrange finite element, so the number of degrees of freedom (DOFs) in the temperature solve is equal to the number of nodes for each mesh level. Mesh Level | Num. Facets | Num. Nodes | Num. Elements
---|---|---|---
1 | $1144$ | $793$ | $2782$
2 | $3120$ | $3081$ | $14339$
3 | $5616$ | $7553$ | $38766$
4 | $9360$ | $13910$ | $72956$
5 | $16224$ | $32097$ | $175422$
6 | $25896$ | $59241$ | $328211$
7 | $34632$ | $90051$ | $505258$

Although this configuration of spheres does not conform to the classical
definition of a “cavity,” i.e. a convex, closed region, the characteristics
listed previously make it an interesting test problem for assessing the
effectiveness of the low-rank solver’s capabilities. The difference between a
closed and open cavity is reflected in the treatment of the view factor matrix
entries as follows. First, we define the (scaled) $i$th row sum

$\displaystyle s_{i}:=\frac{1}{A_{i}}\sum_{j}F_{ij}$ (32)

For the “truth” view factor matrix, we of course have $0\leq s_{i}\leq 1$, but
because of quadrature error and floating point tolerances, we sometimes obtain
row sums slightly outside this range. Therefore, we define the corresponding
“clamped” row sum

$\displaystyle c_{i}:=\left\\{\begin{array}[]{cc}0,&s_{i}<0\\\ 1,&s_{i}>1\\\
s_{i},&\text{otherwise}\\\ \end{array}\right.$ (36)

Then, in the case of a closed cavity, the $i$th row of the view factor matrix
is scaled by $c_{i}$, i.e. if $c_{i}\neq 0$,

$\displaystyle F_{ij}\leftarrow\frac{F_{ij}}{c_{i}}\quad\forall i,j$ (37)

In an open cavity, on the other hand, we do not scale the view factor matrix
rows, but instead define an “ambient” temperature, $T_{\infty}$, and include
an additional nonlinear boundary flux term in the residual (6) which is given
by

$\displaystyle
B_{N}(T_{h}):=\sum_{i}\int_{\Gamma_{i}}\left(1-c_{i}\right)\left(T_{h}^{4}-T_{\infty}^{4}\right)\varphi_{N}\,\mathrm{d}{S}$
(38)

to represent thermal energy which is “radiated to ambient.”

While the most “realistic” approach would be to treat this model as an open
cavity problem, with almost all the thermal energy being radiated to ambient,
and little heat actually being transferred between the spheres, we found that
by instead treating this model as an artificially (mathematically) closed
cavity using the approach described above, the solutions had a more pronounced
transient behavior and hence provided a better test case for our
implementation. That is, in the open cavity case, the surrounding spheres’
temperature remains nearly unchanged, while in the closed cavity case, the
surrounding spheres are heated non-uniformly. This choice allowed us to better
classify the accuracy of the different low-rank approaches relative to the
Direct solve approach.

The initial temperature of all the spheres is set to 300K, except for the
central sphere, which is set to 1000K. There is no internal heat generation,
and the radiation cavity is defined to be the union of the surfaces of all the
spheres. The final solution time is set to 1000s, and we use a uniform time
step of $\Delta t=25$s, so that a total of 40 time steps are computed. We use
the same time discretization for all mesh levels, though as mentioned
elsewhere, a more realistic mesh refinement study would refine in both time
and space simultaneously to reduce both temporal and spatial discretization
errors at approximately the same rate. We treat this model as a closed cavity,
which means that all thermal energy radiated from the central sphere is
eventually absorbed (and possibly re-emitted) by the other spheres.

Figure 6: Surface temperature of central spheres at time $t=1000s$.

At steady state, all spheres should reach the same equilibrium temperature. A
representative solution at time $t=1000s$ on the Level 7 mesh is shown in Fig.
6, where the color bar is set to a maximum of 380K in order to better
highlight the temperature variation on the non-central spheres; this causes
the central sphere to appear as a single color which is outside the maximum
temperature range. Also, note that we are only showing the central seven
spheres in Fig. 6 since the outer spheres are both too far away to show in the
same figure with a reasonable level of detail, and the temperature fields on
the outer spheres are more uniform and hence less interesting than the inner
spheres.

(a) $t=0$s

(b) $t=200$s

(c) $t=400$s

(d) $t=600$s

(e) $t=800$s

(f) $t=1000$s

Figure 7: Surface temperature of sphere 2 for various simulation times.

More detailed images of sphere 2 (the sphere closest to the central sphere) at
evenly-spaced time steps are shown in Fig. 7. This figure clearly shows that
the side of sphere 2 which is closest to the central sphere is preferentially
heated relative to the far side, which remains several degrees cooler. In Fig.
8, we show the internal details of the temperature field in sphere 3, by
creating a clip plane whose normal is orthogonal to the line segment which
joins the centers of spheres 1 and 3. This figure shows how heat is
transferred internally in the sphere (via conduction) after the surface is
heated via radiation. The temperature field contours within the sphere are
curved rather than straight, a result of the non-uniform radiation heat flux
being applied to the surface.

(a) $t=0$s

(b) $t=200$s

(c) $t=400$s

(d) $t=600$s

(e) $t=800$s

(f) $t=1000$s

Figure 8: Clip plane showing internal temperature of sphere 3 for various
simulation times.

###  6.2 Accuracy comparison with Direct solve method

When possible, we compare the performance of the low-rank solver to the so-
called “Direct solve” method in which the full view factor and reflection
matrices are explicitly constructed, and the reflection matrix inverse is
explicitly computed. That is, the $J^{\text{cav}}$ contribution from (25) is a
true dense matrix, and the Direct solver thus has no way of taking advantage
of e.g. the zero rows arising from the isolated facets mentioned previously.

In our comparisons below in which the low-rank truncation tolerance
$\epsilon_{\text{rel}}$ is varied, the Direct solution method is taken as the
“truth” solution, and the relative error in the low-rank solution methods is
computed via the discrete space-time norm

$\displaystyle
e(T):=\max_{t_{n}}\frac{\|T-T_{\text{Direct}}\|}{\|T_{\text{Direct}}\|}$ (39)

where $T_{\text{Direct}}$ is the solution computed by the Direct solve method,
the maximum is taken over all time steps $t_{n}$,
$n=1,\ldots,n_{\text{steps}}$, and $\|\cdot\|$ is the discrete
$\ell^{2}$-norm. In addition to comparing the temperature fields computed by
the Direct and low-rank approaches, we also provide detailed comparisons of
the overall runtime and memory usage of the different approaches, in order to
quantify the effectiveness of the low-rank solver approach.

In our comparisons we employed the sequence of meshes from Level 1 to 5 shown
in Fig. 5(b). On each mesh considered, we also varied the truncation tolerance
$\epsilon_{\text{rel}}$ (cf. (31)) from a maximum value of 0.75 to a minimum
value of $10^{-6}$. The results of these calculations are shown in Fig. 9(a),
and provide good evidence that the accuracy of the low-rank solver depends
continuously on the truncation accuracy $\epsilon_{\text{rel}}$, but does not
depend strongly on the mesh refinement level. In practice, this means that one
can choose a single value of $\epsilon_{\text{rel}}$ for an entire mesh
refinement study, and expect a consistent level of agreement with the
corresponding Direct solve at that mesh refinement level.

The relative error comparisons stop at Level 5 since the Direct solver
required more RAM than was available on the compute node in order to
successfully solve the problem on the Level 6 and 7 meshes. The low-rank
solver, on the other hand, had relatively modest peak memory requirements, as
shown in Fig. 9(b). For the low-rank solver, the peak memory usage increases
quickly for decreasing $\epsilon_{\text{rel}}$ on fine grids, but for the
largest $\epsilon_{\text{rel}}$ values (coarsest truncation accuracy), we
observe that the peak RAM usage for the Level 7 mesh is only about 1.87 times
larger than for the Level 1 mesh, despite the Level 7 mesh having over 30
times as many facets in the cavity.

(a) Relative error $e(T)$

(b) Peak RAM usage

Figure 9: Plots comparing (a) relative error $e(T)$ defined in (39) and (b)
peak RAM usage vs. $\epsilon_{\text{rel}}$ for all mesh refinement levels.

###  6.3 The choice of $\epsilon_{\text{rel}}$

Due to the importance of the choice of $\epsilon_{\text{rel}}$, we make some
further comments about this tolerance here. Recall that this tolerance refers
to the truncation accuracy of the data-sparse blocks within the block low-rank
approximation, whereas the dense blocks are represented exactly. This means
that in general we expect a higher-accuracy approximation of $F$ and $C$ than
the level specified by $\epsilon_{\text{rel}}$, given that we have a mix of
($\epsilon_{\text{rel}}$ approximate) low-rank and (exact) dense blocks. We
quantify this effect in Table 2, which shows that the approximation error in
$F$ is typically an order of magnitude lower than $\epsilon_{\text{rel}}$.

Table 2: Relative error in the low-rank approximation, $F_{\epsilon_{\text{rel}}}$, with respect to the “true” view factor matrix, $F$, for different $\epsilon_{\text{rel}}$ values. $\|\cdot\|$ represents the Frobenius norm. $\epsilon_{\text{rel}}$ | $\|F-F_{\epsilon_{\text{rel}}}\|/\|F\|$
---|---
$1.0\text{\times}{10}^{-1}$ | $9.529\text{\times}{10}^{-03}$
$1.0\text{\times}{10}^{-2}$ | $9.668\text{\times}{10}^{-04}$
$1.0\text{\times}{10}^{-3}$ | $9.416\text{\times}{10}^{-05}$
$1.0\text{\times}{10}^{-4}$ | $9.399\text{\times}{10}^{-06}$
$1.0\text{\times}{10}^{-5}$ | $9.372\text{\times}{10}^{-07}$
$1.0\text{\times}{10}^{-6}$ | $9.415\text{\times}{10}^{-08}$

A second point we note is that we are ultimately only interested in the
accuracy of the solution field, $T$, and the cavity radiation terms are simply
a “stepping stone” to $T$. This is an important point because the solution
fields in heat transfer are well known to be smooth, which tends to compensate
for approximation introduced in the representation of $F$ and $C$. Taken
together, these two effects explain why in general we observe
$e(T)\ll\epsilon_{\text{rel}}$, and hence this provides justification for
choosing relatively large values of $\epsilon_{\text{rel}}$ such as
$\epsilon_{\text{rel}}=0.1$.

This is welcome since these relatively large values of $\epsilon_{\text{rel}}$
are in the region in which we observe good scalability and high speed from the
block low-rank approach in our numerical results. Nevertheless, it is
important to keep in mind that in general the behavior of the block low-rank
approximation with respect to truncation tolerance will be problem-dependent,
and hence one should ideally justify a particular choice of
$\epsilon_{\text{rel}}$ by first running some test cases targeted at the model
problem of interest.

###  6.4 Performance vs. $\epsilon_{\text{rel}}$

The performance of the low-rank solver can roughly be categorized into two
main parts: (i) the time spent in computing and applying the low-rank LU
factorization, and (ii) everything else including finite element assembly,
GMRES iterations, file I/O, and other tasks. Here we will mainly be interested
in the time spent in part (i), since part (ii) represents tasks that must also
be performed in the Direct solver. Within part (i), it is useful to further
separate the time spent (a) building the view factor matrix, (b) computing the
low-rank LU factorization, and (c) applying the low-rank LU factorization.
Items (a) and (b) happen only once per simulation, and thus their costs are
amortized across all time steps, while the time required for item (c)
increases in direct proportion to the number of time steps computed.

The time spent in part (i), i.e. building the view factor matrix, computing
the low-rank factorization, and applying the LU factorization is shown, for a
range of $\epsilon_{\text{rel}}$ values, in Fig. 10. For coarse meshes, we
find that as the truncation tolerance $\epsilon_{\text{rel}}$ is decreased,
the percentage of time spent building the view factor matrix (Fig. 10(a))
plateaus or even decreases, while the percentage of time spent computing the
LU factorization (Fig. 10(b)) increases steadily. In general, decreasing the
truncation tolerance $\epsilon_{\text{rel}}$ has a more pronounced effect on
the LU factorization time for the finer (Level 6 and 7) grids; the LU
factorization time exceeds 50% of the total runtime in the finest case.

(a) Build view factor matrix

(b) Build LU factorization

(c) Apply LU factorization

Figure 10: Time required to (a) build the view factor matrix, (b) build the
low-rank LU factorization of the cavity Jacobian, and (c) apply the low-rank
LU factorization (as a percentage of the total runtime) vs.
$\epsilon_{\text{rel}}$ for each mesh level.

Our test problem only computes 40 time steps, so we would of course expect the
percentage of time spent building the view factor matrix and LU factorization
to be lower in a model that performed hundreds or even thousands of time
steps. Regarding the time spent applying the LU factorization, for the coarse
mesh levels tested, the percentage of time spent applying the factorization
increases with decreasing $\epsilon_{\text{rel}}$, but on finer meshes it
actually decreases with decreasing $\epsilon_{\text{rel}}$. This indicates the
relatively larger expense of building the low-rank LU factorization vs.
applying it for fine grids. In problems with more time steps, we would expect
the application time of the low-rank factorization to eventually dominate the
creation time.

###  6.5 Performance vs. Direct solver

The previous plots have focused on comparing the low-rank solver’s performance
for different mesh refinement levels and different $\epsilon_{\text{rel}}$
values, but it is of course also necessary to compare the performance of the
low-rank solver to the Direct solve method. In Fig. 11, we therefore switch
the $x$-axis of the plots from $\epsilon_{\text{rel}}$ to the number of cavity
facets, a quantity which is the same for both the Direct and low-rank solver
at all mesh levels. In Fig. 11(a), we compare the full simulation time (also
referred to as “Alive Time”) for the Direct solver and the low-rank solver
with three different $\epsilon_{\text{rel}}$ values of $10^{-1}$, $10^{-2}$,
and $10^{-3}$.

(a) Alive time

(b) Peak RAM usage

Figure 11: Comparison of (a) total simulation time and (b) peak RAM usage for
the Direct and low-rank solvers with
$\epsilon_{\text{rel}}=10^{-1},10^{-2},10^{-3}$ vs. number of radiation cavity
facets.

As mentioned previously, we were unable to run the Direct solver on the Level
6 and 7 meshes due to the RAM required to form the explicit cavity Jacobian
inverse, and this rapid memory growth is depicted in Fig. 11(b). Moreover, the
dashed line in Fig. 11 indicates extrapolated Direct solve simulation times
for the Level 6 and 7 meshes, based on a curve fit of the preceding data
points. For the Level 7 mesh with $\epsilon_{\text{rel}}=10^{-1}$, the low-
rank solver’s simulation time is approximately 15.6 minutes while the
estimated time for the direct solver (on a computer with $\geq$ 40 GB RAM) is
approximately 4.9 hours.

To conclude this discussion of results, we present in Fig. 12 a plot of the
speed-up of the low-rank solver with respect to the direct solver for all
meshes and $\epsilon_{\text{rel}}$ values tested. Speed-up here is defined
simply as the ratio of the Direct solver time to the low-rank solver time, and
in the case of the Level 6 and 7 meshes (dashed lines in the Figure) where the
Direct solve time is not available, we use an estimated time based on the
curve fit discussed previously.

Figure 12: Speed-up relative to Direct solve vs. $\epsilon_{\text{rel}}$ for
different mesh refinement levels. Speed-up for the Level 6 and 7 cases (dashed
lines) are based on extrapolated Direct Solve solution times since the compute
node had insufficient memory to complete the Direct Solve in those cases.

For the coarsest meshes (Levels 1 and 2), the problem size is too small to
realize a large speed-up with respect to the Direct solver; the time spent in
“other” parts of the code such as I/O and sparse matrix finite element
assembly are comparable in both solvers in these cases. On the finest meshes
(Levels 6 and 7), the (estimated) speed-up over the Direct solver is
approximately 17.3$\times$ and 21.3$\times$, respectively, but we note that
the speed-up decreases quickly with decreasing $\epsilon_{\text{rel}}$ (i.e.
increasing accuracy), so a judicious choice of $\epsilon_{\text{rel}}$ is
important for ensuring good performance.

##  7 Conclusions

In this work we presented a new approach to the finite element-based
computation of nonlinear transient heat transfer with cavity radiation.
Conventional finite element formulations for this problem employ dense matrix-
based approaches for the cavity terms, which scale poorly as the number of
cavity facets becomes large (e.g. $>$10000$$). To address the poor scaling of
the dense matrix formulation, we developed a novel hierarchical low-rank
approximation of the cavity terms. As noted in Section 1, this method has
similarities to [8], but also some important differences, including a low-rank
LU-factorization and back/forward substitution framework for efficiently
applying cavity terms at many nonlinear iterations and time steps, and an
efficient ACA method for constructing low-rank blocks.

Our numerical results demonstrated the accuracy, efficiency, and scalability
of the block low-rank framework that we proposed. Using the dense matrix
approach as the reference, we were able to obtain highly accurate results in
all cases, and with a speed-up of more than 20$\times$ for the test cases with
larger $n_{\rm facets}$. Perhaps most importantly, we demonstrated that for a
given $n_{\rm facets}$ the block low-rank approach has much lower memory
requirements than the dense approach. This means that, for a given hardware
configuration, the block low-rank approach can be applied to much larger
models (with higher $n_{\rm facets}$) than is possible with the conventional
“dense matrix” approach.

Based on our findings, it is clear that the computational advantage of the
block low-rank approach would only increase further for models with larger
$n_{\rm facets}$ than those considered here. Hence the methodology proposed in
this work is an enabler for solving large-scale cavity radiation problems that
arise regularly in industrial and scientific applications, without loss of
fidelity in the numerical approximation.

## References

  * [1] D. A. McQuarrie, J. D. Simon, Physical Chemistry: A molecular approach, University Science Books, Sausalito, California, 1997.
  * [2] P. B. Chinoy, D. A. Kaminski, S. K. Ghandhi, Effects of thermal radiation on momentum, heat, and mass transfer in a horizontal chemical vapor deposition reactor, Numerical Heat Transfer, Part A: Applications 19 (1) (1991) 85–100, https://doi.org/10.1080/10407789108944839.
  * [3] P. von Zedtwitz, A. Steinfeld, Steam-gasification of coal in a fluidized-bed/packed-bed reactor exposed to concentrated thermal radiation — Modeling and experimental validation, Industrial & Engineering Chemistry Research 44 (11) (2005) 3852–3861, https://doi.org/10.1021/ie050138w.
  * [4] A. Z’Graggen, A. Steinfeld, A two-phase reactor model for the steam-gasification of carbonaceous materials under concentrated thermal radiation, Chemical Engineering and Processing: Process Intensification 47 (4) (2008) 655–662, https://doi.org/10.1016/j.cep.2006.12.003.
  * [5] T. S. Barman, P. H. Hauschildt, F. Allard, Irradiated planets, The Astrophysical Journal 556 (2) (2001) 885–895, https://doi.org/10.1086/321610.
  * [6] M. Delbo, M. Mueller, J. P. Emery, B. Rozitis, M. T. Capria, Asteroid thermophysical modeling, in: P. Michel, F. E. Demeo, W. F. Bottke (Eds.), Asteroids, Vol. IV, The University of Arizona Press, 2005, pp. 107–128, https://doi.org/10.2458/azu_uapress_9780816532131-ch006.
  * [7] P. O. Hayne, O. Aharonson, N. Schörghofer, Micro cold traps on the moon, Nature Astronomy 5 (2) (2021) 169–175, https://doi.org/10.1038/s41550-020-1198-9.
  * [8] S. F. Potter, S. Bertone, N. Schörghofer, E. Mazarico, Fast hierarchical low-rank view factor matrices for thermal irradiance on planetary surfaces, ArXiv e-print, https://arxiv.org/abs/2209.07632 (2022).  arXiv:2209.07632.
  * [9] W. Hackbusch, Hierarchical Matrices: Algorithms and Analysis, Springer, Berlin, 2015\.
  * [10] S. A. Sauter, C. Schwab, Boundary Element Methods, Springer, Berlin, 2011.
  * [11] M. Bebendorf, Approximation of boundary element matrices, Numerische Mathematik 86 (4) (2000) 565–589, https://doi.org/10.1007/PL00005410.
  * [12] S. Börm, L. Grasedyck, W. Hackbusch, Hierarchical Matrices, Lecture Note No. 21, https://tinyurl.com/ycywrvjv (Jun. 2006).
  * [13] A. M. Franklin, An implementation of surface-to-surface, blackbody radiation heat transfer in a MOOSE application, Master’s thesis, Texas A&M University, https://hdl.handle.net/1969.1/192275 (Aug. 2020).
  * [14] J. P. Holman, Heat Transfer 6th Edition, McGraw-Hill Book Company, Singapore, 1986\.
  * [15] A. J. Chapman, Heat Transfer, Fourth Edition, Macmillan, New York, 1984.
  * [16] Abaqus/CAE User’s Manual, Version 6.11, Online (2011).
  * [17] R. S. Dembo, S. C. Eisenstat, T. Steihaug, Inexact Newton Methods, SIAM Journal on Numerical Analysis 19 (2) (1982) 400–408, https://doi.org/10.1137/0719025.
  * [18] D. A. Knoll, D. E. Keyes, Jacobian-free Newton–Krylov methods: A survey of approaches and applications, Journal of Computational Physics 193 (2) (2004) 357–397, https://doi.org/10.1016/j.jcp.2003.08.010.
  * [19] R. A. Dunlap, The Golden Ratio and Fibonacci Numbers, World Scientific, New Jersey, 1997.