# A Single Binary May Host Recurrent Thermonuclear Supernovae

Kaela J. Lee  Department of Astronomy and Theoretical Astrophysics Center,
University of California, Berkeley, CA 94720, USA  Ken J. Shen Department of
Astronomy and Theoretical Astrophysics Center, University of California,
Berkeley, CA 94720, USA Kaela J. Lee<EMAIL_ADDRESS>
###### Abstract

The most commonly accepted progenitor system for Type Iax supernovae (SNe Iax)
is the partial deflagration of a near-Chandrasekhar-mass white dwarf (WD)
accreting from a non-degenerate helium donor star, leaving a bound remnant
following the explosion. In this paper, we investigate whether the WD remant
can undergo multiple SNe during the system’s lifetime. We use Modules for
Experiments in Astrophysics (MESA) to evolve various single-degenerate
binaries to determine which could plausibly undergo multiple SNe Iax due to
multiple helium accretion phases. We also investigate the possibility for a
subsequent Type Ia SN after the formation of a double WD system. Our work
concludes that close binaries with relatively high-mass donors produce the
highest probability for several thermonuclear SNe.

††software: matplotlib (Hunter, 2007), MESA (Paxton et al., 2011, 2013, 2015,
2018, 2019)

##  1 Introduction

Type Iax supernovae (SNe Iax; also referred to as 2002cx-like based on their
prototypical member) are explosive transients similar in many ways to Type Ia
supernovae (SNe Ia). Both are believed to be thermonuclear explosions of white
dwarfs (WDs) and share spectral similarities such as a lack of hydrogen lines
and the presence of intermediate-mass and iron-group elements (Foley et al.,
2013; Jha, 2017). It is estimated that SNe Iax occur at approximately 5-30%
the rate of SNe Ia (Foley et al., 2013; Liu et al., 2015). Like SNe Ia, there
exists a correlation between luminosity and light curve shape for SNe Iax;
however, the correlation for SNe Iax is offset with increased scatter (Foley
et al., 2013). In comparison to SNe Ia, characteristics of SNe Iax include
lower peak luminosities, lower maximum velocities, lower ejecta masses (Foley
et al., 2013; Takaro et al., 2020), and unique late-time spectra (Li et al.,
2003; Foley et al., 2016). While SNe Ia are found in both early- and late-type
galaxies, all SNe Iax have been found in late-type galaxies with the exception
of SN 2008ge (Foley et al., 2010; Takaro et al., 2020). The presence of SNe
Iax in young host galaxies suggests that progenitor systems must produce and
explode WDs with short delay times between star formation and supernova. The
only thermonuclear explosion for which a potential companion star was observed
before the explosion is SN 2012Z, whose companion was a luminous blue star
with a mass of $\sim 2\,M_{\odot}$ (McCully et al., 2014).

Although the SNe Iax class is relatively new, there has already been
significant progress in narrowing down potential progenitor scenarios. The
current leading contender involves the partial deflagration of a WD with a
mass near the Chandrasekhar limit ($M_{\textrm{Ch}}$) accreting from a helium
donor star through Roche lobe overflow (but see, e.g., Kashyap et al. 2018 and
Bobrick et al. 2022 for alternative scenarios). Iben & Tutukov (1994) first
suggested that a degenerate star accreting from a helium donor could
potentially reach $M_{\textrm{Ch}}$ and produce a supernova; subsequent
calculations by Yoon & Langer (2003) showed that this evolution was indeed
plausible. The helium star donor is necessary to explain the short delay times
in a Chandrasekhar-mass scenario since the stable mass transfer rate of helium
is higher than for hydrogen (Ruiter et al., 2009; Jha, 2017). Deflagration
rather than detonation is thought to explain the lower kinetic energies and
greater chemical mixing of the ejecta (Foley et al., 2013). We also note that
thermonuclear SNe require the ignition of carbon in the core; by contrast,
off-center ignition results in stable burning and the formation of a massive
O/Ne WD or a neutron star (Shen et al., 2012; Schwab et al., 2016). The
deflagration is considered “partial” if the destruction of the WD is
incomplete. In fact, a study by Jha (2017) concluded that since the final
velocities of the ejecta are typically lower than the escape velocity from the
surface of the WD, some WD material likely remains bound. Additionally,
studies of three-dimensional hydrodynamic simulations of the deflagration
model also suggest a bound remnant is left behind (Jordan et al., 2012; Fink
et al., 2014; Long et al., 2014; Lach et al., 2022). This leads to the natural
question of whether the surviving remnant could begin a subsequent accretion
phase through Roche lobe overflow and undergo multiple SNe Iax.

In section 2, we discuss the methods used in this investigation to determine
the viability of multiple SNe Iax. In section 3, we discuss the results of our
simulations and describe the evolution of several systems in detail. We
conclude in section 4.

##  2 Description of models and calculations

We model the evolution of a variety of single-degenerate binary systems using
Modules for Experiments in Stellar Astrophysics (MESA) (Paxton et al., 2011,
2013, 2015, 2018, 2019). The donor star for each of our systems starts as a
helium main sequence star. The white dwarf (WD) accretor is modeled as a point
source. When the WD accretes up to the Chandrasekhar limit
($M_{\textrm{Ch}}$), we simulate a Type Iax event by decreasing the WD mass
and adjusting the binary separation based on the angular momentum loss due to
the ejected mass. We do not alter the mass or structure of the helium star,
which is a reasonable assumption for compact donors, but is less so for helium
giants due to the impact of the SN ejecta (Pan et al., 2012, 2013; Liu et al.,
2013; Bauer et al., 2019). We examine which systems begin mass accretion again
and bring the remnant back up to $M_{\textrm{Ch}}$. Additionally, each binary
eventually becomes a double WD system with the potential to undergo a SN Ia by
in-spiralling due to gravitational wave emission (Guillochon et al., 2010; Dan
et al., 2012). We vary two main initial parameters: the helium star’s initial
mass and the initial period of the system. We keep the initial mass of the
white dwarf fixed at $1.0\,M_{\odot}$ for this investigation. Notably, studies
by Wong & Schwab (2019) and Wong et al. (2021) performed similar simulations
using MESA to investigate thermonuclear explosions of helium star and WD
binaries. However, our investigation differs by focusing on the post-explosion
evolution of these systems.

We begin by using MESA to create helium main sequence star models with masses
of $1.0$, $1.5$, and $2.0\,M_{\odot}$. Details for using MESA to create helium
stars are located in the Appendix. We set the chemical composition to be 98%
4He and 2% 14N, by mass; we assume nitrogen is the endpoint of the CNO
isotopes after CNO burning. For the initial masses we consider, the helium
star expands after it exhausts helium in the core, potentially leading to mass
transfer due to Roche lobe overflow.

After performing the same stellar evolution with four nuclear reaction
networks within MESA (basic, approx21, co_burn, and mesa_201), we determine
that the mass fractions for the most abundant isotopes (4He, 12C, 14N, and
16O) are consistent across all networks within a few percent. Additionally,
these various networks do not noticeably alter the stellar evolution of the
models. We thus choose the basic network for simplicity.

We then evolve the previously created helium star models in a binary system
with a $1.0\,M_{\odot}$ point mass that is implied to be a C/O WD. It has been
suggested that more massive WDs that have undergone previous carbon-burning
may have a hybrid structure with an O/Ne mantle and a C/O core (Denissenkov et
al., 2013), which could still allow for central carbon ignition and an
eventual SN Iax (e.g., Meng & Podsiadlowski 2014; Kromer et al. 2015; Willcox
et al. 2016); however, for simplicity and because WDs with such structures may
not exist (Lecoanet et al., 2016; Brooks et al., 2017; Schwab & Garaud, 2019),
we implicitly assume an initially lower-mass C/O WD accretor. For each helium
star mass, we vary the initial period of the system. During our process, some
MESA settings differed between systems in order to allow them to run more
smoothly. Details are located in the Appendix.

In general, Roche lobe overflow can occur due to stellar expansion or
shrinking binary separation due to loss of angular momentum. In our
investigation, we turn off angular momentum loss from magnetic braking since
our helium main sequence stars do not have large convective zones at their
surface. We do allow angular momentum loss due to gravitational wave
radiation. However, since the main sequence lifetime for helium stars of this
mass is short (on the order $10^{7}$ years), the probability of encountering
systems with small initial binary separations such that mass transfer begins
during the main sequence phase is low.

The smallest initial period for each donor star mass is set such that the
Roche lobe radius is slightly larger than the helium star’s radius. Additional
systems were constructed with periods increasing by factors of approximately
2. The WD is approximated as a point mass with perfectly conservative mass
transfer, which allows transfer to occur both below and above the thermally
stable WD accretion rates (Brooks et al., 2016). The limitation of simulating
the WD as a point source is the inability to determine whether carbon ignition
begins centrally or off-center, which determines whether a supernova occurs.
By assuming conservative mass transfer, we maximize the total mass and orbital
angular momentum of the system. The results in this work can thus be regarded
as maximizing the number of SNe Iax. Although the studies by Wong & Schwab
(2019) and Wong et al. (2021) consider the stellar interior of the WD and
effects of non-conservative mass transfer, they do not investigate the
evolution of systems beyond a single thermonuclear explosion. We leave a more
in-depth analysis of recurrent thermonuclear SNe to future work.

If the WD reaches $M_{\textrm{Ch}}=1.4\,M_{\odot}$, we assume a SN Iax occurs
and examine how the dynamics of the system change due to the sudden mass
ejection from the WD. We assume no mass is stripped from the helium donor star
during the explosion. We also assume that the ejecta does not impart a recoil
kick on the bound remnant. Kicks up to several hundred ${\rm km}\,{\rm
s^{-1}}$ have been found by some WD deflagration studies for certain
conditions (Jordan et al., 2012; Lach et al., 2022) but not by others (Kromer
et al., 2013; Fink et al., 2014; Lach et al., 2022). The actual recoil kick
realized in nature may have a complicated dependence on initial conditions and
explosion properties (Lach et al., 2022); for simplicity, we neglect it in
this study.

We use the results from Hills (1983) for the ratio of pre- and post-SN semi-
major axes and eccentricities due to a mass ejection of $\Delta M$. Since
there is ample time for circularization prior to the explosion, our initial
binary systems have no eccentricity, and so we modify the Hills (1983) result
by setting the initial eccentricity $e_{0}=0$. The ratio of semi-major axes is
then:

$\frac{a}{a_{0}}=\frac{1}{2}\left(\frac{M_{0}-\Delta M}{\frac{M_{0}}{2}-\Delta
M}\right),$ (1)

where $a$ is the post-SN semi-major axis, $a_{0}$ is the pre-SN semi-major
axis, $M_{0}$ is the initial total mass of the system, and $\Delta M$ is the
mass ejected from the WD due to the SN Iax.

The final eccentricity from Hills (1983) is then:

$e=\left[1-\frac{\left(1-\frac{2\Delta M}{M_{0}}\right)}{\left(1-\frac{\Delta
M}{M_{0}}\right)^{2}}\right]^{\frac{1}{2}}.$ (2)

We now assume that after the SN Iax occurs, tidal circularization happens
rapidly because the helium star donor is nearly filling its Roche lobe. We
calculate a new semi-major axis assuming that the orbital angular momentum is
conserved as the eccentricity dissipates to zero via tides. The orbital
angular momentum is given by:

$L^{2}=\frac{G(M_{0}-\Delta M)\,\mu^{3}}{2E}(e^{2}-1)^{2},$ (3)

where $M_{0}-\Delta M$ is the total mass of the system after mass ejection,
$\mu$ is the reduced mass, and $E$ is the energy of the system.

We then find the new binary separation $a_{\rm recirc}$ is given by:

$\frac{a_{\rm recirc}}{a_{0}}=\frac{M_{0}}{M_{0}-\Delta M}.$ (4)

Using Equation 4, we reinitialize our MESA system with a WD mass of
$1.2\,M_{\odot}$, which simulates a mass ejection of $\Delta
M=0.2\,M_{\odot}$, and with a new semi-major axis $a_{\rm recirc}$. We repeat
this process every time the WD reaches $M_{\textrm{Ch}}$. The relatively small
$\Delta M=0.2\,M_{\odot}$ is on the low end of observationally inferred SN Iax
ejecta masses but matches those of low-luminosity SNe Iax such as SN 2008ha
(Foley et al., 2009). We note that such low ejecta masses are not found by
some theoretical studies (Jordan et al., 2012; Long et al., 2014) but are
produced by others (Fink et al., 2014; Lach et al., 2022). We choose this
ejecta mass to maximize the possibility of recurrent SNe Iax and leave an
exploration of a range of ejecta masses to future work.

Eventually, the systems reach a natural endpoint in MESA. The most common
endpoint is when the helium star ceases helium-burning, evolves into a second
white dwarf, and cools. Gravitational wave radiation will cause the resulting
double WD systems to merge and possibly result in a SN Ia. The merging times
of such systems are:

$t_{merge}=\frac{5c^{5}}{256G^{3}}\frac{a^{4}}{M_{1}M_{2}M_{tot}}.$ (5)

We count a SN Ia as occurring if two conditions are met. The first is if this
merging time is shorter than a Hubble time. The second is if the mass of the
primary WD is greater than $\sim 0.85\,M_{\odot}$ and the secondary WD is also
suitably massive (Dan et al., 2012; Guillochon et al., 2010; Shen et al.,
2017, 2018), which are always the case for the systems we consider.

Another possible endpoint is when the helium star expands into a giant and
causes the mass transfer rate to increase rapidly, which MESA is unable to
handle self-consistently, preventing us from simulating the full evolution of
these systems. Physically, when the mass transfer rate becomes very high, the
white dwarf is unable to stably accrete, causing the material to accumulate in
the circumbinary environment and yielding a common envelope (Iben & Livio,
1993; Ivanova et al., 2013). Frictional forces from the common envelope cause
the two stars to spiral towards each other. The difference in orbital energy
is used to eject the envelope. Thus, for these systems in our simulation, we
assume a common envelope with the mass of the unburned helium envelope of the
donor star is formed and ejected, which therefore alters the binary
separation. The remaining carbon core of the donor star is assumed to become a
degenerate carbon-oxygen (CO) WD. Due to the assumption in our simulation that
mass transfer onto the WD is stable at any accretion rate handled by MESA, we
expect to see fewer common envelopes in our simulation than in reality.

We assume that the entire helium envelope is ejected and find the binding
energy of that envelope:

$E_{bind}\sim-\frac{GM_{\rm He\,core}M_{\rm He\,env}}{R_{\rm He\,star}}.$ (6)

Figure 1: Summary of simulation results for various combinations of helium
star mass and initial period. Fixed parameters are an initial WD mass of
$1.0\,M_{\odot}$ and ejected mass $\delta m=0.2\,M_{\odot}$. Text to the right
of each symbol represents the final mass of the WD point mass in solar masses,
the final mass of the WD that the helium star evolves into in solar masses,
and the merging time of the double WD binary in years. Symbols are summarized
in the legend.

We then assume that in order to unbind the envelope, orbital energy is lost.
Thus, we can find the relationship between the initial and final binary
separations:

$\frac{E_{\rm bind}}{\alpha}=\frac{GM_{\rm He\,star}M_{\rm
WD}}{a_{i}}-\frac{GM_{\rm He\,core}M_{\rm WD}}{a_{f}},$ (7)

where $\alpha$ is an efficiency parameter that we assume to be 1, and we
assume the commonly used structure parameter, $\lambda$, is also 1.

After the common envelope ejection, the binary becomes a double WD system, and
we calculate the expected merging time using Equation 5.

The final possible endpoint we consider is the evolution of the helium star
through carbon-burning to become an oxygen/neon/magnesium (O/Ne) WD. The
evolution towards an O/Ne WD is signalled by the ignition of off-center
convective carbon-burning in the interior of the helium star. The calculation
of such models in MESA becomes very time-intensive, so we do not further
evolve these systems. However, since in each of our simulations at the stage
where carbon-burning occurs, mass transfer has already decreased to negligible
amounts, we assume that these carbon-burning stars in binary systems will
become O/Ne WDs with the same mass.

##  3 Results

We find that systems of recurrent SNe Iax are plausible. In particular, in
systems with greater helium star masses and smaller initial periods, the
potential for several SNe Iax increases.

However, our estimates are optimistic. We made several approximations that
increase the likelihood of multiple SNe Iax. First, we assume only
$0.2\,M_{\odot}$ of material is lost from the WD during each SNe Iax. Such a
small ejecta mass may be the case for very faint SNe Iax such as SN 2008ha
(Foley et al., 2009), but typical SNe Iax ejecta masses are likely several
times larger (Foley et al., 2013). Second, we assume that no mass is stripped
from the helium donor during the SN event. Third, we assume that all mass
transferred off the helium star remains on the WD. In nature, the range of
accretion rates leading to steady and stable helium-burning is only a factor
of a few (Brooks et al., 2016); outside of this regime, mass is lost from the
system via a combination of helium novae, optically thick winds, and binary-
driven mass loss.

Figure 1 demonstrates a summary of our findings. The data points represent
each MESA simulation with their symbol meanings explained in the legend. The
text to the right of each data point represents the final mass of the WD point
mass in solar masses, the final mass of the WD that results from the evolution
of the helium star in solar masses, and the merging time of the subsequent
double WD binary in years.

###  3.1 Systems with $1.0\,M_{\odot}$ helium star donors

Let us examine the possible scenarios for a $1.0\,M_{\odot}$ helium star. None
of these systems undergo a SNe Iax because the WD does not become massive
enough to explode prior to the cessation of mass transfer and the evolution of
the helium star to degeneracy. However, for systems with an initial period of
$\leq$ 0.5 days, the two resulting degenerate stars will merge within a Hubble
time due to the loss of orbital angular momentum by gravitational wave
emission, and possibly yield a SN Ia. Systems with an initial period $\geq$
800 days go through a common envelope phase. Since the ejected mass is only a
few hundredths of a solar mass in both cases we explored, the effect of the
common envelope ejection on the merging timescale is minimal.

Figure 2: Temperature vs. density profile at the time when mass transfer rate
is the highest for a system with a $1.0\,M_{\odot}$ helium star with an
initial period of 0.105 days. Zones where the net nuclear energy (produced
from nuclear reactions minus neutrino losses) is greater than
$1\,\mathrm{erg\,g^{-1}\,s^{-1}}$ are shown in orange. The dashed line shows
conditions for which the radiation pressure equals the ideal gas pressure,
assuming a chemical composition of fully ionized helium. Labelled black dots
represent mass fractions contained within a certain zone. Figure 3:
Temperature vs. density profile at the time when mass transfer rate is the
highest for a system with a $1.0\,M_{\odot}$ helium star with an initial
period of 400 days. In this profile, labelled mass fractions show a dense core
and massive envelope, which demonstrates that the star is in the giant phase.

For a $1.0\,M_{\odot}$ helium star, when the initial period is small, such as
in the case of an initial period of 0.105 days, mass transfer begins as soon
as the helium star begins to expand. Thus, the bulk of mass transfer occurs
along the Hertzsprung gap before the star becomes a giant. This is evident in
the temperature vs. density profile of the star at the time when the mass
transfer rate is highest, which is shown in Figure 2. In this profile, we see
that the star is no longer on the main sequence because the burning region is
outside the core. However, the star is not yet a giant because much of the
mass is confined within a relatively compact inner region of the star. The
density at an enclosed mass fraction of 0.99 is $3~{}\text{g/cm}^{3}$. In
contrast, Figure 3 demonstrates the temperature vs. density profile of a
$1.0\,M_{\odot}$ helium star that is a giant at the time when the mass
transfer rate is highest. This system began with a much larger initial period
of 400 days, and therefore requires the star to expand significantly before it
fills its Roche Lobe. The temperature vs. density profile shows a more massive
low density envelope. The density of this star at a mass fraction of 0.99 is
$4\times 10^{-8}~{}\text{g/cm}^{3}$.

Figure 4: Mass transfer rates for systems with an initial helium star mass of
$1.0\,M_{\odot}$ and differing initial periods. Dotted lines are the upper and
lower bounds of accretion that lead to steady helium-burning on the WD’s
surface according to Brooks et al. (2016).

Figure 4 demonstrates the mass transfer rate vs. WD mass for several systems
with an initial $1.0\,M_{\odot}$ helium star. We see that none of the WDs
reach the Chandrasekhar-mass. The black dotted lines represent the upper and
lower boundaries of accretion required for steady and stable helium-burning on
the WD according to Figure 3 from Brooks et al. (2016). All of the systems
shown in Figure 4 transfer a significant amount of mass outside these
boundaries, which demonstrates the limitations of our assumption that all mass
is stably transferred onto the WD point mass. The mass transfer rates run away
above rates of a few times $10^{-5}M_{\odot}$/yr for systems with an initial
period $\geq$ 800 days.

###  3.2 Systems with $1.5\,M_{\odot}$ helium star donors

Next, we examine systems with $1.5\,M_{\odot}$ helium stars. Systems with an
initial period of approximately $\leq$ 0.45 days can go through two sequential
SNe Iax. Within this group, systems with an initial period of approximately
$\leq$ 0.139 days can also merge within a Hubble time, which allows for the
potential for a SN Ia following the previous two SNe Iax. For systems with
periods between 1 and 400 days, we only find potential for one SN Iax based on
the previously described constraints. Notably, for systems with periods of 800
days, zero SNe are expected. Additionally, a common envelope phase occurs and
ejects approximately $0.5\,M_{\odot}$ of material, or approximately 34% of the
helium star’s original mass. This brings the system into closer orbit and
reduces the merging time slightly. We also expect the common envelope phase
will likely occur for systems with initial periods $\geq$ 800 days.

Figure 5: Temperature vs. density profile at the time when mass transfer rate
is the highest for a system with a $1.5\,M_{\odot}$ helium star with an
initial period of 800 days. The mass transfer rate for this system becomes
high enough to enter a common envelope phase. In this profile, labelled mass
fractions show a dense core and massive envelope, which demonstrates that the
star is in the giant phase. In contrast to Figure 3, a greater fraction of the
mass resides in this low density envelope, which is ejected as a common
envelope.

Figure 5 demonstrates the temperature vs. density profile for the system with
an initial helium star mass of $1.5\,M_{\odot}$ and an initial period of 800
days. This profile corresponds to the time when the mass transfer rate is
highest, which is directly prior to the end of MESA’s ability to calculate the
remaining binary evolution. We notice that at least 60% of the star’s mass is
in a compact region at the core of the star at a density $\geq 8\times
10^{5}~{}\text{g/cm}^{3}$ and a massive envelope contains at least 30% of the
mass at a density $\leq 10^{-5}~{}\text{g/cm}^{3}$. This convective envelope
contributes to runaway mass transfer and the formation of a common envelope
(Ge et al., 2010).

Figure 6: Mass transfer rates for systems with helium star mass of
$1.5\,M_{\odot}$ and differing initial periods. The mass transfer evolution
prior to the first and second SN Iax are shown as solid and dashed lines,
respectively, but any further evolution is neglected for visual simplicity.

Figure 6 demonstrates the mass transfer rate vs. WD mass for several systems
with a $1.5\,M_{\odot}$ helium star. The solid lines represent the mass
transfer prior to the first SN and the dashed lines represent mass transfer
prior to the second SN. If a second SN occurs, all further mass transfer is
neglected for visual simplicity. This figure demonstrates that most mass
transfer for these systems occurs above the rate at which transfer is stable.
The mass transfer rate for the system with an initial period of 800 days runs
away above a few times $10^{-4}M_{\odot}$/yr.

###  3.3 Systems with $2.0\,M_{\odot}$ helium star donors

Figure 7: Temperature vs. density profile for a system with a $2.0\,M_{\odot}$
helium star with an initial period of 2 days. Outer and inner burning
locations (shown in orange) correspond to helium-burning and carbon-burning
regions, respectively. Figure 8: Mass transfer rates for systems with helium
star mass of $2.0\,M_{\odot}$ and differing initial periods. The mass transfer
evolution beyond the first SN Iax is neglected for visual simplicity.

Finally, consider systems with $2.0\,M_{\odot}$ helium stars. Notably, for the
entire range of initial periods that we explored, the systems undergo at least
three Type Iax SNe. Systems with an initial period of 0.121 days have the
potential for four SNe Iax, and systems with an initial period $\leq$ 0.5 days
have the potential to merge within a Hubble time. Systems with an initial
period of 512 days will undergo a common envelope phase and eject
$0.15\,M_{\odot}$ of material. Additionally, helium stars that will eventually
evolve into WDs with a mass $\geq 1.04\,M_{\odot}$ ignite core carbon-burning.
The inward propagation of the carbon-burning wave by diffusion results in very
small simulation timesteps, and therefore we halt these calculations before
their natural endpoints. However, we can assume that carbon-burning will cause
the star to become an oxygen/neon/magnesium (O/Ne) WD. The temperature vs.
density profile for a star becoming an O/Ne WD is shown in Figure 7. The outer
and inner burning regions correspond to helium-burning and convective carbon-
burning regions, respectively. Much of the mass of this star is at a high
density. The density at an enclosed mass fraction of 0.99 is $4\times
10^{3}~{}\text{g/cm}^{3}$.

Figure 8 demonstrates the mass transfer rate vs. WD mass for several systems
with a $2.0\,M_{\odot}$ helium star. Only the mass transfer prior to the first
SN is shown for visual simplicity. From this figure, we can determine that
none of the mass transfer from $2.0\,M_{\odot}$ helium stars occurs in the
stable region prior to the first SN, which highlights the extreme optimism of
our assumptions. In reality, since mass transfer occurs above the stable
boundary, we expect the formation of an optically thick wind or a common
envelope, which we did not account for in systems with a period $<512$ days.

##  4 Conclusions

In this investigation, we used MESA to explore the potential for recurrent
thermonuclear SNe from binary systems containing a helium star donor and WD
accretor. Using our optimistic assumptions, we found that the possibility for
several SNe Iax exists in systems with high-mass helium stars and small
initial periods. We also find that systems with small initial periods are
likely to produce SNe Ia by forming double WD systems and merging within a
Hubble time, suggesting the possibility for both single- and double-degenerate
thermonuclear SNe produced by the same system.

However, we also find that the systems most likely to create recurrent SNe Iax
have the highest degree of error in assumptions in our simulation. Each of
these systems transfers mass at rates above the WD stable accretion rate set
by steady and stable helium-burning on the surface of the WD. Thus, our
assumption of conservative mass transfer onto the WD is very optimistic.
Likely, these systems would undergo a common envelope phase and eject much of
the unaccreted material, which would greatly decrease the chances for repeated
SNe Iax.

Our simulation is also optimistic due to the treatment of the WD as a point
mass, and resulting inability to determine whether carbon ignition begins
centrally or off-center. In cases of off-center carbon ignition, the WD burns
non-explosively and is converted into an O/Ne WD (Brooks et al., 2016) instead
of leading to a SN.

The possibility that binaries can host multiple thermonuclear SNe has obvious
implications for predictions of SN rates and galactic chemical enrichment.
This possibility also suggests polluted helium stars in WD binaries may exist
and may be observable throughout the galaxy. Furthermore, thermonuclear ash
from previous SNe that remains inside a surviving WD that subsequently
undergoes another SN may alter the resulting nucleosynthesis in interesting
ways. However, given the small parameter space that leads to multiple SNe
(restricted to short orbital periods and massive helium star donors) and the
overly optimistic assumptions employed in this work, it appears unlikely that
systems with multiple SNe play a significant role in nature. Still, given the
exciting ramifications that are possible, we recommend that future studies be
performed with more accurate physics regarding the stability and retention of
accreted helium.

This project was initiated with Saurabh Jha and Todd Thompson during a Scialog
Time Domain Astrophysics workshop. We further thank Todd Thompson for helpful
comments that improved this manuscript. Support for this work was provided by
NASA/ESA Hubble Space Telescope program #15918 and by the National Science
Foundation under Grant No. NSF PHY-1748958. The following inlist_project was
used to create helium star models in MESA:
&star_job

! do not begin with a pre-main sequence model
create_pre_main_sequence_model = .false.

! save a model once helium star is created
save_model_when_terminate = .false.
load_saved_model = .false.
save_model_number = 10
save_model_filename = ’He1.0.mod’

! display on-screen plots
pgstar_flag = .true.

! change chemical composition of the star
relax_initial_to_xaccrete = .true.

/ !end of star_job namelist
&controls

! set initial mass
initial_mass = 1.0 ! in Msun units

! options for energy conservation (see MESA V, Section 3)
use_dedt_form_of_energy_eqn = .true.

! set chemical composition of the star
accrete_given_mass_fractions = .true.
num_accretion_species = 2
accretion_species_id(1) = ’he4’
accretion_species_xa(1) = 0.98
accretion_species_id(2) = ’n14’
accretion_species_xa(2) = 0.02

/ ! end of controls namelist
When our binary systems encountered difficulties in MESA, a few settings were
toggled to allow the evolution to occur more smoothly. These included the
allowed variation in stellar structure between models (varcontrol_target), the
method of numerical calculation (min_mdot_for_implicit), and the efficiency of
energy transport (okay_to_reduce_gradT_excess). When some systems encountered
unstable periods, we briefly adjusted the convective mixing length
(mixing_length_alpha) in order to allow the system to progress. The following
inlist_project is an example of how we evolved our binaries:
&binary_job

inlist_names(1) = ’inlist1’
inlist_names(2) = ’inlist2’

evolve_both_stars = .false.

/ ! end of binary_job namelist
&binary_controls

terminate_if_initial_overflow = .false.

m2 = 1.0d0 ! point mass in Msun; reduced by 0.2 in each sequential run to
account for mass loss from SN
initial_period_in_days = 1
initial_separation_in_Rsuns = 100 ! use initial separation calculated after
SNe, rather than initial period

! angular momentum calculations
do_jdot_mb = .false.
do_jdot_gr = .true.
do_jdot_ml = .false.
do_jdot_ls = .false.

! set whether calculation is implicit or explicit
min_mdot_for_implicit = 1d99 ! default is 1d-16; when mdot in Msun/secyer >
min_mdot_for_implicit, then calculation is explicit

/ ! end of binary_controls namelist
The following inlist is an example of how we evolved our helium stars in the
binary systems:
&star_job

! begin from a previously created helium model
load_saved_model = .true.
saved_model_name = ’He1.0.mod’

! if the helium star reaches a mass such that the point mass is 1.4 Msun, then
terminate the run and save the helium star model
save_model_when_terminate = .true.
required_termination_code_string = ’star_mass_min_limit’
save_model_filename = ’He1.0_WD1.0_per1_SN1.mod’

/ ! end of star_job namelist
&controls

! minimum mass for helium star such that point mass reaches maximum of 1.4
Msun
star_mass_min_limit = 0.6d0

okay_to_reduce_gradT_excess = .false. ! turned to .true. for the entire run if
energy transport is causing run to evolve too slowly

varcontrol_target = 1d-4 ! increased to 3d-4 if varcontrol is the limiter of
evolution between steps

mixing_length_alpha = 2 ! increased to 4 or 8 for brief periods to move past
structural instability of helium star

/ ! end of controls namelist

## References

  * Bauer et al. (2019) Bauer, E. B., White, C. J., & Bildsten, L. 2019, ApJ, 887, 68, doi: 10.3847/1538-4357/ab4ea4
  * Bobrick et al. (2022) Bobrick, A., Zenati, Y., Perets, H. B., Davies, M. B., & Church, R. 2022, MNRAS, 510, 3758, doi: 10.1093/mnras/stab357410.48550/arXiv.2104.03415
  * Brooks et al. (2016) Brooks, J., Bildsten, L., Schwab, J., & Paxton, B. 2016, ApJ, 821, 28, doi: 10.3847/0004-637X/821/1/28
  * Brooks et al. (2017) Brooks, J., Schwab, J., Bildsten, L., Quataert, E., & Paxton, B. 2017, ApJ, 834, L9
  * Dan et al. (2012) Dan, M., Rosswog, S., Guillochon, J., & Ramirez-Ruiz, E. 2012, MNRAS, 422, 2417, doi: 10.1111/j.1365-2966.2012.20794.x
  * Denissenkov et al. (2013) Denissenkov, P. A., Herwig, F., Truran, J. W., & Paxton, B. 2013, ApJ, 772, 37
  * Fink et al. (2014) Fink, M., Kromer, M., Seitenzahl, I. R., et al. 2014, MNRAS, 438, 1762, doi: 10.1093/mnras/stt2315
  * Foley et al. (2016) Foley, R. J., Jha, S. W., Pan, Y.-C., et al. 2016, MNRAS, 461, 433, doi: 10.1093/mnras/stw1320
  * Foley et al. (2009) Foley, R. J., Chornock, R., Filippenko, A. V., et al. 2009, AJ, 138, 376, doi: 10.1088/0004-6256/138/2/376
  * Foley et al. (2010) Foley, R. J., Rest, A., Stritzinger, M., et al. 2010, AJ, 140, 1321, doi: 10.1088/0004-6256/140/5/1321
  * Foley et al. (2013) Foley, R. J., Challis, P. J., Chornock, R., et al. 2013, ApJ, 767, 57, doi: 10.1088/0004-637X/767/1/57
  * Ge et al. (2010) Ge, H., Hjellming, M. S., Webbink, R. F., Chen, X., & Han, Z. 2010, ApJ, 717, 724, doi: 10.1088/0004-637X/717/2/724
  * Guillochon et al. (2010) Guillochon, J., Dan, M., Ramirez-Ruiz, E., & Rosswog, S. 2010, ApJ, 709, L64, doi: 10.1088/2041-8205/709/1/L64
  * Hills (1983) Hills, J. G. 1983, ApJ, 267, 322, doi: 10.1086/160871
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iben & Livio (1993) Iben, Icko, J., & Livio, M. 1993, PASP, 105, 1373, doi: 10.1086/133321
  * Iben & Tutukov (1994) Iben, Icko, J., & Tutukov, A. V. 1994, ApJ, 431, 264, doi: 10.1086/174484
  * Ivanova et al. (2013) Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59, doi: 10.1007/s00159-013-0059-2
  * Jha (2017) Jha, S. W. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 375, doi: 10.1007/978-3-319-21846-5_42
  * Jordan et al. (2012) Jordan, IV, G. C., Perets, H. B., Fisher, R. T., & van Rossum, D. R. 2012, ApJ, 761, L23
  * Kashyap et al. (2018) Kashyap, R., Haque, T., Lorén-Aguilar, P., García-Berro, E., & Fisher, R. 2018, ApJ, 869, 140, doi: 10.3847/1538-4357/aaedb710.48550/arXiv.1811.00013
  * Kromer et al. (2013) Kromer, M., Fink, M., Stanishev, V., et al. 2013, MNRAS, 429, 2287
  * Kromer et al. (2015) Kromer, M., Ohlmann, S. T., Pakmor, R., et al. 2015, MNRAS, 450, 3045
  * Lach et al. (2022) Lach, F., Callan, F. P., Bubeck, D., et al. 2022, A&A, 658, A179, doi: 10.1051/0004-6361/202141453
  * Lecoanet et al. (2016) Lecoanet, D., Schwab, J., Quataert, E., et al. 2016, ApJ, 832, 71
  * Li et al. (2003) Li, W., Filippenko, A. V., Chornock, R., et al. 2003, PASP, 115, 453, doi: 10.1086/374200
  * Liu et al. (2015) Liu, Z.-W., Stancliffe, R. J., Abate, C., & Wang, B. 2015, ApJ, 808, 138, doi: 10.1088/0004-637X/808/2/138
  * Liu et al. (2013) Liu, Z.-W., Pakmor, R., Seitenzahl, I. R., et al. 2013, ApJ, 774, 37, doi: 10.1088/0004-637X/774/1/37
  * Long et al. (2014) Long, M., Jordan, George C., I., van Rossum, D. R., et al. 2014, ApJ, 789, 103, doi: 10.1088/0004-637X/789/2/103
  * McCully et al. (2014) McCully, C., Jha, S. W., Foley, R. J., et al. 2014, Nature, 512, 54, doi: 10.1038/nature13615
  * Meng & Podsiadlowski (2014) Meng, X., & Podsiadlowski, P. 2014, ApJ, 789, L45
  * Pan et al. (2012) Pan, K.-C., Ricker, P. M., & Taam, R. E. 2012, ApJ, 750, 151, doi: 10.1088/0004-637X/750/2/151
  * Pan et al. (2013) —. 2013, ApJ, 773, 49, doi: 10.1088/0004-637X/773/1/49
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Ruiter et al. (2009) Ruiter, A. J., Belczynski, K., & Fryer, C. 2009, ApJ, 699, 2026, doi: 10.1088/0004-637X/699/2/2026
  * Schwab et al. (2016) Schwab, J., Quataert, E., & Kasen, D. 2016, MNRAS, 463, 3461, doi: 10.1093/mnras/stw2249
  * Schwab & Garaud (2019) Schwab, J., & Garaud, P. 2019, ApJ, 876, 10
  * Shen et al. (2012) Shen, K. J., Bildsten, L., Kasen, D., & Quataert, E. 2012, ApJ, 748, 35, doi: 10.1088/0004-637X/748/1/35
  * Shen et al. (2018) Shen, K. J., Kasen, D., Miles, B. J., & Townsley, D. M. 2018, ApJ, 854, 52, doi: 10.3847/1538-4357/aaa8de
  * Shen et al. (2017) Shen, K. J., Toonen, S., & Graur, O. 2017, ApJ, 851, L50, doi: 10.3847/2041-8213/aaa015
  * Takaro et al. (2020) Takaro, T., Foley, R. J., McCully, C., et al. 2020, MNRAS, 493, 986, doi: 10.1093/mnras/staa294
  * Willcox et al. (2016) Willcox, D. E., Townsley, D. M., Calder, A. C., Denissenkov, P. A., & Herwig, F. 2016, ApJ, 832, 13
  * Wong & Schwab (2019) Wong, T. L. S., & Schwab, J. 2019, ApJ, 878, 100, doi: 10.3847/1538-4357/ab1b49
  * Wong et al. (2021) Wong, T. L. S., Schwab, J., & Götberg, Y. 2021, ApJ, 922, 241, doi: 10.3847/1538-4357/ac27ae
  * Yoon & Langer (2003) Yoon, S. C., & Langer, N. 2003, A&A, 412, L53, doi: 10.1051/0004-6361:20034607


# Accurate Benchmark Computations of the Polarizability Tensor for
Characterising Small Conducting Inclusions

A.A.S. Amad1, P.D. Ledger2, T. Betcke3, and D. Praetorius4 Manuscript received
2 June 2021. Corresponding author: A.A.S. Amad (email: alan.amad@gmail.com).
1Zienkiewicz Centre for Computational Engineering, College of Engineering,
Swansea University Bay Campus, Wales, UK  2School of Computing & Mathematics,
Keele University, Staffordshire, UK  3Centre for Inverse Problems, University
College London, Gower Street, London, UK  4TU Wien, Institute for Analysis and
Scientific Computing, Wiedner Hauptstr. 8-10 / E101 / 4, 1040 Wien, Austria

###### Abstract

The characterisation of small conducting inclusions with low conductivities in
an otherwise uniform background from low-frequency electrical field
measurements has important applications in medical imaging using electrical
impedance tomography as well as in geological imaging using electrical
resistivity tomography. It is known that such objects can be characterised by
a Póyla-Szegö polarizability tensor. The coefficients of this tensor for
different shapes and contrasts involves the solution of a transmission
problem, which in general must be done numerically. With a dictionary
classification algorithm in mind, we provide a series of accurate tensor
characterisations of objects that can be used as benchmarks for developers of
software tools for solving the transmission problem. To this end, we apply
adaptive the boundary element method to achieve our accurate solutions.

###### Index Terms:

Boundary element method; Adaptive mesh; Calderón’s problem; Conductivity
problem.

This work has been submitted to the IEEE for possible publication. Copyright
will be transferred without notice, after which this version will no longer be
accessible.

##  I Introduction

The characterisation of conducting objects with low conductivities in an
otherwise uniform background using low frequency electric fields has a wide
range of important applications including medical imaging and geology using
electrical impedance tomography (EIT) and electrical resistivity tomography,
respectively, e.g. [1, 2, 3] and references therein. In particular, being able
to characterise small objects is important for the related inverse problem
where the task is to be able to determine the shape and material contrasts of
hidden conducting inclusions from measurements of the induced voltage at the
boundary from induced current distributions [1]. This inverse problem was
first studied by Calderón [4]. There is a huge literature on the EIT inverse
problem and many different computational approaches have been proposed.
Solution approaches include total variational regularisation [5], topological
derivative [6, 7], machine learning [8], among others; please see [9, 2, 1, 3]
for an overview and the challenges in obtaining good solutions.

The characterisation of small three-dimensional objects also has important
applications in hydrodynamics, electrostatics, magnetostatics as well as in
models of dilute composition such as biological cell suspensions [10]. They
also appear in electromagnetic scattering [11, 12, 13] and it has been
proposed that electro-sensing fish use them for sensing their pray [14, 15].
Across this spectrum of different applications, it has been established that
asymptotic expansions of the perturbed field for small objects in the form
$B_{\alpha}=\alpha B+{\bm{z}}$, where $\alpha$ denotes the object size, $B$ a
unit sized object placed at the origin and ${\bm{z}}$ is the translation, has
a leading order term, which separates object position from its shape and
material contrast information. In EIT this takes the form

$\displaystyle(u_{\alpha}-u_{0})(\bm{x})=$
$\displaystyle\nabla_{\bm{x}}G(\bm{x},\bm{z})\cdot[{\mathcal{T}}(k,B)\nabla_{\bm{z}}u_{0}(\bm{z})]+O(\alpha^{4})$
(1)

as $\alpha\to 0$ where $u_{\alpha}$ denotes the field in the presence of the
object, $u_{0}$ denotes the background field,
$G(\bm{x},\bm{z})\mathrel{\mathop{\mathchar
58\relax}}=1/{(4\pi|\bm{x}-\bm{z}|)}$ is the free space Laplace Green’s
function, $k$ denotes the (conductivity) contrast,
${\mathcal{T}}(k,B)=\sum_{i,j=1}^{3}{\mathcal{T}}_{ij}{\bm{e}}_{i}\otimes{\bm{e}}_{j}$
is the appropriate polarizability tensor (also sometimes known as a
polarization tensor) for this problem, which coincides with the Pólya-Szegö
tensor (PST) and ${\bm{e}}_{i}$ the orthonormal coordinate directions. The PST
carries important geometric information, such as the volume of the inclusion
and the contrast [10]. For the EIT problem, Ammari and Kang [10] have obtained
a complete asymptotic expansion of the perturbed field leading to improved
object characterisations using generalised polarizability tensors. Different
polarizability tensor descriptions have been obtained for the eddy current
problem relevant for metal detection [16, 17, 18, 19, 20, 21] and also in
elasticity [10].

In the majority of EIT approaches, the hidden object’s position and its shape
is inferred by finding information about the conductivity distribution in the
imaging region. However, an alternative approach is offered by combining
dictionaries of object characterisations using polarizability tensors with
classification techniques, such as those in machine learning [22]. For such an
approach to be effective, a description of accurate polarizability tensor
coefficients is required. The purpose of this paper is to provide novel
solutions to a series of benchmark polarizability object characterisations so
that others can benchmark their solutions against ours. To achieve this, the
coefficients of the PST are computed by post processing the approximate
solutions to a transmission problem, which is governed by Laplace’s equation
set on an unbounded region surrounding the (scaled) object of interest [10].
Inspired by works [23, 14, 24], the transmission problem will be solved using
a boundary integral approach using the BEM++ library [25], leading to
approximate tensor coefficients for different objects. An adaptive boundary
integral scheme, which has already been shown to be effective for a capacitive
problem [26], is novelly applied to the computation of the PST to generate a
set of accurate benchmark solutions for different geometries.

The work proceeds as follows: In Section II, the mathematical model and
boundary integral/layer potential formulae to compute the PST are introduced.
Then, in Section III, a brief description of how BEM++ can be applied to
compute the PST is presented and, in Section IV, the application of an
existing adaptive algorithm to the computation of the PST is described.
Section V describes our benchmark geometries and presents numerical results
for the accurate computation of the PST using the adaptive algorithm. The
paper ends with some concluding remarks in Section VI.

##  II Small Object Characterisation

We recall a series of alternative, but equivalent formulations for computing
the coefficients of the real symmetric rank 2 PST
${\mathcal{T}}\mathrel{\mathop{\mathchar 58\relax}}={\mathcal{T}}(k,B)$ for an
object $B_{\alpha}=\alpha B+{\bm{z}}$.

###  II-A Transmission Problem Formulation

Given the solutions $\vartheta_{i}$, $i=1,2,3$, to the transmission problem

$\displaystyle\nabla^{2}\vartheta_{i}$ $\displaystyle=0$
$\displaystyle\hbox{in $B\bigcup B^{c}$},$ (2a)
$\displaystyle\left[\vartheta_{i}\right]_{\Gamma}$ $\displaystyle=0$
$\displaystyle\hbox{on $\Gamma$},$ (2b)
$\displaystyle{\mbox{\boldmath$n$}}\cdot\nabla\vartheta_{i}\Big{|}_{+}-{\mbox{\boldmath$n$}}\cdot
k\nabla\vartheta_{i}\Big{|}_{-}$
$\displaystyle={\mbox{\boldmath$n$}}\cdot\nabla_{\bm{\xi}}\xi_{i}$
$\displaystyle\hbox{on $\Gamma$},$ (2c) $\displaystyle\vartheta_{i}$
$\displaystyle\to 0$ $\displaystyle\hbox{as $|\bm{\xi}|\to\infty$},$ (2d)

where $B^{c}\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{R}^{3}\setminus\overline{B}$,
$\Gamma\mathrel{\mathop{\mathchar 58\relax}}=\partial B$,
$\left[\cdot\right]_{\Gamma}=\cdot|_{+}-\cdot|_{-}$ denotes the jump,
$\big{|}_{-}$ the interior trace and $\big{|}_{+}$ the exterior trace,
$\xi_{j}={\bm{e}}_{j}\cdot{\bm{\xi}}$ denotes the $j$-th the component of
${\bm{\xi}}$ measured from an origin in $B$ and differentiation is with
respect to ${\bm{\xi}}$, the coefficients of the PST follow from [10]

$\displaystyle{\mathcal{T}}_{ij}\mathrel{\mathop{\mathchar
58\relax}}=\alpha^{3}\left(\left(k-1\right)|B|\delta_{ij}+\left(k-1\right)^{2}\int_{\Gamma}{\mbox{\boldmath$n$}}\cdot\nabla\vartheta_{i}\Big{|}_{-}\xi_{j}\,\operatorname{d\\!}\bm{\xi}\right).$
(3)

In the above, $|B|$ is the volume of the object and $\delta_{ij}$ is the
Kronecker delta. By integration by parts, (3) can be shown to be equivalent to

$\displaystyle{\mathcal{T}}_{ij}\mathrel{\mathop{\mathchar 58\relax}}=$
$\displaystyle\alpha^{3}\left(\left(k-1\right)|B|\delta_{ij}-\left(k-1\right)^{2}\left[\int_{B^{c}}\nabla\vartheta_{i}\cdot\nabla\vartheta_{j}\operatorname{d\\!}{\bm{\xi}}\right.\right.$
$\displaystyle\left.\left.+\int_{B}k\nabla\vartheta_{i}\cdot\nabla\vartheta_{j}\operatorname{d\\!}{\bm{\xi}}\right]\right),$
(4)

which makes the symmetry of the PST clear. Its independence of ${\bm{z}}$ can
similarly be established [10].

###  II-B Boundary Integral (BI) Formulation

By following [27, equation 56-58] it is possible to rewrite (3), as follows

$\displaystyle{\mathcal{T}}_{ij}=\alpha^{3}\left(\dfrac{1}{r}|B|\delta_{ij}+\left(1-\dfrac{1}{r}\right)\int_{\Gamma}\psi_{i}n_{j}\,\operatorname{d\\!}\bm{\xi}\right),$
(5)

where $r=\dfrac{1}{k}$ and $\psi_{i}\mathrel{\mathop{\mathchar
58\relax}}=-k\vartheta_{i}$ can be expressed as [27]

$\psi_{i}(\bm{\xi})\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{r-1}\mathcal{S}\left[\left(\lambda
I+\mathcal{K}^{\ast}\right)^{-1}\left({\mbox{\boldmath$n$}}_{\bm{\nu}}\cdot\nabla_{\bm{\nu}}\nu_{i}\right)\right](\bm{\xi}),\quad\forall\bm{\xi}\in\mathbb{R}^{3}.$
(6)

Here, $\mathcal{S}$ is the single layer potential operator [10],

$\mathcal{S}[f]({\mbox{\boldmath$\xi$}})\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{|\mathscr{A}|}\int_{\Gamma}\frac{1}{|\bm{\xi}-\bm{\nu}|}f({\mbox{\boldmath$\nu$}})\,\operatorname{d\\!}{\mbox{\boldmath$\nu$}},$
(7)

$I$ is the identity operator and $\mathcal{K}^{*}$ is the adjoint double layer
boundary operator [10],

$\mathcal{K}^{\ast}[f](\bm{\xi})\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{|\mathscr{A}|}\mbox{
p.v.}\int_{\Gamma}\frac{(\bm{\xi}-\bm{\nu})\cdot{\mbox{\boldmath$n$}}_{\bm{\xi}}}{|\bm{\xi}-\bm{\nu}|^{3}}f(\bm{\nu})\,\operatorname{d\\!}\bm{\nu}.$
(8)

In addition, $|\mathscr{A}|$ denotes the area of the unit sphere in
$\mathbb{R}^{3}$, p.v. denotes the Cauchy principal value and
$|\bm{\xi}-\bm{\nu}|$ is the distance between $\bm{\xi}$ and $\bm{\nu}$, for
further details see [10].

###  II-C Layer Potential (LP) Formulation

In [10], the coefficients of ${\mathcal{T}}$ are shown to be equivalently
obtained as follows

${\mathcal{T}}_{ij}=\alpha^{3}\int_{\Gamma}\xi_{j}\phi_{i}\,\operatorname{d\\!}\bm{\xi},$
(9)

where $\phi_{i}$ is given by

$\phi_{i}(\bm{\xi})\mathrel{\mathop{\mathchar 58\relax}}=\left(\lambda
I+\mathcal{K}^{*}\right)^{-1}\left({\mbox{\boldmath$n$}}_{\bm{\nu}}\cdot\nabla_{\bm{\nu}}\nu_{i}\right)(\bm{\xi}),\quad\bm{\xi}\in\Gamma,$
(10)

and $\lambda\mathrel{\mathop{\mathchar 58\relax}}=\dfrac{k+1}{2(k-1)}$. Note
that, for exact computations, we have
$\psi_{i}(\bm{\xi})=\dfrac{1}{r-1}\mathcal{S}\left[\phi_{i}\right]({\bm{\xi}})$.

###  II-D Weighted (Alternative) Formulation

By taking a weighted average of (5) and (9) we get

$\displaystyle{\mathcal{T}}_{ij}=\alpha^{3}\left(\dfrac{\beta}{r}|B|\delta_{ij}+\int_{\Gamma}\left(\beta\left(1-\dfrac{1}{r}\right)\psi_{j}n_{i}+\left(1-\beta\right)\phi_{j}\xi_{i}\right)\,\operatorname{d\\!}{\mbox{\boldmath$\xi$}}\right),$
(11)

with $0\leq\beta\leq 1$. Note that, for $\beta=0$ we obtain the formula (9),
whereas for $\beta=1$ we have (5).

##  III Computation of the Pólya-Szegö Tensor Coefficients

In this work, we use the BEM++ library [25] to compute the coefficients of the
PST. BEM++ is a C++ library with Python bindings for the solution of boundary
integral equations using Galerkin discretisations [25, 28]. Given the
similarity of the transmission problem (2) and the boundary integral solution
formulations (6) and (10) to the capacitive problem already successfully
considered in [26], we will focus on a dual mesh discretisation although a
primal discretisation is also possible and may also offer further efficencies.
In practice, given a surface mesh discretisation of $B$ of size $h$ we solve
for $\phi_{i}^{h}$ and obtain
$\psi_{i}^{h}({\bm{\xi}})=\dfrac{1}{r-1}\mathcal{S}\left[\phi_{i}^{h}\right]({\bm{\xi}})$.
Following the numerical solution, the discrete approximation
${\mathcal{T}}_{ij}^{h}$ to ${\mathcal{T}}_{ij}$ then follows from a simple
postprocessing of the solutions using (11).

Although the coefficients of the true tensor ${\mathcal{T}}$ are symmetric,
the coefficients of its numerical approximation ${\mathcal{T}}^{h}$ using
boundary elements are not guaranteed to have this property. Thus, to ensure
the symmetry of the computed PST, we will symmetrise it as follows:

${\mathcal{T}}^{h,sym}_{ij}\mathrel{\mathop{\mathchar
58\relax}}=\dfrac{1}{2}\left({\mathcal{T}}^{h}_{ij}+{\mathcal{T}}^{h}_{ji}\right).$
(12)

Henceforth, when presenting results, we will only consider the symmetrised
tensor and drop the superscript $sym$. In the next section we illustrate the
convergence behaviour of the numerical scheme for uniform grid refinement for
a PST with smooth analytical solution.

###  III-A Numerical Convergence Study for Ellipsoidal Objects

The convergence of the approximate tensor $\mathcal{T}^{h}$ to the true PST
$\mathcal{T}$ obtained on a sequence of uniformly meshes for ellipsoidal
objects is investigated. In these experiments, we choose $B_{\alpha}=\alpha B$
and define $B$ by

$\dfrac{\xi_{1}^{2}}{a^{2}}+\dfrac{\xi_{2}^{2}}{b^{2}}+\dfrac{\xi_{3}^{2}}{c^{2}}=1,\quad
0<c\leq b\leq a.$ (13)

We choose $\alpha=0.01\,\mathrm{m}$ and $k=10$ and consider the cases
$a=b=c=1$ so that $B$ degenerates to a sphere and consider the ellipsoid
defined by $a=1$, $b=0.7$, $c=0.5$. In Figures 1(a) and 1(b) we show the
behaviour of $\|\mathcal{T}^{Exact}-\mathcal{T}^{h}\|/\|\mathcal{T}^{Exact}\|$
with uniform mesh refinement against the number of degrees of freedom
$\\#(ndof)$ for different values of $0\leq\beta\leq 1$ and both primal and
dual mesh discretisations. Here,
$\|\mathcal{T}^{Exact}\|=\left(\sum_{i,j=1}^{3}({\mathcal{T}}_{ij}^{Exact)})^{2}\right)^{1/2}$
is the Frobenius norm and ${\mathcal{T}}^{Exact}$ is the exact solution for
these objects [10, page 83]. Specifically, in Figure 1(a), the primal mesh is
initially chosen as surface mesh of $232$ surface triangles and then uniformly
refined $6$ times until a mesh of $31\,616$ surface triangles is obtained.
While in Figure 1(a), the primal mesh is initially chosen as surface mesh of
$48$ surface triangles and then uniformly refined $7$ times until a mesh of
$65\,566$ surface triangles is obtained. The corresponding dual meshes are
obtained from the primal meshes using the procedure described in [26]. Except
for a couple of exceptions (eg $\beta=0.4$ shown in 1(a) for the dual mesh),
the convergence behaviour for both primal and dual meshes is similar and
practically identical for different values of $\beta$. We attribute the
exceptions to some super-convergence behaviour in exceptional cases. The fixed
convergence rate of the PST under uniform mesh refinement is not surprising
given the results of Appendix A, which establishes that the convergence
behaviour of the coefficients of the PST is the same as the convergence of the
solutions to the transmission problem under uniform mesh refinement.

(a) Sphere

(b) Ellipsoid

Figure 1: Relative error between the study of the polarization tensor
${\mathcal{T}}^{h}$ (11) computed on primal mesh $\mathcal{M}$ and dual mesh
$\mathcal{M}^{dual}$ for $\beta=0.0,\ldots,1.0$ and compared with the exact
solution in [10, page 83] (a) sphere with radius $a=b=c=1$ and (b) ellipsoid
with radii $a=1$, $b=0.7$, $c=0.5$.

##  IV Adaptive Computation of PST Coefficients using BEM++

Apart from limited smooth geometries, such as the ellipsoid considered above,
the PST does not have an analytical solution and instead numerical
approximations are essential to find its coefficients. Furthermore, for
realistic objects with edges and corners, refinement of the mesh towards these
features is expected to be desirable to accurately capture the solution to the
transmission problem (2), and, hence, the PST coefficents. To efficiently and
automatically guide this mesh refinement, we employ an existing adaptive mesh
method based on [29, 30, 26] to enhance the numerical computation
${\mathcal{T}}^{h}$ for objects with edges and corners. Alternatively, other
adaptive schemes could be used e.g. [31, 32]. In particular, these schemes
allow elements to be selected for refinement according to the magnitude of the
contributions to the error indicator with those with the largest contributions
being refined. The solution on the refined mesh leading to a better
approximation for an economical increase in computational effort.

The approach used here closely follows the successful scheme presented for a
related problem in [26] for which an existing code was available that could be
suitably modified. The basic procedure could be summarised as follows: we
solve for $\phi_{i}^{h}$ on a dual grid discretisation of ${\mathcal{M}}_{0}$
and use a Zienkiewicz-Zhu type error estimator
$\eta_{i}\approx\|\phi_{i}^{h}-\phi_{i}\|_{E}$, where $\|\cdot\|_{E}$ denotes
the energy norm. This is obtained by interpolating solutions on the primal
grid; see [26] for further details. We consider two different adaptive
schemes, which are driven by flagging elements for refinement according to the
elemental contributions to $\eta=\max\eta_{i}$ or
$\eta=\sum_{i=1}^{3}\eta_{i}$ leading to a new mesh ${\mathcal{M}}_{1}$. Then,
by computing $\phi_{i}^{h}$ on ${\mathcal{M}}_{1}$, and repeating the error
estimation process, a new grid ${\mathcal{M}}_{2}$ is obtained. By repeating
this process, the sequence of grids ${\mathcal{M}}_{\ell}$, $\ell=1,\ldots$ is
generated. Note that the Döfler criterion [26] is used to mark the elements
that will be refined. The refinements are controlled by $\theta$, where
$\theta=1$ coincides with uniform mesh-refinement. To maintain a balance
between the number of refinements and the control of the elements to be
refined, the Döfler parameter is chosen as $0.4\leq\theta\leq 0.6$.

##  V Numerical Benchmarks of PST Coefficients using Adaptive BEM

In this section, we present benchmark computations of the PST obtained using
the aforementioned adaptive scheme. We first define the benchmark geometries
that will be considered. Then, we consider the performance in detail of the
L-shape benchmark geometry before presenting benchmark solutions for all
geometries.

###  V-A Benchmark Geometries

TABLE I: Benchmark objects.  L-shape | Cube
---|---
|
Tetrahedron | Key
|

####  V-A1 L-shape

For the L-shape object, we consider $B_{\alpha}\mathrel{\mathop{\mathchar
58\relax}}=\alpha B$ with $\alpha=0.01\,\mathrm{m}$, $B=B_{1}\bigcup B_{2}$,
where

$\displaystyle B_{1}=$ $\displaystyle[0.0,7.8]\times[0.0,2.0]\times[0.0,1.5],$
$\displaystyle B_{2}=$ $\displaystyle[0.0,2.2]\times[2.0,5.6]\times[0.0,1.5].$

The reflectional symmetries of this object mean that ${\mathcal{T}}$ has 4
independent coefficients.

####  V-A2 Cube

The second object considered is where $B_{\alpha}\mathrel{\mathop{\mathchar
58\relax}}=\alpha B$ is a cube, with $\alpha=0.01\,\mathrm{m}$ and $B$ is a
unitary cube. The rotational and reflectional symmetries of this object mean
that ${\mathcal{T}}$ is a multiple of identity and has 1 independent
coefficient.

####  V-A3 Tetrahedron

The third object, we consider $B_{\alpha}\mathrel{\mathop{\mathchar
58\relax}}=\alpha B$ to be a tetrahedron, with $\alpha=0.01\,\mathrm{m}$ and
the vertices of $B$ are chosen to be at the locations

${\mbox{\boldmath$v$}}_{1}=\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix},\
{\mbox{\boldmath$v$}}_{2}=\begin{pmatrix}7\\\ 0\\\ 0\end{pmatrix},\
{\mbox{\boldmath$v$}}_{3}=\begin{pmatrix}5.5\\\ 4.6\\\ 0.0\end{pmatrix},\
{\mbox{\boldmath$v$}}_{4}=\begin{pmatrix}3.3\\\ 2.0\\\ 5.0\end{pmatrix}.$ (14)

This object has no symmetries and, hence, ${\mathcal{T}}$ has 6 independent
coefficients.

####  V-A4 Key

The last object, we will consider is a key
$B_{\alpha}\mathrel{\mathop{\mathchar 58\relax}}=\alpha B$, as illustrated in
Figure 2, with $\alpha=0.001\,\mathrm{m}$, which has additional geometrical
complexities. We have divided the key in three parts, bow, shoulder and blade,

$B=\mathrm{bow}\bigcup\mathrm{shoulder}\bigcup\mathrm{blade},$

where

$\mathrm{bow}=[-7.0,7.0]\times[3.0,15.0]\times[-1.25,1.25],$
$\mathrm{shoulder}=[-4.0,4.0]\times[0.0,4.0]\times[-1.25,1.25],\qquad$

and the maximum dimensions of the blade are

$\mathrm{blade}=[-3.25,3.25]\times[-19.0,0.0]\times[-1.25,1.25].$

To enforce that the shape is more realistic, some cuts have been made in the
blade, as presented in Figure 2. This object has no symmetries and, hence,
${\mathcal{T}}$ has 6 independent coefficients.

Figure 2: Key: Sketch of a simplified key showing the physical object
$B_{\alpha}$.

###  V-B Detailed study of the L-Shape Geometry

As an illustration, we provide a detailed convergence study of the L-shape
geometry. In this section, we provide a convergence study as well as images of
the resulting adaptive surface meshes for this object. Throughout this
section, the surface of $B$ is discretised with an initial mesh
${\mathcal{M}}_{0}$ of $64$ surface triangles and in absence of an analytical
solution, we consider $\mathcal{T}^{h}_{\mathrm{Fixed}}$ obtained on a fine
mesh ${\mathcal{M}}_{\rm{Fixed}}$ with $262\,144$ surface triangles to be
indistinguishable from the exact.

####  V-B1 Convergence Study for Varying $\beta$

We consider the accuracy of ${\mathcal{T}}^{h}$ obtained from (11) when
weights $\beta=0.0,0.1,\ldots,1.0$ are considered. For a contrast $k=10$, the
convergence study of the relative error of the PST given by

$\mathcal{E}\mathrel{\mathop{\mathchar
58\relax}}=\dfrac{\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}-\mathcal{T}^{h}\right\|}{\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}\right\|},$
(15)

is presented in Figures 3-4. In Figure 3, the convergence behaviour obtained
when the error estimator $\eta=\sum_{i=1}^{3}\eta_{i}$ is used to drive the
adaptive process for the cases of $\theta=0.4$, $\theta=0.5$ and $\theta=0.6$.
Analogously, Figure 4 shows the corresponding convergence behaviour obtained
using the error estimator $\eta=\max\eta_{i}$.

(a) $\theta=0.4$

(b) $\theta=0.5$

(c) $\theta=0.6$

Figure 3: L-shape: Convergence study of
${\mathcal{E}}=\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}-\mathcal{T}^{h}\right\|/\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}\right\|$,
using the $\eta=\sum_{i=1}^{3}\eta_{i}$ estimator for $\beta=0.0,\ldots,1.0$
(a) $\theta=0.4$ (b) $\theta=0.5$ (c) $\theta=0.6$.

(a) $\theta=0.4$

(b) $\theta=0.5$

(c) $\theta=0.6$

Figure 4: L-shape: Convergence study of
${\mathcal{E}}=\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}-\mathcal{T}^{h}\right\|/\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}\right\|$,
using the $\eta=\max\eta_{i}$ estimator for $\beta=0.0,\ldots,1.0$ (a)
$\theta=0.4$ (b) $\theta=0.5$ (c) $\theta=0.6$.

Comparing the results shown in Figures 3 and 4, we first observe that, unlike
the results presented in Figures 1(a) and 1(b), changing $\beta$ can improve
the convergence rate when combined with the adaptive scheme. In particular, it
can be observed that for both error estimators $\eta=\sum_{i=1}^{3}\eta_{i}$
and $\eta=\max\eta_{i}$, the best performance is obtained by choosing
$\beta=0.4$, which also achieves small error values. Similar convergence rates
were found for $\beta=0.5,0.6$, but with larger errors. We have repeated this
study for a range of other objects and have found that in all cases the best
performance for the parameters considered was $0.3\leq\beta\leq 0.4$, for
$k=10$.

####  V-B2 Convergence Study for Varying Contrasts $k$

We have undertaken adaptive mesh approximations for different contrasts
$k=0.1,0.2,1.5,10,15,100,150$ and summarise the $\beta$ for which the best
results were obtained using $\theta=0.6$ in Table II.

TABLE II: L-shape: Contrast $k$ and best $\beta$.  Contrast $k$ | Best weight $\beta$ | Measurement option
---|---|---
$0.1$ | $0.0$ | max
$0.2$ | $0.4$ | max
$1.5$ | $0.6$ | max
$10.0$ | $0.4$ | max
$15.0$ | $0.3$ | sum
$100.0$ | $0.1$ | sum
$150.0$ | $0.1$ | sum

If adaptive mesh approximations for different contrasts
$k=0.1,0.2,1.5,10,15,100,150$ are compared with uniform mesh refinements, we
find that the adaptive algorithm offers a superior performance over uniform
refinement for $k=0.1,0.2,1.5,10,15$ with the best performance for
$k=0.1,0.2$. For $k=100,150$, the adaptive approach does not show an
improvement over uniform refinement. We illustrate the best performing cases
in Figure 5(a) and show how error estimators $\eta=\sum_{i=1}^{3}\eta_{i}$ and
$\eta=\max\eta_{i}$ converge for this case and different values of $\beta$ in
Figure 5(b). In these experiments, a maximum of $25\,000$ surface triangles
was generated by the adaptive algorithm while the uniform refinement led to
$65\,536$ surface triangles.

(a) Comparison of the relative error for $k=0.2$

(b) ZZ-error estimator for $k=0.2$

Figure 5: L-shape: (a) Comparative study between uniform refinement and
adaptive mesh solutions for
${\mathcal{E}}=\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}-\mathcal{T}^{h}\right\|/\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}\right\|$
for $\beta=0.4$, $\theta=0.6$ and $k=0.2$ (b) Convergence study of the ZZ-
error estimator, using $\eta=\sum_{i=1}^{3}\eta_{i}$ and $\eta=\max\eta_{i}$
for $\beta=0.0,\ldots,1.0$, $\theta=0.6$ and $k=0.2$.

####  V-B3 Illustrations of Adaptive Meshes

For the case of $k=10$, the initial grid ${\mathcal{M}}_{0}$ and exemplar
adaptive grids ${\mathcal{M}}_{\ell}$, $\ell=2,4,6$ obtained for $\beta=0.6$,
$\theta=0.6$ using $\eta=\max\eta_{i}$ are shown in Figures 6(a)-6(d),
respectively, which have $832$ to $6\,022$ surface triangles. Note that the
meshes generated by the algorithm are more refined along the edges and
corners, this trend is expected, since the solution $\phi_{i,h}$ is less
regular at these locations and elements in the vicinity of edges and corners
contribute most to the error estimate $\eta_{i}$. The distribution of the
refined elements are similar for other objects.

(a) Max - $\ell=0$

(b) Max - $\ell=2$

(c) Max - $\ell=4$

(d) Max - $\ell=6$

Figure 6: L-shape: Adaptive mesh generated using the ZZ-type error estimator
using max direction measurement $\eta_{\ell}(\vartheta_{i}^{dual},T)$ for
$\mathcal{M}^{dual}_{\ell}$ (a) $\ell=0$ (b) $\ell=2$ (c) $\ell=4$ (d)
$\ell=6$.

###  V-C Reference Polarization Tensors and Benchmark Solutions

The benchmark solutions obtained using the adaptive mesh algorithm for each of
the benchmark geometries in Section V-A are now presented for the case of
$k=10$. We state ${\mathcal{E}}$ for each of the tensors resulting from the
final adaptive mesh. In addition, by introducing the splitting

${\mathcal{T}}={\mathcal{O}}+{\mathcal{D}},$ (16)

where ${\mathcal{D}}$ is the diagonal part of ${\mathcal{T}}$ and
${\mathcal{O}}$ the off-diagonal part, we provide a measure of the error in
the off-diagonal part of the tensor as

$\mathcal{E}^{\mathrm{off}}\mathrel{\mathop{\mathchar
58\relax}}=\dfrac{\left\|\mathcal{O}^{h}_{\mathrm{Fixed}}-\mathcal{O}^{h}\right\|}{\left\|\mathcal{T}^{h}_{\mathrm{Fixed}}\right\|}.$
(17)

for $\mathcal{T}^{h}$ obtained on the final adaptive mesh. The results are
summarised in Table III for each benchmark object, using $\eta=\max\eta_{i}$
and $\theta=0.6$.

TABLE III: Relative error and off-diagonal error for the best convergence result.  Object | Best $\beta$ | Maximum $\\#ndof$ | Number of elements on $\mathcal{M}$ | $\mathcal{E}$ | $\mathcal{E}^{\mathrm{off}}$
---|---|---|---|---|---
L-shape | $0.4$ | $11\,681$ | $23\,358$ | $0.0213\%$ | 0.0019%
Cube | $0.4$ | $11\,522$ | $23\,040$ | $0.0013\%$ | 0.0001%
Tetrahedron | $0.5$ | $12\,078$ | $24\,152$ | $0.0170\%$ | 0.0140%
Key | $0.3$ | $17\,179$ | $34\,354$ | $0.0132\%$ | 0.0051%

In the following sections, we present for each benchmark object the best
$\mathcal{T}^{h}$ solution obtained from the adaptive mesh method with the PST
coefficients arranged as a $3\times 3$ matrix in each case. Note that the
solutions are restricted to a mesh with a maximum of $25\,000$ surface
triangles, the only exception being the key object that the mesh can contain
up to $40\,000$ surface triangles. We believe that at least the first four
significant figures in our benchmark tensor characterisations are accurate.

####  V-C1 L-shape: Best Approximation

The best approximation for the L-shape object is obtained for $\beta=0.4$,
$\theta=0.6$ and $\eta=\max\eta_{i}$ and has a final adapted mesh with
$23\,358$ surface triangles. This approximate tensor $\mathcal{T}^{h}$ is
shown in (21)

$\displaystyle\mathcal{T}^{h}=10^{-4}$
$\displaystyle\left(\begin{array}[]{cc}\phantom{-}1.54683606\phantom{\times
10^{-6}}&-0.15091553\phantom{\times 10^{-6}}\\\ -0.15091553\phantom{\times
10^{-6}}&\phantom{-}1.16202896\phantom{\times 10^{-6}}\\\ -6.99198691\times
10^{-6}&\phantom{-}3.33569676\times 10^{-6}\end{array}\right.$ (21)
$\displaystyle\qquad\qquad\qquad\left.\begin{array}[]{c}-6.99198691\times
10^{-6}\\\ \phantom{-}3.33569676\times 10^{-6}\\\
\phantom{-}0.57033981\phantom{\times 10^{-6}}\end{array}\right).$ (25)

####  V-C2 Cube: Best Approximation

For the cube object the best approximation is obtained by choosing
$\beta=0.4$, $\theta=0.6$ and $\eta=\max\eta_{i}$ on a mesh with $23\,040$
surface triangles. This approximate tensor $\mathcal{T}^{h}$ is shown in (29)

$\displaystyle\mathcal{T}^{h}=10^{-6}$
$\displaystyle\left(\begin{array}[]{cc}\phantom{-}2.51110996\phantom{\times
10^{-6}}&\phantom{-}1.92945670\times 10^{-7}\\\ \phantom{-}1.92945670\times
10^{-7}&\phantom{-}2.51111340\phantom{\times 10^{-6}}\\\ -1.59571932\times
10^{-6}&-4.11272795\times 10^{-6}\end{array}\right.$ (29)
$\displaystyle\qquad\qquad\left.\begin{array}[]{c}-1.59571932\times 10^{-6}\\\
-4.11272795\times 10^{-6}\\\ \phantom{-}2.51110887\phantom{\times
10^{-6}}\end{array}\right).$ (33)

####  V-C3 Tetrahedron: Best Approximation

For the tetrahedron the best approximation is obtained by choosing
$\beta=0.5$, $\theta=0.6$ and $\eta=\max\eta_{i}$ on a mesh with $24\,152$
surface triangles. This approximate tensor $\mathcal{T}^{h}$ is shown in (34)

$\mathcal{T}^{h}=10^{-5}\begin{pmatrix}\phantom{-}9.30682676&1.12847255&-0.76359289\\\
\phantom{-}1.12847255&6.83952305&\phantom{-}0.43032199\\\
-0.76359289&0.43032199&\phantom{-}7.80618516\end{pmatrix}.$ (34)

####  V-C4 Key: Best Approximation

For the key the best approximation is obtained by choosing $\beta=0.3$,
$\theta=0.6$ and $\eta=\max\eta_{i}$ on a mesh with $34\,354$ surface
triangles. This approximate tensor $\mathcal{T}^{h}$ is shown in (38).

$\displaystyle\mathcal{T}^{h}=10^{-6}$
$\displaystyle\left(\begin{array}[]{cc}\phantom{-}2.66099087\phantom{\times
10^{-7}}&-2.40425167\times 10^{-2}\\\ -2.40425167\times
10^{-2}&\phantom{-}4.49014335\phantom{\times 10^{-5}}\\\
\phantom{-}1.94050498\times 10^{-5}&-3.78250864\times
10^{-6}\end{array}\right.$ (38)
$\displaystyle\qquad\qquad\left.\begin{array}[]{c}\phantom{-}1.94050498\times
10^{-5}\\\ -3.78250864\times 10^{-6}\\\ \phantom{-}0.96008445\phantom{\times
10^{-7}}\end{array}\right).$ (42)

##  VI Conclusions

In this work, we propose a series of benchmark solutions for the PST for
different objects. These benchmark objects were defined to compare the
accuracy of the tensor coefficients. Since there is no exact solution of
${\mathcal{T}}$ for these benchmark objects, we have assumed that the tensor
coefficients computed using a fine mesh ${\mathcal{M}}_{\rm{Fixed}}$ is
indistinguishable from the exact solution.

Several convergence studies have been presented to confirm that the adaptive
boundary element algorithm reduces the computational cost and obtains an
accurate approximation ${\mathcal{T}}^{h}$. In general, the results showed
greatest improvements for small contrasts $k$ and for larger values of $k$
were no worse than uniform refinement. We highlight the cases with contrasts
$0.1\leq k\leq 1.5$, for these cases we have obtained the best convergence
rate using less degree of freedom for the adaptive boundary element algorithm.
Finally, benchmark solutions for the objects are presented, as well as the
best approximation of relative error $\mathcal{E}$ (15) and the off-diagonal
error $\mathcal{E}^{\mathrm{off}}$ (17) are presented for each benchmark
object.

##  Appendix A Convergence of the Pólya-Szegö Tensor

In this section, a convergence result for the numerical approximation of
${\mathcal{T}}^{h}$ to ${\mathcal{T}}$ under uniform mesh refinement is
presented.

In line with the solution to related boundary integral approximations of
scalar Laplace transmission problems (e.g. [33, pg. 199, Proposition 4.1.31]),
we conjecture that the approximate solution $\phi_{i}^{h}$ to $\phi_{i}$
satisfies the following a-priori type error estimate under uniform mesh
refinement

$\left\|\phi_{i}-\phi_{i}^{h}\right\|_{L^{2}(\Gamma)}\leq
Ch^{s}\|\phi_{i}\|_{H^{s}(\Gamma)},0\leq s\leq 1.$ (43)

It follows that we can establish in Proof 1 that the convergence of the
coefficients of the PST also have the same rate

$\left\|{\mathcal{T}}-{\mathcal{T}}^{h}\right\|\leq Ch^{s},\quad 0\leq s\leq
1.$ (44)

###### Proof 1

Consider the coefficients of the exact PST ${\mathcal{T}}$, using the LP form,
which are given by

${\mathcal{T}}_{ij}=\int_{\Gamma}\xi_{j}\phi_{i}\,\operatorname{d\\!}\bm{\xi}.$
(45)

Let ${\mathcal{T}}^{h}_{ij}$ be a tensor computed numerically written as

${\mathcal{T}}^{h}_{ij}=\int_{\Gamma}\xi_{j}\phi_{i}^{h}\,\operatorname{d\\!}\bm{\xi},$
(46)

where $\phi_{i}^{h}$ is an approximate solution computed using (10) and the
effects of approximate numerical integration are ignored. Thus,

$\displaystyle\left\|{\mathcal{T}}-{\mathcal{T}}^{h}\right\|=$
$\displaystyle\left(\sum^{3}_{i,j=1}\left|{\mathcal{T}}_{ij}-{\mathcal{T}}^{h}_{ij}\right|^{2}\right)^{1/2}$
$\displaystyle=$
$\displaystyle\left(\sum^{3}_{i,j=1}\left|\int_{\Gamma}\xi_{j}\left(\phi_{i}-\phi_{i}^{h}\right)\,\operatorname{d\\!}\bm{\xi}\right|^{2}\right)^{1/2}$
$\displaystyle\leq$
$\displaystyle\left(\sum^{3}_{i=1}\left|\int_{\Gamma}\left(\phi_{i}-\phi_{i}^{h}\right)\,\operatorname{d\\!}\bm{\xi}\right|^{2}\right)^{1/2}\left(\sum^{3}_{j=1}\left|\int_{\Gamma}\xi_{j}\,\operatorname{d\\!}\bm{\xi}\right|^{2}\right)^{1/2},$

where the Cauchy-Schwarz inequality has been applied. This can be written as

$\displaystyle\leq$
$\displaystyle\left(\sum_{i=1}^{3}\left\|\phi_{i}-\phi_{i}^{h}\right\|_{L^{2}(\Gamma)}^{2}\right)^{1/2}\left(\sum_{j=1}^{3}\|\xi_{j}\|_{L^{2}(\Gamma)}^{2}\right)^{1/2}.$
(47)

Note that, $\|\xi_{j}\|_{L^{2}(\Gamma)}\leq C$ with $C$ independent of $h$ and
then we have

$\displaystyle\left\|{\mathcal{T}}-{\mathcal{T}}^{h}\right\|\leq$
$\displaystyle
C\left(\sum^{3}_{i=1}\left\|\phi_{i}-\phi_{i}^{h}\right\|_{L^{2}(\Gamma)}^{2}\right)^{1/2}.$
(48)

Substituting the error estimate in (43) into (48), we have

$\displaystyle\left\|{\mathcal{T}}-{\mathcal{T}}^{h}\right\|\leq$
$\displaystyle Ch^{s}\max_{i=1,2,3}\|\phi_{i}\|_{H^{s}(\Gamma)},$ (49)

with $C$ independent of $h$. Note that,

$\displaystyle\|\phi_{i}\|_{H^{s}(\Gamma)}\leq$ $\displaystyle C,$ (50)

where the constant $C$ depends of the shape of $B$, but not on $h$. Hence
leading to the quoted result.

## Acknowledgment

A.A.S. Amad and P.D. Ledger gratefully acknowledges the financial support
received from EPSRC in the form of grant EP/R002134/2 and T. Betcke gratefully
acknowledges the financial support received from EPSRC in the form of grant
EP/R002274/1. D. Praetorius acknowledges the support of the Austria Science
Fund FWF through grant P27005 and through the special research program Taming
complexity in partial differential systems (grant F65)

## References

  * [1] J. L. Mueller and S. Siltanen, _Linear and Nonlinear Inverse Problems with Practical Applications_. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2012.
  * [2] W. R. B. Lionheart, “EIT reconstruction algorithms: pitfalls, challenges and recent developments,” _Physiological Measurement_ , vol. 25, pp. 125–142, 2004.
  * [3] G. J. Saulnier, R. S. Blue, J. C. Newell, D. Isaacson, and P. M. Edic, “Electrical impedance tomography,” _IEEE Signal Processing Magazine_ , vol. 18, no. 6, pp. 31–43, 2001.
  * [4] A. Calderón, “On an inverse boundary value problem,” _Computational and Applied Mathematics_ , vol. 25, no. 2-3, pp. 133–138, 2006, reprinted from the Seminar on Numerical Analysis and its Applications to Continuum Physics, Sociedade Brasileira de Matemática, Rio de Janeiro, 1980.
  * [5] E. T. Chung, T. F. Chan, and X. C. Tai, “Electrical impedance tomography using level set representation and total variational regularization,” _Journal of Computational Physics_ , vol. 205, no. 1, pp. 357–372, 2005.
  * [6] M. Hintermüller, A. Laurain, and A. A. Novotny, “Second-order topological expansion for electrical impedance tomography,” _Advances in Computational Mathematics_ , vol. 36, no. 2, pp. 235–265, 2012.
  * [7] A. Ferreira and A. A. Novotny, “A new non-iterative reconstruction method for the electrical impedance tomography problem,” _Inverse Problems_ , vol. 33, no. 3, p. 035005, 2017.
  * [8] T. Rymarczyk, G. Kłosowski, and E. Kozłowski, “A non-destructive system based on electrical tomography and machine learning to analyze the moisture of buildings,” _Sensors_ , vol. 18, pp. 1–21, 2018.
  * [9] M. Cheney, D. Issacson, and J. Newell, “Electrical impedance tomography,” _SIAM Review_ , vol. 41, no. 1, pp. 85–101, 1999.
  * [10] H. Ammari and H. Kang, _Polarization and Moment Tensors: With Applications to Inverse Problems_ , ser. Applied Mathematical Sciences. Springer-Verlag, New York, 2007.
  * [11] P. D. Ledger and W. R. B. Lionheart, “The perturbation of electromagnetics fields at distances that are large compared with the object’s size,” _IMA Journal of Applied Mathematics_ , vol. 80, pp. 865–892, 2015.
  * [12] F. Watson, “Better imaging for landmine detection: An exploration of 3D full-wave inversion for ground penetrating radar,” Ph.D. dissertation, School of Mathematics, The University of Manchester, 2015.
  * [13] F. Watson and W. Lionheart, “SVD analysis of GPR full-wave inversion,” in _Proceedings of the 15th International Conference on Ground Penetrating Radar_ , June 2014, pp. 484–490.
  * [14] T. A. Khairuddin and W. R. B. Lionheart, “Characterization of objects by electrosensing fish based on the first order polarization tensor,” _Bioinspiration & Biomimetics_, vol. 11, no. 5, p. 055004, sep 2016.
  * [15] H. Ammari, T. Boulier, J. Garnier, and H. Wang, “Shape recognition and classification in electro-sensing,” _Proceedings of the National Academy of Sciences_ , vol. 111, no. 32, pp. 11 652–11 657, 2014.
  * [16] P. D. Ledger and W. R. B. Lionheart, “Characterising the shape and material properties of hidden targets from magnetic induction data,” _IMA Journal of Applied Mathematics_ , vol. 80, pp. 1776–1798, 2015.
  * [17] ——, “Understanding the magnetic polarizability tensor,” _IEEE Transactions on Magnetics_ , vol. 52, p. 6201216, 2016.
  * [18] ——, “Generalised magnetic polarizability tensors,” _Mathematical Methods in the Applied Sciences_ , vol. 41, pp. 3175–3196, 2018.
  * [19] Y. Zhao, W. Yin, C. Ktistis, D. Butterfield, and A. J. Peyton, “On the low-frequency electromagnetic responses of in-line metal detectors to metal contaminants,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 63, pp. 3181–3189, 2014.
  * [20] O. A. A. Rehim, J. L. Davidson, L. A. Marsh, M. D. O’Toole, D. W. Armitage, and A. J. Peyton, “Measurement system for determining the magnetic polarizability tensor of small metallic objects,” in _IEEE Sensor Application Symposium_ , 2015.
  * [21] Y. Zhao, W. Yin, C. Ktistis, D. Butterworth, and A. J. Peyton, “Determining the electromagnetic polarizability tensors of metal objects during in-line scanning,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 65, pp. 1172–1181, 2016.
  * [22] C. M. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006.
  * [23] T. A. Khairuddin and W. R. B. Lionheart, “Computing the first order polarization tensor: Welcome BEM++!” _Discovering Mathematics (Menemui Matematik)_ , vol. 35, no. 2, pp. 15–20, Dec. 2013.
  * [24] M. Lu, Q. Zhao, P. Hu, W. Yin, and A. J. Peyton, “Prediction of the asymptotical magnetic polarization tensors for cylindrical samples using the boundary element method,” in _2015 IEEE Sensors Applications Symposium (SAS)_ , 2015, pp. 1–4.
  * [25] W. Śmigaj, T. Betcke, S. Arridge, J. Phillips, and M. Schweiger, “Solving boundary integral problems with BEM++,” _ACM Transactions on Mathematical Software_ , vol. 41, no. 2, pp. 6:1–6:40, 2015.
  * [26] T. Betcke, A. Haberl, and D. Praetorius, “Adaptive boundary element methods for the computation of the electrostatic capacity on complex polyhedra,” _Journal of Computational Physics_ , vol. 397, p. 108837, 2019.
  * [27] D. J. Cedio-Fengya, S. Moskow, and M. S. Vogelius, “Identification of conductivity imperfections of small diameter by boundary measurements. continuous dependence and computational reconstruction,” _Inverse Problems_ , vol. 14, pp. 553–595, 1998.
  * [28] T. Betcke, M. W. Scroggs, and W. Śmigaj, “Product algebras for Galerkin discretisations of boundary integral operators and their applications,” _ACM Transactions on Mathematical Software_ , vol. 46, no. 2020, 2020.
  * [29] M. Feischl, T. Führer, M. Karkulik, and D. Praetorius, “ZZ-type a posteriori error estimators for adaptive boundary element methods on a curve,” _Engineering Analysis with Boundary Elements_ , vol. 38, pp. 49 – 60, 2014.
  * [30] A. Bespalov, T. Betcke, A. Haberl, and D. Praetorius, “Adaptive bem with optimal convergence rates for the Helmholtz equation,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 346, pp. 260 – 287, 2019\.
  * [31] C. T. M. Choi, “Adaptive finite element analysis of unbounded electromagnetic problems,” _IEEE Transactions on Magnetics_ , vol. 36, no. 4, pp. 705–707, 2000.
  * [32] S. Polstyanko and Jin-Fa Lee, “Adaptive finite element electrostatic solver,” _IEEE Transactions on Magnetics_ , vol. 37, no. 5, pp. 3120–3124, 2001.
  * [33] S. A. Sauter and C. Schwab, _Boundary Element Methods_ , ser. Springer Series in Computational Mathematics. Berlin: Springer, 2011.