# Toroidal Dipolar Supersolid with a Rotating Weak Link

M. Nilsson Tengstrand  P. Stürmer<EMAIL_ADDRESS>J. Ribbing
S.M. Reimann  Division of Mathematical Physics and NanoLund, Lund University,
Box 118, SE-221 00 Lund, Sweden

###### Abstract

Ring-shaped superfluids with weak links provide a perfect environment for
studying persistent currents and dynamic stirring protocols. Here, we
investigate the effects of a weak-link system on dipolar supersolids. By
calculating the ground state energy at fixed angular momenta, we find that
metastable persistent currents may exist in the supersolid phase near the
superfluid transition point. When stirring the weak link rapidly enough, we
show that vortices can enter the supersolid. These vortex entries cause phase
slips, emitting solitonic excitations that interfere with the crystalline
structure of the supersolid, leading to a continuous melting and
recrystallization of the droplets. Finally, we examine the release of vortex-
carrying supersolids from the trap, observing that the released density
exhibits a discrete structure associated with the density modulation and a
central hole resulting from the vortex core.

##  I Introduction

A many-body quantum state with off-diagonal long-range order exhibits
dissipationless flow, known as superfluidity Gross (1957); Yang (1962). This
leads to remarkable phenomena, such as quantized vortices and persistent
currents, observed in electron systems due to Bose-condensed Cooper pairs,
leading to SQUIDS and their applications in closed superconducting loops
Jaklevic  _et al._ (1964); Sato and Packard (2011). Similarly, persistent
currents occur in neutral systems like the early experiments in superfluid
liquid helium Allum  _et al._ (1977); Avenel and Varoquaux (1985); Hoskinson
_et al._ (2006) or Bose-condensed ultracold alkali atoms Leggett (1999);
Mueller (2002); Ryu  _et al._ (2007); Moulder  _et al._ (2012); Beattie  _et
al._ (2013). These persistent currents are metastable, have integer values of
angular momentum per particle in units of $\hbar$ Bloch (1973), and are
insensitive to perturbations and disorder. Including a small potential
barrier, a so-called ‘weak link’, allows for states of different integer
angular momentum to mix, effectively creating a two-level system Anderson  _et
al._ (2003); Solenov and Mozyrsky (2010a, b); Ryu  _et al._ (2013); Cominotti
_et al._ (2014); Aghamalyan  _et al._ (2015); Amico  _et al._ (2021, 2022).
Stirring the weak link dynamically Wright  _et al._ (2013a, b); Piazza  _et
al._ (2013); Ramanathan  _et al._ (2011); Mathey  _et al._ (2014) can induce
hysteresis between different persistent current states Eckel  _et al._
(2014a). As with any hysteresis process, the different states are well-
separated by a potential barrier. This potential barrier, however, can be
overcome with solitonic excitations which provide a steady-state solution and
allow for non-quantized persistent currents Muñoz Mateo  _et al._ (2015).
If the many-body quantum state exhibits off-diagonal and diagonal long-range
order simultaneously it results in the counter-intuitive state of
supersolidity Chester (1970); Leggett (1970), where superfluidity and
crystalline structure coexist. While long sought after in 4He Kim and Chan
(2004, 2012); Balibar (2010); Boninsegni and Prokof’ev (2012); Kim and Chan
(2012); Chan  _et al._ (2013), supersolidity has been observed in spin-orbit
or light-coupled condensates Henkel  _et al._ (2010); Lin  _et al._ (2011); Li
_et al._ (2016); Léonard  _et al._ (2017), where, however, it occurs due to
external factors and not due to the spontaneous symmetry breaking of both
translational and gauge symmetry. The spontaneous breaking of these symmetries
can be induced in ultracold dipolar Bose gases of dysprosium (164Dy) and
erbium (164Er) Tanzi  _et al._ (2019); Böttcher  _et al._ (2019); Chomaz  _et
al._ (2019), following the stabilization of a dipolar system against collapse
when reducing its scattering length Kadau  _et al._ (2016). The stabilization
is due to a beyond-mean-field correction term that is usually negligible and
occurs due to the interplay of different density dependencies of the mean-
field and beyond-mean-field terms Lee  _et al._ (1957); Lima and Pelster
(2011, 2012); Petrov (2015); Wächtler and Santos (2016). Only when tuning the
system’s scattering length to particular values via Feshbach resonances do
these quantum fluctuations have a significant effect. The increase of the
effective dipolar interaction strength
$\varepsilon_{\mathrm{dd}}=a_{\mathrm{dd}}/a$, where $a_{\mathrm{dd}}$ is the
dipolar length and $a$ the $s$-wave scattering length, beyond a critical value
accompanies a softening of the roton mode in the exctiation spectrum Pomeau
and Rica (1994); Saccani  _et al._ (2012); Macia  _et al._ (2012); Roccuzzo
and Ancilotto (2019); Natale  _et al._ (2019). The finite momentum value
associated with the roton minimum relates to the periodicity of the
crystalline structure within the supersolid or the isolated droplets Santos
_et al._ (2003).

Figure 1: (Color online) Density isosurfaces for a toroidal supersolid with
$\varepsilon_{\mathrm{dd}}=1.366$ and our wide link system, where the
transparent and solid isosurfaces are taken at density values $n=7.5\times
10^{13}$ cm-3 and $n=2\times 10^{14}$ cm-3, respectively. The grey-shaded area
is the contour of the toroidal trap with a weak link, and the red-shaded area
is the difference from a purely toroidal trap.

In previous work, we found a hysteretic cycle in the angular momentum states
in a toroidal supersolid of dipolar dysprosium atoms and demonstrated the
existence of persistent currents in these systems Nilsson Tengstrand  _et al._
(2021). We now turn to a similar system, however, with a weak link (see Fig.
1). The weak link causes density modulations to form around it, close to the
regime where the roton softens, similar to a dark solitary wave in a dipolar
superfluid Edmonds  _et al._ (2016). This weak link creates an avoided level-
crossing of angular momentum states, depending on the relative dipolar
strength and the weak link strength. When rotating the weak link, the droplets
move along while the superfluid stays stationary. However, as with purely
contact-interacting superfluids, there is a local flow within the weak link
Dubessy  _et al._ (2012); Woo and Son (2012); Mathey  _et al._ (2014). If the
flow velocity in this area approaches the local speed of sound it becomes
energetically favorable for a supercurrent to exist. With the emergence of the
persistent current, a phase slip and a corresponding vortex entry occur. We
differentiate between two different regimes. In the well-established
superfluid phonon-dominated regime, a solitonic wave-like excitation
accompanies the phase slip. The solitonic excitation decays over time as it
creates phonon excitations in the system. Further, the exact type of solitonic
excitation depends on the chemical potential. In a regime where the roton mode
softens, however, a roton excitation accompanies the emerging dark solitary
wave. This roton leads to a density modulation and a still strongly connected
supersolid. By bringing the system past the critical point into the supersolid
phase, the emitted roton and dark solitary excitation lead to complicated
dynamics between the droplets, resulting in a roton-mitigated melting and
recrystallization of neighboring crystal sites, reminiscent of the Lindemann
melting criterion Lindemann (1910); Khrapak (2020).

##  II Model and Methods

We model our dipolar condensate using a nonlocal extended Gross-Pitaevskii
equation

$i\hbar\partial_{t}\psi(\mathbf{r},t)=H_{\mathrm{GP}}\psi(\mathbf{r},t),$ (1)

where

$\displaystyle H_{\mathrm{GP}}=$
$\displaystyle-\frac{\hbar^{2}\nabla^{2}}{2m}+V_{\mathrm{trap}}(\mathbf{r})+g|\psi(\mathbf{r},t)|^{2}$
(2) $\displaystyle+\int
d\mathbf{r}^{\prime}V_{\mathrm{dd}}(\mathbf{r}-\mathbf{r}^{\prime})|\psi(\mathbf{r}^{\prime},t)|^{2}+\gamma|\psi(\mathbf{r},t)|^{3}.$

Here the order parameter $\psi$ is normalized to the total number of particles
$N$, $m$ is the mass of a single particle, and $g=4\pi\hbar^{2}a/m$ is the
contact interaction coupling constant with the $s$-wave scattering length $a$.
The nonlocal dipolar potential is
$V_{\mathrm{dd}}(\mathbf{r})=\frac{\mu_{0}\mu^{2}}{4\pi}\frac{1-\cos^{2}\beta}{|\mathbf{r}|^{3}}$,
where $\mu_{0}$ is the permeability in a vacuum, $\mu$ the magnetic moment,
and $\beta$ the angle between $\mathbf{r}$ and the dipole alignment direction,
which is taken to be along the $z$-axis. The final term is the beyond-mean-
field term in the local density approximation Lima and Pelster (2011, 2012);
Wächtler and Santos (2016), where
$\gamma=\frac{128\sqrt{\pi}}{3}\frac{\hbar^{2}a^{5/2}}{m}(1+\frac{3}{2}\varepsilon_{\mathrm{dd}}^{2})$,
and $a_{\mathrm{dd}}=m\mu_{0}\mu^{2}/(12\pi\hbar^{2})$. The ring-shaped trap
with a weak link is modeled according to

$V_{\mathrm{trap}}(\mathbf{r})=\frac{m\omega^{2}(\theta)}{2}\left[(\rho-\rho_{0})^{2}+\lambda^{2}z^{2}\right],$
(3)

where the dependence on the azimuthal angle $\theta$ of
$\omega(\theta)=\omega_{0}[1+Ae^{-\theta^{2}/(2w^{2})}]$ causes a constriction
around $\theta=0$. Here $\rho_{0}$ is the radius of the ring, $\omega_{0}$ the
axial trapping frequency, $\lambda$ the ratio of transversal and axial
trapping frequencies, and $A$ and $w$ are dimensionless numbers that
characterize the strength and width of the link, respectively. To study a
rotating link we consider a frame that moves along with it, adding to the
right-hand side of Eq. (1) the term $-\Omega(t)L_{z}\psi$, where $\Omega(t)$
is the time-dependent rotating frequency, and
$L_{z}=i\hbar(y\partial_{x}-x\partial_{y})$. To obtain the ground state at
some imposed angular momentum $L_{0}$ we minimize the quantity
$\tilde{E}=E+C(L-L_{0})^{2}$, where $E$ is the energy corresponding to Eq.
(1), $L=\int d\mathbf{r}\psi^{\ast}L_{z}\psi$, and $C>0$ a number that, when
sufficiently large, ensures that the energetic minimum occurs at $L\approx
L_{0}$ Komineas  _et al._ (2005). We solve the extended Gross-Pitaevskii
equation in imaginary and real-time by using the split-step Fourier method.

##  III Results

To exemplify our findings we consider a system consisting of 164Dy atoms and
fix the parameters $N=1\times 10^{5}$, $\rho_{0}=5$ $\mu$m, $\omega_{0}=100$
Hz, and $\lambda=1.1$. The amplitude of the weak link is set to $A=0.5$ with a
width equal to $w=0.3$ for a wide link or $w=0.03$ for a thin link.
For these parameters, we find that in the absence of a weak link the ground
state in the supersolid phase contains ten localized high-density sites.
Adding a weak link this amount reduces by one for the wide link and remains
unchanged for the thin link (see Fig. 2) where the integrated density along
the azimuthal angle $n_{\theta}(\theta)=\iint d\rho dz|\psi|^{2}$ is displayed
for some values of $\varepsilon_{\mathrm{dd}}$. The emergence of the
supersolid phase in the presence of a weak link is perhaps not as evident as
for the system without, since in the latter the supersolid transition in the
ring occurs precisely when azimuthal symmetry breaks. This spontaneous
symmetry breaking can not be taken as the transition point for the weak-link
system since the trap itself breaks azimuthal symmetry.

Figure 2: (Color online) Densities $n_{\theta}(\theta)=\iint d\rho
dz|\psi|^{2}$ as a function of the azimuthal angle $\theta$ when
$\varepsilon_{\mathrm{dd}}=0.7,1.3,1.354,1.366$ as indicated by the legend.
The top panel shows densities for the wide link ($w=0.3$) and the bottom panel
for the thin link ($w=0.03$). The insets in the top left corner of each panel
show density isosurfaces for the case $\varepsilon_{\mathrm{dd}}=1.366$, where
the transparent and solid isosurfaces are taken at density values $n=7.5\times
10^{13}$ cm-3 and $n=2\times 10^{14}$ cm-3, respectively.

We can see that when the dipolar interaction is weak compared to the contact
interaction the density is nearly homogeneous away from the link. This
behavior changes with increased relative interaction strength, and the atoms
start to cluster at the edges of the weak link, as shown here for
$\varepsilon_{\mathrm{dd}}=1.3$. This can be understood since the diminished
density in the link results in a decreased side-by-side repulsion in its
vicinity. For stronger interactions, the density acquires an oscillatory
behavior where the amplitude diminishes away from the link. For even larger
values of $\varepsilon_{\mathrm{dd}}$, the overall behavior stays the same as
the oscillation amplitude increases and clear regular localizations in the
density appear.
To gain some insight regarding the rotational properties of the weak-link
system, we consider the ground-state energy as a function of the angular
momentum for different values of $\varepsilon_{\mathrm{dd}}$, see Fig. 3 for
results with a wide link (we find the same qualitative behavior for $w=0.3$
and $w=0.03$ and choose to display results only for the former here). In the
superfluid phase, the energy forms dampened intersecting parabolas, similar to
a mixed Bose-Bose system in the bubble phase Stürmer  _et al._ (2022). Such
dampened behavior builds the basis of the neutral atomic qubit occurring for
repulsive Bose gases in the superfluid phase Anderson  _et al._ (2003);
Solenov and Mozyrsky (2010a, b); Ryu  _et al._ (2013); Cominotti  _et al._
(2014); Aghamalyan  _et al._ (2015); Amico  _et al._ (2021, 2022). The
dampened behavior is in contrast to a superfluid in an azimuthally symmetric
torus, where $E(L)$ is concave with a V-shaped minimum at $L/(N\hbar)=1$ Bloch
(1973). As we increase the relative interaction strength and approach the
supersolid phase, the intersection of the parabolas sharpens and moves towards
$L/(N\hbar)\sim 0.5$, and eventually, the energetic behavior becomes that of
two intersecting parabolas. This is similar to the toroidal dipolar supersolid
without a link Nilsson Tengstrand  _et al._ (2021), and has the interpretation
of rigid motion of the localized high-density regions together with
counterflow in the low-density regions for angular momenta smaller than
$L/(N\hbar)=0.5$, and rigid motion plus a singly quantized vortex nucleated in
the center for angular momenta larger than $L/(N\hbar)=0.5$ (and smaller than
$L/(N\hbar)=1.5$, beyond which multiply-quantized vortices become
energetically favorable). Since local energetic minima exist at finite $L>0$
for small enough $\varepsilon_{\mathrm{dd}}$, the supersolid may consequently
support persistent currents when a weak link is present.

Figure 3: (Color online) Ground-state energy as a function of angular momentum
for the wide ($w=0.3$) link when
$\varepsilon_{\mathrm{dd}}=1.3,1.34,1.36,1.366,1.375$, as indicated by the
legend. The energy unit is $E_{0}=\hbar^{2}/(m\xi^{2})$, where $\xi=1$
$\mu\mathrm{m}$.

An interesting question is whether it is possible to introduce vorticity to a
supersolid through the rotation of a weak link. To this end, we consider a
protocol where the link’s rotation frequency $\Omega$ is ramped linearly from
$\Omega=0$ to $\Omega=\Omega_{\mathrm{max}}$ over 250 ms and then kept at this
maximal value for another 250 ms. Fig. 4 shows the angular momentum as a
function of time $L(t)$ for both links for four different values of
$\Omega_{\mathrm{max}}$ 111See the Supplemental Material for movies of the
full time-evolution of the corresponding densities and phases. Before the
sharp jumps in angular momentum, all systems oscillate regularly in
$L/(N\hbar)$. This oscillation originates from the excitation of a mode
corresponding to the swaying back and forth of the localized high-density
regions of the supersolid. When the maximal rotation frequency is sufficiently
high, $L(t)$ is characterized by a quick increase in the angular momentum,
which corresponds to the entering of vortices into the condensate. We see that
the behavior for the different links is quite different and that the critical
rotation frequencies also differ. To estimate this frequency (where an
$s$-times quantized vortex is energetically favorable) it is possible to use
the analytical expression derived in Ref. Nilsson Tengstrand  _et al._ (2021)
for a torus without a weak link

$\Omega_{\mathrm{crit}}=\frac{(2s-1)\hbar}{2m\rho_{0}^{2}},$ (4)

which was derived for a system where the radius of the torus is much larger
than the length scale set by the trapping frequency, i.e.
$\rho_{0}\gg\sqrt{\hbar/(m\omega_{0})}$. For the parameters in this study this
translates to $\Omega_{\mathrm{crit}}/\omega_{0}\approx 0.012$ for $s=1$,
$\Omega_{\mathrm{crit}}/\omega_{0}\approx 0.037$ for $s=2$, and
$\Omega_{\mathrm{crit}}/\omega_{0}\approx 0.062$ for $s=3$. Although these
values serve as decent estimates, an inspection of the results in Fig. 4 shows
that these estimates are not exact. In particular, for the thin link when
$\Omega_{\mathrm{max}}=0.025$, two vortices enter in quick succession around
$t=250$ ms followed by the almost immediate exit of one of these, leaving only
a single vortex remaining in the system. Following the quick increase in
angular momentum associated with vorticity, $L(t)$ becomes strongly
oscillatory. As the vortices enter, a solitonic excitation appears that moves
more slowly than the link. This results in an interference between the
solitonic excitation and the localized structure of the supersolid, and gives
rise to additional oscillations on top of the already excited mode.

Figure 4: (Color online) Angular momentum as a function of time when the link
is set to rotate for $\varepsilon_{\mathrm{dd}}=1.366$ for the wide link (top
panel) and the thin link (bottom panel). The rotation frequency is ramped
linearly from $\Omega=0$ to $\Omega=\Omega_{\mathrm{max}}$ over 250 ms and
then kept at this maximal value for another 250 ms.

To understand this mechanism, let us investigate the entry of solitonic
excitations due to the phase slip associated with the stirring protocol. In
such a case, solitonic excitations provide a pathway for non-quantized
currents. Firstly, we recover the effects known for a Bose-gas in the phonon
regime (either purely contact interacting or at low relative dipolar
interaction strengths, here $\varepsilon_{\mathrm{dd}}=0.7$) in a torus with a
weak link Muñoz Mateo  _et al._ (2015). In this case, a vortex-antivortex pair
creates the counterflow in the region of the weak link, similar to the
Meissner effect Woo and Son (2012); Dubessy  _et al._ (2012). The resulting
phase jump increases with rotational speed until the flow velocity inside the
weak link is equal to the local speed of sound. At this point a phase slip
occurs, creating a solitonic vortex that moves along the torus Abad  _et al._
(2011); Yakimenko  _et al._ (2015). This solitonic excitation decays by
creating phononic excitations itself, moving closer to the outer edge of the
superfluid in the process.
For a system close to the roton-instability (here for
$\varepsilon_{\mathrm{dd}}=1.3$, see Fig. 5) this process is equal up to the
point of the emission of solitonic excitations. Here we find a dark solitary
wave moving through the condensate instead of a solitonic excitation.
Similarly to the above, the solitary wave creates quasi-particle excitations
which move with the local speed of sound. However, the excited quasi-particles
are rotons. This triggers the roton instability, such that one observes a
periodic density modulation preceding the dark solitary wave with the local
speed of sound. The density modulation itself disappears and reemerges
periodically as the solitary wave and rotonic excitation move through the
system.

Figure 5: (Color online) (top panel) Density distribution snapshots at given
times with isosurfaces taken at $n=7.5\times 10^{13}\,\mathrm{cm}^{-3}$ and
$n=2.0\times 10^{14}\,\mathrm{cm}^{-3}$. Note the varying amount of droplets
throughout. (bottom panel) Time evolution of the integrated angular density
$n_{\theta}(\theta)=\iint d\rho\,dz|\psi|^{2}$ for
$\varepsilon_{\mathrm{dd}}=1.3$.

Equipped with these insights, let us now turn to the supersolid phase (see
Fig. 6). As the phase slip occurs, similar to the superfluid case in the roton
regime, the dark solitary wave emerges. However, here the solitary wave
becomes untraceable as it encounters the first droplet. As before, the
solitary wave excites rotons, which move at the local speed of sound.
Additionally, the density modulation associated with the solitary wave creates
a new droplet between the solitary wave and the weak link. After the solitary
wave decays, the initial kick and propagation of the excited rotons determine
the dynamics and propagation of the system. These rotons interact with the
existing droplets, leading to a continuous melting and recrystallization of
droplets. If two or more droplets move too close to each other, they may melt
into a large, lower-density area (see the first two droplets below the weak
link at around $400$ ms in Fig. 6). This melting is similar to the mean
displacement argument in Lindemann’s melting criterion Lindemann (1910), where
one, based on a system’s temperature, estimates the solid-liquid phase
transition by calculating the mean displacement necessary for two crystal
sites to touch. On the other hand, if two individual droplets move too far
apart, a new droplet forms in between them (see the first two droplets below
the weak link at around $550$ ms in Fig. 6). The propagating roton-excitation
further mitigate this behavior when interacting with the already existing
structure.

Figure 6: (Color online) (top panel) Density distribution snapshots at given
times with isosurfaces taken at $n=7.5\times 10^{13}\,\mathrm{cm}^{-3}$ and
$n=2.0\times 10^{14}\,\mathrm{cm}^{-3}$. Note the varying amount of droplets
throughout. (bottom panel) Time evolution of the integrated angular density
$n_{\theta}(\theta)=\iint d\rho\,dz|\psi|^{2}$ for
$\varepsilon_{\mathrm{dd}}=1.366$. White dots mark the positions of each local
maxima with $n_{\theta}(\theta)>4000\,(\mu\mathrm{m})^{-1}$.

All three cases outlined above for the wide link are equivalent to the tight
link. In the supersolid case, however, the dark solitary wave dissipates much
faster.

Figure 7: (Color online) Densities in the $z=0$ plane 5 ms after the radial
trap has been turned off for the thin link ($w=0.03$) for different values of
the vortex quantization $s=0,1,3,5$, as indicated by the figure. The relative
interaction strength is $\varepsilon_{\mathrm{dd}}=1.366$ and the initial
systems prior to release are generated as ground states in the rotating frame
at rotation frequencies $\Omega/\omega_{0}=0,0.02,0.07,0.11$, respectively.
Note that the color scale is saturated and that the peak density at the center
of the top left panel is $n\approx 570$ $(\mu\mathrm{m})^{-3}$. The insets in
the bottom right corner of each panel show a zoom-in of the central region.

Finally, we look at the release of a supersolid from the weak-link trap for
configurations with different amounts of initial vorticity, created by
obtaining ground states in a rotating frame. Fig. 7 shows densities in the
$z=0$ plane 5 ms after the radial trap has been turned off (thus also removing
the link). When there is no initial vortex ($s=0$), a high-density region
forms at the center due to interference of the dispersing condensate (see the
top left panel in Fig. 7). Such a central density peak is prevented when
$s\neq 0$ due to the existence of a vortex core there, and the size of this
core increases with increasing $s$, similarly to what occurs for a regular BEC
made out of 23Na atoms Murray  _et al._ (2013). The coexistence of
superfluidity and a crystalline structure is thus clearly established for the
vortical supersolid in this type of time-of-flight setup, exhibiting both a
central hole and a discrete localized structure.

##  IV Conclusions and Outlook

In conclusion, we have studied dipolar supersolids in ring-shaped traps with a
weak link. Through a study of the energetic behavior, we found that persistent
currents may be supported in these systems. Simulating the real-time dynamics
of a rotating link, we saw that the entry of vortices resulted in the creation
of solitonic excitations traveling along the ring, interacting with the
supersolid structure. When released from the radial trap, we observed a
characteristic density pattern with a discrete structure together with a
central hole, providing a useful experimental signature for the coexistence of
superfluidity and spontaneous density modulation. Superfluid Bose-Einstein
condensates made out of 23Na in ring-shaped traps with a weak link have been
realized experimentally Wright  _et al._ (2013c); Eckel  _et al._ (2014b), and
the parameters used throughout this paper are of similar orders of magnitude
as in these experiments. Unlike the pure superfluid, however, the high-density
regions in the supersolid phase set a boundary on the lifetimes due to the
increased impact of three-body recombination losses. Finally, it would be
interesting to study the behavior of an asymmetric binary Bose mixture in the
droplet-superfluid phase in the presence of a weak link. These systems are
analogous to the dipolar supersolids studied here in many ways Petrov (2015);
Tengstrand and Reimann (2022), and may help pave the way in increasing our
understanding of beyond-mean-field physics in the context of atomtronics.

###### Acknowledgements.

This work was financially supported by the Knut and Alice Wallenberg
Foundation, the Swedish Research Council, and NanoLund. We are grateful for
many fruitful discussions with K. Mukherjee, T. Arnone Cardinale, T. Bland, as
well as A. Recati and S. Stringari at the initial stage of the project.

## References

  * Gross (1957) Eugene P. Gross, “Unified theory of interacting bosons,” Phys. Rev. 106, 161–162 (1957).
  * Yang (1962) C. N. Yang, “Concept of off-diagonal long-range order and the quantum phases of liquid he and of superconductors,” Rev. Mod. Phys. 34, 694–704 (1962).
  * Jaklevic  _et al._ (1964) R. C. Jaklevic, John Lambe, A. H. Silver, and J. E. Mercereau, “Quantum interference effects in josephson tunneling,” Phys. Rev. Lett. 12, 159–160 (1964).
  * Sato and Packard (2011) Y Sato and R E Packard, “Superfluid helium quantum interference devices: physics and applications,” Reports on Progress in Physics 75, 016401 (2011).
  * Allum  _et al._ (1977) DR Allum, Peter VE McClintock, A Phillips, and RM Bowley, “The breakdown of superfluidity in liquid 4he: an experimental test of landau’s theory,” Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 284, 179–224 (1977).
  * Avenel and Varoquaux (1985) O. Avenel and E. Varoquaux, “Observation of singly quantized dissipation events obeying the josephson frequency relation in the critical flow of superfluid ${}^{4}\mathrm{He}$ through an aperture,” Phys. Rev. Lett. 55, 2704–2707 (1985).
  * Hoskinson  _et al._ (2006) Emile Hoskinson, Yuki Sato, Inseob Hahn, and RE Packard, “Transition from phase slips to the josephson effect in a superfluid 4he weak link,” Nature Physics 2, 23–26 (2006).
  * Leggett (1999) A. J. Leggett, “Superfluidity,” Rev. Mod. Phys. 71, S318–S323 (1999).
  * Mueller (2002) Erich J. Mueller, “Superfluidity and mean-field energy loops: Hysteretic behavior in bose-einstein condensates,” Phys. Rev. A 66, 063603 (2002).
  * Ryu  _et al._ (2007) C. Ryu, M. F. Andersen, P. Cladé, Vasant Natarajan, K. Helmerson, and W. D. Phillips, “Observation of persistent flow of a bose-einstein condensate in a toroidal trap,” Phys. Rev. Lett. 99, 260401 (2007).
  * Moulder  _et al._ (2012) Stuart Moulder, Scott Beattie, Robert P. Smith, Naaman Tammuz, and Zoran Hadzibabic, “Quantized supercurrent decay in an annular bose-einstein condensate,” Phys. Rev. A 86, 013629 (2012).
  * Beattie  _et al._ (2013) Scott Beattie, Stuart Moulder, Richard J. Fletcher, and Zoran Hadzibabic, “Persistent currents in spinor condensates,” Phys. Rev. Lett. 110, 025301 (2013).
  * Bloch (1973) F. Bloch, “Superfluidity in a ring,” Phys. Rev. A 7, 2187–2191 (1973).
  * Anderson  _et al._ (2003) B. P. Anderson, K. Dholakia, and E. M. Wright, “Atomic-phase interference devices based on ring-shaped bose-einstein condensates: Two-ring case,” Phys. Rev. A 67, 033601 (2003).
  * Solenov and Mozyrsky (2010a) Dmitry Solenov and Dmitry Mozyrsky, “Metastable states and macroscopic quantum tunneling in a cold-atom josephson ring,” Phys. Rev. Lett. 104, 150405 (2010a).
  * Solenov and Mozyrsky (2010b) Dmitry Solenov and Dmitry Mozyrsky, “Macroscopic two-state systems in trapped atomic condensates,” Phys. Rev. A 82, 061601 (2010b).
  * Ryu  _et al._ (2013) C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, “Experimental realization of josephson junctions for an atom squid,” Phys. Rev. Lett. 111, 205301 (2013).
  * Cominotti  _et al._ (2014) Marco Cominotti, Davide Rossini, Matteo Rizzi, Frank Hekking, and Anna Minguzzi, “Optimal persistent currents for interacting bosons on a ring with a gauge field,” Phys. Rev. Lett. 113, 025301 (2014).
  * Aghamalyan  _et al._ (2015) Davit Aghamalyan, Marco Cominotti, Matteo Rizzi, Davide Rossini, Frank Hekking, Anna Minguzzi, Leong-Chuan Kwek, and Luigi Amico, “Coherent superposition of current flows in an atomtronic quantum interference device,” New Journal of Physics 17, 045023 (2015).
  * Amico  _et al._ (2021) L. Amico, M. Boshier, G. Birkl, A. Minguzzi, C. Miniatura, L.-C. Kwek, D. Aghamalyan, V. Ahufinger, D. Anderson, N. Andrei, A. S. Arnold, M. Baker, T. A. Bell, T. Bland, J. P. Brantut, D. Cassettari, W. J. Chetcuti, F. Chevy, R. Citro, S. De Palo, R. Dumke, M. Edwards, R. Folman, J. Fortagh, S. A. Gardiner, B. M. Garraway, G. Gauthier, A. Günther, T. Haug, C. Hufnagel, M. Keil, P. Ireland, M. Lebrat, W. Li, L. Longchambon, J. Mompart, O. Morsch, P. Naldesi, T. W. Neely, M. Olshanii, E. Orignac, S. Pandey, A. Pérez-Obiol, H. Perrin, L. Piroli, J. Polo, A. L. Pritchard, N. P. Proukakis, C. Rylands, H. Rubinsztein-Dunlop, F. Scazza, S. Stringari, F. Tosto, A. Trombettoni, N. Victorin, W. von Klitzing, D. Wilkowski, K. Xhani, and A. Yakimenko, “Roadmap on Atomtronics: State of the art and perspective,” AVS Quantum Science 3 (2021), 10.1116/5.0026178, 039201.
  * Amico  _et al._ (2022) Luigi Amico, Dana Anderson, Malcolm Boshier, Jean-Philippe Brantut, Leong-Chuan Kwek, Anna Minguzzi, and Wolf von Klitzing, “Colloquium: Atomtronic circuits: From many-body physics to quantum technologies,” Rev. Mod. Phys. 94, 041001 (2022).
  * Wright  _et al._ (2013a) K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, “Driving phase slips in a superfluid atom circuit with a rotating weak link,” Phys. Rev. Lett. 110, 025302 (2013a).
  * Wright  _et al._ (2013b) K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, “Threshold for creating excitations in a stirred superfluid ring,” Phys. Rev. A 88, 063633 (2013b).
  * Piazza  _et al._ (2013) F Piazza, L A Collins, and A Smerzi, “Critical velocity for a toroidal bose–einstein condensate flowing through a barrier,” Journal of Physics B: Atomic, Molecular and Optical Physics 46, 095302 (2013).
  * Ramanathan  _et al._ (2011) A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, “Superflow in a toroidal bose-einstein condensate: An atom circuit with a tunable weak link,” Phys. Rev. Lett. 106, 130401 (2011).
  * Mathey  _et al._ (2014) Amy C. Mathey, Charles W. Clark, and L. Mathey, “Decay of a superfluid current of ultracold atoms in a toroidal trap,” Phys. Rev. A 90, 023604 (2014).
  * Eckel  _et al._ (2014a) Stephen Eckel, Jeffrey G Lee, Fred Jendrzejewski, Noel Murray, Charles W Clark, Christopher J Lobb, William D Phillips, Mark Edwards, and Gretchen K Campbell, “Hysteresis in a quantized superfluid ‘atomtronic’circuit,” Nature 506, 200–203 (2014a).
  * Muñoz Mateo  _et al._ (2015) A. Muñoz Mateo, A. Gallemí, M. Guilleumas, and R. Mayol, “Persistent currents supported by solitary waves in toroidal bose-einstein condensates,” Phys. Rev. A 91, 063625 (2015).
  * Chester (1970) G. V. Chester, “Speculations on bose-einstein condensation and quantum crystals,” Phys. Rev. A 2, 256–258 (1970).
  * Leggett (1970) A. J. Leggett, “Can a solid be ”superfluid”?” Phys. Rev. Lett. 25, 1543–1546 (1970).
  * Kim and Chan (2004) Eunseong Kim and Moses Hung-Wai Chan, “Probable observation of a supersolid helium phase,” Nature 427, 225–227 (2004).
  * Kim and Chan (2012) Duk Y. Kim and Moses H. W. Chan, “Absence of supersolidity in solid helium in porous vycor glass,” Phys. Rev. Lett. 109, 155301 (2012).
  * Balibar (2010) Sebastien Balibar, “The enigma of supersolidity,” Nature 464, 176–182 (2010).
  * Boninsegni and Prokof’ev (2012) Massimo Boninsegni and Nikolay V. Prokof’ev, “Colloquium: Supersolids: What and where are they?” Rev. Mod. Phys. 84, 759–776 (2012).
  * Chan  _et al._ (2013) Moses Hung-Wai Chan, RB Hallock, and L Reatto, “Overview on solid 4 he and the issue of supersolidity,” Journal of Low Temperature Physics 172, 317–363 (2013).
  * Henkel  _et al._ (2010) N. Henkel, R. Nath, and T. Pohl, “Three-dimensional roton excitations and supersolid formation in rydberg-excited bose-einstein condensates,” Phys. Rev. Lett. 104, 195302 (2010).
  * Lin  _et al._ (2011) Y-J Lin, K Jiménez-García, and Ian B Spielman, “Spin-orbit-coupled bose–einstein condensates,” Nature 471, 83–86 (2011).
  * Li  _et al._ (2016) Junru Li, Wujie Huang, Boris Shteynas, Sean Burchesky, Furkan Çağrı Top, Edward Su, Jeongwon Lee, Alan O. Jamison, and Wolfgang Ketterle, “Spin-orbit coupling and spin textures in optical superlattices,” Phys. Rev. Lett. 117, 185301 (2016).
  * Léonard  _et al._ (2017) Julian Léonard, Andrea Morales, Philip Zupancic, Tilman Esslinger, and Tobias Donner, “Supersolid formation in a quantum gas breaking a continuous translational symmetry,” Nature 543, 87–90 (2017).
  * Tanzi  _et al._ (2019) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, “Observation of a dipolar quantum gas with metastable supersolid properties,” Phys. Rev. Lett. 122, 130405 (2019).
  * Böttcher  _et al._ (2019) Fabian Böttcher, Jan-Niklas Schmidt, Matthias Wenzel, Jens Hertkorn, Mingyang Guo, Tim Langen, and Tilman Pfau, “Transient supersolid properties in an array of dipolar quantum droplets,” Phys. Rev. X 9, 011051 (2019).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, “Long-lived and transient supersolid behaviors in dipolar quantum gases,” Phys. Rev. X 9, 021012 (2019).
  * Kadau  _et al._ (2016) Holger Kadau, Matthias Schmitt, Matthias Wenzel, Clarissa Wink, Thomas Maier, Igor Ferrier-Barbut, and Tilman Pfau, “Observing the rosensweig instability of a quantum ferrofluid,” Nature 530, 194–197 (2016).
  * Lee  _et al._ (1957) T. D. Lee, Kerson Huang, and C. N. Yang, “Eigenvalues and eigenfunctions of a bose system of hard spheres and its low-temperature properties,” Phys. Rev. 106, 1135–1145 (1957).
  * Lima and Pelster (2011) Aristeu R. P. Lima and Axel Pelster, “Quantum fluctuations in dipolar bose gases,” Phys. Rev. A 84, 041604 (2011).
  * Lima and Pelster (2012) A. R. P. Lima and A. Pelster, “Beyond mean-field low-lying excitations of dipolar bose gases,” Phys. Rev. A 86, 063609 (2012).
  * Petrov (2015) D. S. Petrov, “Quantum mechanical stabilization of a collapsing bose-bose mixture,” Phys. Rev. Lett. 115, 155302 (2015).
  * Wächtler and Santos (2016) F. Wächtler and L. Santos, “Quantum filaments in dipolar bose-einstein condensates,” Phys. Rev. A 93, 061603 (2016).
  * Pomeau and Rica (1994) Yves Pomeau and Sergio Rica, “Dynamics of a model of supersolid,” Phys. Rev. Lett. 72, 2426–2429 (1994).
  * Saccani  _et al._ (2012) S. Saccani, S. Moroni, and M. Boninsegni, “Excitation spectrum of a supersolid,” Phys. Rev. Lett. 108, 175301 (2012).
  * Macia  _et al._ (2012) A. Macia, D. Hufnagl, F. Mazzanti, J. Boronat, and R. E. Zillich, “Excitations and stripe phase formation in a two-dimensional dipolar bose gas with tilted polarization,” Phys. Rev. Lett. 109, 235307 (2012).
  * Roccuzzo and Ancilotto (2019) Santo Maria Roccuzzo and Francesco Ancilotto, “Supersolid behavior of a dipolar bose-einstein condensate confined in a tube,” Phys. Rev. A 99, 041601 (2019).
  * Natale  _et al._ (2019) G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, “Excitation spectrum of a trapped dipolar supersolid and its experimental evidence,” Phys. Rev. Lett. 123, 050402 (2019).
  * Santos  _et al._ (2003) L. Santos, G. V. Shlyapnikov, and M. Lewenstein, “Roton-maxon spectrum and stability of trapped dipolar bose-einstein condensates,” Phys. Rev. Lett. 90, 250403 (2003).
  * Nilsson Tengstrand  _et al._ (2021) M. Nilsson Tengstrand, D. Boholm, R. Sachdeva, J. Bengtsson, and S. M. Reimann, “Persistent currents in toroidal dipolar supersolids,” Phys. Rev. A 103, 013313 (2021).
  * Edmonds  _et al._ (2016) M. J. Edmonds, T. Bland, D. H. J. O’Dell, and N. G. Parker, “Exploring the stability and dynamics of dipolar matter-wave dark solitons,” Phys. Rev. A 93, 063617 (2016).
  * Dubessy  _et al._ (2012) R. Dubessy, T. Liennard, P. Pedri, and H. Perrin, “Critical rotation of an annular superfluid bose-einstein condensate,” Phys. Rev. A 86, 011602 (2012).
  * Woo and Son (2012) S. J. Woo and Young-Woo Son, “Vortex dynamics in an annular bose-einstein condensate,” Phys. Rev. A 86, 011604 (2012).
  * Lindemann (1910) F. Lindemann, “The calculation of molecular vibration frequencies,” Z. Phys. 11, 609–614 (1910).
  * Khrapak (2020) Sergey A. Khrapak, “Lindemann melting criterion in two dimensions,” Phys. Rev. Res. 2, 012040 (2020).
  * Komineas  _et al._ (2005) S. Komineas, N. R. Cooper, and N. Papanicolaou, “Single vortex states in a confined bose-einstein condensate,” Phys. Rev. A 72, 053624 (2005).
  * Stürmer  _et al._ (2022) P. Stürmer, M. Nilsson Tengstrand, and S. M. Reimann, “Mixed bubbles in a one-dimensional bose-bose mixture,” Phys. Rev. Res. 4, 043182 (2022).
  * Note (1) See the Supplemental Material for movies of the full time-evolution of the corresponding densities and phases.
  * Abad  _et al._ (2011) M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M. Jezek, “Phase slippage and self-trapping in a self-induced bosonic josephson junction,” Phys. Rev. A 84, 035601 (2011).
  * Yakimenko  _et al._ (2015) A. I. Yakimenko, K. O. Isaieva, S. I. Vilchinskii, and E. A. Ostrovskaya, “Vortex excitation in a stirred toroidal bose-einstein condensate,” Phys. Rev. A 91, 023607 (2015).
  * Murray  _et al._ (2013) Noel Murray, Michael Krygier, Mark Edwards, K. C. Wright, G. K. Campbell, and Charles W. Clark, “Probing the circulation of ring-shaped bose-einstein condensates,” Phys. Rev. A 88, 053615 (2013).
  * Wright  _et al._ (2013c) K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, “Driving phase slips in a superfluid atom circuit with a rotating weak link,” Phys. Rev. Lett. 110, 025302 (2013c).
  * Eckel  _et al._ (2014b) S. Eckel, J. G. Lee, F. Jendrzejewski, N. Murray, C. W. Clark, C. J. Lobb, W. D. Phillips, M. Edwards, and G. K. Campbell, “Hysteresis in a quantized superfluid ’atomtronic’ circuit,” Nature 506, 200–203 (2014b).
  * Tengstrand and Reimann (2022) M. Nilsson Tengstrand and S. M. Reimann, “Droplet-superfluid compounds in binary bosonic mixtures,” Phys. Rev. A 105, 033319 (2022).


# Precise Damage Shaping in Self-Sensing Composites Using Electrical Impedance
Tomography and Genetic Algorithms

Hashim Hassan  Tyler N. Tallman  School of Aeronautics and Astronautics,
Purdue University, West Lafayette, IN, 47907, United States

###### Abstract

Fiber-reinforced composites with nanofiller-modified polymer matrices have
immense potential to improve the safety of high-risk engineering structures.
These materials are intrinsically self-sensing because their electrical
conductivity is affected by deformations and damage. This property, known as
piezoresistivity, has been extensively leveraged for conductivity-based damage
detection via electrical resistance change methods and tomographic imaging
techniques such as electrical impedance tomography (EIT). Although these
techniques are very effective at detecting the presence of damage, they suffer
from an inability to provide precise information about damage shape, size, or
mechanism. This is particularly detrimental for laminated composites which can
suffer from complex failure modes, such as delaminations, that are difficult
to detect. To that end, we herein propose a new technique for precisely
determining damage shape and size in self-sensing composites. Our technique
makes use of a genetic algorithm (GA) integrated with realistic physics-based
damage models to recover precise damage shape from conductivity changes imaged
via EIT. We experimentally validate this technique on carbon nanofiber
(CNF)-modified glass fiber-reinforced polymer (GFRP) laminates by considering
two specific damage mechanisms: through-holes and delaminations. Our results
show that this novel technique can accurately reconstruct multiple through-
holes with radii as small as 1.19 mm and delaminations caused by low velocity
impacts. These findings illustrate that coupling piezoresistivity with
conductivity-based spatial imaging techniques and physics-based inversion
strategies can enable damage shaping capabilities in self-sensing composite
structures.

###### keywords:

structural health monitoring , self-sensing composites , electrical impedance
tomography , damage shaping , genetic algorithms

††journal: arXiv

[table]capposition=top

##  1 Introduction

Fiber-reinforced composites are ubiquitous in aerospace, automotive,
mechanical, and light-weight civil structures due to their high strength-to-
weight ratios, good fatigue performance, and design flexibility. However, it
is well known that composite structures are susceptible to matrix cracking,
delamination, and fiber breakage. If left unchecked, these damage modes can
accumulate and cause catastrophic failure. For structures in safety-critical
applications such as aircraft and bridges, this may result in significant loss
of life. Therefore, damage must be inspected and preventive measures taken
before it reaches a critical level.

Traditional damage tolerant design approaches that rely on intermittent
inspections via nondestructive evaluation (NDE) often prove insufficient for
composites due to the complexity of their damage initiation and growth
mechanisms. Structural health monitoring (SHM) [1, 2] is a promising
complement to traditional NDE as it allows for continuous, real-time
monitoring of damage initiation and growth in composite structures.
Conventional SHM methodologies rely on the use of external sensors to monitor
the condition of the structure. Three common categorizations of conventional
SHM are vibration-based SHM, guided wave-based SHM, and embedded sensor-based
SHM.

Vibration-based SHM [3, 4] monitors changes in the modal parameters of the
structure such as natural frequencies, mode shapes, and curvature changes.
These methods have been used to detect delaminations [5, 6] and matrix cracks
[7] in composite plates and also in circular cylindrical composite shells [8].
Guided wave-based SHM [9, 10, 11] uses ultrasonic wave propagation to detect
the presence of damage. The waves are dispersed or attenuated when they
encounter a flaw or a crack. The most common types of waves used for detecting
cracks and delaminations include Lamb waves [12, 13] and Rayleigh waves [14,
15]. Embedded sensor-based SHM uses discrete sensors integrated within the
structure to monitor its condition. The most common types of sensors used for
composites are strain gauges [16, 17, 18], fiber-optic sensors [19, 20, 21,
22], and piezoelectric sensors [23, 24].

Although conventional SHM methodologies are more widely used, they have
certain limitations. Vibration-based SHM, for example, is often not sensitive
enough to detect small-sized defects. Both guided wave-based and embedded
sensor-based SHM require physical integration of the sensors and actuators
within the structure which can sometimes weaken the material. Another
overarching limitation is that these methods only provide localized, discrete,
and often limited information about the condition of the structure.
Consequently, recent research focus has shifted toward the development of
composites with intrinsic self-sensing properties to alleviate the limitations
of conventional SHM.

Self-sensing composites possess one or more properties that allow them to
‘sense’ changes to their mechanical, thermal, or chemical state. That is, the
material itself functions as a sensor. In this area, nanofiller-modified
composites have received considerable attention [25, 26, 27, 28, 29]. These
materials exhibit changes in their electrical behavior due to deformations and
damage. This self-sensing property, known as piezoresistivity, has been
extensively studied for SHM in combination with conductivity-based spatial
imaging techniques such as electrical impedance tomography (EIT) because this
combination allows for the mapping of mechanical effects through conductivity
changes.

Beyond just allowing for spatial mapping of mechanics in self-sensing
composites, EIT has recently emerged as a potential SHM modality because it is
low-cost, non-invasive, and has nearly real-time imaging capabilities. Broadly
speaking, and as summarized in [30], EIT has been used for SHM in
intrinsically and extrinsically self-sensing composites. Intrinsically self-
sensing composites incorporate conductive nanofillers in the polymer matrix.
For such materials, EIT has been used for detecting through-holes [31], impact
damage [32], matrix cracking [33], delaminations [34], and UV radiation damage
[35, 36] in composites incorporating carbon nanotubes (CNTs) [33], carbon
nanofibers (CNFs) [37, 38, 39, 40], and carbon black (CB) [31]. Although most
research has focused on thin plates, there has been some work in using EIT for
damage detection in non-planar structures as well [41, 42]. Extrinsically
self-sensing composites do not incorporate nanofillers but make use of
externally applied sensing media. In this area, EIT has been used for
detecting strain and damage in composites using nanofiller-based thin films
[43, 44] and paints [45, 46]. Aside from composites, EIT has been used
extensively for detecting damage due to holes [47, 48], cracks [49], and
moisture flows [50] in self-sensing cementitious materials [51, 52].

From the preceding discussion, it is clear that EIT can effectively detect and
localize numerous damage modes in self-sensing materials. However, a critical
limitation is that it provides little-to-no information about precise damage
shape, size, or mechanism. This information is crucial for accurate
prognostication of structural health. For example, if the precise size and
location of a delamination can be ascertained, structural failure can be
predicted and averted. Herein, we recognize the opportunity to address this
limitation and advance the state-of-the-art in SHM of self-sensing composites.
We propose a novel methodology for precisely determining damage shape and size
from electrical measurements in self-sensing composites. Our unique approach
uses a genetic algorithm (GA) to inversely compute damage geometry based on
imaged conductivity changes and boundary voltage measurements. We
experimentally investigate the potential of this approach on two damage
mechanisms—through-holes and delaminations—in CNF-modified GFRP laminates.

The remainder of this manuscript is organized as follows. First, we
mathematically formulate the EIT forward and inverse problems. Second, we
present and formulate the GA-enabled damage shaping problem. Third, we
describe the composite manufacturing procedures. Fourth, we discuss the
electrical, through-hole, and impact testing procedures. Fifth, we present the
EIT results along with the results from our GA-enabled damage shaping
methodology and discuss their significance. Finally, we end with a brief
summary, conclusions, and possible paths for future work.

##  2 Electrical Impedance Tomography

###  2.1 Forward Problem

As discussed earlier, EIT is a non-invasive method of imaging the internal
conductivity distribution of a domain. EIT works by minimizing the difference
between a vector of experimentally measured voltages and a vector of voltages
computed numerically. The experimental voltage collection procedure is
illustrated in Figure 1. The domain to be imaged is lined with electrodes.
Current is injected between the first pair of electrodes and voltage
differences are measured between the remaining electrode pairs. The current
injection is then moved to the next electrode pair and voltage differences are
again measured between the remaining electrode pairs. This process is repeated
until all electrode pairs have received one current injection and a vector of
$L(L-3)$ voltage differences has been obtained, where $L$ is the total number
of electrodes.

Figure 1: Illustration of first and second current injection and voltage
measurement schemes. The periphery of the domain is instrumented with
electrodes, indicated by the red rectangles. Current is injected between the
first pair of electrodes and voltage differences are measured between the
remaining electrode pairs. The current injection is then moved to the next
electrode pair and voltage differences are again measured. This process
repeats until all adjacent electrode pairs have received a current injection.

The numerical voltages are computed by solving the EIT forward problem. To
formulate the forward problem, consider Laplace’s equation for steady-state
diffusion in the absence of internal current sources. This is shown in
equation (1), where $\sigma$ is the internal conductivity of the domain and
$\phi$ is the domain potential. We enforce two complete electrode model
boundary conditions on equation (1). The first boundary condition, shown in
equation (2), assumes a voltage drop across the electrodes due to the contact
impedance between the domain and the electrodes. The second boundary
condition, shown in equation (3), enforces conservation of charge. In these
equations, $z_{l}$ is the contact impedance between the $l$th electrode and
the domain, $\bm{n}$ is an outward pointing normal vector, $E_{l}$ is the
length of the $l$th electrode, and $V_{l}$ is the voltage of the $l$th
electrode.

$\nabla\cdot\sigma\nabla\phi=0$ (1)

$\phi+z_{l}\sigma\nabla\phi\cdot\bm{n}=V_{l}$ (2)

$\sum_{l=1}^{L}\int_{E_{l}}\sigma\nabla\phi\cdot\bm{n}\hskip
3.0pt\mathrm{d}S_{l}=0$ (3)

Equations (1) to (3) are often solved via the finite element (FE) method, as
shown in equation (4). In this equation, $\bm{\Phi}$ is a vector of domain
potentials, $\bm{V}$ is a vector of electrode voltages and $\bm{I}$ is a
vector of injected currents. The matrices $\bm{A_{M}}$, $\bm{A_{Z}}$,
$\bm{A_{W}}$, and $\bm{A_{D}}$ are formed as shown in equations (5) to (8),
where $N_{i}$ is the $i$th FE basis function. For an explicit solution to
these equations for linear elements, interested readers are directed to [30].

$\begin{bmatrix}\bm{A_{M}}+\bm{A_{Z}}&\bm{A_{W}}\\\
\bm{A_{W}^{T}}&\bm{A_{D}}\\\ \end{bmatrix}\begin{bmatrix}\bm{\Phi}\\\
\bm{V}\end{bmatrix}=\begin{bmatrix}\bm{0}\\\ \bm{I}\end{bmatrix}$ (4)

$A^{e}_{M\hskip 3.0ptij}=\int_{\Omega_{e}}\frac{\partial N_{i}}{\partial
x_{k}}\sigma_{kl}\frac{\partial N_{j}}{\partial x_{l}}\hskip
3.0pt\mathrm{d}\Omega_{e}$ (5)

$A_{Z\hskip 3.0ptij}=\sum_{l=1}^{L}\int_{E_{l}}\frac{1}{z_{l}}N_{i}N_{j}\hskip
3.0pt\mathrm{d}S_{l}$ (6)

$A_{W\hskip 3.0ptli}=-\int_{E_{l}}\frac{1}{z_{l}}N_{i}\hskip
3.0pt\mathrm{d}S_{l}$ (7)

$A_{D}=\text{diag}\bigg{(}\frac{E_{l}}{z_{l}}\bigg{)}$ (8)

###  2.2 Inverse Problem

The goal of the EIT inverse problem is to recover the conductivity
distribution of the domain interior using the measured and numerically
computed boundary voltages. Herein, we will make use of a one-step
linearization scheme to solve the inverse problem. This involves collecting a
set of pre-damage voltages, a set of post-damage voltages, and then finding
the conductivity change distribution that minimizes the difference between the
two sets of voltages. We begin by defining the vector $\delta\bm{V}$ as the
difference between voltages collected at times $t_{2}$ and $t_{1}$, as shown
in equation (9).

$\delta\bm{V}=\bm{V}(\sigma_{2},t_{2})-\bm{V}(\sigma_{1},t_{1})$ (9)

Next, we define an equivalent numerical vector, $\bm{W}(\delta\bm{\sigma})$,
as shown in equation (10). Here, $\bm{F}(\cdot)$ represent the voltages
computed by solving the forward problem at the conductivity in the argument,
$\bm{\sigma_{0}}$ is the baseline or undamaged conductivity, and
$\delta\bm{\sigma}$ is the conductivity change vector we seek. The boldfaced
symbols indicate vector quantities that have been discretized via finite
elements.

$\bm{W}(\delta\bm{\sigma})=\bm{F}(\bm{\sigma_{0}}+\delta\bm{\sigma})-\bm{F}(\bm{\sigma_{0}})$
(10)

We then linearize the first term on the right-hand side of equation (10) using
a Taylor series expansion about $\bm{\sigma_{0}}$ and retain only the linear
terms, as shown in equation (11).

$\bm{F}(\bm{\sigma_{0}}+\delta\bm{\sigma})\approx\bm{F}(\bm{\sigma_{0}})+\frac{\partial\bm{F}(\bm{\sigma_{0}})}{\partial\bm{\sigma}}\delta\bm{\sigma}$
(11)

Substituting equation (11) into equation (10) and defining $\bm{J}$ =
$\partial\bm{F}(\bm{\sigma_{0}})/\partial\bm{\sigma}$ as the sensitivity
matrix yields equation (12). The one-step linearized inverse problem then
seeks the conductivity change, $\delta\bm{\sigma}^{*}$, that minimizes the
difference between $\bm{W}(\delta\bm{\sigma})$ and $\delta\bm{V}$ in the
constrained linear least squares sense, as shown in equation (13).

$\bm{W}(\delta\bm{\sigma})\approx\bm{J}\delta\bm{\sigma}$ (12)

$\delta\bm{\sigma}^{*}=\underset{\bm{-\sigma_{0}}\leq\delta\bm{\sigma}\leq
0.01\bm{\sigma_{0}}}{\text{min}}\frac{1}{2}\Bigg{(}\Bigg{|}\Bigg{|}\begin{bmatrix}\bm{J}\\\
\alpha\bm{L}\\\ \end{bmatrix}\delta\bm{\sigma}-\begin{bmatrix}\delta\bm{V}\\\
\bm{0}\\\ \end{bmatrix}\Bigg{|}\Bigg{|}^{2}_{2}\Bigg{)}$ (13)

A regularization term, $\bm{L}$, has been included in equation (13) due to the
ill-posed nature of the EIT inverse problem. Here, we will use the discrete
Laplace operator as the regularization term. The amount of regularization is
controlled by the scalar parameter $\alpha$. Also, note that the search for
$\delta\bm{\sigma}^{*}$ is constrained. The constraints used here are based on
realistic assumptions about the conductivity change. It is well understood
that damage such as cracks, through-holes, and delaminations manifests as a
decrease in conductivity. The maximum decrease that is physically possible has
a magnitude equal to the baseline conductivity and we expect the maximum
increase to be relatively small. Therefore, we use 100% change for the lower
bound and 1% change for the upper bound. This is reflected in the limits shown
in equation (13).

##  3 GA-Enabled Damage Shaping Problem

As discussed earlier, EIT has been successfully used to detect various kinds
of damage in self-sensing composites. However, due to the ill-posed nature of
the inverse problem, the spatial resolution of EIT is somewhat indistinct
inasmuch as precise geometric features of conductivity artifacts cannot be
sharply reconstructed. Considerable research has aimed at developing
methodologies for precise artifact shape reconstruction from EIT images [53,
54, 55]. While these approaches are undoubtedly powerful, they are unaware of
the underlying physics of the material and consequently cannot discern between
various damage modes. Therefore, in order to accurately determine damage shape
and size from EIT-imaged conductivity changes, specific damage mechanisms must
be integrated with suitable inversion strategies.

The goal of the damage shaping problem presented here is to recover the
geometry of specific damage modes from EIT-imaged conductivity changes and
boundary voltages by building some knowledge about the damage physics into the
shaping algorithm. However, recovering mechanics from electrical measurements
is an ill-posed, multi-modal inverse problem. This is because there is no
well-defined mathematical relationship between damage geometry and
conductivity change. In the absence of a suitable inversion algorithm and
realistic constraints on the search space, the search may converge to a
locally optimum, non-physical solution. In order to address this problem we
need a global search algorithm and constraints based on realistic damage
physics.

Genetic algorithms (GAs) [56, 57] are a family of global search algorithms
inspired by natural evolution. GAs are widely used for optimizing multi-modal
functions by generating and evolving a population of solutions dispersed
inside a search space. The population of solutions evolves through a process
analogous to natural selection. Herein, we do not aim to develop a new
algorithm but to instead use a well-established GA to solve the multi-modal
damage shaping problem. To that end, we will make use of a GA developed by WA
Crossley and used by Raghavan et al. [58] for jet engine turbine blade NDE.

In addition to a global search algorithm, we require geometric models that
realistically describe how different damage modes affect material
conductivity. These models will need to be integrated with the GA so that only
physically viable solutions are generated and selected. We can now formulate
the GA-enabled damage shaping problem. We begin by defining the vector
$\delta\bm{F}(s_{GA})$ as shown in equation (14), where
$\bm{F}(\bm{\sigma}(s_{GA}))$ are the voltages predicted by solving the EIT
forward problem for a domain containing damage described by the GA-generated
parameter, $s_{GA}$. This parameter describes the shape and location of a
specific damage mode based on a geometric damage model.

$\delta\bm{F}(s_{GA})=\bm{F}(\bm{\sigma}(s_{GA}))-\bm{F}(\bm{\sigma_{0}})$
(14)

The goal of the damage shaping problem then is to seek the optimum shape
parameter, $s_{GA}^{*}$, that minimizes the $l_{1}$-norm of the difference
between $\delta\bm{F}(s_{GA})$, given by equation (14), and $\delta\bm{V}$,
given by equation (9). This can be cast as the optimization problem shown in
equation (15). In this work, we will develop geometric models for two damage
modes—through-holes and delaminations—and integrate them with the GA to solve
equation (15). These models are described later in §6.

$s_{GA}^{*}=\underset{s_{\text{min}}\leq s_{GA}\leq
s_{\text{max}}}{\arg\min}(\underbrace{|\delta\bm{V}-\delta\bm{F}(s_{GA})|_{1}}_{f})$
(15)

In the above equation, the parenthetical term is the fitness function, $f$.
During a single search, as the population of solutions evolves, $f$ decreases
and finally attains a stable value. This occurs when the candidates in the
population have achieved a certain level of genetic similarity. This is
quantified as the bit-string affinity (BSA) and is a pre-specified value. When
the BSA is achieved, further evolution ceases. The BSA is one of two
convergence criteria that we will use herein. The second convergence criterion
is the maximum number of generations. That is, after the population has
evolved a pre-specified number of times, further evolution ceases. At this
point the optimum solution, $s^{*}_{GA}$, and minimum fitness function value,
$f^{*}$, are recorded and a second search initiates using updated values of
$s_{\text{min}}$ and $s_{\text{max}}$. This is repeated until the following
global convergence criterion is satisfied, where $f_{n}^{*}$ is the minimum
fitness function value obtained after the $n$th search.

$|f_{n}^{*}-f_{n+1}^{*}|\leq 1\times 10^{-3}$ (16)

At this point, no further searches are performed using the GA and the final
solution attained is treated as the converged solution.

##  4 Composite Manufacturing

To experimentally validate the proposed GA-enabled damage shaping methodology,
GFRP laminates with 1.0 wt.% CNF-modified polymer matrices were manufactured
based on the procedure described by Tallman et al. [38]. The CNF-modified
matrix was prepared by mixing epoxy resin (Fibre Glast 2000) with the
appropriate amounts of Pyrograf III PR-24-XT-HHT CNFs (Applied Sciences),
surfactant (Triton X-100), and acetone. Acetone aids mixing by lowering the
viscosity and the surfactant facilitates dispersion of the CNFs by modifying
the surface chemistry of the mixture. An epoxy-to-acetone volume ratio of 2:1
and a surfactant-to-CNF weight ratio of 0.76:1 were used. The modified epoxy
was then mixed in a planetary centrifuge for 3 minutes and sonicated in a bath
sonicator for 1 hour for every 10 g of the mixture. The bath sonicator
operated at a power of 35 W and a frequency of 45 kHz. After sonication, the
mixture was stirred on a magnetic hot plate stirrer for 24 hours at a
temperature of 60 ∘C and a stirring speed of 600 rpm. The epoxy was then
allowed to cool to room temperature and hardener was added using an epoxy-to-
hardener weight ratio of 100:27. Air release agent (BYK A-501) was also added
using an air release agent-to-total mixture weight ratio of 0.001:1. The
mixture was then stirred by hand for 5 minutes and degassed at room
temperature for 20 minutes in a vacuum chamber.

Unidirectional glass fiber sheets (Fibre Glast Saertex) were impregnated with
the degassed CNF-modified epoxy using hand lay-up to produce two 10” $\times$
10” laminates with stacking sequences of $[0/90]_{\text{s}}$ and thicknesses
of 3 mm. We will refer to these as laminates 1 and 2. Both laminates were
vacuum bagged and cured in an oven for 5 hours at 60 ∘C. After curing, two
square specimens measuring 3.25” $\times$ 3.25” were cut from each laminate.
using a water-cooled tile saw. This resulted in a total of four plate-like
specimens, one of which is used for later through-hole tests and three of
which are used for later impact tests. Additionally, ten smaller specimens
measuring 0.25” $\times$ 0.25” were cut out of the remaining material from the
10” $\times$ 10” laminates for the purpose of measuring the average in-plane
and through-thickness conductivity of each plate.

##  5 Experimental Procedures

###  5.1 Electrical Testing

The conductivities of CNF/GFRP were measured to estimate the baseline
conductivity for EIT and also for the damage shaping algorithm. For this,
conductivity measurements were collected from the ten 0.25” $\times$ 0.25”
specimens cut from each laminate. The conductivities of each square were
measured in the principal ($x$, $y$, $z$) coordinate system, illustrated in
Figure 2. To measure the conductivity in the $x$-direction, electrodes were
applied by painting colloidal silver patches on the faces of the square with
normal vectors pointing in the $\pm x$-directions. The resistance, $R$,
between the electrodes was measured using a multi-meter and the conductivity
was calculated as $\sigma_{x}=l_{x}/RA_{e}$, where $l_{x}$ is the length of
the square in the $x$-direction and $A_{e}$ is the area of an electrode. The
conductivities in the $y$\- and $z$-directions were measured similarly. The
measurement scheme is illustrated in Figure 2 and the mean conductivities
measured from both laminates are listed in Table 1.

Figure 2: Illustration of conductivity measurement scheme and principal coordinate system. The ellipses represent the fiber cross-sections. Top: Schematic of a sub-section of the laminate showing the principal coordinate system and the layer orientations. Bottom: Conductivity measurement in the $x$-(left), $y$-(middle), and $z$-directions (right). The filled red rectangles represent electrodes on the front faces and the outlined red rectangles represent electrodes on the back faces. Table 1: Mean and standard deviations of conductivities in the $x$-, $y$-, and $z$-directions measured from laminates 1 and 2. Laminate |  $\sigma_{x}$ [S/m] |  $\sigma_{y}$ [S/m] |  $\sigma_{z}$ [S/m]
---|---|---|---
1 | $0.024\pm 0.003$ | $0.023\pm 0.004$ | $(3.5\pm 1.8)\times 10^{-4}$
2 | $0.054\pm 0.003$ | $0.057\pm 0.004$ | $(5.6\pm 2.2)\times 10^{-4}$

###  5.2 Through-Hole and Impact Testing

The four 3.25” $\times$ 3.25” specimens were used for through-hole and impact
testing as follows. For the two specimens cut from laminate 1, one was used
for through-hole testing and the other was impacted with 25 J. For the two
specimens cut from laminate 2, one was impacted with 23 J and the other was
impacted with 28 J. Note, however, that all specimens tested were made using
the same procedure described above and used the same CNF weight fraction. In
order to collect EIT data, electrodes were applied to each specimen by
painting evenly spaced colloidal silver patches on each edge. Each specimen
was then adhered to an acrylic base and additional colloidal silver patches
were painted on the acrylic to act as extended electrode tabs. A
representative specimen with electrodes painted is shown in Figure 3. Each
specimen was connected to a current source (Keithley 6221) and a DAQ system
(National Instruments PXIe-6368) to measure the electrode voltages.

(a)

(b)

Figure 3: Representative CNF/GFRP square laminate used for through-hole and
impact testing. Left: Top-down view of laminate without electrodes. Right:
Laminate with electrodes painted and adhered to acrylic base.

For the through-hole testing specimen, a current magnitude of 0.25 mA was
used. One set of voltages was collected from the specimen in its undamaged
state. Next, three circular holes of radii 1.19 mm, 2.38 mm, and 3.18 mm were
successively drilled and voltages were collected after drilling each new hole.

For the three impact testing specimens, the energies (23 J, 25 J, and 28 J)
were chosen deliberately to induce a measurable delamination without
completely perforating the specimens. A current magnitude of 0.25 mA was used
for the 23 J and 28 J impact specimens while a current magnitude of 0.1 mA was
used for the 25 J impact specimen. One set of voltages was collected from each
specimen in its undamaged state. The specimens were then mounted on a 6”
$\times$ 4” aluminum plate with a centrally located 2” $\times$ 2” window and
impacted in a drop tower (CEAST 9340) using a steel hemispherical striking
head with a diameter of 15.8 mm. Post-impact voltages were then collected from
each specimen. All three specimens were then destructively evaluated in order
to determine the true size and shape of the delaminations. This was done by
making multiple cuts through each specimen using a water-cooled tile saw and
examining the delamination length at the cross-section of each cut using a
Zeiss Axioskop 2 MAT microscope. Several images along the length of the
delamination at each cross-section were taken and were used to reconstruct the
full delamination shape.

##  6 Results and Discussion

###  6.1 Through-Holes

The mean conductivities measured from laminate 1, shown in Table 1, were used
as the initial estimate for the baseline conductivity of the through-hole
testing specimen. The initial estimate was then iteratively updated to find
the optimal baseline conductivity that minimized the difference between the
pre-damage voltages and the forward operator predicted voltages. Using this
approach, the optimum in-plane baseline conductivities were found to be
$\sigma_{x0}$ = $\sigma_{y0}$ = 0.025 S/m. EIT was then performed using this
baseline on a mesh consisting of 5,670 linear triangular elements and the EIT-
imaged conductivity changes are shown in the left column in Figure 4. It can
be seen from Figure 4 that EIT is able to successfully detect the presence of
the holes. However, it provides little-to-no information about the precise
shape and size of each hole. This information must be indirectly inferred from
the magnitude of the observed conductivity change. That is, a larger
conductivity change corresponds to a bigger hole.

In order to reconstruct the precise hole size using the damage shaping
methodology, a geometric model for through holes must be integrated with the
GA. In this case the model consisted of a circular hole inside the FE mesh and
the damage shape parameter was specified as $s_{GA}$ = $[x_{c},y_{c},r]$,
where $x_{c}$ and $y_{c}$ are the $x$\- and $y$-coordinates of the center of
the hole, respectively, and $r$ is the radius of the hole. The GA initiates
the search by generating a population between the bounds $s_{\text{min}}$ and
$s_{\text{max}}$, which are the lower and upper bounds, respectively, on the
location and size of the hole. These are specified based on the size and
location of the artifacts observed in the EIT images. For example, for the
case of one hole, $s_{\text{min}}$ was specified as [55, 55, 0.5] mm and
$s_{\text{max}}$ was [62, 62, 5] mm. An adaptive meshing algorithm [59] was
used to integrate this geometric model with the GA. Additionally, a population
size of 50, a BSA of 99%, and a maximum of 30 generations were used for the
GA. The results are shown in the right column in Figure 4.

Figure 4: EIT and GA-enabled damage shaping results from through-hole testing.
The black circles indicate the true hole sizes and locations. Left column:
EIT-imaged conductivity change. We can see that EIT is able to detect the
presence of the holes but the exact shape and location are not clear. Right
column: GA-enabled damage shaping results. It can be seen that by integrating
a realistic damage model with the GA, the shape and size of the holes can be
determined much more precisely than standard EIT.

In the first case, the GA is able to reconstruct a hole with a radius of 1.33
mm. In the second case, the GA is able to reconstruct one hole with a radius
of 1.34 mm and a second hole with a radius of 2.79 mm. And in the third case,
the GA is able to reconstruct three holes with radii of 1.34 mm, 2.79 mm, and
3.43 mm. We can see that these results are considerably more accurate and
provide much more precise information about the damage shape and size than
standard EIT, which often suffers from background noise and stray artifacts.
The average percent errors in the GA predicted hole sizes and locations are
listed in Table 2.

Table 2: Average percent error in GA-predicted solutions relative to actual hole sizes and locations. Here, $e_{x_{c}}$ and $e_{y_{c}}$ are the average percent errors in the $x$\- and $y$-coordinates of the center of the hole and $e_{r}$ is the average percent error in the hole radius. Hole number |  $e_{x_{c}}$ [%] |  $e_{y_{c}}$ [%] |  $e_{r}$ [%]
---|---|---|---
1 | 1.67 | 3.33 | 11.76
2 | 3.33 | 5.77 | 17.22
3 | 9.52 | 3.70 | 7.86

The BSA convergence for the first search of each case and the fitness function
convergence for all searches of each case are shown in Figure 5. For one hole,
we observe good convergence with one search. For two holes, we require a
higher number of generations to reach the BSA stopping criterion. In this
case, successive searches slightly improve the solution. For three holes, the
maximum number of generations is reached before the BSA stopping criterion is
met. In this case, successive searches significantly improve the solution.
This is because as the number of holes increases, the number of variables in
$s_{GA}$ also increases and the GA requires more computations to achieve a
genetically similar population.

Figure 5: Convergence plots for through-hole reconstruction. Left: BSA
convergence for first search case of each case. Right: Fitness function
convergence for all searches of each case.

###  6.2 Delaminations

The optimum baseline conductivities of the impact testing specimens were
estimated using a similar approach as the through-hole testing specimen. The
conductivities of the 23 J and 28 J impact specimens (from laminate 2) were
found to be $\sigma_{x0}$ = $\sigma_{y0}$ = 0.055 S/m, and
$\sigma_{z0}=5\times 10^{-4}$ S/m, and the conductivities of the 25 J impact
specimen (from laminate 1) were found to be $\sigma_{x0}$ = $\sigma_{y0}$ =
0.02 S/m, and $\sigma_{z0}=3\times 10^{-4}$ S/m. EIT was performed on a mesh
of 5,670 linear triangular elements and the conductivity change results from
each impact are shown in the top row in Figure 9. Note that EIT endeavored to
reconstruct the in-plane conductivity, $\sigma_{x0}$ = $\sigma_{y0}$ =
$\sigma$, which is isotropic – even though each lamina is mildly electrically
anisotropic, the symmetric layup renders the in-plane conductivity of the
entire laminate as very nearly isotropic. We can immediately observe that EIT
is able to detect each of the impacts and as the impact energy increases the
intensity of the conductivity change at the impact location also increases.
This indicates that increasing the impact energy is increasing the amount of
damage within the laminate. However, from the EIT image alone, the size and
shape of the induced delamination cannot be obviously determined.

To compute the delamination shape and size, the GA was integrated with a
geometric model for delaminations as follows. FE models of each laminate were
constructed using $\sigma_{x0}$, $\sigma_{y0}$, and $\sigma_{z0}$. These
models represent homogenized or ‘effective’ laminates that have the same bulk
electrical behavior as the physical laminate. Each FE model consisted of one
layer of 2,700 linear quadrilateral elements sandwiched between two layers of
2,700 linear hexahedral elements. The GA was then used to generate elliptical
conductivity artifacts described by the shape parameter $s_{GA}$ =
$[x_{c},y_{c},r_{x},r_{y}]$ inside the middle layer. A conductivity of
$\sigma_{x}$ = $\sigma_{y}$ = $\sigma_{z}$ = $\sigma_{d}$ = $1\times 10^{-6}$
S/m was assigned to the elements within the elliptical artifact. This concept
is based on the fact that a delamination can be modeled as a failure at the
interface between two layers of the laminate. In this case, these are the two
layers of hexahedral elements that represent the homogenized composite. The
individual layers of the physical laminate were not explicitly modeled for the
inversion. These were not modeled because the EIT electrodes spanned the
thickness of the laminate thereby making it impossible to extract through-
thickness information. As such, the GA-produced delamination was inserted at
the mid-plane of the homogenized models. A schematic of the delamination
modeling procedure is shown in Figure 6. For the GA, a population size of 20,
a BSA of 65%, and a maximum of 20 generations were used.

Figure 6: Schematic of delamination model integrated with the GA. The
experimentally measured bulk conductivities are used to build a homogenized FE
model with three layers of hexahedral elements. The middle layer is then
collapsed to infinitesimal thickness and the GA is used to generate elliptical
conductivity artifacts within this layer of elements.

As described earlier, optical microscopy was used to destructively evaluate
the specimens. Figure 7 shows an optical micrograph of the cross-section of
the 23 J laminate. A delamination is clearly visible between the third
($0^{\circ}$) and fourth ($90^{\circ}$) layers. In fact, it was observed that
the delamination always occurred between the third and fourth layers for each
impact. The full shape of the delamination was then reconstructed as follows.
For each cut made through a specimen, the delamination length was recorded and
images were taken at regular intervals along the cross-section. A new cut was
then made, the delamination length recorded again, and more images taken. In
other words, the delamination length was measured with the aid of an optical
microscope for each slice. This was repeated until a delamination was no
longer observed using the microscope. The recorded delamination lengths and
cut thicknesses were then used to reconstruct the in-plane shape of the
delamination. This process is illustrated for a single cut in Figure 8.

Figure 7: Optical microscopy image of cross-section of 23 J impact specimen. A
delamination is clearly visible between the third and fourth layers, as
indicated by the white arrows. The fibers in the $90^{\circ}$ layers are
pointing out of the page. Figure 8: Illustration of actual delamination shape
reconstruction using optical microscopy. Left: A cut is made through the
laminate to expose the cross-section. Middle: Optical microscopy is used to
measure the delamination length at the cross-section. Right: The cut thickness
and delamination length are used to piece together the delamination. This
procedure is repeated until the full delamination shape is obtained.

The GA-generated delamination shapes along with the actual delamination shapes
reconstructed using optical microscopy are shown in the second and third rows,
respectively, in Figure 9. We can see immediately that the GA-generated
solutions agree very well with the actual delamination shapes. As the impact
energy increases, the size of the delamination also increases, which is what
we observe from both the GA results and optical microscopy. There are
nonetheless some differences between the GA-generated solutions and the actual
delamination shapes. This may be attributed to the fact that an impact does
not cause a pure delamination but a combination of delamination, matrix
cracking, and fiber breakage. However, the proposed model does not account for
these additional damage modes. As such, the algorithm tries to compensate for
the additional damage by increasing the size of the predicted delamination.
Regardless, we can see from Figure 9 that the GA-predicted solutions are a
significant improvement over standard EIT in terms of damage shape.

Figure 9: EIT, delamination shaping, and optical microscopy results from
impact testing specimens. First (top) row: EIT-imaged conductivity change. The
magnitude of the conductivity change increases as the impact energy increases.
Second row: GA-enabled delamination shaping results. Third row: Actual
delamination shapes reconstructed using optical microscopy. Bottom row: Photos
of post-impacted specimens. The silver circles indicate the actual impact
locations. Note that damage caused by the impacts is barely visible on the
surface of the specimens.

Lastly, we examine the convergence of the delamination shaping algorithm. The
fitness function convergence for the first search and the minimum fitness
function value for all searches are shown in Figure 10. We can see that the
fitness function attains its minimum value for a relatively small number of
generations. Additionally, the change in the minimum fitness function value
between successive searches is not significant. This indicates that successive
searches may not necessarily produce better quality solutions for the
delamination shaping problem.

Figure 10: Fitness function convergence plots for GA-enabled delamination
shaping algorithm. Left: Fitness function convergence for first search of each
impact case. Right: Fitness function convergence for all searches of each
case.

##  7 Summary and Conclusions

In this article we have presented a new technique for more accurately
determining the shape and size of damage in self-sensing composites via EIT.
Our technique uses a GA to inversely determine the shape of particular damage
modes from EIT-imaged conductivity changes and boundary voltages. This is a
significant advancement to state-of-the-art conductivity-based SHM which seeks
only to spatially localize damage and does not provide any information about
the underlying mechanics. The technique proposed here addresses this
limitation and can transform SHM and NDE in self-sensing structures from mere
damage detection and localization to much more robust damage characterization.

We began by manufacturing self-sensing CNF-modified GFRP laminates. The
laminates were damaged by drilling multiple through-holes and impacting with
three different energies to cause a delamination. EIT was performed on each
laminate to image the damage-induced conductivity change. Two physics-based
geometric models for damage were then developed—one for through-holes and one
for delaminations—and each was integrated with a GA. The GA-integrated damage
shaping algorithms were then used to reconstruct the shape and size of the
holes and delaminations based on the observed conductivity changes and
boundary voltages. The inversely reconstructed through-holes showed good
agreement with the actual hole sizes. Furthermore, the algorithm was able to
accurately reconstruct multiple holes. The delamination shaping results also
showed good agreement with the actual delamination shapes determined
destructively using optical microscopy. The work presented here has potential
to improve the safety of self-sensing structures in safety-critical
applications by providing precise information about damage condition. However,
some limitations of this approach must be recognized and addressed in future
work. First, the use of a GA adds significant computational burden to EIT. We
anticipate that this can be overcome by using more efficient high-performance
computing architecture. Second, although this approach seems to show much
potential, this exploratory study was limited to through-holes and
delaminations. Additional damage models should be developed to capture more
realistic and complex damage scenarios. Appropriate damage models will have to
developed for these cases too. Finally, machine learning and artificial neural
networks have shown promise in solving damage inversion and categorization
problems [60]. These should also be explored in future work.

## Acknowledgements

The authors gratefully thank Professor WA Crossley at Purdue University for
the use of his GA, and Professor W Chen and graduate research assistant
Nesredin Kedir at the Impact Science Lab at Purdue University for their
assistance with optical microscopy. The authors also thank the Purdue
University Graduate School for providing financial support through the
Bilsland Dissertation Fellowship.

## References

  * [1] H. Sohn, C. R. Farrar, F. M. Hemez, D. D. Shunk, D. W. Stinemates, B. R. Nadler, J. J. Czarnecki, A review of structural health monitoring literature: 1996–2001, Los Alamos National Laboratory, USA 1 (2003).
  * [2] J. P. Lynch, An overview of wireless structural health monitoring for civil structures, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365 (1851) (2007) 345–372.
  * [3] E. P. Carden, P. Fanning, Vibration based condition monitoring: a review, Structural health monitoring 3 (4) (2004) 355–377.
  * [4] D. Montalvao, N. M. M. Maia, A. M. R. Ribeiro, A review of vibration-based structural health monitoring with special emphasis on composite materials, Shock and vibration digest 38 (4) (2006) 295–324.
  * [5] Y. Zou, L. Tong, G. P. Steven, Vibration-based model-dependent damage (delamination) identification and health monitoring for composite structures—a review, Journal of Sound and vibration 230 (2) (2000) 357–378.
  * [6] S. D. Valdes, C. Soutis, Delamination detection in composite laminates from variations of their modal characteristics, Journal of sound and vibration 228 (1) (1999) 1–9.
  * [7] S. S. Kessler, S. M. Spearing, M. J. Atalla, C. E. Cesnik, C. Soutis, Damage detection in composite materials using frequency response methods, Composites Part B: Engineering 33 (1) (2002) 87–95.
  * [8] K.-H. Ip, P.-C. Tse, Locating damage in circular cylindrical composite shells based on frequency sensitivities and mode shapes, European Journal of Mechanics-A/Solids 21 (4) (2002) 615–628.
  * [9] M. Mitra, S. Gopalakrishnan, Guided wave based structural health monitoring: A review, Smart Materials and Structures 25 (5) (2016) 053001.
  * [10] A. Raghavan, Guided-wave structural health monitoring, Ph.D. thesis, University of Michigan Ann Arbor (2007).
  * [11] A. J. Croxford, P. D. Wilcox, B. W. Drinkwater, G. Konstantinidis, Strategies for guided-wave structural health monitoring, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463 (2087) (2007) 2961–2981.
  * [12] Z. Su, L. Ye, Y. Lu, Guided lamb waves for identification of damage in composite structures: A review, Journal of sound and vibration 295 (3-5) (2006) 753–780.
  * [13] C. A. Paget, S. Grondel, K. Levin, C. Delebarre, Damage assessment in composites by lamb waves and wavelet coefficients, Smart materials and Structures 12 (3) (2003) 393.
  * [14] S. K. Chakrapani, V. Dayal, The interaction of rayleigh waves with delaminations in composite laminates, The Journal of the Acoustical Society of America 135 (5) (2014) 2646–2653.
  * [15] J. M. Hughes, J. Vidler, C.-T. Ng, A. Khanna, M. Mohabuth, L. F. Rose, A. Kotousov, Comparative evaluation of in situ stress monitoring with rayleigh waves, Structural Health Monitoring 18 (1) (2019) 205–215.
  * [16] A. Ajovalasit, Advances in strain gauge measurement on composite materials, Strain 47 (4) (2011) 313–325.
  * [17] S. Zike, L. P. Mikkelsen, Correction of gauge factor for strain gauges used in polymer composite testing, Experimental Mechanics 54 (3) (2014) 393–403.
  * [18] K. Zhou, Z. Y. Wu, Strain gauge placement optimization for structural performance assessment, Engineering Structures 141 (2017) 184–197.
  * [19] C. K. Y. Leung, Z. Yang, Y. Xu, P. Tong, S. K. L. Lee, Delamination detection in laminate composites with an embedded fiber optical interferometric sensor, Sensors and Actuators A: Physical 119 (2) (2005) 336–344.
  * [20] H. Lee, Y. Ochi, T. Matsui, Y. Matsumoto, Y. Tanaka, H. Nakamura, Y. Mizuno, K. Nakamura, Distributed strain measurement and possible breakage detection of optical-fiber-embedded composite structure using slope-assisted brillouin optical correlation-domain reflectometry, Applied Physics Express 11 (7) (2018) 072501.
  * [21] M. Ramakrishnan, G. Rajan, Y. Semenova, G. Farrell, Overview of fiber optic sensor technologies for strain/temperature sensing applications in composite materials, Sensors 16 (1) (2016) 99.
  * [22] J. Leng, A. Asundi, Structural health monitoring of smart composite materials by using efpi and fbg sensors, Sensors and Actuators A: Physical 103 (3) (2003) 330–340.
  * [23] V. Giurgiutiu, G. Santoni-Bottai, Structural health monitoring of composite structures with piezoelectric-wafer active sensors, AIAA journal 49 (3) (2011) 565–581.
  * [24] H. Kalhori, L. Ye, S. Mustapha, Inverse estimation of impact force on a composite panel using a single piezoelectric sensor, Journal of Intelligent Material Systems and Structures 28 (6) (2017) 799–810.
  * [25] Y. Li, X. Huang, L. Zeng, R. Li, H. Tian, X. Fu, Y. Wang, W.-H. Zhong, A review of the electrical and mechanical properties of carbon nanofiller-reinforced polymer composites, Journal of Materials Science 54 (2) (2019) 1036–1076.
  * [26] P. Costa, J. Nunes-Pereira, J. Oliveira, J. Silva, J. A. Moreira, S. Carabineiro, J. Buijnsters, S. Lanceros-Mendez, High-performance graphene-based carbon nanofiller/polymer composites for piezoresistive sensor applications, Composites Science and Technology 153 (2017) 241–252.
  * [27] S. Stassi, V. Cauda, G. Canavese, C. F. Pirri, Flexible tactile sensing based on piezoresistive composites: A review, Sensors 14 (3) (2014) 5296–5332.
  * [28] L. Groo, J. Nasser, L. Zhang, K. Steinke, D. Inman, H. Sodano, Laser induced graphene in fiberglass-reinforced composites for strain and damage sensing, Composites Science and Technology 199 (2020) 108367.
  * [29] L. Groo, J. Nasser, D. Inman, H. Sodano, Fatigue damage tracking and life prediction of fiberglass composites using a laser induced graphene interlayer, Composites Part B: Engineering 218 (2021) 108935.
  * [30] T. N. Tallman, D. J. Smyl, Structural health and condition monitoring via electrical impedance tomography in self-sensing materials: a review, Smart Materials and Structures 29 (12) (2020) 123001.
  * [31] T. N. Tallman, S. Gungor, K. Wang, C. E. Bakis, Damage detection via electrical impedance tomography in glass fiber/epoxy laminates with carbon black filler, Structural Health Monitoring 14 (1) (2015) 100–109.
  * [32] G. J. Gallo, E. T. Thostenson, Spatial damage detection in electrically anisotropic fiber-reinforced composites using carbon nanotube networks, Composite Structures 141 (2016) 14–23.
  * [33] H. Dai, G. J. Gallo, T. Schumacher, E. T. Thostenson, A novel methodology for spatial damage detection and imaging using a distributed carbon nanotube-based composite sensor combined with electrical impedance tomography, Journal of Nondestructive Evaluation 35 (2) (2016) 1–15.
  * [34] M. Haingartner, S. Gschoßmann, M. Cichocki, M. Schagerl, Improved current injection pattern for the detection of delaminations in carbon fiber reinforced polymer plates using electrical impedance tomography, Structural Health Monitoring (2020) 1475921720972308.
  * [35] S. Rastogi, D. Bartolo, S. Gurses, C. Kronawitter, V. La Saponara, Structural health monitoring of irradiated high-density polyethylene samples with electrical resistance tomography, Journal of Materials Science 56 (31) (2021) 17824–17842.
  * [36] M. Clausi, E. Toto, S. Botti, S. Laurenzi, V. La Saponara, M. G. Santonicola, Direct effects of UV irradiation on graphene-based nanocomposite films revealed by electrical resistance tomography, Composites Science and Technology 183 (2019) 107823.
  * [37] T. N. Tallman, S. Gungor, G. M. Koo, C. E. Bakis, On the inverse determination of displacements, strains, and stresses in a carbon nanofiber/polyurethane nanocomposite from conductivity data obtained via electrical impedance tomography, Journal of Intelligent Material Systems and Structures 28 (18) (2017) 2617–2629.
  * [38] T. N. Tallman, S. Gungor, K. W. Wang, C. E. Bakis, Damage detection and conductivity evolution in carbon nanofiber epoxy via electrical impedance tomography, Smart Materials and Structures 23 (4) (2014) 045034.
  * [39] H. Hassan, T. N. Tallman, Failure prediction in self-sensing nanocomposites via genetic algorithm-enabled piezoresistive inversion, Structural Health Monitoring 19 (3) (2020) 765–780.
  * [40] H. Hassan, T. N. Tallman, A comparison of metaheuristic algorithms for solving the piezoresistive inverse problem in self-sensing materials, IEEE Sensors Journal 21 (1) (2020) 659–666.
  * [41] A. J. Thomas, J. J. Kim, T. N. Tallman, C. E. Bakis, Damage detection in self-sensing composite tubes via electrical impedance tomography, Composites Part B: Engineering 177 (2019) 107276.
  * [42] J. Jauhiainen, M. Pour-Ghaz, T. Valkonen, A. Seppänen, Nonplanar sensing skins for structural health monitoring based on electrical resistance tomography, Computer-Aided Civil and Infrastructure Engineering (2021).
  * [43] T.-C. Hou, K. J. Loh, J. P. Lynch, Spatial conductivity mapping of carbon nanotube composite thin films by electrical impedance tomography for sensing applications, Nanotechnology 18 (31) (2007) 315501.
  * [44] T.-C. Hou, K. J. Loh, J. P. Lynch, Electrical impedance tomography of carbon nanotube composite materials, in: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 2007, Vol. 6529, International Society for Optics and Photonics, 2007, p. 652926.
  * [45] W. Lestari, B. Pinto, V. La Saponara, J. Yasui, K. J. Loh, Sensing uniaxial tensile damage in fiber-reinforced polymer composites using electrical resistance tomography, Smart Materials and Structures 25 (8) (2016) 085016.
  * [46] B. R. Loyola, V. La Saponara, K. J. Loh, T. M. Briggs, G. O’Bryan, J. L. Skinner, Spatial sensing using electrical impedance tomography, IEEE Sensors Journal 13 (6) (2013) 2357–2367.
  * [47] S. Gupta, J. G. Gonzalez, K. J. Loh, Self-sensing concrete enabled by nano-engineered cement-aggregate interfaces, Structural Health Monitoring 16 (3) (2017) 309–323.
  * [48] S. Nayak, S. Das, Spatial damage sensing ability of metallic particulate-reinforced cementitious composites: Insights from electrical resistance tomography, Materials & Design 175 (2019) 107817.
  * [49] D. Smyl, M. Pour-Ghaz, A. Seppänen, Detection and reconstruction of complex structural cracking patterns with electrical imaging, NDT & E International 99 (2018) 123–133.
  * [50] M. Hallaji, A. Seppänen, M. Pour-Ghaz, Electrical resistance tomography to monitor unsaturated moisture flow in cementitious materials, Cement and Concrete Research 69 (2015) 10–18.
  * [51] M. Hallaji, A. Seppänen, M. Pour-Ghaz, Electrical impedance tomography-based sensing skin for quantitative imaging of damage in concrete, Smart Materials and Structures 23 (8) (2014) 085001.
  * [52] M. Hallaji, M. Pour-Ghaz, A new sensing skin for qualitative damage detection in concrete elements: Rapid difference imaging with electrical resistance tomography, NDT & E International 68 (2014) 13–21.
  * [53] D. Liu, D. Smyl, D. Gu, J. Du, Shape-driven difference electrical impedance tomography, IEEE Transactions on Medical Imaging 39 (12) (2020) 3801–3812.
  * [54] D. Liu, D. Gu, D. Smyl, J. Deng, J. Du, Shape reconstruction using boolean operations in electrical impedance tomography, IEEE transactions on medical imaging 39 (9) (2020) 2954–2964.
  * [55] W. Fan, C. Wang, A new damage estimation method for carbon fiber reinforced polymer based on electrical impedance tomography, Review of Scientific Instruments 92 (2) (2021) 025102.
  * [56] D. E. Golberg, Genetic algorithms in search, optimization, and machine learning, Addion wesley 1989 (102) (1989) 36.
  * [57] M. Mitchell, An introduction to genetic algorithms, MIT press, 1998.
  * [58] S. Raghavan, P. Imbrie, W. A. Crossley, Spectral analysis of R-lines and vibronic sidebands in the emission spectrum of ruby using genetic algorithms, Applied spectroscopy 62 (7) (2008) 759–765.
  * [59] D. Engwirda, Locally optimal Delaunay-refinement and optimisation-based mesh generation, Ph.D. thesis, University of Sydney (2014).
  * [60] G. Lin, L. Zhao, T. N. Tallman, Real-time precise damage characterization in self-sensing materials via neural network-aided electrical impedance tomography: A computational study, ES Materials & Manufacturing (2020).