# Relativistic first principles theory of Yu–Shiba–Rusinov states applied to
an Mn adatom and Mn dimers on Nb(110)

Bendegúz Nyári1 András Lászlóffy1,2 László Szunyogh1,3 Gábor Csire2,4 Kyungwha
Park5 Balázs Ujfalussy1,2 1Department of Theoretical Physics, Budapest
University of Technology and Economics, Budafoki út 8., HU-1111 Budapest,
Hungary  2Wigner Research Centre for Physics, Institute for Solid State
Physics and Optics, H-1525 Budapest, Hungary  3MTA-BME Condensed Matter
Research Group, Budapest University of Technology and Economics, Budafoki út
8., HU-1111 Budapest, Hungary  4Catalan Institute of Nanoscience and
Nanotechnology (ICN2), CSIC, BIST, Campus UAB, Bellaterra, Barcelona, 08193,
Spain  5Department of Physics, Virginia Tech, Blacksburg, Virginia 24061, USA

###### Abstract

We present a fully relativistic first principles based theoretical approach
for the calculation of the spectral properties of magnetic impurities on the
surface of a superconducting substrate, providing a material specific
framework for the investigation of the Yu–Shiba–Rusinov (YSR) states. By using
a suitable orbital decomposition of the local densities of states we discuss
in great details the formation of the YSR states for an Mn adatom and for two
kinds of Mn dimers placed on the Nb(110) surface and compare our results to
recent experimental findings. In case of the adatom we find that the spin-
orbit coupling slightly shifts some of the YSR peaks and also the local spin-
polarization on the Nb atoms have marginal effects to the their positions.
Moreover, by scaling the exchange field on the Mn site we could explain the
lack of the $d_{x^{2}-y^{2}}$-like YSR state in the spectrum. While our
results for a close packed ferromagnetic dimer are in satisfactory agreement
with the experimentally observed splitting of the YSR states, in case of an
antiferromagnetic dimer we find that the spin-orbit coupling is not
sufficiently large to explain the splitting of the YSR states seen in the
experiment. Changing the relative orientation of the magnetic moments in this
dimer induces splitting of the YSR states and also shifts their energy,
leading even to the formation of a zero bias peak in case of the deepest YSR
state.

##  I Introduction

Local magnetic moments inside or at surfaces of metals in the superconducting
state lead to the formation of localized bound states within the
superconducting gap, referred to as Yu-Shiba-Rusinov (YSR) states 1 ; 2 ; 3 .
In the simplest picture, the magnetic moment is viewed as a classical spin
which is exchange-coupled to the Cooper pairs, and a single particle-hole
symmetric pair of YSR states is predicted within the gap 3 ; 4 ; 5 ; 6 . With
the high resolution of scanning tunneling microscopy (STM) and spectroscopy
(STS) techniques available nowadays, experiments are able to resolve several
pairs of YSR resonances within the superconducting gap 9 ; 10 ; 11 . This has
been explained due to the crystal field of the substrate which lifts the
degeneracy of the adatom’s $d$ level and, consequently, also the degeneracy of
the YSR states 6 ; 9 .

While these studies revealed the mechanism of the formation of the YSR states
in the magnetic impurity-superconductor systems in many details, a materials
specific description is still missing. In the present paper, we adopt and
further develop the solution of the Bogoliubov–de Gennes (BdG) equations with
band theoretical methods treating, therefore, magnetism, relativistic (spin-
orbit coupling) effects, superconductivity and the self-consistent electronic
structure of the host and the impurity on the same footing 100 ; 101 . Such a
solution incorporates the full orbital character of the electronic structure
of real impurities and is not only able to address both the general and the
specific properties of multi-orbital YSR states, but also allows a
straightforward comparison to local experimental probes having direct access
to the local density of states (LDOS).

Understanding inhomogeneous superconductors at the ab initio  level is rather
challenging, even for conventional electron-phonon coupled systems.
Previously, some of the present authors demonstrated that by solving the
Dirac-Bogoliubov–de Gennes (DBdG) equations via Multiple Scattering Theory
(MST) is a very powerful tool to explain induced superconductivity101 ; 104 ,
induced triplet pairing in superconductorscsire-relativistic2018  or the
effects of Fermi surface anisotropy 107 . Almost twenty years ago, the
Embedded Cluster Method (ECM) has been employed with the MST to describe
finite magnetic clusters of atoms in the normal state lazarovits-2002 .
Recently, the non-relativistic BdG-MST and the ECM have been combined and used
to study bulk impurity-superconductor systems saunderson-2021 .

In this work we extend this method by merging the fully relativistic DBdG-MST
equations with the ECM to perform a first principles study of the YSR states
induced by magnetic impurities at superconducting surfaces. In the next
section we briefly describe the elements of the theory and our computational
implementation. Then we apply the newly developed code to an Mn adatom and Mn
dimers on top of the (110) surface of bcc Niobium and we compare or results
with the experimental observationsbeck-natcomm2021 . By decomposing the local
density of states (LDOS) according to atomic orbitals we perform a thorough
analysis of the YSR states and also investigate the effect of the spin-orbit
coupling, of the local spin-polarization in the substrate and of the strength
of the exchange field at the Mn site. We also consider two kinds of close
packed Mn dimers, one with ferromagnetic (FM) and the other with
antiferromagnetic (AFM) coupling. We analyze the hybridization and
corresponding splitting of the YSR states in the FM coupled dimer. Our
calculations show that this effect is strongly suppressed in the AFM dimer
even in the presence of SOC. Finally, we studied the effect of the
noncollinearity of the spin moments in the AFM dimer and, in case of the
deepest YSR state, we show the formation of a zero bias peak (ZBP) at a
certain tilting angle.

##  II Methods

###  II.1 The Dirac-Bogoliubov–de Gennes equation

Following Ref. csire-relativistic2018,  we seek for the solutions of the
Dirac–Bogoliubov–de Gennes (DBdG) equation,

$\left[\varepsilon-\left(\begin{array}[]{cc}H_{\text{D}}&\Delta_{\text{eff}}\eta\\\
\Delta^{\ast}_{\text{eff}}\eta^{T}&-H^{\ast}_{\text{D}}\\\
\end{array}\right)\right]\left(\begin{array}[]{c}\Psi_{e}(\varepsilon)\\\
\Psi_{h}(\varepsilon)\end{array}\right)=0\;,$ (1)

where $\varepsilon$ is the energy relative to the Fermi energy $E_{F}$,
$\Psi_{e}(\varepsilon)$ and $\Psi_{h}(\varepsilon)$ are the four-component
electron and hole part of the wavefunction, respectively, $\Delta_{\rm eff}$
is the effective pair interaction, while

$H_{D}=c\vec{\alpha}\vec{p}+(\beta-\mathbb{I}_{4})c^{2}/2+\left(V_{\text{eff}}-E_{F}\right)(\vec{r})\mathbb{I}_{4}+\vec{\Sigma}\vec{B}_{\text{eff}}(\vec{r})\,,$
(2)

is the Kohn–Sham–Dirac Hamiltonian in Rydberg units ($\hbar=1$, $m=1/2$,
$e^{2}=2$, $c=274.072$) with the effective atomic potential
$V_{\text{eff}}(\vec{r})$ and effective exchange field
$\vec{B}_{\text{eff}}(\vec{r})$. The $4\times 4$ matrices in Eqs. (1) and (2)
are defined as

$\vec{\alpha}=\left(\begin{array}[]{cc}0&\vec{\sigma}\\\
\vec{\sigma}&0\end{array}\right)\,,\quad\beta=\left(\begin{array}[]{cc}\mathbb{I}_{2}&0\\\
0&-\mathbb{I}_{2}\end{array}\right)\,,$
$\vec{\Sigma}=\left(\begin{array}[]{cc}\vec{\sigma}&0\\\
0&\vec{\sigma}\end{array}\right)\,,\quad\eta=\left(\begin{array}[]{cc}i\sigma_{y}&0\\\
0&i\sigma_{y}\end{array}\right)\,,$

where $\vec{\sigma}$ stand for the Pauli matrices. As described by for example
in Ref. 101, , first $V_{\text{eff}}(r)$ and $\vec{B}_{\text{eff}}(r)$ are
calculated self-consistently in the normal state, then the solution of (1)
proceeds by either fitting the $\Delta_{\rm eff}$ parameter to the
experimental value of the energy gap in the bulk superconductor, or, by
assuming a suitable approximate energy functional and completing the self-
consistency cycle in the superconducting state 107 . In the present paper we
choose the simpler first route by fitting $\Delta_{\rm eff}$ to the Nb bulk
superconducting gap. This seems to be a reasonable approximation, because, on
the one hand, a magnetic adatom is not expected to support any singlet pairing
interaction, therefore $\Delta_{\rm eff}=0$ can be assumed on the impurity
sites, while, on the other hand, the pairing in the superconducting host is
not going to be influenced by a few impurity atoms.

###  II.2 Embedded Cluster Green’s function technique

In order to describe an impurity on a superconducting surface, it is desirable
to avoid the supercell approximation used frequently in normal-state
calculations. Namely, the superconducting coherence length is usually very
large compared to the lattice constant, thus it is expected that the
interference between supercells might seriously affect the obtained results.
The Screened Korringa-Kohn-Rostoker (SKKR) method used to solve the DBdG
equations for layered systems in Ref. csire-relativistic2018,  provides an
excellent basis to achieve this goal. First, the Green’s function of a semi-
infinite surface is calculated. The second step is the description of the
surface impurity via the Green’s function embedding technique lazarovits-2002
. This procedure trivially leads to the Green’s function of the impurity-
superconductor system,

$\displaystyle G^{\rm
ab}(\varepsilon,\vec{r}+\vec{R}_{i},\vec{r}^{\;\prime}+\vec{R}_{j})=\sum_{QQ^{\prime}}Z^{{\rm
a},R}_{i,Q}(\varepsilon,\vec{r})\tau^{{\rm
ab},ij}_{QQ^{\prime}}(\varepsilon)Z^{{\rm
b},L}_{j,Q^{\prime}}(\varepsilon,\vec{r}^{\;\prime})$
$\displaystyle-\delta_{ij}\left\\{\theta(r^{\prime}-r)\sum_{Q}Z^{{\rm
a},R}_{i,Q}(\varepsilon,\vec{r})J_{i,Q}^{{\rm
b},L}(\varepsilon,\vec{r}^{\;\prime})\right.$
$\displaystyle\quad\quad+\left.\theta(r-r^{\prime})\sum_{Q}J^{{\rm
a},R}_{i,Q}(\varepsilon,\vec{r})Z_{i,Q}^{{\rm
b},L}(\varepsilon,\vec{r}^{\;\prime})\right\\}$ (3)

where ${\rm a,b}\in\left\\{{\rm e},{\rm h}\right\\}$, $Q$ and $Q^{\prime}$ are
total angular momentum indices, $i$ and $j$ denote sites with position vectors
$\vec{R}_{i}$ and $\vec{R}_{j}$, respectively, while $Z^{{\rm a},R/L}_{i,Q}$
and $J^{{\rm a},R/L}_{i,Q}$ are in order properly normalized regular and
irregular right-hand-side(R) and left-hand-side(L) solutions of the DBdG
equation (1) csire-relativistic2018 . In this equation the scattering due to
the impurity is taken into account by replacing the inverse scattering
matrices ${\bm{t}}_{\text{host}}^{-1}(\varepsilon)$ of the host atoms by those
of the impurity or a of a cluster of impurities
${\bm{t}}_{\text{clus}}^{-1}(\varepsilon)$ to obtain the matrix of the
scattering path operator (SPO) at the impurity sites,

${\bm{\tau}}_{\text{clus}}(\varepsilon)={\bm{\tau}}_{\text{host}}(\varepsilon)\left\\{{\bm{I}}-\left[{\bm{t}}_{\text{host}}^{-1}(\varepsilon)-{\bm{t}}_{\text{clus}}^{-1}(\varepsilon)\right]{\bm{\tau}}_{\text{host}}(\varepsilon)\right\\}^{-1}.\;$
(4)

Here the bold symbols denote matrices with both angular momentum and site
indices, whereas the latter ones are restricted to the sites of the impurity
cluster. For the details of the embedded cluster method see for example Ref.
lazarovits-2002, . Note that the solution of the DBdG equation instead of the
Dirac equation implies a doubled matrix dimension for the SPO matrix due to
the appearance of the electron-hole index,

${\bm{\tau}}(\varepsilon)=\left(\begin{array}[]{cc}{\bm{\tau}}^{\rm
ee}(\varepsilon)&{\bm{\tau}}^{\rm eh}(\varepsilon)\\\ {\bm{\tau}}^{\rm
he}(\varepsilon)&{\bm{\tau}}^{\rm hh}(\varepsilon)\end{array}\right)\;.$ (5)

Once the SPO for the impurity system is obtained, the Green’s function can be
calculated according to Eq. (3). It should be noted that the embedding process
via Eq. (4) requires the calculation of the site off-diagonal elements of the
SPO. This is the most time consuming step of the calculation, however, it
enables the embedding of a significantly larger cluster of impurities as
compared to the ordinary size of the supercells in supercell methods. Once the
Green’s function is calculated, the local density of states (LDOS) averaged
inside the atomic cell $V_{i}$ can be obtained as usual, albeit now including
new component indices for the electrons and the holes,

$n_{i}^{\rm
ab}(\varepsilon)=-\dfrac{1}{\pi}\text{Im}\int_{V_{i}}d^{3}r\,\text{Tr}\,G^{\rm
ab}(\varepsilon,\vec{r}+\vec{R}_{i},\vec{r}+\vec{R}_{i})\;,$ (6)

where the Tr denotes the trace in the four-component Hilbert space. In this
work we will present electron and hole LDOS’s, $n_{i}^{\rm ee}(\varepsilon)$
and $n_{i}^{\rm hh}(\varepsilon)$, respectively, and the sum of them termed as
the total LDOS.

##  III Results

In the present work we apply the theory described above to an Mn adatom and
various Mn dimers on the Nb(110) surface, studied in detail recently by STM
experiments and also theoretically using a tight-binding model beck-
natcomm2021 . The self-consistent calculations in the normal state were
performed by employing the local density approximation according to Vosko et
al. vosko-1980  and the atomic sphere approximation (ASA) for the description
of the single-site effective potentials and exchange fields, while the
scattering of the electrons was treated using an orbital momentum cutoff of
$\ell_{\text{max}}=2$. For the energy integration we used a semi-circular
energy contour on the upper half plane with 16 energy points. For the self-
consistent calculation of the Nb(110) surface layers, forming two-dimensional
(2D) centered rectangular lattices with $C_{2v}$ point-group symmetry, we
sampled 253 $k$ points in the irreducible wedge of the 2D Brillouin zone (BZ).
By neglecting geometrical relaxions, both in the absence and the presence of
impurities, the lattice positions were considered in an ideal parent bulk
geometry with a lattice constant of $a=6.237\,a_{0}$ straumanis-1970 . The
monomer and the dimer atoms were then placed in hollow positions above the
surface layer of Nb(110), i.e. in the lattice positions of the first vacuum
(empty sphere) layer.

To take into account charge- and spin-density oscillations caused by the
impurity in the substrate, beyond the magnetic impurities we inlcuded
substrate atoms and also empty spheres representing the vacuum in the self-
consistently treated embedded cluster within a radius of $r=1.21\,a_{\rm 2D}$
($a_{\rm 2D}=\sqrt{2}a$), often referred to as an extended Impurity Cluster
(IC). By solving the DBdG equations for the extended impurity we calculated
the LDOS for each embedded atom according to Eq. (6). In this final step we
increased the size of the IC to $3a_{\rm 2D}$ (including e.g. 118 atoms in the
case of a single impurity), to explore the spatial extent of the YSR states.
During the calculations in the superconducting state, we assumed an effective
pair potential of $\Delta_{\rm eff}=1.51$ meV beck-natcomm2021  on the Nb
atoms and $\Delta_{\rm eff}=0$ on the Mn impurities.

###  III.1 Mn monomer on Nb(110)

In the case of the Mn monomer, we considered a perpendicular direction of the
magnetization and performed the calculations of the self-consistent normal
state and of the LDOS in the superconducting gap as outlined above. Motivated
by the simple tight-binding analysis beck-natcomm2021  we decomposed both the
electron and hole LDOS according to orbitals related to the irreducible
representations of the $C_{2v}$ point group. As the $s$ and $p$ orbitals
provide negligible contributions to the LDOS in the SC gap, in Fig. 1 we only
show the calculated $d$-like electron and hole contributions to the LDOS for
the Mn impurity, both in the absence and in the presence of spin-orbit
coupling. The former one is possible by manipulating the Dirac equation as
proposed originally by Ebert at al. ebert . In both cases, four pairs of sharp
peaks in the LDOS are seen that clearly can be associated with Yu-Shiba-
Rusinov states induced by the magnetic impurity. It is obvious that the
electron-hole symmetry of the DBdG equations is well reflected in our results.
One can also observe that the states below the Fermi level are mostly hole-
like, while the symmetric states with positive energy are mostly electron-
like.

Figure 1: The LDOS of the Mn adatom on the Nb(110) surface decomposed
according to canonical $d$ orbitals without (top) and with (bottom) spin-orbit
coupling. The solid lines correspond to electron densities and the dashed
lines to hole densities. The total DOS summed for all orbitals and for the
electrons and holes is plotted with dash-dotted line. For a better visibility,
the curves are shifted with respect to each other.

To get a deeper insight into the symmetry of the YSR states let us first
consider the case without SOC. From the orbital resolution of the LDOS we
conclude that three YSR states can be assigned with single $d$ orbitals:
$d_{yz}$ at $\pm$0.98 meV, $d_{xz}$ at $\pm$1.1 meV and $d_{xy}$ at $\pm$1.2
meV, while the peak at $\pm$1.3 meV having dominant $d_{z^{2}}$ character also
has a small $d_{x^{2}-y^{2}}$ contribution. Since in the absence of SOC the
Hamilton operator is diagonal in spin space, it is not surprising to recognize
that the observed orbital decomposition of the YSR states reflects the (one-
dimensional) irreducible representations of the simple $C_{2v}$ point group:
the first three $d$ orbitals correspond to different irreducible
representations, while both the $d_{z^{2}}$ and the $d_{x^{2}-y^{2}}$ orbitals
belong to the total symmetric $A_{1}$ representation, therefore, they can
hybridize. According to this reasoning there should exist another
$d_{z^{2}}-d_{x^{2}-y^{2}}$ peak with $d_{x^{2}-y^{2}}$ dominance. Most
possibly, this peak can not be seen as it might form outside the coherence
region $(E_{F}-\Delta,E_{F}+\Delta)$, see later.

By switching on the SOC, the position of some of the peaks apparently changes:
the $d_{yz}$ and the $d_{xy}$ peaks move closer to $E_{F}$ by about 0.2 meV
and 0.1 meV, respectively, while the $d_{xz}$ and $d_{z^{2}}-d_{x^{2}-y^{2}}$
peaks roughly keep their positions in energy. As a result, the order of the
$d_{xy}$ and $d_{xz}$ peaks becomes reversed. It is also expected that the SOC
mixes both the spin and orbital characters of the YSR states. However, such
mixing of the orbital characters of the YSR states is hardly visible in the
lower panel of Fig. 1: only a minor peak in the $d_{xz}$ LDOS can be inferred
at the position of the large $d_{yz}$ peak at $\pm 0.75$ meV due to the spin-
orbit coupling between states with the orbital quantum numbers $m_{\ell}=\pm
1$. We repeated the DBdG simulations in the presence of SOC also with the
effective potentials and fields obtained from the self-consistent calculations
without SOC and found almost the same positions for the YSR states as without
SOC in both steps of the calculations (upper panel of Fig. 1). From this
observation we conclude that the shift of the YSR states due to the SOC (lower
panel of Fig. 1) originates from its indirect effect in changing the self-
consistent potentials and fields.

As mentioned already, the magnetic impurity causes a spatial perturbation in
the charge and magnetization density in the normal state, which manifests
itself also in the formation of local spin-polarization on the atoms around
the impurity. Within the local density approximation (LDA) of density
functional theory (DFT), this local spin-polarization is represented by the
local exchange field being zero for a non-magnetic atom. In order to
investigate the role that the local spin-polarization in the nonmagnetic host
plays in the formation of the YSR states we performed a series of calculations
including an increasing number of neighbor shells in the IC with local
exchange fields switched on. In practice, in all calculations we used the same
IC containing four neighbor shells and we switching off the selfconsistently
calculated exchange field at the same site.

We start from a zero shell neighborhood meaning that at all sites in the IC
but at the Mn impurity the exchange field is switched off and then we turn on
the exchange fields determined from the normal-state self-consistent
calculations shell by shell. The results are shown in Fig. 2, where we plot
the positions of the YSR peaks as a function of atomic shell number for which
the exchange fields were included in the calculations. It can be seen that
even in the zero shell calculation we get a rather accurate position of the
YSR peaks, which changes only a little as we include the first shell. There is
an even smaller change with the inclusion of the farther nearest-neighbor
shells. As can be seen from the figure, the above findings apply irrespective
of whether the SOC is switched on or off, which is not surprising, since the
formation of the local spin-polarization is usually insensitive of the SOC.

Figure 2:  The energy of the YSR states below $E_{F}$ of the Mn adatom on
Nb(110) as a function of the number of atomic shells around the adatom formed
by host atoms with nonzero exchange field included in the calculations. The
solid and dashed lines stand for the cases with and without spin-orbit
coupling, respectively.

When comparing our results with the experimental findings in Ref. beck-
natcomm2021, , one can easily see that regarding both the number of peaks and
their positions they agree favorably. The only peak which appears to be off in
position is the one closest to the center of the gap, being shifted
considerably towards the gap edge in our calculations. In order to determine
the orbital character of the YSR states the authors of Ref. beck-natcomm2021,
relied on the study of the spatial symmetry of the scanning tunneling
spectroscopy (STS) pictures. As discussed in context of Fig. 1, in the current
formalism it is straightforward to decompose the LDOS at the Mn impurity
according to orbitals corresponding to the irreducible representations of the
point group of the system, in the case of Nb(110), of the $C_{2v}$ group. The
question then arises whether this decomposition can be correlated with the
spatial shape of the STS pictures. To answer this question we solved the DBdG
equations for a large number of atoms in the first vacuum layer around the
impurity site and calculated the LDOS on each site in this IC. Due to the
atomic sphere approximation used in our method, a single averaged value of the
LDOS can be attributed to each atomic cell. Since such a visualization of the
LDOS would be rather course grained, we applied some interpolation to the data
within the atomic cells.

We presented the obtained LDOS maps in Fig. 3 at the energies of the YSR peaks
seen in Fig. 1 below the Fermi energy. These patterns reflect the spatial
distribution, in particular, the symmetry of the corresponding YSR states. For
example, the LDOS map of the $d_{xy}$ YSR peak at $-1.01$ meV displays zero
values at the $y-z$ ($x=0$) and the $x-z$ ($y=0$) planes, thus it can safely
be associated with the $-\beta$ state of $d_{xy}$ character in Fig. 1. of Ref.
beck-natcomm2021, . Similarly, the $d_{yz}$ state at $-0.58$ meV, the $d_{xz}$
state at $-0.96$ meV and the $d_{z^{2}}$ state at $-1.28$ meV can be directly
related to the $-\delta$, $-\gamma$ and $-\alpha$ states, respectively,
detected experimentally and identified also in terms of tight-binding
simulationsbeck-natcomm2021 . Thus, our calculations reproduce all the YSR
states of a Mn adatom on Nb(110) with the correct order in energy and with the
correct orbital character seen from STS experimentsbeck-natcomm2021 .

Figure 3: The color map of the LDOS in arbitrary units in the vacuum layer
containing the Mn adatom at the energies corresponding to the four YSR states
identified in Fig.1. The LDOS in a 2D unit cell is associated with the
averaged LDOS from the corresponding atomic sphere. The LDOS of the Mn atom is
replaced by the average of the LDOS from the first nearest neighbor vacuum
cells. To get a smooth picture in space an interpolation scheme was used. The
position of the Mn atom is marked by a small red circle. The white bar denotes
a lattice spacing along the [100] direction ($a_{\rm 2D}$).

As proposed originally by Shiba 2 , the YSR states are known to vary in energy
as a function of the atomic spin ($S$). In order to study this relationship
within a first principles multi-orbital scenario, we performed a computer
experiment, where we artificially scaled the exchange field of the Mn atom and
calculated the LDOS in the superconducting state as a function of the scaling
factor, $B/B_{0}$. Because of the missing self-consistency, we simultaneously
set the exchange fields to zero for all the atoms in the extended impurity
cluster. Note that a more sophisticated scaling could be obtained by the so-
called ”fixed moment search”, where the size of the impurity’s spin-moment is
constrained, while the effective potentials and fields, as well as the induced
exchange fields are recalculated self-consistently.

Figure 4: The energy of all the five $d$-like YSR states for the Mn adatom on
Nb(110) as a function of the scaling of the exchange field, $B/B_{0}$, on the
Mn atom. The calculations have been performed without SOC.

Fig. 4 shows the energy of the YSR states for $0<B/B_{0}<2$ using the
simplified scaling picture described above. One can see that the peaks enter
into the energy range of the superconducting gap at rather low fields and as
we increase the exchange field they cross the Fermi energy and then leave the
gap. This trend is compatible with the original prediction of Shiba, even
though a direct comparison is quite difficult. One important feature of the
figure that the scattering chanel with $x^{2}-y^{2}$ orbital character appears
only in a rather narrow range of the exchange field, and at the actual self-
consistent value ($B/B_{0}=1.0$) it does not provide any contribution to the
density of states within the gap. This is in perfect agreement with the
experimental findings of Beck et al. beck-natcomm2021 , who could not find a
peak with $x^{2}-y^{2}$ character.

###  III.2 Mn dimers on Nb(110)

In the case of two impurities forming a dimer, three adjacent positions were
reported in Ref. beck-natcomm2021,  that can be characterized by the direction
of the vectors connecting the two Mn atoms, namely, the nearest-neighbor (NN)
position by [100], the second NN neighbor position by [1$\bar{1}0$], and the
third nearest NN position by [1$\bar{1}$1]. Here we consider the [1$\bar{1}0$]
the [1$\bar{1}$1] dimers only. Based both on experiments on short Mn chains on
Nb(110)schneider-2021  and on ab-initio calculations laszloffy-2021  it is
suggested that the coupling of the Mn atoms in the [1$\bar{1}$0] dimer is
ferromagnetic (FM), while in the [1$\bar{1}1$] dimer it is antiferromagnetic
(AFM).

Figure 5: LDOS in the superconducting state calculated without and with SOC
for the ferromagnetic [1$\bar{1}$0] and the antiferromagnetic [1$\bar{1}$1] Mn
dimers on Nb(110). Solid and dashed lines stand for the electron and hole part
of the LDOS.

The LDOS for these two dimers is plotted in Fig. 5 both without and with SOC
taken into account in the calculations. Let’s consider first the results
without SOC. The LDOS for the FM [1$\bar{1}$0] dimer shows a complex structure
formed by eight peaks each for negative and positive energies upon electron-
hole symmetry, which is twice the number of the peaks as seen in the case of a
single impurity. The number of these states and their positions are again in
reasonable agreement with experimentbeck-natcomm2021 , though - similarly to
the single impurity case - the peaks appear somewhat shifted towards the gap
edges. On the contrary, in the case of the AFM [1$\bar{1}$1] dimer, the YSR
states appear in a very similar way as for a single impurity density of states
on Fig. 1 and the doubling of the number of YSR states can not be observed.
This is in contrast with the experiment where six peaks were detectedbeck-
natcomm2021 .

As well-known, if the localized YSR states of the two impurities overlap, they
can hybridize leading to the formation of symmetric or bonding and anti-
symmetric or anti-bonding pairs, where the symmetric/antisymmetric one has an
enhanced/zero density at the mirror plane between the two impurities. However,
at least within the a non-relativistic theory3 , such a hybridization and
splitting of the YSR states occurs only in FM-coupled dimers, while the YSR
states remain degenerate in an AFM-coupled dimer as reported also
experimentally flatte ; morr .

Figure 6: Upper panel: The decomposition of the LDOS below $E_{F}$ to $d$-like
orbitals for the FM [1$\bar{1}$0] dimer. The classification of the split YSR
peaks is also displayed as obtained from the simulated STS pictures. Lower
panels: Color map of the simulated STS pictures (see Fig. 3) corresponding to
the eight YSR peaks marked in the upper panel. The left/right four entries
refer to the antisymmetric/symmetric YSR states of the dimer. The positions of
the Mn atoms are marked by small red circles. The white bar denotes a lattice
spacing along the [100] direction ($a_{\rm 2D}$).

To analyze the hybridizations and the splitting of the YSR states for the FM
[1$\bar{1}$0] dimer in more details, in the top panel of Fig. 6 we plot the
decomposition of the LDOS into $d$-like contributions near the lower edge of
the SC gap. In addition, we plotted the color maps of the LDOS in the vacuum
layer containing the Mn dimer, see also Fig. 3 for the single impurity. Based
on these figures, the peaks at -1.22 meV and -1.32 meV with large $xy$
contributions could easily be identified as the symmetric and antisymmetric
components of the split $xy$ YSR state, denoted by $xy_{s}$ and $xy_{a}$,
respectively. Similarly, the $xz_{a}$ and $xz_{s}$ states can be
straightforwardly recognized because of their overwhelming $xy$ contributions
and the shapes in the lowest two panels in Fig. 6. The peak at -0.59 meV has
dominant $yz$ character and from the corresponding entry in Fig. 6 it can be
identified as the bonding $yz_{s}$ state. At -1.27 meV one can see a peak with
dominant $z^{2}$ contribution, somewhat hybridized with $yz$ states.
Nevertheless, based on the respective color map in Fig. 6 we can safely
associate it with the antisymmetric $z^{2}_{a}$ state. There are two peaks
remaining at -1.36 meV and -1.14 meV, both of them containing almost equal
$yz$ and $z^{2}$ orbital contributions. From the corresponding color maps it
is clear that they display an antisymmetric and a symmetric hybridized state,
referred to as the $yz_{a}$ and $z_{s}^{2}$ states, respectively.

The above analysis demonstrates that it is possible to classify the split
multi-orbital YSR states of the FM dimer derived from our first principles
theory and they show similar patterns as the experimentsbeck-natcomm2021 . In
contrast to the experiments, in our calculations the symmetric states are
higher in energy than the corresponding antisymmetric ones. However, note that
for their positive-energy counterparts containing mostly electron states, this
relationship is reserved. It is also worthwhile to mention that the splitting
of the $xy$ ($\beta$) states could not be resolved in the experiment.

Inferring the LDOS’s in Fig. 5 with SOC turned on, it is obvious that for both
kind of dimers the SOC only slightly alters the positions of the YSR states.
In case of the FM coupled [1$\bar{1}$0] dimer it seems that additional
splittings happen due to the SOC, but a more detailed investigation shows that
this feature happens due to non-uniform shifts of the peak positions with
different orbital characters. In case of the AFM [1$\bar{1}$1] dimer the
changes caused by the SOC are hardly visible. In sharp contrast with the
theoretical explanation in Ref. beck-natcomm2021, , we conclude that the spin-
orbit coupling in this system is not strong enough to explain the splitting of
the YSR states found experimentally for the AFM dimer. This conclusion is
consistent with our previous findings for the Mn adatom, where we also found a
rather weak effect of SOC, and explains the similarity between the LDOS of the
single impurity and the AFM coupled dimer as it was put forward earlier in the
context of a non-relativistic theory3 .

The theoretical study of Lászlóffy et al.laszloffy-2021  suggests that a
slightly non-collinear and/or canted configuration can exist for the dimers in
the normal state. The possibility of the formation of such tilted ground
states and the lack of split YSR states in the AFM dimer compared to
experiments motivated us to study the effect of the non-collinear spin-
configurations to the YSR states. In these calculations we kept fixed the
effective potentials and fields obtained selfconsistently in the collinear
state and rotated the direction of one of the Mn moments continuously between
the AFM ($\theta=180^{\circ}$) and the FM ($\theta=0^{\circ}$) states. In Fig.
7 we show the change of the energy of the YSR states for the [1$\bar{1}$1]
dimer as a function of the opening angle $\theta$.

As also shown based on tight-binding calculations in Ref. beck-natcomm2021, ,
a deviation from the collinear magnetic configuration lifts the degeneracy of
the YSR states present in the collinear AFM dimer. In Table 1 we presented the
splittings calculated for the relative angles $170^{\circ}$ and $160^{\circ}$,
together with the those found experimentallybeck-natcomm2021 . Apparently, for
these angles the calculated values are in the range of the measured ones, with
exception of the splitting of the $\gamma$ state. Nevertheless, such big
deviations from the AFM spin configuration seem to be unrealistic based on the
spin-model simulations of Ref. laszloffy-2021, , where for the NN
[1$\bar{1}$1] dimer a $\theta\simeq 179.5^{\circ}$ was predicted.

Another important feature that can be seen in Fig. 7 that one line of the
split $d_{yz}$ $(\delta)$ states moves rapidly towards the Fermi level, while
the energies of the other YSR states are shifted more moderately when going
from the AFM to the FM configuration. At roughly $\theta=40^{\circ}$ this peak
can be found right in the middle of the energy range, namely at zero bias.
Such a behavior of the upper lying $yz$ peak has been also found in Ref. beck-
natcomm2021,  using a multi-orbital tight-binding approach, though, because of
the much deeper position of the $\delta$-states, already at about
$\theta=140^{\circ}-150^{\circ}$ depending on whether the SOC was turned off
or on. Zero bias peaks (ZBP) have the property that they are necessarily have
an equal electron-like and hole-like component in the LDOS. Such zero bias
peaks in extended chains are of primary interest nowadays because they are
supposed to realize topological quantum bitsbeenakker-2011 ; alicea-2012 ;
elliott-2015 . Quite interestingly, theoretical calculations on the same
platform as in this work reveal that ZBP’s can be formed in single magnetic
impurities on top of superconducting Pb(110) by tilting the magnetic
orientation of the impurity with respect to the surface kyungwha-2021 . The
present calculations demonstrate that ZBP’s can be created with the
manipulation of the relative angle in dimers. It is therefore expected that
non-collinearity, e.g. spin-spiral states, can induce topological ZBP’s in
long chains.

Figure 7: The energy of the YSR states in the [1$\bar{1}$1] Mn dimer on Nb(110) as a function of the relative angle $\theta$ between the two Mn moments. Note that the FM and the AFM states correspond to $\theta=0$ and $\theta=180^{\circ}$, respectively. Table 1: Calculated energy splittings of the YSR states in the [1$\bar{1}$1] Mn dimer on Nb(110) for Mn moments with a relative angle $\theta$. The last line contains the experimental values from Ref. beck-natcomm2021, .  $\theta$ |  $z^{2}$ ($\alpha$) |  $xy$ ($\beta$) |  $xz$ ($\gamma$) |  $yz$ ($\delta$)
---|---|---|---|---
170∘ | 37 | 70 | 41 | 19
160∘ | 84 | 40 | 39 | 53
exp | 130 | 400 | - | 110

##  IV Summary and conclusions

In order to provide a material specific description of magnetic impurities on
superconducting surfaces, in the present paper we generalised the Embedding
Method of Multiple Scattering Theory to solve the Dirac-Bogoliubov-de Gennes
equations. By using the new method, first-principles band structure,
relativistic effects and superconductivity is treated on the same footing. In
particular, the multi-orbital treatment of the in-gap states is a priori
inherently involved in the theory.

For a single Mn impurity placed on Nb(110) we found a series of Yu-Shiba-
Rusinov states in the superconducting gap, in excellent agreement with
experiments both concerning the energy and the orbital decomposition of these
states. Moreover, we were able to prove that the orbital character of the YSR
states inferred from the LDOS of the Mn impurity is directly related to the
spatial shape and symmetry of the corresponding simulated STS pictures. Our
calculations demonstrate that the spin-orbit coupling has quite a little
effect in shifting the energy of these states. We also showed that the local
spin-polarization on the magnetic impurity beyond the first coordination shell
influence negligibly the energy of the YSR states. In line with earlier
theoretical predictions, we also demonstrated that the positions of the YSR
states are very sensitive on the size of the magnetic moments, which can be
used to explain the missing in-gap state with $d_{x^{2}-y^{2}}$ orbital
character in case of the Mn adatom on Nb(110).

In case of the close packed [1$\bar{1}$0] ferromagnetic Mn dimer, our
calculations demonstrated the well-established splitting of the YSR states due
to hybridization. By analyzing the symmetry of the simulated STS maps we
managed to identify the orbital character of the split states and found a
satisfactory agreement with the experiment. Quite interestingly, however, we
obtained a reversed order for the energies of the bonding and antibonding
states as compared to the experiment.

In good agreement with the prediction of non-relativistic theories we found no
splitting of the YSR states for the antiferromagnetic [1$\bar{1}$1] Mn dimer
when the SOC was switched off in the calculations. By including the SOC, we
observed very small changes in the energy of the YSR states, in sharp
contradiction with the experiment and related tight-binding calculations
including SOC which showed a sizeable splitting of the YSR states also in this
case. Our first principles calculations thus indicate that the SOC is not
sufficiently strong to explain every detail of the experimental findings for
the AFM coupled [1$\bar{1}$1] Mn dimer.

As a possibility to resolve this contradiction, we investigated the effect of
the relative magnetic orientation in the dimer and found splittings of the YSR
states being comparable with the measured ones for about $20^{\circ}$
deviation from the AFM arrangement. Since such spin-configurations seem to be
unrealistic to form the ground state of the Mn dimer, further investigations
are needed to improve the agreement between theory and experiment. It is a
subject of future work to take into account e.g. the relaxation of the lateral
position of the impurities, where an increased hybridization with the
substrate is expected and a more accurate description of the dimer YSR states
can be obtained.

Remarkably, at about $40^{\circ}$ to the parallel alignment of the moments the
energy of the deepest YSR state reached the Fermi level, forming thus a zero
bias peak. This hints the possibility of creating zero bias peaks in longer
chains upon non-collinear arrangements of the spin.

###### Acknowledgements.

Among the authors B.Ny., A.L., L.Sz. and B.U. are grateful for the financial
support by the Ministry of Innovation and Technology and the Hungarian
National Scientific Research Fund under project No. K131938 and within the
Quantum Information National Laboratory of Hungary. B.Ny. also greatful for
the support of ÚNKP-20-3-II-BME-249 New National Excellence Program of the
Ministry for Innovation and Technology of Hungary. G.Cs. gratefully
acknowledges support from the European Union’s Horizon 2020 research and
innovation program under the Marie Sklodowska-Curie Grant Agreement No.
754510. This work was supported by Spanish MINECO (the Severo Ochoa Centers of
Excellence Program under Grant No. SEV- 2017-0706), Spanish MICIU, AEI and EU
FEDER (Grant No. PGC2018-096955-B-C43), and Generalitat de Catalunya (Grant
No. 2017SGR1506 and the CERCA Program). The work was also supported by the
European Union MaX Center of Excellence (EU-H2020 Grant No. 824143).

## References

  * (1) L. Yu, Acta Phys. Sin. 21, 75 (1965).
  * (2) H. Shiba, Prog. Theor. Phys. 40, 435 (1968).
  * (3) A. I. Rusinov, Zh. Eksp. Teor. Fiz. Pisma Red. 9, 146 (1968) [JETP Lett. 9, 85 (1969)].
  * (4) A. L. Fetter, Phys. Rev. 140, A1921 (1965).
  * (5) A. Yazdani, B. A. Jones, C. P. Lutz, M. F. Crommie, and D. M. Eigler, Science 275, 1767 (1997).
  * (6) A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
  * (7) S.-H. Ji, T. Zhang, Y.-S. Fu, X. Chen, X.-C. Ma, J. Li,W.-H. Duan, J.-F. Jia, and Q.-K. Xue, Phys. Rev. Lett. 100, 226801 (2008).
  * (8) M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B.W. Heinrich and K. J. Franke, Phys. Rev. Lett. 115, 087001 (2015).
  * (9) N. Hatter, B.W. Heinrich, M. Ruby, J. I. Pascual, and K. J. Franke, Nat. Commun. 6, 8988 (2015).
  * (10) G. Csire, B. Újfalussy, J. Cserti, and B. L. Gyorffy. Phys. Rev. B, 91, 165142, (2015).
  * (11) G Csire, A. Deák, B. Nyári, H. Ebert, J. Annett, and B. Ujfalussy. Phys. Rev. B, 97, 024514 (2018).
  * (12) B. Lazarovits, L. Szunyogh, and P. Weinberger. Phys. Rev. B, 65, 104441, (2002).
  * (13) G. Csire, B. Ujfalussy, and J. F. Annett. The European Physical Journal B, 91(10),217, (2018).
  * (14) G. Csire, J. Cserti, I Tütto and B. Ujfalussy. Phys.Rev. B, 94:104511, 2016.
  * (15) Gábor Csire, Stephan Schönecker, and Balázs Újfalussy. First-principles approach to thin superconducting slabs and heterostructures. Phys. Rev. B, 94:140502, Oct 2016.
  * (16) S. K. Ghosh, G. Csire, P. Whittlesea, J. F. Annett, M. Gradhand, B. Ujfalussy and J. Quintanilla. Phys. Rev. B, 101, 100506, (2020).
  * (17) T. G. Saunderson, J. F. Annett, B. Újfalussy, G. Csire, and M. Gradhand, Phys. Rev. B 101, 064510, 2020.
  * (18) T. G. Saunderson, J. F. Annett, G. Csire, and M. Gradhand, arXiv:2107.00237 (2021).
  * (19) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * (20) M. E. Straumanis, and S. Zyszczynski, J. Appl. Crystallogr., 3, 1 (1970)
  * (21) P. Beck, L. Schneider, L. Rózsa, K. Palotás, A. Lászlóffy, L. Szunyogh and R. Wiesendanger, Nature Communications, 12(1), 1-9 (2021).
  * (22) H. Ebert and H. Freyer and A. Vernes and G.-Y. Guo, Phys. Rev. B, 53 7721 (1996)
  * (23) Schneider, L., Beck, P., Wiebe, J. and Wiesendanger, R. Sci. Adv. 7, eabd7302 (2021).
  * (24) A. Lászlóffy, K. Palotás, L. Rózsa, and L. Szunyogh, Nanomaterials 11, 1933 (2021).
  * (25) M. E. Flatté and D. E. Reynolds, Phys. Rev. B 61, 14810–14814 (2000).
  * (26) D. K. Morr and N. A. Stavropoulos, Phys. Rev. B 67 020502 (2003).
  * (27) C. W. J. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113–136 (2011).
  * (28) J. Alicea, Reports Prog. Phys. 75, 76501 (2012).
  * (29) S. R. Elliott and M. Franz, Rev. Mod. Phys. 87, 137–163 (2015).
  * (30) K. Park, B. Nyári, A. Lászlóffy, L. Szunyogh, and B. Újfalussy, submitted for publication.


22footnotetext<EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding Author: Ranjini
Bandyopadhyay; Phone: +91 9480836210; Email<EMAIL_ADDRESS>
# Dichotomous behaviors of stress and dielectric relaxations in dense
suspensions of swollen thermoreversible hydrogel microparticles

Chandeshwar Misra  Soft Condensed Matter Group, Raman Research Institute, C.
V. Raman Avenue, Sadashivanagar, Bangalore 560 080, INDIA Paramesh Gadige
Department of Physics, Sri Sathya Sai Institute of Higher Learning, Prasanthi
Nilayam-515134, INDIA Ranjini Bandyopadhyay  Soft Condensed Matter Group,
Raman Research Institute, C. V. Raman Avenue, Sadashivanagar, Bangalore 560
080, INDIA

(August 27, 2024)

###### Abstract

Hypothesis: While the mechanical disruption of microscopic structures in
complex fluids by large shear flows has been studied extensively, the effects
of applied strains on the dielectric properties of macromolecular aggregates
has received far less attention. Simultaneous rheology and dielectric
experiments can be employed to study the dynamics of sheared colloidal
suspensions over spatiotemporal scales spanning several decades.

####

Experiments: Using a precision impedance analyzer, we study the dielectric
behavior of strongly sheared aqueous suspensions of thermoreversible hydrogel
poly($N$-isopropylacrylamide) (PNIPAM) particles at different temperatures. We
also perform stress relaxation experiments to uncover the influence of large
deformations on the bulk mechanical moduli of these suspensions.

####

Findings: The real parts of the complex dielectric permittivities of all the
sheared PNIPAM suspensions exhibit distinct relaxation processes in the low
and high frequency regimes. At a temperature below the lower consolute
solution temperature (LCST), both real and imaginary parts of the
permittivities of highly dense PNIPAM suspensions decrease with increase in
applied oscillatory strain amplitudes. Simultaneously, we note a counter-
intuitive slowdown of the dielectric relaxation dynamics. Contrary to our
rheo-dielectric findings, our bulk rheology experiments, performed under
identical conditions, reveal shear-thinning dynamics with increasing strain
amplitudes. We propose the shear-induced rupture of fragile clusters of
swollen PNIPAM particles to explain our observations.

####

Keywords : Rheo-dielectric; Thermoresponsive hydrogels; Segmental motion;
Counterion polarization; Interfacial polarization; Dense suspensions;
Rheology; Shear-thinning.

####

Abbreviations : LCST, Lower consolute solution temperature; PNIPAM,
Poly($N$-isopropylacrylamide); ERFs, Electrorheological fluids; DLS, Dynamic
light scattering; DSC, Differential scanning calorimetry; LAOS, Large
amplitude oscillatory strains; RPM, Revolutions per minute; SI, Supporting
Information; EP, Electrode polarization.

##  1 Introduction

Soft materials such as colloidal suspensions, gels, glasses, polymers and
biological macromolecules can deform easily under the influence of thermal
stresses at room temperature. When subjected to large external shear stresses
and strains, they display fascinating flow properties over a hierarchy of
spatiotemporal scales [1]. The application of appropriate shear profiles can
cause build-up or disruption of fragile suspension structures, in phenomena
that are commonly referred to as shear-induced thickening or thinning
respectively. It would be of interest to systematically study how a material’s
microscopic properties dictate its bulk flow behavior. Dielectric
spectroscopy, which probes how an externally applied oscillatory electric
field interacts with a material’s permanent or field-induced dipoles at
different oscillatory frequencies, is a widely used technique to study complex
metastable phenomena such as the dynamical slowing down process in glass-
forming liquids [2, 3]. The rheo-dielectric technique, which involves the
study of the dielectric response of a material in the presence of shear
forces, has been used to directly determine how shear-induced structural
changes alter charge distributions and polarization fluctuations in nematic
liquid crystals, polymers, electrorheological fluids (ERFs), viscoplastic
droplets and elastomers [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Simultaneous
measurements of dielectric and viscoelastic properties can therefore help us
to correlate material properties over length and time scales spanning several
decades.

####

Thermoresponsive poly($N$-isopropylacrylamide) (PNIPAM) hydrogel particles
when dispersed in water show swelling and deswelling behaviors respectively
below and above the lower consolute solution temperature (LCST, also referred
to as the volume phase transition temperature, VPTT) of $\approx$ 34∘C [14,
15]. Below the LCST, the hydrophilic amide groups of PNIPAM form hydrogen
bonds with water molecules such that individual PNIPAM particles absorb water
and swell significantly. Due to increased entropy above the LCST, water
molecules reorient around the non-polar regions of the PNIPAM macromolecules,
thereby restricting hydrogen bond formation with the constituent amide groups.
PNIPAM particles, therefore, exhibit hydrophobic behavior above their LCST and
shrink by expelling water [16]. As a consequence of this thermoreversible
swelling behavior, the volume fraction of PNIPAM particles in aqueous
suspensions can be controlled simply by altering the temperature of the medium
at a fixed particle concentration [17]. This offers an advantage over many
temperature-insensitive hard colloidal particles such as poly(methyl
methacrylate) (PMMA) and polystyrene (PS). Dense aqueous suspensions of PNIPAM
particles have therefore been recognized as excellent model colloidal systems
in the study of glass transition dynamics [18, 19]. Various protocols, for
example, synthesizing the particles in different ionic liquids or changing the
concentration and nature of crosslinkers during synthesis have been
implemented to alter the LCST of these particles [20]. Besides being
thermoresponsive, PNIPAM particles are deformable and compressible in aqueous
suspensions at high particle volume fractions and a fixed temperature. Given
their compressibility and deformability, these particles can be packed well
above the random close packing fraction, $\phi_{rcp}$ = 0.64, of monodisperse
hard spheres. The unique characteristics of PNIPAM particles have led to their
use as drug delivery agents [21, 22] and biosensors [23].

####

In the existing scientific literature, researchers have independently
investigated the temperature dependent phase behavior of aqueous PNIPAM
suspensions using rheological and dielectric measurements [24, 25, 26, 27].
Romeo et al. have reported $via$ bulk rheology experiments that PNIPAM
particles in dense aqueous suspensions undergo a glass-liquid-gel transition
with increase in temperature across the LCST due to changes in the Flory-
Huggins parameter [28]. Using UV-visible spectroscopy and thermodynamic
perturbation theory, Wu et al. [27] have investigated the temperature
dependent phase behavior of aqueous PNIPAM suspensions. While PNIPAM particles
were seen to form crystalline structures below the LCST, the suspension
liquefied at a temperature near the LCST and phase separated at higher
temperatures. Su et al. [29] have systematically studied the dielectric
behavior of suspensions of PNIPAM particles of different crosslinker
distributions (dense core (DC), loose core (LC) and homogeneously (HOMO)
cross-linked particles). While the LCST does not change very much, the
dielectric parameters of the aqueous PNIPAM suspensions near the LCST display
significant changes with change in the distribution of crosslinkers.
Dielectric spectroscopy experiments performed by Young et al. [30] have
confirmed that dense PNIPAM suspensions exhibit a colloidal crystal to liquid
transition upon increasing temperature below the LCST. While most of the
previous experimental studies focused exclusively on the rheological or
dielectric properties of aqueous PNIPAM suspensions [29, 30, 31], rheo-
dielectric studies investigating the dielectric response of colloidal hydrogel
suspensions under applied strains have never been reported to the best of our
knowledge.

####

We attempt to bridge this gap by performing rheo-dielectric experiments to
study the dielectric response of PNIPAM colloidal suspensions under large
applied deformations for different particle effective volume fractions and
solvent temperatures. We synthesize PNIPAM particles using a free radical
precipitation polymerization procedure and characterize the thermoresponsive
behavior and LCST of these particles in aqueous suspensions using dynamic
light scattering (DLS) and differential scanning calorimetry (DSC). In the
rheo-dielectric measurements, we apply large amplitude oscillatory strains
(LAOS) to the sample while simultaneously measuring their dielectric responses
using a precision impedance analyzer. We analyze the real and imaginary parts
of the dielectric permittivities of aqueous PNIPAM suspensions of different
effective volume fractions, subjected to large applied strain amplitudes, at
temperatures below, near and above the LCST over a wide frequency range of the
applied electrical voltage. For all the samples investigated in this study,
the real parts of the dielectric permittivities display two distinct
dielectric relaxation processes in the low and high frequency regimes. We
extract the dielectric relaxation strengths and relaxation times at different
strain amplitudes by fitting the experimental data to the Cole-Davidson
relaxation model [32]. In earlier studies, a low frequency dielectric
relaxation process was attributed to segmental motion of the polymers and the
motion of counterions along the polymer chain over the entire temperature
range [30]. The high frequency dielectric relaxation was understood to
originate from counterion fluctuations and interfacial polarization below and
above the LCST respectively [30].

####

For densely-packed aqueous PNIPAM suspensions below the LCST, we report that
both the real and imaginary contributions of the dielectric permittivities
decrease with increase in strain amplitudes. In contrast, the dielectric
responses of densely-packed PNIPAM suspensions at temperatures near and above
the LCST, and of loosely-packed suspensions at all experimental temperatures,
are insensitive to the applied strain. When the applied strain amplitude is
increased systematically, we observe significant increases of the dielectric
relaxation times of densely-packed suspensions at a temperature below the
LCST. This indicates the presence of a counter-intuitive dynamical slowing
down process at very short time scales. Interestingly, we also observe a
simultaneous speeding up of the stress relaxation process in bulk rheological
measurements. Therefore, while our dielectric relaxation data indicates
dynamical slowing down at the shortest accessible length scales, our bulk
rheological data points to simultaneous shear thinning rheology at macroscopic
length scales. We attribute the observed length scale dependent response
functions in densely-packed suspensions below the LCST to the shear-induced
formation and rupture of clusters of swollen PNIPAM particles. We propose that
while the restricted motion of the entangled PNIPAM chains under shear
increases the dielectric relaxation time, the shear-induced disruption of the
swollen PNIPAM microstructures accelerates the stress relaxation in the bulk.
Our study, therefore, provides insights into the dynamics of PNIPAM particles
in aqueous suspensions over several decades of length and time scales. Such
length scale dependent rheology of sheared densely-packed PNIPAM suspensions
can have far-reaching implications in the use of colloidal hydrogels in
electrical energy storage devices such as flow batteries [33] and flow
capacitors [34].

##  2 Material and methods

###  2.1 Synthesis of poly($N$-isopropylacrylamide) (PNIPAM) particles

PNIPAM particles were synthesized by following the one-pot free radical
precipitation polymerization method. All the chemicals were purchased from
Sigma-Aldrich and used as received without further purification. In the
polymerization reaction, 7.0g $N$-isopropylacrylamide (NIPAM) (>99%), 0.875g
$N,N\textquoteright$-methylenebisacrylamide (MBA) (>99.5%) and 0.03g sodium
dodecyl sulphate (SDS) were dissolved in 470ml Milli-Q water (Millipore
Corporation, resistivity 18.2 M$\Omega$ cm) in a three-necked round bottomed
(RB) flask attached with a reflux condenser, a magnetic stirrer with heating
(Heidolph), a platinum sensor and a nitrogen gas (N2) inlet/outlet. The
solution was stirred at 600 RPM and purged with N2 gas for 30min to remove the
oxygen dissolved in water. The free radical precipitation polymerization
reaction was initiated by the addition of 0.28g of potassium persulphate (KPS)
(99.9%) dissolved in 30ml Milli-Q water after heating the mixture to 70∘C. The
reaction was allowed to proceed for 4h with a constant stirring speed of 600
RPM. After 4h, the suspension was cooled down to room temperature. Four
successive centrifugations and re-dispersions at a rotational speed of 20,000
RPM for 60min were used to purify the suspension. After centrifugation, the
supernatant was filtered out and the remaining sample was dried by evaporating
the water. A fine powder was prepared by grinding the dried particles using a
mortar and pestle.

###  2.2 Sample preparation and characterization of PNIPAM particles

Aqueous PNIPAM suspensions of desired concentrations/effective volume
fractions were prepared by stirring PNIPAM powder in Milli-Q water. PNIPAM
particles are soft and can deform when packed above the random close packing
volume fraction of monodisperse hard spheres. Due to their deformability, the
volume fraction $\phi$ is not the appropriate parameter to quantify the
packing of these particles in suspension. Therefore, we use a modified
parameter called the effective volume fraction, $\phi_{eff}$, as a
quantitative measure of the extent of particle packing in suspension. The
protocol for computing $\phi_{eff}$ for suspensions of PNIPAM particles in
water is discussed in detail in section 1 of Supporting Information (SI). The
average sizes and the thermoreversible swelling-deswelling transition of
PNIPAM particles in a dilute aqueous suspension were recorded using a
Brookhaven Instruments Corporation (BIC) BI-200SM dynamic light scattering
(DLS) setup attached with a 150mW solid-state laser (NdYVO4, Coherent Inc.,
Spectra Physics) having an emission wavelength of 532nm. We observe that the
average hydrodynamic diameter of PNIPAM decreases abruptly when temperature is
raised (Fig. 1(a)), thereby signalling a swelling-deswelling transition of the
hydrogel particles. We note that the thermoreversible PNIPAM particles swell
maximally below 20∘C and collapse fully above 45∘C. The lower consolute
solution temperature (LCST) of the aqueous PNIPAM suspensions was quantified
using differential scanning calorimetry (DSC) (Mettler Toledo, DSC 3). The
temperature corresponding to the endothermic peak ($\approx$ 34.3∘C) in Fig.
S2 of SI represents the LCST of the PNIPAM particles in their aqueous
suspensions.

Fig. 1: (a) DLS measurements of the hydrodynamic diameters of thermoreversible
PNIPAM particles during heating and cooling temperature cycles. The different
background colors represent the three different regimes showing the distinct
temperature dependent responses of the thermoreversible PNIPAM; from left to
right: below LCST, near LCST and above LCST, respectively. Measurements were
made at 20∘C, 33∘C and 55∘C, which are marked with vertical lines. (b)
Schematic representation of a rheo-dielectric setup. (c) Dependence of the
real part of the relative dielectric permittivity
($\varepsilon_{r}^{{}^{\prime}}$) of an aqueous PNIPAM suspension of volume
fraction $\phi_{eff}$ = 1.9 on the frequency of the AC electrical voltage for
various applied strain amplitudes at 20∘C (below the LCST). The solid lines
are fits to the real part of Eqn. 1. The inset shows the imaginary part of the
dielectric permittivity ($\varepsilon_{r}^{{}^{\prime\prime}}$) of an aqueous
PNIPAM suspension of volume fraction $\phi_{eff}$ = 1.9 for the same applied
strain amplitudes at 20∘C.

###  2.3 Rheo-dielectric measurements

Dielectric measurements of aqueous PNIPAM suspensions were performed using a
precision impedance analyzer (Wayne Kerr Electronics, 6500B series). An
alternating electrical voltage of amplitude 750mV over a frequency range
between 20Hz and 5MHz was applied, while a stress-controlled Anton Paar MCR
702 rheometer was used to simultaneously apply oscillatory shears to the
sample. A parallel plate sample geometry (PP-25/DI/TI) with a diameter of 25mm
and a measuring gap of 0.5mm was used. The sample was sandwiched between two
plates that simultaneously served as capacitor electrodes and shearing
geometry. Fig. 1(b) displays a schematic diagram of the rheo-dielectric setup.
Oscillatory shear strains of amplitudes between 10% to 2000% at a fixed
angular frequency of 10rad/sec were applied to study the dielectric properties
of aqueous PNIPAM suspensions under flow. The impedance analyzer was connected
to the electrodes (both top and bottom plates) with shielded wires. The
measured capacitance ($C$) and loss factor ($D$) of aqueous PNIPAM suspensions
were used to estimate the real and imaginary parts of the relative dielectric
permittivities ($\varepsilon_{r}^{{}^{\prime}}$ and
$\varepsilon_{r}^{{}^{\prime\prime}}$) using the equations:
$\varepsilon_{r}^{{}^{\prime}}=C(d/A\varepsilon_{0})$ and
$\varepsilon_{r}^{{}^{\prime\prime}}=D\varepsilon_{r}^{{}^{\prime}}$, where
$d$, $A$ and $\varepsilon_{0}$ are the gap between the two electrode plates
(0.5mm in these experiments), the cross-sectional area of the top plate and
the permittivity of free space respectively. Representative plots of the real
and imaginary parts of the dielectric permittivities, measured for a PNIPAM
suspension of effective volume fraction 1.9 below the LCST under different
oscillatory strains, are displayed in Fig. 1(c). In observations that are
consistent with previous studies [29, 30], we note that the real part of the
dielectric permittivity, $\varepsilon_{r}^{{}^{\prime}}$, exhibits two
distinct relaxation processes in the low and high frequency regimes. As
discussed earlier in the context of PNIPAM suspensions [30], the relaxation
observed at frequencies below approximately 1kHz (the low frequency
relaxation) can be attributed to the local segmental motion of the polymer and
also to the motion of counterions along the polymer chain over the entire
temperature range. Fig. S3(a) of SI represents a schematic representation of
the segmental motion of a PNIPAM chain.

####

The relaxation process that we observe at frequencies greater than 1kHz (high
frequency dielectric relaxation) in Fig. 1(c) can be attributed to the
fluctuations of counterions (due to the motion of Na+ counterions around the
SO${}_{4}^{-}$ groups dissociated from sodium dodecyl sulphate (SDS) used
during PNIPAM synthesis) below the LCST and to interfacial polarization
(schematic illustration displayed in Fig. S3(b)) above the LCST [30]. Since
PNIPAM particles swell fully by absorbing water at a temperature below the
LCST, a clear interface between PNIPAM particles and water does not exist.
Therefore, the main contribution to the fast dielectric response below the
LCST arises from the electric field-induced displacement and rearrangement of
counterions bound to the charged colloidal particles [30]. In contrast, PNIPAM
particles collapse and expel most of the absorbed water above the LCST,
thereby revealing a clear particle-water interface. Hence, at a temperature
above the LCST, interfacial polarization (the separation of counterion charges
at the interface) dominates the high frequency relaxation process. The
acquired dielectric data has tiny fluctuations presumably arising due to the
presence of metallic electrical contacts [7]. In the present work, the raw
dielectric permittivity data have been smoothed by averaging over 10 adjacent
data points (Fig. S4 of SI).

The acquired frequency dependent dielectric permittivity data was fitted to
the Cole-Davidson relaxation model by incorporating two relaxation terms that
correspond to the low and high frequency regimes (the slow and fast relaxation
processes respectively), while also accounting for electrode polarization [29,
30]:

$\varepsilon_{r}^{*}=\varepsilon_{h}+\frac{\varepsilon_{l}-\varepsilon_{m}}{1+(j\omega\uptau_{s})^{\beta_{s}}}+\frac{\varepsilon_{m}-\varepsilon_{h}}{1+(j\omega\uptau_{f})^{\beta_{f}}}+A\omega^{-m}$
(1)

Here, $\varepsilon_{l},\varepsilon_{m}$ and $\varepsilon_{h}$ are the
dielectric permittivities in the low, intermediate and high frequency regimes
respectively, $\uptau_{s}$ and $\uptau_{f}$ are respectively the slow and fast
dielectric relaxation times, while $\beta_{s}$ and $\beta_{f}$ are the
corresponding exponents which vary between 0 to 1. It is to be noted that
$\Delta\varepsilon_{l}$ (=$\varepsilon_{l}-\varepsilon_{m}$) and
$\Delta\varepsilon_{h}$ (=$\varepsilon_{m}-\varepsilon_{h}$) are the
dielectric relaxation strengths in the low and high frequency regimes
respectively. The term $A\omega^{-m}$, where $A$ and $m$ are adjustable
parameters, is introduced to account for electrode polarization (EP) which
occurs at lower frequencies due to the accumulation of free charges at the
electrode interface. These charged species presumably originate from the
electrolyte (potassium persulfate KPS) that we employed as a reaction
initiator during the synthesis process.

Fig. 2: Three dimensional representations at 20∘C (below the LCST) of the (a)
real ($\varepsilon_{r}^{{}^{\prime}}$) and (b) imaginary
($\varepsilon_{r}^{{}^{\prime\prime}}$) parts of the relative dielectric
permittivities of aqueous PNIPAM suspensions of effective volume fraction
$\phi_{eff}$ = 1.9, (c) $\varepsilon_{r}^{{}^{\prime}}$ and (d)
$\varepsilon_{r}^{{}^{\prime\prime}}$ for effective volume fraction
$\phi_{eff}$ = 1.6, and (e) $\varepsilon_{r}^{{}^{\prime}}$ and (f)
$\varepsilon_{r}^{{}^{\prime\prime}}$ for effective volume fraction
$\phi_{eff}$ = 1.3, over a range of AC electrical voltage frequencies and
applied oscillatory strains. The color bars represent the values of
$\varepsilon_{r}^{{}^{\prime}}$ and $\varepsilon_{r}^{{}^{\prime\prime}}$ with
their magnitudes increasing from bottom to top.

##  3 Results and Discussion

Fig. 2 shows three-dimensional representations of the real and imaginary parts
of the dielectric permittivities of aqueous poly($N$-isopropylacrylamide)
(PNIPAM) suspensions at 20∘C (below the LCST) versus frequency of the applied
alternating electrical voltage and the applied strain amplitude. At very low
strains and at a temperature below the LCST, the elasticities of PNIPAM
suspensions of the higher effective volume fractions ($\phi_{eff}$ = 1.9 and
1.6) are comparable to each other, but are very different from the sample of
the lowest effective volume fraction ($\phi_{eff}$ = 1.3) used here (Figs.
S5(a) of SI). We designate the samples of $\phi_{eff}$ = 1.9 and 1.6 as
densely-packed suspensions and the one of $\phi_{eff}$ = 1.3 as a loosely-
packed suspension. Figs. S6 and S7 of SI show three-dimensional
representations of the dielectric permittivities as functions of strain
amplitudes and AC electrical voltage frequencies for the densely-packed and
loosely-packed aqueous PNIPAM suspensions at temperatures 33∘C (near the LCST)
and 55∘C (above the LCST) respectively. From Figs. 2(a-d), we note that both
real and imaginary dielectric permittivities decrease with increasing applied
strain amplitudes for the densely-packed PNIPAM suspensions ($\phi_{eff}$ =
1.9 and 1.6) below the LCST. The observed decreases in
$\varepsilon_{r}^{{}^{\prime}}$ and $\varepsilon_{r}^{{}^{\prime\prime}}$ due
to an increase in applied strains can be attributed to enhanced structural
deformations of the swollen PNIPAM particles in aqueous suspensions at high
effective volume fractions. The color bars in Figs. 2, S6 and S7 of SI
represent the values of $\varepsilon_{r}^{{}^{\prime}}$ and
$\varepsilon_{r}^{{}^{\prime\prime}}$ with their magnitudes increasing from
bottom to top. We observe that for the loosely-packed PNIPAM suspension, the
dielectric permittivities are insensitive to the applied strain amplitude at
all temperatures (Figs. (e-f) of 2, S6 and S7). For the densely-packed PNIPAM
suspensions, in contrast, the real and imaginary permittivities are
independent of the applied shear only at temperatures near and above the LCST
(Figs. (a-d) of S6 and S7 of SI). For significant modifications in the
dielectric properties of homogeneous polymeric systems, the applied
deformation must be greater than the equilibrium molecular motion (shear rate
> relaxation time-1) [5, 7]. Small deformations can affect the dielectric
properties of glassy systems due to the presence of dynamically heterogeneous
structures [5, 7]. Our results indicate that densely-packed PNIPAM suspensions
at a temperature below the LCST are possibly in a disordered and metastable
(glassy) state constituted by PNIPAM particle assemblies. For the loosely-
packed aqueous PNIPAM suspension at all temperatures and the densely-packed
suspensions at temperatures near or above the LCST, the observed insensitivity
of the dielectric responses to applied deformations presumably arises because
of an absence of self-assembled particle aggregates. Smaller particle
assemblies have much higher relaxation rates and therefore require much larger
shears to induce appreciable changes in the dielectric properties, hence the
observed insensitivity of the measured dielectric responses to applied shears
in the last scenario.

Fig. 3: (a) Low frequency ($\Delta\varepsilon_{l}$) and (b) high frequency
($\Delta\varepsilon_{h}$) dielectric relaxation strengths and (c) slow
($\uptau_{s}$) and (d) fast ($\uptau_{f}$) relaxation times as functions of
applied oscillatory strain amplitudes for aqueous suspensions of PNIPAM of
volume fraction $\phi_{eff}$ = 1.9 at temperatures below, near and above the
LCST. These values are obtained by fitting the data in Figs. 2(a), S6(a) and
S7(a) to the real part of Eqn. 1.

####

The real parts of the dielectric permittivities of aqueous PNIPAM suspensions
show two distinct dielectric relaxation processes in the low and high
frequency regimes, as can be clearly seen in (a), (c) and (e) of Figs. 2, S6
and S7. We compute the relaxation times and dielectric strengths (determined
by the number density of molecules, magnitudes of the moments and orientation
fluctuations of the dipoles [4, 5]) of all the PNIPAM suspensions prepared at
temperatures below, near or above the LCST, by fitting the acquired
$\varepsilon_{r}^{{}^{\prime}}$ data to the real part of Eqn. 1 (fitted data
are shown in Fig. S8 of SI). The dielectric strengths and relaxation times
provide important information about the orientation and dynamics of the PNIPAM
particles under shear flow. It has been reported that the alignment of dipoles
parallel and perpendicular to the electric field respectively increase and
decrease the measured dielectric strengths [4]. The values of the fitted
parameters in the electrode polarization term ($A$ and $m$) of Eqn. 1, which
originates from the mobility of the free ions and contributes to the low
frequency dielectric permittivities, are listed in Tables 2-4 of SI. In Figs.
3 (a-b) and S9 (a-b) of SI, we show the low- and high frequency dielectric
relaxation strengths ($\Delta\varepsilon_{l}$ and $\Delta\varepsilon_{h}$
respectively) of densely-packed PNIPAM suspensions ($\phi_{eff}$ = 1.9 in
Figs. 3(a-b) and $\phi_{eff}$ = 1.6 in Figs. S9(a-b) of SI) at all the
experimental temperatures. Figs. S10 (a-b) display the values of
$\Delta\varepsilon_{l}$ and $\Delta\varepsilon_{h}$ when oscillatory strains
were applied to the loosely-packed PNIPAM suspension ($\phi_{eff}$ = 1.3) at
the same temperatures. For all the PNIPAM suspensions,
$\Delta\varepsilon_{l}$, which is understood to originate from the segmental
motion of PNIPAM chains over the entire temperature range as reported in an
earlier work [30], shows an enhancement with increasing suspension
temperatures. This can be attributed to the increase in polymer density inside
the PNIPAM particles due to particle shrinkage/ collapse above the LCST.
Interestingly, we note that both $\Delta\varepsilon_{l}$ and
$\Delta\varepsilon_{h}$ decrease strongly with increase in applied strain
amplitudes for densely-packed PNIPAM suspensions at a temperature below the
LCST. It has been proposed that the application of shear strain disrupts the
structured polymer particles [4] that are present in densely-packed (jammed)
suspensions and randomizes the directions of the dipoles. Strong shears are
expected to result in the alignment of dipoles parallel to the flow direction.
In our experimental geometry, the flow direction is perpendicular to the
electric field and results in the observed decrease in net polarization and
dielectric strengths. In contrast, no significant strain-induced changes in
$\Delta\varepsilon_{l}$ and $\Delta\varepsilon_{h}$ are observed for loosely-
packed PNIPAM suspension at all experimental temperatures and for densely-
packed PNIPAM suspensions at temperatures near and above the LCST. The strain
insensitivity of the dielectric strengths under these conditions arises as the
PNIPAM particles do not self-assemble to form strong structures due to low
effective volume fractions. We therefore conclude that the packing densities
and temperature dependent morphological changes of PNIPAM particles
effectively control the dielectric properties and shear responses of their
aqueous suspensions.

####

Apart from the dielectric strength, the other parameters that we study to
quantify the responses of the charged species are the fast and slow relaxation
times [29, 30], extracted by fitting the observed two-step decay of the real
part of the dielectric permittivities to Eqn. 1. Figs. 3(c-d) and S9(c-d)
display the slow ($\uptau_{s}$) and fast ($\uptau_{f}$) relaxation times of
densely-packed aqueous PNIPAM suspensions at $\phi_{eff}$ = 1.9 and 1.6
respectively as a function of strain amplitudes at all the temperatures
explored in this study. Figs. S10(c-d) represent the slow and fast relaxation
times as a function of strain amplitude for the loosely-packed PNIPAM
suspension (i.e. $\phi_{eff}$ = 1.3) at the same temperatures. We see that
both $\uptau_{s}$ and $\uptau_{f}$ decrease with increase in suspension
temperature for all the samples at a fixed effective volume fraction. We
attribute the observed decrease in time scales with increasing temperatures to
the reduction in effective volume fractions and enhanced mobilities of the
collapsed PNIPAM particles [35]. Interestingly, we see from Figs. 3(c-d) and
Figs. S9(c-d) of SI that $\uptau_{s}$ and $\uptau_{f}$ increase monotonically
with increasing strain amplitudes for the densely-packed PNIPAM suspensions at
a temperature below the LCST, thereby demonstrating a slowdown of the
polarization fluctuation dynamics with increasing strain amplitude. However,
the dielectric relaxation time scales of these samples are insensitive to the
applied shear at temperatures near and above the LCST where the particles have
already collapsed due to expulsion of water. Furthermore, we note that
$\uptau_{s}$ and $\uptau_{f}$ are relatively unaffected by the applied strain
amplitudes at all temperatures when PNIPAM particles are loosely-packed in
aqueous suspension (Figs. S10(c-d) of SI). The time scales extracted from our
rheo-dielectric measurements therefore reveal the distinct dynamical
properties of PNIPAM particles in aqueous suspensions at different
temperatures and under different externally imposed mechanical stresses.

####

Moreover, from Figs. 2, S6 and S7, we observe that the imaginary contributions
to the dielectric permittivities show a single loss peak. Figs. S11(a-c) of SI
show the relaxation times estimated from the peak frequency values of the
dielectric losses ($\uptau=1/2\pi f_{peak}$), where $f_{peak}$ is the
frequency at which the dielectric loss peak is noted. We see that the
relaxation times $\uptau$ are of the order of the fast relaxation processes
estimated from fits of the real part of the dielectric permittivities to the
Cole-Davidson relaxation model (Figs. 3(d) and S9(d) of SI). We therefore
conclude that the dominant relaxation mechanism responsible for the observed
dielectric loss originates from counterion polarization below the LCST. In
contrast, at temperatures above the LCST, the dielectric loss can be
attributed to interfacial polarization generated by the free charged species
present in the suspension and the dynamics of counterions around the
particles. In results that are consistent with data displayed in Figs. 3(d)
and S9(d), we note from Figs. S11(a-b) that the dielectric relaxation times,
$\uptau$, of densely-packed PNIPAM suspensions ($\phi_{eff}$ = 1.9 and 1.6)
increase with increase in applied strain amplitudes at a temperature below the
LCST. However, the relaxation rates are independent of temperature for
densely-packed PNIPAM suspensions at temperatures near and above the LCST and
for the loosely-packed suspension ($\phi_{eff}$ = 1.3) at all temperatures.

Fig. 4: Oscillatory shear stress response vs. time exhibited by PNIPAM
suspensions of effective volume fractions $\phi_{eff}$ = (a) 1.9, (b) 1.6 and
(c) 1.3 at temperature 20∘C (below the LCST) and for different applied
oscillatory strains $\gamma$. Fig. 5: Stress relaxation times, $\uptau_{1}$,
of densely-packed aqueous suspensions of PNIPAM particles at temperatures 20∘C
and 55∘C as a function of applied oscillatory strain amplitude.

####

To complement our rheo-dielectric study that uncovers counter-intuitive
polarization fluctuation dynamics upon applying oscillatory shear
deformations, we next study the influence of applied large amplitude
oscillatory strains on the bulk mechanical moduli of PNIPAM suspensions by
performing independent stress relaxation experiments. In these measurements,
the range of the large amplitude oscillatory strains is comparable to that
applied in our rheo-dielectric experiments. Our experimental design therefore
allows us to study the responses of PNIPAM suspensions over a length scale
range varying over several (7-8) orders of magnitude. Stress relaxation
decays, exhibited by densely-packed (i.e. $\phi_{eff}$ = 1.9 and 1.6) and
loosely-packed (i.e. $\phi_{eff}$ = 1.3) aqueous PNIPAM suspensions at
temperatures below, near and above the LCST, upon the application of
oscillatory strains of pre-determined amplitudes, are shown in Figs. 4, S12,
and S13, respectively. For the densely-packed PNIPAM suspensions at/near the
LCST, power-law stress relaxations are observed (Fig. S12(a-b)). For the
loosely-packed suspension under the same condition, the stress decay flattens
out very rapidly (Fig. S12(c)). Interestingly, the two-step stress relaxation
processes observed for the densely-packed PNIPAM suspensions (Figs. 4(a-b) and
S13(a-b)) below and above the LCST are reminiscent of the relaxation of glassy
materials or entangled polymeric melts [36, 37]. Given the extremely rapid
decay of the short-time stress relaxation, we could not acquire adequate
reliable data to model the fast stress relaxation process satisfactorily. We
therefore focus only on the slow bulk stress relaxation data which is well
described by the Kohlrausch–Williams–Watts (KWW) model [38]. The KWW model
represents the time dependent stress relaxation of the slower mode in glassy
materials as a stretched exponential function:

$\sigma(\uptau)=Aexp[-(\uptau/\uptau_{1})^{\beta}]$ (2)

where, $\sigma(\uptau)$, $A$, $\uptau{{}_{1}}$ and $\beta$ are the time
dependent shear stress, intercept of the slow relaxation decay, time scale
corresponding to the slow mode of the stress relaxation process and stretching
exponent (0 < $\beta$ < 1) respectively. We note that the stress relaxation
times for densely-packed PNIPAM suspensions, which are obtained by fitting the
data plotted in Figs. 4(a-b) to Eqn. 2, decrease with applied shear strain
amplitude below the LCST (Fig. 5), thereby revealing the presence of bulk
shear thinning behavior at low temperatures. Interestingly, therefore, while
the polarization fluctuation in a suspension of densely-packed PNIPAM
particles slows down with an increase in the applied strain below the LCST,
the bulk stress relaxation speeds up simultaneously.

Fig. 6: A schematic illustration of a densely-packed PNIPAM suspension under
an applied strain. The green arrow at the top and the red arrow on the left
represent, respectively, the increase in applied strain amplitude and the
direction of the electric field.

The observed dichotomous dynamics at macroscopic and nanometre length scales
can be understood by considering the shear-induced formation and rupture of
fragile PNIPAM particle clusters in densely-packed suspensions. It has been
reported in earlier work that PNIPAM particles in dense aqueous suspensions
phase separate by forming clusters under shear flow [35, 39]. We attribute the
slowing down of the slow and fast dielectric relaxation processes to the
formation of such microscopic PNIPAM particle clusters. This is illustrated in
Fig. 6 which schematically represents polarization fluctuations in the sample
(left-panel) and breakup of PNIPAM clusters (right-panel) due to stresses
imposed by electric fields and shear strains respectively. Badiger et al. [35]
have observed stress overshoots in aqueous suspensions of PNIPAM particles
undergoing shear flow due to the presence of entangled PNIPAM polymers. Such
shear-induced entanglements of PNIPAM chains in densely-packed suspensions
below the LCST are expected to restrict the segmental motion of PNIPAM chains,
resulting in the observed increase of the slow dielectric relaxation times.
Furthermore, counterion separation is enhanced due to the fragmentation of
particle clusters undergoing shear deformation (Fig. 6, right-hand panel).
This enhances the particle dipole moments and causes the dielectric fast
relaxation processes to slow down. The simultaneous decrease in the dielectric
strengths with increasing shears arises from the randomization of the dipoles
and the distortions and breakups of PNIPAM clusters by the applied shear.
Furthermore, the decrease in the bulk viscosity of PNIPAM suspensions due to
shear induced breakdown of the fragile clusters results in accelerated
particle dynamics and shear thinning at macroscopic length scales.

##  4 Conclusions

The bulk mechanical response and the polarization fluctuations of sheared
complex fluids have independently been studied in the literature using
rheometry and dielectric spectroscopy [24, 25, 26, 27, 28, 29, 30, 31].
However, to the best of our knowledge, simultaneous rheological and dielectric
experiments have never been performed to study the dynamics of sheared
colloidal suspensions over spatiotemporal scales spanning several decades. Our
work attempts to bridge this gap by simultaneously studying the dielectric and
stress relaxation behaviors of aqueous suspensions of colloidal hydrogel
poly($N$-isopropylacrylamide) (PNIPAM) particles under large applied
oscillatory shear strains. We observe two distinct dielectric relaxation
processes in the low and high frequency regimes. The low frequency relaxation
process is understood to originate from the segmental motion of the polymer
and the motion of counterions along the polymers chains [30] over the entire
temperature range. On the other hand, the high frequency relaxation is
believed to arise from counterion fluctuations at a temperature below the LCST
and from interfacial polarization above the LCST. We propose that the observed
decrease of the dielectric relaxation strengths in densely-packed PNIPAM
suspensions upon the application of oscillatory shear strains at a temperature
below the LCST is due to the distortion of dynamically heterogeneous
suspension structures and randomization in the orientations of the constituent
dipoles. Notably, the dielectric responses of loosely-packed PNIPAM
suspensions at all temperatures and densely-packed PNIPAM suspensions at
temperatures near and above the LCST are insensitive to the applied
deformations. The most significant finding of this study is that while the
polarization fluctuations of densely-packed PNIPAM suspensions slow down with
increase in applied oscillatory strain amplitudes at a temperature below the
LCST, the bulk stress relaxation dynamics speed up under identical conditions.
Earlier studies such as those by Su et al. [29] and Yang et al. [30] reported
the dielectric behavior of relatively dilute aqueous PNIPAM suspensions (0.005
wt.% and 5.6 wt.% respectively). The dielectric strengths and time scales
reported in the present study are therefore much higher than for the
relaxation dynamics reported earlier [29, 30]. The dense PNIPAM suspensions in
our experiments are characterised by a significantly enhanced number of
dipoles that restrict the segmental motion of PNIPAM chains, resulting in the
novel dynamics reported in our study.

####

Besides revealing the rich dynamical behavior of dense suspensions of PNIPAM
hydrogels at different length scales, we show here that the self-assembly and
rheological properties of these samples can be fine-tuned by controlling
temperature and the applied strain deformation. These hydrogels are therefore
attractive candidates for engineering energy storage devices, such as
stretchable, soft and flow batteries [33] and flow capacitors [34]. We
strongly believe that the present study can serve as a starting point for
further research into the length scale dependent structural and dynamical
properties of dense colloidal (soft glassy) suspensions and other kinetically
constrained complex systems.

## Acknowledgments

The authors thank K. N. Vasudha for help with DSC measurements and N. V.
Madhusudana for useful discussions. We acknowledge Department of Science and
Technology, Science and Engineering Research Board for funding our research
(grant number EMR/2016/006757).

## References

  * [1] H. Watanabe, Y. Matsumiya, K. Horio, Y. Masubuchi, and T. Uneyama, Rheo-Dielectric Behavior of Soft Matters, pp. 37–87.
  * [2] T. Pakula, D. Vlassopoulos, G. Fytas, and J. Roovers, “Structure and dynamics of melts of multiarm polymer stars,” Macromolecules, vol. 31, no. 25, pp. 8931–8940, 1998.
  * [3] T. Nicolai and G. Floudas, “Dynamics of linear and star poly(oxypropylene) studied by dielectric spectroscopy and rheology,” Macromolecules, vol. 31, no. 8, pp. 2578–2585, 1998.
  * [4] S. Capaccioli, D. Prevosto, A. Best, A. Hanewald, and T. Pakula, “Applications of the rheo-dielectric technique,” Journal of Non-Crystalline Solids, vol. 353, no. 47, pp. 4267–4272, 2007.  Dielectric Relaxation and Related Phenomena.
  * [5] H. Watanabe, T. Sato, Y. Matsumiya, T. Inoue, and K. Osaki, “Rheo-dielectrics : Its applicability,” Nihon Reoroji Gakkaishi, vol. 27, no. 2, pp. 121–125, 1999.
  * [6] Y. Peng, Y. M. Shkel, and G. Kim, “Stress dielectric response in liquid polymers,” Journal of Rheology, vol. 49, no. 1, pp. 297–311, 2005.
  * [7] H. Watanabe, Y. Matsumiya, and T. Inoue, “Rheo-dielectrics in oligomeric and polymeric fluids: a review of recent findings,” J. Phys.: Condens. Matter, vol. 15, pp. S909–S921, mar 2003.
  * [8] J. Knapik-Kowalczuk, K. Jurkiewicz, A. Kocot, and M. Paluch, “Rheo-dielectric studies of the kinetics of shear-induced nematic alignment changes in itraconazole,” Journal of Molecular Liquids, vol. 302, p. 112494, 2020\.
  * [9] K. Negita and H. Kaneko, “Rheodielectric study on shear-induced structural change in the smectic-$a$ phase of $4\text{$-$}n\text{-octyl-}{4}^{{}^{\prime}}\text{-cyanobiphenyl}$ (8cb),” Phys. Rev. E, vol. 80, p. 011705, Jul 2009.
  * [10] K. Horio, T. Uneyama, Y. Matsumiya, Y. Masubuchi, and H. Watanabe, “Rheo-dielectric responses of entangled cis-polyisoprene under uniform steady shear and laos,” Macromolecules, vol. 47, no. 1, pp. 246–255, 2014\.
  * [11] P. Placke, R. Richert, and E. W. Fischer, “Rheodielectric spectroscopy on heterogeneous polymer systems —electrorheological fluids,” Colloid and Polymer Science, vol. 273, pp. 848–856, 1995.
  * [12] B. Nath, M. P. Borthakur, and G. Biswas, “Electric field induced dynamics of viscoplastic droplets in shear flow,” Physics of Fluids, vol. 32, no. 9, p. 092110, 2020.
  * [13] D. Steinhauser, M. Möwes, and M. Klüppel, “Carbon black networking in elastomers monitored by simultaneous rheological and dielectric investigations,” Journal of Physics: Condensed Matter, vol. 28, p. 495103, oct 2016.
  * [14] M. Heskins and J. E. Guillet, “Solution properties of poly(n-isopropylacrylamide),” Journal of Macromolecular Science: Part A \- Chemistry, vol. 2, no. 8, pp. 1441–1455, 1968.
  * [15] Y. Hirokawa and T. Tanaka, “Volume phase transition in a non-ionic gel,” AIP Conference Proceedings, vol. 107, no. 1, pp. 203–208, 1984.
  * [16] H. Schild, “Poly(n-isopropylacrylamide): experiment, theory and application,” Progress in Polymer Science, vol. 17, no. 2, pp. 163–249, 1992.
  * [17] J. Appel, B. Fölker, and J. Sprakel, “Mechanics at the glass-to-gel transition of thermoresponsive microgel suspensions,” Soft Matter, vol. 12, pp. 2515–2522, 2016.
  * [18] P. J. Yunker, K. Chen, M. D. Gratale, M. A. Lohr, T. Still, and A. G. Yodh, “Physics in ordered and disordered colloidal matter composed of poly(n-isopropylacrylamide) microgel particles,” vol. 77, p. 056601, may 2014\.
  * [19] S. K. Behera, D. Saha, P. Gadige, and R. Bandyopadhyay, “Effects of polydispersity on the glass transition dynamics of aqueous suspensions of soft spherical colloidal particles,” Phys. Rev. Materials, vol. 1, p. 055603, Oct 2017.
  * [20] K. Jain, R. Vedarajan, M. Watanabe, M. Ishikiriyama, and N. Matsumi, “Tunable lcst behavior of poly(n-isopropylacrylamide/ionic liquid) copolymers,” Polym. Chem., vol. 6, pp. 6819–6825, 2015.
  * [21] L. chang Dong and A. S. Hoffman, “A novel approach for preparation of ph-sensitive hydrogels for enteric drug delivery,” Journal of Controlled Release, vol. 15, no. 2, pp. 141–152, 1991.
  * [22] C. M. Nolan, C. D. Reyes, J. D. Debord, A. J. García, and L. A. Lyon, “Phase transition behavior, protein adsorption, and cell adhesion resistance of poly(ethylene glycol) cross-linked microgel particles,” Biomacromolecules, vol. 6, no. 4, pp. 2032–2039, 2005.  PMID: 16004442.
  * [23] M. R. Islam, A. Ahiabu, X. Li, and M. J. Serpe, “Poly (n-isopropylacrylamide) microgel-based optical devices for sensing and biosensing,” Sensors, vol. 14, no. 5, pp. 8984–8995, 2014.
  * [24] H. Senff and W. Richtering, “Temperature sensitive microgel suspensions: Colloidal phase behavior and rheology of soft spheres,” The Journal of Chemical Physics, vol. 111, no. 4, pp. 1705–1711, 1999.
  * [25] C. Misra, S. K. Behera, and R. Bandyopadhyay, “Influence of particle size on the thermoresponsive and rheological properties of aqueous poly(n-isopropylacrylamide) colloidal suspensions,” Bulletin of Materials Science, vol. 43, p. 182, Jul 2020.
  * [26] J. Wu, G. Huang, and Z. Hu, “Interparticle potential and the phase behavior of temperature-sensitive microgel dispersions,” Macromolecules, vol. 36, no. 2, pp. 440–448, 2003.
  * [27] J. Wu, B. Zhou, and Z. Hu, “Phase behavior of thermally responsive microgel colloids,” Phys. Rev. Lett., vol. 90, p. 048304, Jan 2003.
  * [28] G. Romeo, A. Fernandez-Nieves, H. M. Wyss, D. Acierno, and D. A. Weitz, “Temperature-controlled transitions between glass, liquid, and gel states in dense p-nipa suspensions,” Advanced Materials, vol. 22, no. 31, pp. 3441–3445, 2010.
  * [29] W. Su, K. Zhao, J. Wei, and T. Ngai, “Dielectric relaxations of poly(n-isopropylacrylamide) microgels near the volume phase transition temperature: impact of cross-linking density distribution on the volume phase transition,” Soft Matter, vol. 10, pp. 8711–8723, 2014.
  * [30] M. Yang, C. Liu, Y. Lian, K. Zhao, D. Zhu, and J. Zhou, “Relaxations and phase transitions during the collapse of a dense pnipam microgel suspension—thorough insight using dielectric spectroscopy,” Soft Matter, vol. 13, pp. 2663–2676, 2017.
  * [31] J. Zhou, J. Wei, T. Ngai, L. Wang, D. Zhu, and J. Shen, “Correlation between dielectric/electric properties and cross-linking/charge density distributions of thermally sensitive spherical pnipam microgels,” Macromolecules, vol. 45, no. 15, pp. 6158–6167, 2012.
  * [32] K. Asami, “Characterization of heterogeneous systems by dielectric spectroscopy,” Progress in Polymer Science, vol. 27, no. 8, pp. 1617–1659, 2002.
  * [33] A. Helal, T. Divoux, and G. H. McKinley, “Simultaneous rheoelectric measurements of strongly conductive complex fluids,” Phys. Rev. Applied, vol. 6, p. 064004, Dec 2016.
  * [34] V. Presser, C. R. Dennison, J. Campos, K. W. Knehr, E. C. Kumbur, and Y. Gogotsi, “The electrochemical flow capacitor: A new concept for rapid energy storage and recovery,” Advanced Energy Materials, vol. 2, no. 7, pp. 895–902, 2012.
  * [35] M. V. Badiger and B. A. Wolf, “Shear induced demixing and rheological behavior of aqueous solutions of poly(n-isopropylacrylamide),” Macromolecular Chemistry and Physics, vol. 204, no. 4, pp. 600–606, 2003.
  * [36] R. Bandyopadhyay, P. H. Mohan, and Y. M. Joshi, “Stress relaxation in aging soft colloidal glasses,” Soft Matter, vol. 6, pp. 1462–1466, 2010.
  * [37] D. Feldman, “The theory of polymer dynamics, by m. doi and s. f. edwards, the clarendon press, oxford university press, new york, 1986, 391 pp. price: $78.50.,” Journal of Polymer Science Part C: Polymer Letters, vol. 27, no. 7, pp. 239–240, 1989.
  * [38] G. Williams and D. C. Watts, “Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function,” Trans. Faraday Soc., vol. 66, pp. 80–85, 1970.
  * [39] M. Stieger and W. Richtering, “Shear-Induced phase separation in aqueous polymer solutions: Temperature-Sensitive microgels and linear polymer chains,” Macromolecules, vol. 36, pp. 8811–8818, Nov. 2003.