\section{Conclusions

We investigated, via large scale simulations, the role that surface tension and mobility models play in continuum descriptions of grain growth.
This was accomplished by applying in the context of very large numbers of grains the new, extremely simple and streamlined threshold dynamics algorithm recently introduced in \cite{SalvadorEsedogluSimplified} together with the algorithm previously introduced in \cite{SelimFelix} to perform 2D and 3D simulations.
Grain statistics, such as the GSD and the MDF, for three distinct choices of surface tension / mobility models were compared to each other and to available experimental data in existing literature.
The 2D simulations did not reveal a strong dependence of the GSD on the surface tension / mobility model used.
However, in 3D, a clear distinction emerged between the stationary GSDs of the model with anisotropic reduced mobilities and the other two models considered that had isotropic reduced mobilities.
The asymptotic behavior of the MDF was also studied, both in 2D and 3D, using 2D and 3D crystallography, respectively.
In 3D, it was observed that as long as the orientation texture of the initial data is random, even an initial MDF that is very far from the Mackenzie distribution eventually enters a very close vicinity of it.
\section{Introduction

The shapes and sizes of grains constituting polycrystalline materials, such as most metals and ceramics, influence important physical attributes of these materials, for example their conductivity and yield strength.
Therefore, there has been much interest in modeling and simulation of grain boundary motion during common manufacturing processes, such as annealing.
In particular, the evolution of grain size distribution (GSD), the grain boundary character distribution (GBCD), as well as a simpler version of the latter, the misorientation distribution function (MDF), have been a central focus of research in materials science.

In this paper, we report simulation results for grain boundary motion at the mesoscale, using a recent, very concise and flexible numerical algorithm developed for the continuum model of grain boundaries that was introduced by Mullins \cite{MullinsGrainBoundaryMotion} and earlier authors.
According to this well known model, interfaces in the grain network evolve via steepest descent for an energy that penalizes surface areas of the interfaces, often weighted by a factor, known as the surface tension, that depends on the misorientation across each interface.
The resulting dynamics is described by a system of partial differential equations (PDEs), coupled to each other along junctions by what are known as Herring angle conditions \cite{HerringAngleCondition}, and allows specifying in addition a mobility factor for each interface.
The results we present are the first large scale simulations using the new algorithm proposed in \cite{SalvadorEsedogluSimplified} for this PDE system.

We investigate several questions that have repeatedly come up in previous literature concerning certain statistics of the grain network.
In particular, we study, in both two and three dimensions:
\begin{itemize}
\item The effect of surface tension and mobility model used on statistics such as the GSD.
\item The asymptotic behavior of the MDF and its dependence on initial configuration.
\end{itemize}
Specifically, we simulate in both 2D and 3D Mullins' mesoscale model with
\begin{enumerate}
\item Read-Shockley surface tensions \cite{ReadShockley} (and its extension \cite{HolmMisorientationAngles} to 3D anisotropy) along with isotropic mobilities,
\item Read-Shockley surface tensions with ``reciprocal'' mobilities resulting in isotropic {\it reduced} mobilities but asymmetric junction angles, and
\item Isotropic surface tensions and mobilities, resulting in isotropic reduced mobilities as well as symmetric (all $120^\circ$) junction angles.
\end{enumerate}
We compare to experimental data whenever it is available in existing literature.
Some of our results and conclusions differ substantially from earlier studies.
Highlights of our findings are:
\begin{itemize}
\item In 2D simulations, the stationary GSD that is reached appears to not depend on the surface tension and mobility models used: it is essentially the same in all three surface tension / mobility combinations (i), (ii), (iii) used in our simulations.
This is in agreement with a variety of previous studies, e.g. \cite{upmanyu,Hallberg2014}.
\item However, in 3D, the stationary GSD reached is substantially different when Read-Shockley surface tensions along with isotropic mobilities (i) are used, compared to (ii) and (iii) that result in isotropic reduced mobilities.
This finding appears to be distinct from earlier studies, e.g. \cite{GruberMCPart1}.
\item Moreover, the different GSD reached in 3D by models (i) and (ii) that differ only in their choice of mobilities constitutes a contrary point to previous studies that suggest mobilities do not play a significant role in the long time behavior of the GSD.
\item The GSD reached by (ii) and (iii) are almost identical in both 2D and 3D, even though junction angles in these simulations are very different, suggesting that perhaps the isotropy vs. anisotropy of reduced mobilities is what is significant.
\item The GSD reached in our 3D simulations with surface tension / mobility combination (i) is not less consistent with experimental measurements \cite{rowenhurst,zhang,groeber}, which have substantial variation themselves.
\item Statistics of grain {\em shapes}, such as that of isoperimetric ratios (eccentricities), also display considerable dependence on the choice of surface tension / mobility model.
\item Models (i) \& (ii) allow us to study the asymptotic behavior of the MDF.
Confirming earlier studies, the two drastically different choices of mobility in (i) \& (ii) lead to almost identical behavior, both in 2D and 3D.
Distinct from  earlier studies, we exhibit 3D initial data with MDFs perturbed substantially away from the Mackenzie distribution and observe the MDF return to a very close vicinity of the Mackenzie distribution, using either (i) or (ii) as the surface tension / mobility model.
\end{itemize}

The code for the simulations is publicly available and can be found at \url{https://github.com/tiagosalvador/statisticsGBM}.

\section{Acknowledgements
\noindent The authors gratefully acknowledge support from the NSF grant DMS-1719727.

\bibliographystyle{amsalpha}

\section{Our model and algorithm}\label{sec:model

Our starting point is the curvature driven grain boundary motion model that goes back at least to Mullins \cite{MullinsGrainBoundaryMotion}.
According to this model, grain boundaries evolve via steepest descent for the energy
\begin{equation}\label{eq:energy}
E = \sum_{i<j} \sigma_{ij}\text{Area}(\Gamma_{ij})
\end{equation}
where $\Gamma_{i,j}$ denotes the boundary between two adjacent grains $\Sigma_i$ and $\Sigma_j$.
In particular we will neglect the normal dependence of the energy density, but allow each grain boundary $\Gamma_{i,j}$ to have a distinct surface tension $\sigma_{i,j}$.
The very special case $\sigma_{i,j} = 1$ for all $i\not= j$ is often referred to as the {\em isotropic} surface tension model.
We will also be interested in the surface tension model due to Read \& Shockley \cite{ReadShockley} and its extension \cite{HolmMisorientationAngles} to 3D crystallography:
\begin{equation}\label{eq:RSmodel}
\sigma_{ij} = \begin{cases}
\frac{\theta_{ij}}{\theta_*}\left(1-\log\left(\frac{\theta_{ij}}{\theta_*}\right)\right)	& \text{if } \theta_{ij} \in [0,\theta_*],\\
1	& \text{if } \theta_{ij} \geq \theta_*.
\end{cases}
\end{equation}
where $\theta_{i,j}$ is the misorientation angle between $\Sigma_i$ and $\Sigma_j$ and $\theta_*$ is the  Brandon angle \cite{BrandonHighAngle} , often taken to be between $15^\circ$ and $30^\circ$.
Following the convention in materials science literature, we will refer to (\ref{eq:RSmodel}) as an {\em anisotropic} surface tension model, even though there is still no dependence on the normal to the interface.

In two dimensions, each grain $\Sigma_i$ is modeled as a square lattice and is assigned an orientation angle $\theta_i$ of a clockwise rotation about the origin that maps it back to the standard two-dimensional lattice $\mathbb{Z}^2$. Thus, the misorientation between two grains $\Sigma_i$ and $\Sigma_j$ is given by
\begin{equation}\label{eq:misorentation_angle2D}
\theta_{i,j} = \min_{k\in\mathbb{Z}} \abs{\theta_i-\theta_j+k\frac{\pi}{2}}.
\end{equation}
In three dimensions, the orientation of a grain with cubic lattice can be described (nonuniquely) by a matrix $g \in SO(3)$ that corresponds to the rotation required to obtain the lattice of the grain from the standard integer lattice $\mathbb{Z}^3$. The misorientation angle between two grains $\Sigma_i$ and $\Sigma_j$ is defined as
\begin{equation}\label{eq:misorentation_angle3D}
\theta_{i,j} = \min_{r\in\mathcal{O}} \arccos\left(\frac{1}{2}(\text{tr}( rg_ig_j^T)-1)\right),
\end{equation}
where $\mathcal{O}$ denotes the octahedral group (of symmetries of the cube in the three dimensions).

According to Mullins, the dynamics associated with energy \eqref{eq:energy} is given by $L^2$ gradient descent for the interfaces, leading to the normal speed
\begin{equation}\label{eq:normalspeed}
v_{i,j} = \mu_{i,j} \sigma_{i,j} \kappa_{i,j}
\end{equation}
for interface  $\Gamma_{i,j}$, where $\kappa_{i,j}$ denotes the mean curvature of $\Gamma_{i,j}$ and $\mu_{i,j}>0$ is a mobility factor. In addition, a condition known as the \emph{Herring angle condition} \cite{HerringAngleCondition} holds along triple junctions. At a junction formed by the meeting of the three distinct grains $\Sigma_i$, $\Sigma_j$ and $\Sigma_k$, this condition reads
\begin{equation}\label{eq:YoungsLaw}
\sigma_{i,j} n_{i,j} + \sigma_{j,k} n_{j,k} + \sigma_{k,i} n_{k,i} = 0.
\end{equation}
As a consequence the angles formed by normals $n_{i,j}$, $n_{j,k}$ and $n_{k,i}$ to the three interfaces $\Gamma_{i,j}$, $\Gamma_{j,k}$ and $\Gamma_{k,i}$ along the triple junction are determined by their associated surface tensions; this relation is also known as Young's law.

The algorithms used in this study are obtained from a non-local approximation to energy \eqref{eq:energy}
\begin{equation}\label{eq:approximate_energy}
\frac{1}{\sqrt{\delta t}} \sum_{i<j} \sigma_{i,j} \int_{\Sigma_i} G_{\sqrt{\delta t}} \ast \mathbb{1}_{\Sigma_j} \:dx
\end{equation}
where $G_t$ denotes the Gaussian kernel
\[
G_{\sqrt{t}}(x) = \frac{1}{(4\pi t)^{d/2}}\exp\left(-\frac{\abs{x}^2}{4t}\right)
\]
and $\mathbb{1}_\Sigma(x)$ denotes the characteristic function for a set $\Sigma$:
\[
\mathbb{1}_\Sigma(x) = \begin{cases}
1 & \text{if } x\in\Sigma,\\
0 & \text{otherwise}.
\end{cases}
\]
The width $\delta t$ of the Gaussian kernel appearing in \eqref{eq:approximate_energy} ends up playing the role of the time step size for our scheme, described below, that approximates gradient decent of \eqref{eq:energy} in $L^2$ sense, as described by \cite{MullinsGrainBoundaryMotion}. Energy \eqref{eq:approximate_energy} has been proposed in \cite{SelimFelix} and has been shown to converge in a very precise sense (namely, that of Gamma-convergence) to energy \eqref{eq:energy} in the limit $\delta t\to 0^+$. Intuitively, we have
\[
\frac{1}{\sqrt{\delta t}} \int_{\Sigma_i} G_{\sqrt{\delta t}} \ast \mathbb{1}_{\Sigma_j} \:dx \approx \frac{1}{\sqrt{\pi}}\text{Area}(\Gamma_{i,j})
\]
since the function $\frac{1}{\sqrt{\delta t}} \mathbb{1}_{\Sigma_i} G_{\sqrt{\delta t}} \ast \mathbb{1}_{\Sigma_j}$ approximates a delta function concentrating near $\Gamma_{i,j}$ as $\delta t \to 0^+$. The reason for our interest in this specific - perhaps unusual - approximation of Mullins' energy is that it generates exceptionally simple and efficient algorithms for simulating gradient descent dynamics associated with \eqref{eq:energy}. Indeed, it has been shown in \cite{SelimFelix} to lead, in a systematic way, to the correct multiphase, arbitrary surface tension analogue of a very fast algorithm known as \emph{threshold dynamics} that was originally proposed in \cite{MBO92,MBO94} for networks with all equal surface tensions (i.e., $\theta_{i,j} = 1$ for all $i$ and $j$) and mobilities.

We start by recalling Algorithm \ref{alg:EO}, the simplest version of the generalization of threshold dynamics to arbitrary surface tensions given in \cite{SelimFelix}. Its benefits include its unconditional stability (time step size $\delta t$ can be chosen arbitrarily large, constrained only by accuracy considerations), seamless handling of topological changes in any dimension and very low per time step cost $\mathcal{O}(N\log N)$ on a uniform grid of $N$ points, using the FFT to compute the convolutions. Yet, despite its extreme simplicity, it handles automatically and correctly all the essential features of the dynamics, including Herring angle conditions \eqref{eq:YoungsLaw} along junctions, and countless types of topological changes that may occur during the evolution. Note that Algorithm \ref{alg:EO} can be seen as a type of level set method, where the function
\[
-\psi_i^n(x) + \min_{j\neq i} \psi_j^n(x)
\]
plays the role of the level set function depicting the shape of the $i$-th grain at the $(n+1)$-th time step. It turns out that this very simple algorithm results in a misorientation dependent mobility $\mu_{i,j}$, associated with the boundary $\Gamma_{i,j}$ between grains $i$ and $j$ that is given by
\[
\mu_{i,j} = \frac{1}{\sigma_{i,j}}
\]
leading to the normal speed
\[
v_{i,j} = \mu_{i,j}\sigma_{i,j}\kappa_{i,j}(x) = \kappa_{i,j}(x)
\]
at any point $x\in\Gamma_{i,j}$.
In other words, using the Gaussian as the convolution kernel, Algorithm \ref{alg:EO} naturally prefers  equal (isotropic) {\em reduced mobilities}: Each interface in the network moves with the same (equal) multiple of its (mean) curvature, and the only difference from the isotropic model are the junction angles.

\boxalgorithm{\caption{(in \cite{SelimFelix})}\label{alg:EO}}{Given the initial grain shapes $\Sigma_1^0, \ldots, \Sigma_N^0$ and a time step size $\delta t$, obtain the grain shapes $\Sigma_1^{n+1}, \ldots, \Sigma_N^{n+1}$ at the $(n+1)$-th time from the grain shapes $\Sigma_1^n, \ldots, \Sigma_N^n$ at the end of the $n$-th times as follows:
\begin{enumerate}
	\item Compute the convolutions:
\[
\phi_i^n = G_{\sqrt{\delta t}} \ast \bm{1}_{\Sigma_i^k}.
\]
	\item Form the comparison functions:
\[
\psi^n_i = \sum_{j \not= i} \sigma_{i,j} \phi_j^n.
\]
	\item Update the grain shapes:
\[
\Sigma_i^{n+1} = \left\{x: \psi_i^n(x) < \min_{j\neq i} \psi_j^n(x)\right\}.
\]
\end{enumerate}
}

A more general algorithm was in fact given in \cite{SelimFelix} that allows arbitrary choice of surface tensions {\em and} mobilities, but is considerably more complicated and deviates from the simplicity of Algorithm \ref{alg:EO}.
In this paper, we use for the first time on large scale simulations a more recent version, namely Algorithm \ref{alg:ES}, that was proposed in \cite{SalvadorEsedogluSimplified}.
It is simpler and closer in spirit to Algorithm \ref{alg:EO} while allowing for both arbitrary surface tensions and arbitrary mobilities.
This is achieved by merely replacing the convolution kernel used in Algorithm \ref{alg:EO} by the positive sum of two distinct Gaussians, thus preserving the efficiency and simplicity of the original version.
As shown in \cite{SalvadorEsedogluSimplified}, in the special but important and very common case of {\em equal mobilities} (i.e., $\mu_{i,j} = 1$ for all $i\not= j$), a careful choice of the widths of the two Gaussians, denoted by the parameters $\alpha$ and $\beta$ in what follows, guarantees the retention of  the major benefits of the original Algorithm \ref{alg:EO}, namely unconditional stability and convergence to the correct Mullins' energy \eqref{eq:energy}, while allowing its extension to e.g. Read-Shockley surface tensions \eqref{eq:RSmodel}.

\boxalgorithm{\caption{}\label{alg:ES}}{Given the initial grain shapes $\Sigma_1^0, \ldots, \Sigma_N^0$ and a time step size $\delta t$, obtain the grain shapes $\Sigma_1^{n+1}, \ldots, \Sigma_N^{n+1}$ at the $(n+1)$-th time from the grain shapes $\Sigma_1^n, \ldots, \Sigma_N^n$ at the end of the $n$-th times as follows:
\begin{enumerate}
	\item Compute the convolutions:
\[
\phi^n_{1,i} = G_{\sqrt{\alpha\delta t}} * \mathbf{1}_{\Sigma_i^n} \mbox{ and }
\phi^n_{2,i} = G_{\sqrt{\beta\delta t}} * \mathbf{1}_{\Sigma_i^n}.
\]
	\item Form the comparison functions:
\[
\psi^n_i = \sum_{j \not= i} a_{i,j} \phi_{1,j}^n + b_{i,j} \phi_{2,j}^n
\]
where $a_{i,j}$ and $b_{i,j}$ are given by
\[
a_{i,j} = \frac{\sqrt{\pi}\sqrt{\alpha}}{\alpha-\beta}\left(\sigma_{i,j}-\beta\mu_{i,j}^{-1}\right) \quad \text{and} \quad b_{i,j} = \frac{\sqrt{\pi}\sqrt{\beta}}{\alpha-\beta}\left(-\sigma_{i,j}+\alpha\mu_{i,j}^{-1}\right).
\]
	\item Update the grain shapes:
\[
\Sigma_i^{n+1} = \left\{x: \psi_i^n(x) < \min_{j\neq i} \psi_j^n(x)\right\}.
\]
\end{enumerate}
}

\section{Previous Work

The exising literature on the statistics of grain shapes and sizes during annealing is absolutely vast, both from the experimental and the modeling and simulation side.
Our discussion will therefore be necessarily limited to a few representative studies closest to the questions addressed in the present paper.

There are a plethora of numerical algorithms developed for the simulation of Mullins' mesoscale model of grain boundary motion.
They include the Monte Carlo Potts, front tracking, phase field, level set, and threshold dynamics methods.
In principle, all these methods should converge to the same continuum model described by the same system of PDEs, provided that various parameters appearing in them are scaled appropriately, at least until topological changes start to take place (even then, one hopes there would be agreement in the vast majority of situations).
Nevertheless, there isn't always good agreement even among simulations.

Very large, well resolved simulations of normal grain growth were carried out by Mason et. al. \cite{MasonStatisticsGBM} using the front tracking method in both 2D and 3D.
The authors catalogue a large number of grain size and shape statistics,  for the two dimensional model, for the three dimensional model, and for two dimensional slices of the three dimensional model.
However, this work is restricted to the equal surface tension, equal mobility, i.e. the fully isotropic, version of Mullins' model.
That in particular means the exclusion of such statistics as the GBCD or the MDF which are predicated on variation of surface tension between different interfaces in the network.
Furthermore, left unadressed is the question of how sensitively these shape and size statistics depend on the surface tension and mobility model used in conjunction with Mullins' continuum model.

Upmanyu et. al. \cite{upmanyu} study the dependence of grain boundary statistics on the surface tension and mobility models used via Monte Carlo Potts and phase field simulations in 2D.
The GSDs they report, including when the surface tension, the mobility, and both are anisotropic, show little deviation from that of the isotropic model; they all follow a more or less log-normal distribution.
Their investigation of the MDF in the same vein reveals strong dependence on the surface tension model used, but not so much on mobilities:
The network is driven towards a prevalence of interfaces with low misorientations and therefore low surface tensions.
Other statistics, such as the number of neighbors (or faces), are also observed to have a much more pronounced dependence on surface tensions than mobilities.
Similar MC Potts simulations in \cite{HolmMisorientationAngles} also show little difference between the GSDs of isotropic vs. Read-Shockley surface tensions.

Gruber et. al. \cite{GruberMCPart1} carry out MC Potts simulations in 3D, using isotropic as well as Read-Shockley (anisotropic) surface tension models.
Among the grain properties they document are the GSD and the MDF.
They observe no difference between the stationary GSDs reached using isotropic vs. anisotropic surface tension models.
As regards the MDF, this study suggests surface tension anisotropy has a measurable effect on the steady state MDF reached in simulations that start from an initial data with random orientation texture, regardless of the dimension.
When the initial orientation texture is not random, no steady state for the MDF is observed in the case of anisotropic surface tensions: instead, a proliferation of low angle grain boundaries results in a gradual but persistent concentration of the MDF at the origin.

Kim et. al. \cite{Kim2006} carry out phase field simulations in 2D and 3D using isotropic surface tensions and mobilities. In 3D, the stationary GSD they reach (given in terms of the reduced equivalent radii) follows very closely the Hillert distribution.
On the other hand, the Monte-Carlo Potts simulations of Zoellner et. al. \cite{zoellner} in 3D, using isotropic surface tensions and mobilities, and starting from an initial data with random texture, results in a stationary GSD that is quite far from the Hillert distribution.
The more recent very large scale simulations of Miyoshi et. al. \cite{Miyoshi2017} in 3D using the phase field method also yielded a GSD similar to that of \cite{zoellner} that deviates substantially from the Hillert distribution.

The recent 2D simulations reported in \cite{MiyoshiAnisotropic} using a phase field method treat surface tension / mobility models (i) and (iii), as well as a third model that combines Read-Shockley surface tensions with mobilities that have a sigmoidal dependence on the misorientation angle as proposed by Humphreys \cite{Humphreys}.
For this third surface tension / mobility model, the stationary GSD they obtain is significantly different from those with models (i) \& (iii).
The stationary GSDs of models (i) \& (iii) are reported to be very similar.

\section{Large scale simulations

In this section, we use Algorithms \ref{alg:EO} and \ref{alg:ES} to study the evolution of large networks of grains, in both 2D and 3D.
The 2D simulations are initialized with approximately 100,000 well-resolved grains and are repeated 3 times, while the 3D simulations are initialized with 10,000 grains and repeated 10 times.
Recall from the Introduction the following surface tension and mobility models that will be explored in both cases:
\begin{enumerate}
\item Read-Shockley surface tensions with constant (isotropic) mobilities: $\sigma_{i,j}$ given by (\ref{eq:RSmodel}) and $\mu_{i,j} = 1$.
\item Read-Shockley surface tensions and reciprocal mobilities: $\sigma_{i,j}$ given by (\ref{eq:RSmodel}) and $\mu_{i,j} = \sigma_{i,j}^{-1}$.
This results in equal (isotropic) {\em reduced mobilities}.
\item Equal (isotropic) surface tensions and mobilities: $\sigma_{i,j} = \mu_{i,j} = 1$.
This also results in equal (isotropic) reduced mobilities.
\end{enumerate}
Whenever Read-Shockley surface tensions are used, we take the Brandon angle $\theta_*$ in (\ref{eq:RSmodel}) to be $\theta_* = 30^\circ$ to agree with simulations performed in Refs. \cite{GruberMCPart1,HolmMisorientationAngles} and to lie within the experimentally observed range \cite{SuttonBalluffiBook}.

In all simulations, the initial grain configuration is taken to be the Voronoi diagram of a set of points chosen uniformly at random from a periodic box.
\subsection{Two spatial dimensions
In these 2D simulations, the initial grain configuration of around 100,000 well-resolved grains is assigned a random fiber texture: All grain orientations are obtained from a reference configuration by a rotation about the axis normal to the 2D domain, and the angle of the rotation is taken from $[0,2\pi)$ uniformly at random.
We present results at two distinct times for each model: model (i) at $t_i = 1.907 \times 10^{-5}$ and $t_f = 8.345 \times 10^{-5}$, model (ii) at $t_i =  1.550 \times 10^{-5}$ and $t_f = 6.259 \times 10^{-5}$ and model (iii) at $t_i =  1.788 \times 10^{-5}$ and $t_f = 7.033 \times 10^{-5}$. The different stopping times ensure that at the intermediate time $t_i$, around 30\% of all grains remain in each one of our runs, while at the final time $t_f$ around 10\% of all grains remain.
The simulations are performed on a $4096 \times 4096$ grid discretizing $[0,1]^2$.

Algorithm \ref{alg:ES} requires choosing the parameters $\alpha$ and $\beta$, the widths of the two Gaussian kernels.
They need to satisfy
\[
\alpha > \max_{i\neq j} \sigma_{ij} \quad \text{and} \quad \beta < \min_{i\neq j} \sigma_{ij}
\]
to guarantee the properties described in \autoref{sec:model}.
However, from a practical point of view,  we also need the quotient $\alpha/\beta$ to remain moderate for accuracy considerations so that the Gaussians have comparable widths and thus may be sampled well without requiring unduly high spatial resolution.
Since both $\alpha$ and $\beta$ depend on the surface tensions which in turn depend on the misorientation angle, this can be accomplished by imposing a minimum misorientation angle between any two grains to be $1^\circ$.
This reduces the effective size of the surface tension matrix, while keeping the same number of distinct grains, and ultimately results in a small quotient $\alpha/\beta$.

We achieve this in the following way. Once the grain orientations are assigned, we partition the domain $[0,2\pi)$ into equally spaced bins of $1^\circ$ length. Then all grains with orientations in the same bin are reassigned a new grain orientation, the midpoint of the bin. This is equivalent to assigning grain orientations from the uniform distribution on $\{0.5^\circ, 1.5^\circ,\ldots,359.5^\circ\}$.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth] {microstructure2d_t0} & \includegraphics[width=0.28\textwidth]{microstructure2d_t1} & \includegraphics[width=0.28\textwidth]{microstructure2d_t2}
\end{tabular}
\caption{The initial microstructure contains $99,130$ grains (left). At the intermediate time $t_i = 1.907 \times 10^{-5}$, $29,841$ grains remain in full microstructure while at the final time $t_f = 8.345 \times 10^{-5}$, $10,003$ grains remain (right). These microstructures were obtained with model (i).}
\label{fig:microstructures2dzoomin}
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.48\textwidth]{microstructure2d_t2_zoomin} & \includegraphics[width=0.48\textwidth]{microstructure2dconstant_t2_zoomin}
\end{tabular}
\caption{Close up of the final microstructures for model (i) (left) and model (iii) (right). The former produces more elongated grains with narrow shapes, while the latter more rounded grains.}
\label{fig:microstructures2d}
\end{figure}

Figure \ref{fig:reduced_area2d} shows the time evolution of the GSD for the three different models, while Figure \ref{fig:reduced_effective_radii2d} shows the corresponding distributions of reduced effective radii.
For models (i) \& (ii) with Read-Shockley surface tensions, the distributions appear stationary with no significant difference between the two.
The distributions for the isotropic model (iii) are also very similar, except at final time $t_f$, we see some evidence of the two distinct peaks observed previously in simulations of \cite{MasonStatisticsGBM} with the same isotropic model, but using the very different algorithm of front tracking.
This two peak distribution has not been observed in other recent 2D simulations of the {\em isotropic} (constant) surface tension and mobility model (i) in e.g. \cite{Kinderlehrer2006,Kim2006,Miyoshi2018}, whose results look closer to the ones we obtained with Read-Shockley surface tensions.
On the other hand, the recent study \cite{MiyoshiAnisotropic}, which includes simulations with models (i) \& (iii), also has some evidence of the double peak with model (iii); our results are in close agreement with theirs for both of these models.
\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-areas_t0} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-areas_t1} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-areas_t2}\end{tabular}
\caption{GSD (probability density vs. reduced area $\frac{A}{<A>}$) for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:reduced_area2d}
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-radii_t0} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-radii_t1} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-radii_t2}
\end{tabular}
\caption{Probability density vs. reduced effective radii $\frac{\sqrt{A}}{\sqrt{<A>}}$ for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:reduced_effective_radii2d}
\end{figure}

The similarity in the GSDs of the three surface tension / mobility models considered here notwithstanding, it is natural to expect some differences when we move on to statistics of grain {\em shapes}.
After all, it is well known that the microstructure of anisotropic models can be visually quite different than that of the isotropic model (see Figures \ref{fig:microstructures2dzoomin} and \ref{fig:microstructures2d}).
Indeed, the similarities finally cease when we look at the isoperimetric ratio, also called circularity, which is given by $\frac{4\pi A}{P^2}$ where $A$ and $P$ denote the area and perimeter of the grain, respectively.
Figure \ref{fig:isopratio2d} shows that all three surface tension / mobility models lead to peak distributions but with a higher amplitude for the isotropic model. From time $t_i$ to time $t_f$, while $20,000$ disappear, the isoperimetric ratio retains its peaked distribution but with the peak moving towards 1 which indicates that the grains become more rounded.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-isopratio_t0} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-isopratio_t1} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-isopratio_t2}
\end{tabular}
\caption{Probability density vs isoperimetric ratio $\frac{4\pi A}{L^2}$ for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:isopratio2d}
\end{figure}

Finally, Figure \ref{fig:MDF2d} shows the MDF for the two unequal surface tension models: Read-Shockley with equal and reciprocal mobilities.
Due to its random fiber texture, the MDF of the initial data is uniform.
As time evolves, it concentrates at small misorientations corresponding to small surface tensions.
This is in agreement with previous studies, e.g.
\cite{HolmMisorientationAngles,HolmDimensionalEffects,Kinderlehrer2004,ElseyEsedogluSmerekaActaMaterialia}.
There is a slight but noticeable difference between the MDFs of the two models, which may in any case keep concentrating at the origin with further evolution: Model (i) has a more concentrated MDF compared to (ii), perhaps because low angle boundaries move slower in model (i) and hence persist longer.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-MDF_t0} &\includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-MDF_t1} & \includegraphics[width=0.28\textwidth]{simulation2dcomparison-4096x4096-MDF_t2}
\end{tabular}
\caption{Misorientation distribution for models (i) and (ii): initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:MDF2d}
\end{figure}

\subsection{Three spatial dimensions

For the 3D simulations, we initialize the microstructure with approximately 10,000 well-resolved grains, chosen as the approximate Voronoi diagram to points placed uniformly at random in a periodic box that is discretized into a $256 \times 256 \times 256$ uniform grid.
This results in approximately $12$ grid points across each grain in every direction at time $t=0$, which is roughly the same resolution as in our 2D simulations.
Each grain is assigned a fully random orientations from $SO(3)$.
We present results at two distinct times for each model: model (i) at $t_i = 4.578 \times 10^{-3}$ and $t_f = 7.629 \times 10^{-3}$, model (ii) and (iii) at $t_i =  1.831 \times 10^{-3}$ and $t_f = 4.578 \times 10^{-3}$.
Just like in the 2D simulations, the different stopping times ensure that at the intermediate time $t_i$, around 30\% of all grains remain in each one of our runs, while at the final time $t_f$ around 10\% of all grains remain.
As in our 2D simulations, we impose a minimum misorientation angle, this time of $5^\circ$, between the grains once again in order to keep the quotient $\alpha/\beta$ that appears in our algorithm at a moderate value.
Given that the minimum misorientation angle imposed here is larger than usual, a natural question arises: if it were taken to be smaller, would the different statistics change significantly? In order to investigate this, we run our simulations for model (ii) without doing any merging and observed that the statistics remained essentially the same.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth,trim={4cm 1cm 4cm 1cm},clip] {microstructure3d_t0} & \includegraphics[width=0.28\textwidth,trim={4cm 1cm 4cm 1cm},clip]{microstructure3d_t1} & \includegraphics[width=0.28\textwidth,trim={4cm 1cm 4cm 1cm},clip]{microstructure3d_t2}
\end{tabular}
\caption{The initial microstructure contains $9881$ grains (left). At the intermediate time $t_i = 4.578 \times 10^{-3}$, $3031$ grains remain in full microstructure while at the final time $t_f = 7.629 \times 10^{-3}$, 998 grains remain (right). These microstructures were obtained with model (i).}
\label{fig:microstructures3d}
\end{figure}

Unlike in our 2D simulations, we see substantial differences in the eventual GSD reached by surface tension / mobility model (i) vs. models (ii) \& (iii).
Figures \ref{fig:reduced_volume3d} and Figure \ref{fig:reduced_effective_radii3d} show the time evolution of the GSD in terms of the reduced volume and corresponding effective radii distributions, respectively.
We note that while the effective radii distributions of (ii) \& (iii) remain very close to each other throughout the evolution (and agree with previous 3D simulations of the same model via different numerical methods, e.g. \cite{Kim2006,ElseyEsedogluSmerekaRoyal,Miyoshi2018,MasonStatisticsGBM}), that of (i) is more spread out, and its peak is shifted in the direction of smaller grains -- a difference not seen in the 2D simulations with the same algorithms, and at very similar resolutions.
Figure \ref{fig:reduced_effective_radii3dexperimental} suggests the different eventual GSD observed with surface tension / mobility model (i) may not be inconsistent with experimental data:
Two experimental distributions available in existing literature, namely \cite{zhang} and \cite{rowenhurst}, display substantial deviation from the GSD seen in simulations of the isotropic model (iii) in our and previous simulations such as \cite{Kim2006,ElseyEsedogluSmerekaRoyal,Miyoshi2018,MasonStatisticsGBM}, in the same direction as our simulations with model (i), namely in the direction of smaller grains and more spread out distribution.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-volumes_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-volumes_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-volumes_t2}
\end{tabular}
\caption{GSD (probability density vs. reduced volume $\frac{V}{<V>}$) for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:reduced_volume3d}
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_t1-simulations} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_t2-simulations}
\end{tabular}
\caption{Probability density vs. reduced effective radii $\frac{V^{1/3}}{<V^{1/3}>}$ for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:reduced_effective_radii3d}
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{cc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_t1-experiments} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_t2-experiments}
\end{tabular}
\caption{Comparison of the GSD in terms of the reduced effective radii distributions in models (i), (ii) and (iii) at time $t_i$ when approximately 30\% of grains remain (left) and time $t_f$ when approximately 10\% of grains remain (right) with the experimental data from \cite{zhang} and \cite{rowenhurst}.}
\label{fig:reduced_effective_radii3dexperimental}
\end{figure}

Figure \ref{fig:isopratio3d} shows the distribution of the isoperimetric ratio, namely $\frac{36\pi V^2}{S^3}$, where $V$ denotes the volume and $S$ the surface area of the grain.
Remarkably, we see less difference between surface tension / mobility models (ii) and (iii) in 3D than we did in 2D simulations: their stationary isoperimetric ratio distributions are almost identical.
The distribution for model (i), however, is broader, indicating a greater presence of eccentric grains compared to the isotropic model (iii), as in our 2D results show in Figure \ref{fig:isopratio2d}.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_t2}
\end{tabular}
\caption{Probability Density vs Isoperimetric Ratio $\frac{36\pi V^2}{S^3}$ for the three different models: initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and at time $t_f$ when approximately 10\% of grains remains (right).}
\label{fig:isopratio3d}
\end{figure}

We now turn to statistics of grain shapes and sizes in 2D slices of 3D data, often called 3DX statistics in the literature.
To this end, we took a series of cross sections parallel to the faces of the cubic computational domain $[0,1]^3$.
They were spaced roughly five average grain diameters apart to reduce correlations between neighboring cross sections. Figure \ref{fig:reduced_effective_radii3DX} shows the distribution of GSD in terms of the reduced effective radii distributions.
The results are similar to the fully 3D grains: the surface tension / mobility models (ii) and (iii) are close to each other, while the distribution for model (i) is more spread out and its peak is achieved at smaller grains.
In Figure \ref{fig:isopratio3DX} we compare the distribution of the isoperimetric ratio, namely $\frac{4\pi A}{P^2}$, where $A$ denotes the area and $P$ the perimeter of the grain in the 2D slice. All surface tension / mobility models present almost identical peak distributions, with models (ii) and (iii) attaining a higher peak.

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_slices_3_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_slices_3_t1-experiments} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-radii_slices_3_t2-experiments}
\end{tabular}
\caption{Probability density vs reduced effective radii distributions $\frac{\sqrt{A}}{\sqrt{<A>}}$ taken from the cross-sections of constant $z$-value from the initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and $t_f$ when approximately 10\% of grains remain (right). Experimental data from \cite{hull} and \cite{groeber} is included for comparison.}
\label{fig:reduced_effective_radii3DX}
\end{figure}

\begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_slices_3_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_slices_3_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-isopratio_slices_3_t2}
\end{tabular}
\caption{Probability density vs isoperimetric ratio $\frac{4\pi A}{L^2}$ taken from the cross-sections of constant $z$-value from the initial Voronoi data (left), at time $t_i$ when approximately 30\% of grains remain (center) and $t_f$ when approximately 10\% of grains remain (right).}
\label{fig:isopratio3DX}
\end{figure}

Finally, we report asymptotic behavior of the MDF in our 3D simulations with surface tension / mobility models (i) \& (ii).
When the initial data is generated by assigning orientations at random sampled from the uniform distribution on $SO(3)$ to grains from a Voronoi construction, it is well known that the initial MDF is given by the Mackenzie distribution \cite{Mackenzie}.
With this initial data, both models (i) \& (ii) behave very similarly: Their MDFs stay in close vicinity of the Mackenzie distribution; see the first row of plots in Figure \ref{fig:MDF3d}.
This observed behavior is in close agreement with previous simulation results, e.g. \cite{HolmDimensionalEffects, HolmDimensionalEffects, GruberMCPart1}.
The very slight deviation from the Mackenzie distribution (observed in these previous studies as well) appears to be present mainly in model (ii), and not so much in (i).

We also explore the asymptotic behavior of the MDF when the initial configuration deviates substantially from the Mackenzie distribution.
To that end, we devised and implemented an algorithm to modify the randomly assigned initial orientations so that the initial MDF is perturbed substantially away from the Mackenzie distribution.
We leave the shapes of the initial grains unchanged, as the random Voronoi diagram.
To be specific, consider the function
\[
 \text{MDF}_\delta (x) = \sum_{i<j} \text{Area}(\Gamma_{ij}) b_\delta(x-\theta_{ij})
 \]
 for $\delta>0$ with
 \[
 b_\delta(x) = \frac{1}{\delta\sqrt{2\pi}}e^{-\frac{x^2}{2\delta^2}}.
 \]
 Without loss of generality assume that
 \[
 \sum_{i<j} \text{Area}(\Gamma_{ij}) = 1.
 \]
Then $\text{MDF}_\delta$ is a regularized approximation of the MDF of the grain network.
Let $T(x)$ denote the desired initial MDF.
Based on \eqref{eq:misorentation_angle3D}, $\text{MDF}_\delta$ can be treated as a function of $g_1,\ldots,g_N \in SO(3)$, the orientations of each grain, by expressing the misorientation angles $\theta_{ij}$ in terms of these.
We apply a steepest descent procedure on the $L^2$ norm of the difference between $\text{MDF}_\delta$ and $T$ with respect to the orientations $g_i$ of the grains; this updates the initial orientations $g_j$ until the initial MDF closely matches the desired distribution $T(x)$.
As this procedure has no preference for a particular orientation, we expect that the initial orientation distribution remains uniform.
In other words, out approach allows perturbing the MDF without introducing texture.

The second and third rows of Figure \ref{fig:MDF3d} show the evolution of the MDF when it is initially perturbed away from the Mackenzie distribution by this procedure, by two different sized (and quite large) perturbations.
In each case, the MDF evolves back to a very close vicinity of the Mackenzie distribution, which is quite different from earlier results in \cite{GruberMCPart1}, where no steady state MDF is observed in simulations with initial data whose MDFs differ from the Mackenzie distribution.
We expect the difference is due to the different ways of generating non-Mackenzie initial data: The method of \cite{GruberMCPart1} introduces texture (non-uniform orientation distribution), while ours does not. This is illustrated in Figure \ref{fig:ODF3d} where we plot the orientation angle distribution. Here, by orientation angle of a grain we mean its misorientation angle when measured against the sample reference grain, the one that is completely aligned with the axes. In the presence of a random texture the orientation angle distribution will match the Mackenzie distribution.

 \begin{figure}[h]
 \centering
 \begin{tabular}{ccc}
 \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-MDF_t0} &\includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-MDF_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-1-256x256x256-MDF_t2}\\
 \includegraphics[width=0.28\textwidth]{simulation3dcomparison-2-256x256x256-MDF_t0} &\includegraphics[width=0.28\textwidth]{simulation3dcomparison-2-256x256x256-MDF_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-2-256x256x256-MDF_t2}\\
 \includegraphics[width=0.28\textwidth]{simulation3dcomparison-3-256x256x256-MDF_t0} &\includegraphics[width=0.28\textwidth]{simulation3dcomparison-3-256x256x256-MDF_t1} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-3-256x256x256-MDF_t2}
 \end{tabular}
 \caption{Misorientation distribution function for models (i) and (ii) with different initial MDFs at the initial time (left column), at time $t_i$ when approximately 30\% of grains remain (middle column) and time $t_f$ when approximately 10\% of grains remain (right column).}
 \label{fig:MDF3d}
 \end{figure}

 \begin{figure}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.28\textwidth]{simulation3dcomparison-2-256x256x256-ODF_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-3-256x256x256-ODF_t0} & \includegraphics[width=0.28\textwidth]{simulation3dcomparison-4-256x256x256-ODF_t0}
\end{tabular}
\caption{Orientation angle distribution (probability density vs orientation angle) for the initial Voronoi data when the MDF is perturbed with the steepest descent procedure described here (left and center graphs whose corresponding MDFs displayed in the second and third row of Figure \ref{fig:MDF3d} left column, respectively) and with the method of \cite{GruberMCPart1}.}
\label{fig:ODF3d}
\end{figure}


\section{Topolectrical circuits design}

We design a circuit that implements both non-reciprocal and non-Hermitian couplings in a periodical network. Topological effects in electrical circuit networks have initially been considered for the Hofstadter model~\cite{PhysRevX.5.021031,PhysRevLett.114.173902}. The Laplacian formalism, with the complex admittance as its observable, has been introduced in~\cite{leecommphys}. It has established the notion of admittance band structures, which are fully accessible experimentally~\cite{helbigprb}. Topolectrical circuits are a newly emerging platform for synthetic topological matter which allows to implement any tight-binding model of any coupling range or dimensionality. As such, topolectric circuits manage to host nodal knots or links as zero admittance eigenvalues in their three-dimensional variant~\cite{Research2018Yu,lee2019topological}. Topological corner states as implementations of higher order topology can likewise be realized~\cite{quadrupolar,PhysRevB.100.045407}.
Our non-Hermitian non-reciprocal circuit network is described by the model Laplacian introduced in Eq.~(1). An idealized conceptual unit cell of the circuit is shown in Fig.~\ref{fig:circuit_schematic_simplified}a. The unit cell consists of two sublattice nodes, $A$ and $B$. The sublattices are connected in the unit cell by an $LC$-resonating circuit, which represents the Hermitian and reciprocal coupling $\ii \omega\,v(\omega) = \ii \omega C_1 + 1/(\ii \omega L_1) $. The frequency-dependence of its impedance allows us to modulate the intracell coupling and to resolve the fundamentally different parameter regimes of the model on the same circuit board. The $LC$ oscillator is characterized by its resonance frequency $\omega_0^{-1} = \sqrt{L_1 C_1}$, at which any current is blocked from flowing	through the component by an impedance divergence, such that the connectivity $v(\omega_0) = 0$. The coupling $C_2$ in between unit cells acts as the reciprocal intercell connection $\ii\omega \, r = \ii \omega C_2$.
Inevitable resistive losses in the circuit can, to first approximation, be modeled by an effective resistance $R_\text{x}$ connecting the nodes $A$ and $B$ within an unit cell.

The key non-Hermitian component in our circuit is the negative impedance converter with current inversion (INIC), which implements non-reciprocity~\cite{PhysRevLett.122.247702}. Its setup is shown in Fig.~\ref{fig:circuit_schematic_simplified}b. Through an operational amplifier (OpAmp) in a negative feedback configuration, we realize an antisymmetric connection, that implements a positive capacitance $C_3$ in one direction and a negative capacitance with equal magnitude in the opposing direction, leading to the breaking of reciprocity. Solving Kirchhoff's laws for the INIC leads to the reduced two-node Laplacian
\begin{align}
\begin{pmatrix}
I_\text{in} \\
I_\text{out}
\end{pmatrix} = 	\ii \omega C_3 \,
\underbrace{\begin{pmatrix}
	-1 &  1 \\
	- 1 & 1
	\end{pmatrix}}_{= J_\text{INIC}}
\begin{pmatrix}
V_\text{in} \\
V_\text{out}
\end{pmatrix}.
\end{align}
The in- and output currents are related by $I_\text{in} = - I_\text{out}$, meaning, that both currents flow in opposite directions. For this, the OpAmp must introduce a current source or sink in the system, with currents either flowing inwards to $V_0$ or outwards from $V_0$ in both directions. From the viewpoint of $V_\text{in}$, the INIC acts as a negative capacitance $-C_3$, whereas $V_\text{out}$ experiences the positive capacitance $C_3$. The INIC breaks reciprocity, $J_\text{INIC} \neq J_\text{INIC}^\top$, in an antisymmetric fashion as the off-diagonal elements differ by a minus sign. Its diagonal contribution is proportional to $\sigma_z$. To avoid it, we add a grounding term of $C_0 \equiv 2 C_3$ to the node at $V_\text{in}$ in order to change the diagonal terms to a unit matrix contribution.

The circuit in Fig.~\ref{fig:circuit_schematic_simplified}a is mathematically described by a Laplacian matrix, which in momentum space under PBCs, takes the form
\begin{align}
J_{\omega}(k)  = \ii \omega \Bigg[
&\left(C_1 - \frac{1}{\omega^2 L_1}\right) \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix}
+ C_2 \begin{pmatrix} 1 & -\ee^{-\ii k} \\ -\ee^{\ii k} & 1 \end{pmatrix}
+ C_3 \begin{pmatrix} -1 & 1 \\ -1 & 1 \end{pmatrix} \nonumber \\
&- \frac{\ii}{\omega R_\text{x}} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix} -\left( \frac{1}{\omega^2 L_0} + \frac{\ii}{\omega R_0}\right)\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + C_0 \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}
\Bigg]. \label{eq:SC Laplacian}
\end{align}
Additional to the introduced coupling terms, we include grounding terms $1/(\ii \omega \, L_0)$, $1/R_0$ at all nodes as well as $\ii \omega \, C_0$ at sublattice $A$. By omiiting the resistive intracell connection $R_x$ and representing equation \eqref{eq:SC Laplacian} in the Pauli matrix basis, we arrive at the effective circuit Laplacian stated in equation~(1).

The inductive grounding at each node is used as an $A$-$B$ symmetric grounding, which shifts the imaginary part of the admittance in a frequency-dependent fashion. We damp the circuit globally by grounding all nodes with a resistance $R_0$. The damping is needed to stabilize the circuit and eliminate instabilities (self-amplifying energy gains) in the system, that arise due to the INIC as an active element. The resistance $R_0$ shifts the real part of the admittance eigenvalue spectrum globally by $1/R_0$ but leaves the eigenmodes invariant.

\section{Experimental implementation}

A circuit consisting of ten unit cells was realized on a Printed Circuit Board (PCB). The concise values for the circuit components are detailed in Tab.~\ref{tab:Nominal values}. They are chosen to optimize several aspects of the experimental setup: The Lock-In amplifiers, which are employed for the measurement, operate up to an upper frequency bound of 100\,kHz, the choice of specific parameters facilitates a more stable circuit, in the way that it requires less damping to stabilize it and the non-local response is seen more dominantly if the parameters are adjusted to $r = 2 \gamma$.

In the implementation of the INIC, we employ the unity-gain stable operational amplifier LT1363.
In principle, the auxiliary components in the positive and negative feedback loop of the OpAmp can be chosen arbitrarily, provided they are equal. For a real OpAmp with its finite gain bandwidth product, limited stability conditions, and non-zero output impedance, however, a careful stability analysis~\cite{SuppMat} is required to guarantee the stable operation of LT1363 with our capacitive load $C_3$. This does not only concern the frequency range at which the experiment was conducted, but also frequencies in the MHz range to avoid undesired gain resonances in this regime. Considering those aspects we chose $R_\text{a} = \text{0.2}\,\Omega$ and $C_\text{a} = \text{0.94}\,$ $\mu$F.

As the INIC is an active, non-reciprocal circuit element it pumps energy into the system. This can lead to eigenmodes of the dynamical matrix, which are described by a complex energy eigenvalue with negative imaginary part, such that their magnitude increases exponentially with time. Hence, energy accumulates in the system until the OpAmp shows non-linear saturation effects and discontinues to function properly. In order to prevent this self-sustained energy gains, we add resistors $R_0$ from each node to ground, which consume the desired amount of power. The suited value for $R_0$, which differs between OBC and PBC, should be small enough to avoid instabilities but large enough to avoid oversized damping, which localizes the voltage response and likewise decreases the measurement accuracy. A numerical computation of the eigenfrequencies using the circuit Hamiltonian formalism~\cite{SuppMat} requires $R_0 \leq \text{45.0} \,\Omega$ for PBC and $R_0 \leq \text{66.8} \,\Omega$ for OBC to remove all global divergences. The discrepancy between those values roots in $\min(\mathfrak{Re}\,(J_\text{OBC})) \geq \min(\mathfrak{Re}\,(J_\text{PBC}))$ for the PBC and OBC spectrum, where a negative real part of the admittance is prone to cause global instabilities. In the actual experimental implementation, parasitic resistances help to stabilize the system and further increase the necessary $R_0$ to stabilize the system.
For the spectral measurements we chose $R_0 = \text{20}\,\Omega$, while for the non-local voltage response measurements we used $R_0 = \text{120}\,\Omega$.

In order to further eliminate undesired parasitic effects in the experiment and increase its accuracy, we design the circuit boards such that the unintended electromagnetic coupling between different wires and inductors is minimized. We add additional shielding around the inductors, use two circuit boards with five unit cells each, and place them on a metallic mesh. The shielding of the inductors changes their nominal values specified in Tab.+\ref{tab:Nominal values} and furthermore alters their frequency dependence. Moreover, we move the supply lines for the OpAmps out of the PCB board plane. In order to keep parametric disorder to a minimum and thereby preserve an approximate translational invariance, all circuit components were precharacterized by a BK Precision 894 LCR-meter and sorted into groups with only 1\% tolerance.

To perform the spectral measurements, a constant AC current is fed into one node of the board, while lock-in amplifiers with a high dynamical range are used to measure the voltage response of the circuit. The current is realized by a voltage source connected to the PCB through a shunt resistance of $R_\text{S} = \text{12.0}\,\Omega$. The PBC results were acquired by feeding an input current to sublattice $A$ and measuring the voltage at all the other nodes, then repeating the same procedure for sublattice $B$. The spectrum was then calculated through Fourier transformation, which is possible due to circuit periodicity.  For the OBC spectra, the current was fed at each node individually and the voltages of all the other nodes were measured. The results of these measurements are used to populate the Green’s function, as the inverse of the Laplacian,
\begin{align}
\nvec{V} = \begin{pmatrix}
G_{1,m} \, I_m \\
\vdots \\
G_{N,m} \, I_m
\end{pmatrix} = \begin{pmatrix}
G_{1,1} & \cdots & G_{1,N} \\
\vdots & \ddots & \vdots \\
G_{N,1} & \cdots & G_{N,N}
\end{pmatrix}
\cdot \begin{pmatrix}
0 \\ \vdots \\ I_m \\ \vdots \\ 0
\end{pmatrix} \qquad \Leftrightarrow \qquad G_{n,m} = \frac{V_n}{I_m}
\end{align}
where $I_m$ is the excitation current at node $m$, and $V_n$ is the measured voltage at node $n$. Using this construction, the Laplacian $J$, as well as its eigenvalues and eigenvectors can be computed numerically.

\section{Results and Analysis}

\subsection{Full-scale experimental fit of the theoretical model.}\label{subsec:Full-scale experimental fit of the theoretical model}

To fit the experimental results theoretically, we extend our description of the circuit to a full-scale fitted model that takes the most prominent imperfections into account. They occur as parasitic resistances in the inductive couplings of the circuit, and can be modeled as a serial resistor $R_{L_i}$ which modifies the inductances according to
$L_i \quad \rightarrow \quad  L_i + \frac{R_{L_i}}{\ii \omega }$,
as depicted in Fig.~\ref{fig:parasitics}. This results in the full-scale fitted model
\begin{align}\label{eq:full_scale_Laplacian}
J_{\omega,\text{fit}}(k) = \ii \omega \, \Bigg[&\left(C_1 + C_2 + \frac{C_0}{2} - \frac{1}{\omega^2 L_0 - \ii \omega R_{L_0} } - \frac{1}{\omega^2 L_1 - \ii \omega R_{L_1} } -\frac{\ii}{\omega R_0}\right) \, \mathbbm{1} \nonumber \\
-&\left(C_1 -\frac{1}{\omega^2 L_1 - \ii \omega R_{L_1} } +C_2 \cos(k)  \right) \, \sigma_x - \Big(C_2 \sin(k) - \ii\, C_3 \Big) \, \sigma_y \nonumber \\
&+\left(\frac{C_0}{2} -C_3\right)\sigma_z\Bigg],
\end{align}
with capacitances given by $C_1 =$ 300\,nF, $C_2 =$ 94\,nF, $C_3 =$ 47\,nF and $C_0 =$ 94\, nF. The serial resistance to ground, which ensures stability, is chosen to be $R_0 = 20 \, \Omega$ for spectral measurements in accordance with the stability condition evaluated in~\cite{SuppMat}.

In comparison to the inductors, the serial resistances of the capacitors are negligible. The filling material of the selected capacitances for the experimental implementation is chosen such that the capacitance is stable across the measured frequency regime. In contrast, the inductors and their corresponding serial resistances $L_1$, $L_0$, $R_{L_1}$, $R_{L_0}$ vary significantly with the AC driving frequency $f=\omega/2\pi$. We fit the PBC and OBC spectra of $J_{\omega,\text{fit}}(k)$ to the measured eigenvalues in Fig.~\ref{fig:Spectral measurements} (a1-a5), and obtain the fit parameters for $L_1$, $L_0$, $R_{L_1}$, and $R_{L_0}$ (see Tab.~\ref{tab:Fit parameters}).

Incorporating the parasitic resistance $R_{L_1}$ in the intracell connection given by $L_1$, the intracell parametric hopping $v(\omega) = C_1 -\frac{1}{\omega^2 \, L_1}$ is adjusted to
\begin{align}
v(\omega) \rightarrow \left(C_1 - \frac{L_1}{\omega^2 \, L_1^2 + R_{L_1}^2}\right) +\  \ii \  \left(-\frac{1}{\omega^2} \cdot\frac{ \, R_{L_1}}{\omega^2 \, L_1^2 +  \, R_{L_1}^2}\right) \equiv v_R+i v_I,
\end{align}
where $v_R$ and $v_I$ label its real and imaginary part. In a theoretical analysis, we are able to connect the limit of periodic and open boundary conditions for $J_{\omega,\text{fit}}(k)$, which differ by a connection with the boundary coupling of $\ii \omega C_2$. To find the full OBC admittance spectrum $j_\text{OBC}$ from the PBC spectrum $j(k)$, we recast the boundary coupling as a pumping of an effective imaginary flux. As detailed in section~\cite{SuppMat}, the OBC spectrum is reached at the smallest, in general finite $\kappa$, where $J(k + i \kappa)$, as an analytic continuation of momentum to the complex plane, recovers the Bloch state limit of plane waves.

The magnitude of the Bloch factor $\abs{\alpha}$ for the OBC modes relates to the imaginary momentum $\kappa$ by $\abs{\alpha} = e^{- \kappa}$. For the concrete model, a computation (see~\cite{SuppMat} with $v_R\leftrightarrow v$ and $v_I\leftrightarrow \beta$) yields
\begin{align}\label{eq:OBC Bloch factor}
\abs{\alpha} = \sqrt{\abs{\frac{v_R+ \ii \, v_I + \gamma}{v_R + \ii \, v_I - \gamma}}} = \sqrt[4]{\frac{(v_R+\gamma)^2+v_I^2}{(v_R-\gamma)^2+v_I^2}}.
\end{align}
The corresponding localization length of the bulk eigenmodes is given by
\begin{align}
\xi = \frac{1}{\kappa} = -\frac{1}{\ln(\abs{\alpha})} = \left[\frac{1}{4} \, \ln(\frac{(v_R-\gamma)^2+v_I^2}{(v_R+\gamma)^2+v_I^2})\right]^{-1},
\end{align}
which translates to the OBC bulk modes through $e^{-x/\xi} = e^{-\kappa x}$, where $x$ labels the circuit unit cells from left to right. The inverse localization length $\kappa$ is plotted in Fig.~\ref{fig:Phase diagram} as a function of $\omega$. The OBC eigenvalues are correspondingly given by all analytically continued PBC eigenvalues for different momenta $k$,
\begin{align}
j_{\text{OBC}} = \bigcup_k \,  j(k + i \kappa),
\end{align}
where $\kappa = \ln\sqrt[4]{[(v_R-\gamma)^2+v_I^2]/[(v_R+\gamma)^2+v_I^2]}$. For $v_R \neq 0$, the nonzero $\kappa$ yields boundary localized modes for OBC, with the intrinsic localization length of $\xi$.

In Fig.~\ref{fig:Spectral measurements}a, we plot the PBC and OBC admittance spectrum obtained from experimental measurements at differentfrequencies $\omega$. By continuous curves, we show the theoretical band structure obtained from the fitted model including all parasitic effects and circuit fit parameters chosen as detailed in Tab.~\ref{tab:Fit parameters}. For $v_R \neq 0$, we observe an extensive breakdown of bulk-boundary correspondence, as the PBC eigenvalues lying on closed loops coalesce to open arcs for OBC. The corresponding eigenvectors are calculated from the measurement data and displayed in Fig.~\ref{fig:Spectral measurements}b, showcasing the localization of all bulk voltage eigenmodes at either one boundary of the OBC circuit depending on the sign of $v_R(\omega)$. For $v_R <0$ as in Fig.~\ref{fig:Spectral measurements}\,(b1),(b2), all modes are localized at the left edge due to $\kappa > 0 $ ($\xi>0$) as shown in Fig.~\ref{fig:Phase diagram}. The localization length switches its sign as a function of frequency with $v_R >0$ as seen in Fig.~\ref{fig:Phase diagram} and converts the localization of the OBC modes to the right edge in Fig.~\ref{fig:Spectral measurements}\,(b4),(b5).

In our experimental frequency sweep, we cross one transition point, where all modes are delocalized for both PBCs and OBCs as depicted in Fig.~\ref{fig:Spectral measurements}\,(b3). It is found at $v_R= 0$ at
\begin{align}
\tilde{\omega}_0 = \sqrt{\frac{1}{L_1 \, C_1}- \frac{R_{L_1}}{L_1^2}},
\end{align}
yielding a transition frequency of $\tilde{f}_0 \approx \ $91.3\,kHz. At this point, the localization length $\xi$ for the bulk modes diverges, which is visualized for $\kappa \rightarrow 0$ in Fig.~\ref{fig:Phase diagram}. Bulk-boundary correspondence is restored, leading to a perfect matching of PBC and OBC eigenvalues in Fig.~\ref{fig:Spectral measurements}\,(a3) except for topological SSH eigenvalues. Consequently, the PBC spectrum lies on arcs and follow an extended reciprocity condition, stating that the PBC spectrum is reciprocal around a symmetry point $k_s$ in the Brillouin zone, $j(k_s+k) = j_s(k_s-k)$ (see~\cite{SuppMat}). In the fitted model, the PBC spectrum is symmetric around $k_s = -\arctan(\gamma/v_I) \approx \ $0.4$\pi$ at $\tilde{f}_0$.

For $v_I \neq 0$, due to the incorporation of parasitic resistances, the exceptional points in the OBC spectrum disappear. The strongest localization of the OBC bulk modes is found at the respective global maximum (left edge) and minimum (right edge) of the imaginary momentum $\kappa$, as depicted in Fig.~\ref{fig:Phase diagram}. Due to the small parametric adjustment of $v$ by the incorporation of $R_{L_1}$ as in $v_I(\omega) \ll v_R(\omega)$ and their respective non-trivial frequency dependencies, the frequencies of maximum localization of the OBC modes are marginally adjusted from $\abs{v} = \abs{\gamma}$ to $f_l \approx \ $84.5\,kHz and $f_r \approx \ $100.0\,kHz for the left and right edge localization respectively. Experimental data close to those frequencies is shown in Fig.~\ref{fig:Spectral measurements}\,(a2),(b2) as well as Fig.~\ref{fig:Spectral measurements}\,(a5),(b5), where we observe a strong localization of OBC modes.

In Fig.~\ref{fig:Spectral measurements}, we further notice the emergence of Su-Schrieffer-Heeger type topological modes, whose eigenvalues reside between the two bulk admittance bands. The topological modes are localized at the boundary of the sample. To find the regime, where topological states exist, we resort to the full-scale fitted model. In analogy to the argument given in~\cite{SuppMat}, a similarity transformation involving the Bloch factor $\alpha$ acting on the OBC Laplacian of $J_\text{fit}$ yields a fully reciprocal model. It presents an intracell hopping of
\begin{align}
t_0 = \sqrt{(v_R + \ii \,v_I)^2-\gamma^2}
\end{align}
and an intercell hopping of $r$. Two zero-admittance solutions corresponding to boundary-localized topological modes exist, for $\abs{t_0 } < \abs{r}$, which translates to
\begin{align}
\abs{v_R} < \sqrt{\sqrt{r^4-(2 \gamma v_I)^2}-v_I^2+\gamma^2}.
\end{align}
This leads to a frequency regime of
\begin{align}
f_L < f < f_U,
\end{align}
where we expect zero-admittance solutions in the OBC spectrum and the corresponding topological modes to exist. The zero eigenvalues are shifted according to the frequency-dependent unit matrix contribution in the fitted model. For the given parameters, the topological transition frequencies are given by $f_L \approx \ $78.8\,kHz and $f_U \approx \ $113.0\,kHz. The computed values agree with the measurement, as topological modes occur in Fig.~\ref{fig:Spectral measurements}\,(a2)-(a5), but not in Fig.~\ref{fig:Spectral measurements}\,(a1). Topological modes hence coexist with the localization of all bulk modes for OBC, \textit{i.e.\ }\/ the non-Hermitian skin effect. Furthermore, the emergence of topological modes is parametrically distorted away from the PBC complex gap-closing point (PBC Exceptional point) due to the non-Hermitian contribution.

We visualize the imaginary flux pumping causing the localization of the OBC modes in an experimental measurement (Fig.~\ref{fig:Spectral measurements}c) supported by a numerical computation (Fig.~\ref{fig:Spectral flow}). The spectral flow of the complex admittance spectrum from PBC to OBC is experimentally achieved through the attenuation of the boundary coupling, that modulates the boundary conditions in Fig.~\ref{fig:Spectral measurements}c.
For PBC the coupling is implemented by the capacitance $C_2$. To realize the imaginary flux pumping, we decrease the capacitance by adding additional capacitors in parallel and, eventually, open the connection completely to achieve OBC.

In Tab.~\ref{tab:alpha}, we show the theoretically expected magnitude of the OBC Bloch factor $\abs{\alpha(\omega)}$, as computed in equation \eqref{eq:OBC Bloch factor} while using the parameters given in Tab.~\ref{tab:Fit parameters} for the experimentally studied frequencies. The Bloch factor can also be obtained from $\abs{\alpha(\omega)} = e{-\kappa}$ with $\kappa$ being plotted in Fig.~\ref{fig:Phase diagram}. In comparison, we calculate the mean decay of the OBC bulk eigenmodes, as they are obtained from the measurement data in Fig.~\ref{fig:Spectral measurements}b. The experimental values agree with the theoretical data within the scope of the computed standard deviation.
This confirms the accuracy of the theoretical fitted model.
Discrepancies might root in additional parasitic effects and couplings in the realistic setting, which are not included in the fitted model as parasitic resistances of the capacitances, point-like resistive defects due to solder joints or capacitive couplings between conducting wires.

\subsection{Non-local voltage response}

Consider a linear periodic circuit array with periodic boundary conditions and a single current input at a specific circuit node. Then, the voltage response to the current feed must be localized around the input unit cell, \textit{i.e.\ }\/ the voltage profile decreases with distance from this unit cell.
In the limit of a dissipationless circuit, the response can also be delocalized but never increase in a different unit cell.
Otherwise, the system cannot be stable as the voltage accumulates infinitely in a periodic circuit.
Now consider the same circuit with open boundary conditions and a current input in the bulk of the system. If the circuit exhibited bulk-boundary correspondence, the eigenstates of the system's bulk would only change perturbatively and not effect the principal voltage response in comparison to the PBC case. By implication, a non-local voltage response in a periodic circuit is an inevitable signature of BBC breakdown and localized modes of extensive cardinality.

To compute the voltage response to an external current excitation $\nvec{I}$, we resort to the Greens function $G_J = J^{-1}$ of a system with $N$ unit cells. In its spectral representation, it is given by
\begin{align}
G_J = \sum_{n } \frac{1}{j_n} \cdot \frac{\nvec{V}_n \, \nvec{U}_n^\dagger}{\nvec{U}_n^\dagger \, \nvec{V}_n},
\end{align}
where $\nvec{V}_n$ ($\nvec{U}_n$) are the right (left) eigenvectors of the Laplacian matrix to the eigenvalue $j_n$. The voltage response $\nvec{V}$ is then given by $\nvec{V} = G_J\,\nvec{I}$. The eigenstates, which correspond to the smallest eigenvalues are excited the most.
Assume that one eigenvector, say $\nvec{V}_1$, is dominantly excited by the current profile $\nvec{I} = (0,\,\cdots,\,0,\,I_a,\,0,\,\cdots,\,0)^\top$ with only one nonzero contribution, \textit{i.e.\ }\/ one input current at node $a$. Then, the voltage response behaves as
\begin{align}\label{eq:non-local voltage response approx}
\nvec{V} \sim \frac{1}{j_1} \cdot \frac{U_{1,a}^*}{(\nvec{U}_1^\dagger \, \nvec{V}_1)} \, \nvec{V}_1 ,
\end{align}
where $U_{1,a} \sim \abs{\alpha}^{-a/2}$ assuming $\abs{\alpha}>1$ (right edge localization w.l.o.g.) as the left and right eigenvectors decay with $\abs{\alpha}^{-1}$ into the bulk. The normalization then follows to be $(\nvec{U}_1^\dagger \, \nvec{V}_1) \sim 2 \cdot N \, \abs{\alpha}^{-N}$ in the circuit system with $N$ unit cells. The voltage profile is approximated by
\begin{align}
\nvec{V} \sim \frac{1}{2 N \, j_1} \cdot \abs{\alpha}^{N- \frac{a}{2}} \cdot \nvec{V}_1.
\end{align}
From this formula, we make two observations. First, the voltage response increases with the nodal distance of the current feed to the localized voltage profile at the right edge. We can therefore associate an amplification factor proportional to $\abs{\alpha}$ to each unit cell in the chain. As we move the excitation further away from the mode localization, the magnitude of the response gets exponentially amplified with $\abs{\alpha}$. Second, the normalization factor $(\nvec{U}_1^\dagger \, \nvec{V}_1)$ decreases as we create a system with larger $N$. A feed on the edge opposite to the localized response is therefore also exponentially increased with the number of unit cells $N$. This makes sense, because the number of intermediate amplification cells is increased in the system. In the limit of $N\rightarrow \infty$, the normalization of left and right eigenmodes vanishes.

Similar to Sec.~\ref{subsec:Full-scale experimental fit of the theoretical model}, we find the mean localization of the OBC bulk eigenmodes, that were obtained from the experimental data to $\abs{\alpha}_\text{nl} = ($2.2 $\pm$ 0.3$)$. It agrees with the theoretically computed value in Tab.~\ref{tab:alpha} within the scope of the standard deviation, but the localization is weakened due to non-linear effects in the operational amplifiers, that become noticeable in the regime of small damping. In Fig.~\ref{fig:non-linear response}a, we show a logarithmic plot of the non-local voltage response with the circuit being excited at node 1 and node 7. It agrees with the expected amplification profile due to the localization of the OBC modes with the Bloch factor of $\abs{\alpha}_\text{nl}$, as shown by a straight line in Fig.~\ref{fig:non-linear response}. Fig.~\ref{fig:non-linear response}b depicts the amplification factor of the circuit with respect to the voltage level at the feed node. The maximum voltage decreases exponentially as a function of the nodal distance to the current feed. The expectation in Fig.~\ref{fig:non-linear response}b using  $\abs{\alpha}_\text{nl}$ as the unit-wise amplification factor agrees with the measurement data. Fig.~\ref{fig:non-linear response} highlights two statements on the non-local response. First, the voltage response retains the decay profile of the localized bulk modes with $\abs{\alpha}_\text{nl}$ and second, $\abs{\alpha}_\text{nl}$ also accurately models the amplification factor as a contribution from each unit cell.

Note, that if we decrease $R_0$ and thereby increase the damping in the system, we need to take the excitation of all modes into account and the approximation formulated in equation \eqref{eq:non-local voltage response approx} breaks down. With larger damping, we introduce interference between all modes in the system leading to a voltage response, that is more and more localized at the point of the current feed.

Fig.~\ref{fig:Phase diagram} displays the localization of OBC modes, where $\kappa = 0$ means full delocalization. We notice, that the maximum of $\abs{\kappa}$ is larger for right edge localized modes ($\kappa <0$) than for left edge localized ones ($\kappa <0$). This asymmetry is caused by the serial resistances, which exhibit a different frequency dependence than capacitive and inductive components. Their effect is therefore weakened for larger frequencies leading to a stronger localized voltage profile at the right edge than at left edge at the respective frequencies.

{\bf Caption Table 1:} Nominal values for the components used for the circuit implementation. (${}^*$) The potentiometer implementing $R_0$ was adjusted to $20\,\Omega$ for the measurements of the spectra and to $120\,\Omega$ for the measurement of the non-local voltage response.

{\bf Caption Table 2:} Frequency-dependent fitted inductances and their corresponding serial resistances $L_1, L_0, R_{L_1}, R_{L_0}$ at each experimentally studied frequency. Other experimental parameters are chosen to be $R_0 = 20 \, \Omega$, $C_1 =$ 300\,nF, $C_2 =$ 94\,nF, $C_3 =$ 47\,nF and $C_0 =$ 94\, nF.

{\bf Caption Table 3:} Theoretically calculated and experimentally fitted magnitude of the Bloch factor $\abs{\alpha}$ for the experimentally studied frequencies. Experimental values are obtained as a mean value of the exponential decay profile of all eigenmodes, that have been numerically computed from the measurement data. Their error is given by the standard deviation of the mean value while considering all eigenmodes of the system and therefore quantifies the agreement of the localization of all bulk states.

{\bf Caption Figure 4:} \textbf{(a)} Schematic diagram of one unit cell of the skin circuit including all grounding terms. An $LC$ resonating circuit models the frequency dependent intracell coupling $v(\omega)$. Connectivity between unit cells is modeled by a capacitance $C_2$, which is also depicted in Fig.~1 of the main text. \textbf{(b)} Negative impedance converter with current inversion used in the intracell $A-B$ connection based on an operational amplifier to generate the non-Hermitian, non-reciprocal coupling $\gamma$.

{\bf Caption Figure 5:} Schematic adjustment of the inductive coupling as we incorporate parasitic resistances in the circuit model.

{\bf Caption Figure 6:} \textbf{(a1-5)} Complex admittance spectra for frequencies $f_1=70$\,kHz, $f_2=\ $84.2\,kHz, $f_3= \ $91.5\,kHz, $f_4= \ $95.0\,kHz and $f_5=\ $98.5\,kHz with PBCs in red and OBCs in green. Discrete dots denote measurements, joint lines show the parasitics-corrected theoretical values. \textbf{(b1-5)} Magnitude of voltage modes $\abs{\psi}_n$ for OBCs with $n$ indicating the $n$-th eigenvector of the system, computed from experimental data. The topological modes are separated from the bulk modes by a green line. \textbf{(c)} Spectral flow of the admittance spectrum from PBC to OBC for $f= \ $80\,kHz quantified by $\eta$, which is the boundary bond normalized to its PBC value. Black dots are experimental data points, red and green lines indicate the admittance spectrum of the fitted model for PBC and OBC respectively.

{\bf Caption Figure 7:} Top: Variation of the intracell coupling $v(\omega)$ as a function of frequency. The two OBC exceptional points for the ideal model at $v = \pm \gamma$ and the transition point at $v=0$ with delocalized OBC modes are marked by black lines. Bottom: Localization length $\xi$ of the bulk OBC eigenmodes for the ideal effective model as a dashed red curve and for the full-scale fit model as a blue joint line.

{\bf Caption Figure 8:} Spectral flow of the admittance in the complex plane at the frequency $f = \ $80\,kHz for the effective model (left) and the aberration-corrected fit model (right). The flow is quantified by the parameter $\eta$, which describes the attenuation factor of the 1-20 bond.

{\bf Caption Figure 9:} \textbf{ (a)} Logarithmic plot of the absolute voltage response in a measurement involving an external current feed with AC driving frequency $f = \ $98.5\,kHz at the left edge of the sample at node 1 in red and at node 4 in blue. The voltage increases exponentially to the right side of the circuit due to a unit-wise amplification. The expected voltage profile normalized to the voltage at the feed node is shown as a straight line, where the slope has been numerically calculated as the average decay profile of the OBC bulk eigenmodes.  \textbf{ (b)} Logarithmic plot of the amplification factor of the circuit defined as the ratio of the maximum voltage $V_\text{max}$ to the voltage  $V_\text{feed}$ at the feed node. The voltage response increases with nodal distance to the right edge, which is supplemented by the expected drop in amplification calculated from the mode localization as in \textbf{ (a)}.

\clearpage

\begin{table}
	\centering
	\begin{tabular}{|c|c|c|c|c|c|c|c|}
		\hline
		$C_1$ [nF] & $C_2$ [nF] & $C_3$ [nF] & $L_1$ [$\mu$H] & $L_0$ [$\mu$H] & $R_0$ [$\Omega$] & $R_\text{a}$ [$\Omega$] & $C_\text{a}$ [nF] \\
		\hline
		300 & 94 & 47 & 10 & 33 & 20, $120^*$ & 0.2 & 940 \\
		\hline
	\end{tabular}
	\caption{}
	\label{tab:Nominal values}
\end{table}

\begin{table}
	\centering
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
		Frequency $f$ [kHz] & 70.0 & 80.0 & 84.2 & 91.5 & 95.0 & 98.5 \\
		\hline
		Inductance $L_1$ [$\mu$H] & 10.1 & 10.1 & 10.1 & 10.1 & 10.1 & 10.1 \\
		\hline
		Inductance $L_0$ [$\mu$H] & 33.0 & 32.0 & 31.7 & 30.3 & 30.1 & 30.0 \\
		\hline
		Serial resistance $R_{L_1}$ [m$\Omega$]& 220 & 230 & 250 & 280 & 300 & 320 \\
		\hline
		Serial resistance $R_{L_0}$ [m$\Omega$]	& 660 & 700 & 750 & 780 & 780 & 780 \\
		\hline
	\end{tabular}
	\caption{}
	\label{tab:Fit parameters}
\end{table}

\begin{table}
	\centering
	\begin{tabular}{|c|c|c|c|c|c|}
		\hline
		frequency $f$ [kHz]	& 70.0 & 84.2 & 91.5 & 95.0 & 98.5 \\
		\hline
		Theory (fit): $\abs{\alpha}_\text{theo}$  & 0.80 &  0.42 & 1.02 & 1.60 & 2.53 \\
		\hline
		Experiment:  $\abs{\alpha}_\text{exp}$ & 0.79 $\pm$ 0.03 & 0.42 $\pm$ 0.03 & 1.0 $\pm$ 0.1 &  1.6 $\pm$ 0.2 & 2.7 $\pm$ 0.3 \\
		\hline
	\end{tabular}
	\caption{}
	\label{tab:alpha}
\end{table}

\clearpage

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{SC_figure_comp_unit_cell_INIC.pdf}
	\caption{}
	\label{fig:circuit_schematic_simplified}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.7\linewidth]{SC_Figure_parasitics.pdf}
	\caption{}
	\label{fig:parasitics}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{SC_figure_eigenvalues_appendix_v2.pdf}
	\caption{}
	\label{fig:Spectral measurements}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.5\linewidth]{SC_Figure_phase_diagram.pdf}
	\caption{}
	\label{fig:Phase diagram}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{SC_figure_flow.pdf}
	\caption{}
	\label{fig:Spectral flow}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{SC_figure_amplification.pdf}
	\caption{}
	\label{fig:non-linear response}
\end{figure}

\end{document}