Strain-Induced Alignment in Collagen Gels

D. Vader1, A. Kabla2, D. Weitz1,3, L. Mahadevan1,∗

1 School of Engineering and Applied Sciences, Harvard University, Cambridge,
MA, USA

2 Department of Engineering, University of Cambridge, Cambridge, UK

3 Department of Physics, Harvard University, Cambridge, MA, USA

## Abstract

Collagen is the most abundant extracellular-network-forming protein in animal
biology and is important in both natural and artificial tissues, where it
serves as a material of great mechanical versatility. This versatility arises
from its almost unique ability to remodel under applied loads into anisotropic
and inhomogeneous structures. To explore the origins of this property, we
develop a set of analysis tools and a novel experimental setup that probes the
mechanical response of fibrous networks in a geometry that mimics a typical
deformation profile imposed by cells in vivo. We observe strong fiber
alignment and densification as a function of applied strain for both
uncrosslinked and crosslinked collagenous networks. This alignment is found to
be irreversibly imprinted in uncrosslinked collagen networks, suggesting a
simple mechanism for tissue organization at the microscale. However,
crosslinked networks display similar fiber alignment and the same geometrical
properties as uncrosslinked gels, but with full reversibility. Plasticity is
therefore not required to align fibers. On the contrary, our data show that
this effect is part of the fundamental non-linear properties of fibrous
biological networks.

Citation: Vader D, Kabla A, Weitz D, Mahadevan L (2009) Strain-Induced
Alignment in Collagen Gels. PLoS ONE 4(6): e5902.
doi:10.1371/journal.pone.0005902

Received March 6, 2009; Accepted April 21, 2009; Published June 16, 2009

Copyright: © 2009 Vader et al. This is an open-access article distributed
under the terms of the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided the
original author and source are credited.

Funding: This work was supported by NIH Bioengineering Research Partnership
grant R01 CA085139-01A2 and by NSF IGERT program in biomechanics. The funders
had no role in study design, data collection and analysis, decision to
publish, or preparation of the manuscript.

Competing Interests: The authors have declared that no competing interests
exist.
$\ast$ E-mail<EMAIL_ADDRESS>
## Introduction

Fiber networks, which arise in a range of natural and technological
situations, are prime candidates for a wide spectrum of applications requiring
tunable mechanical, transport and chemical properties [1]. In nature, these
networks self-assemble to generate important structural and functional
elements at various length scales: actin, intermediate filaments and
microtubules are the main components of the cytoskeleton [2]; spectrin confers
versatile qualities to red blood cell membranes [3]; fibrin is an essential
element in hemostasis [4]; collagen is the main component of the extracellular
matrix (ECM) in the animal kingdom [5] and cellulose is used by plants to
build cell walls [6].

The mechanical function of biological fiber networks is essentially two-fold:
(i) at the subcellular (actin, spectrin) and supracellular (collagen, fibrin)
scales, the material offers little resistance and high sensitivity to small
deformations, allowing it to be easily remodeled locally ; (ii) at larger
strains it stiffens strongly to ensure cell and tissue integrity [7]. The non-
linear stiffening, while observed in many biological systems [8, 9], is not
fully understood yet, with theories focusing on one of two broad mechanisms:
(i) microstructural nonlinearities of individual filaments [8], and (ii)
collective non-affine deformations of multiple filaments [10, 11]. To unravel
the relative importance of these mechanisms, a range of experimental tools
have been developed to quantify the network’s mechanical non-linearity in
systematic ways and relate the material micro-structure (network density and
morphology, fiber behavior) to the mesoscopic stress-strain laws. These tools
fall into two broad categories: simple shear in cone-plate or parallel plate
geometries, and uniaxial/biaxial stretch.

Simple shear deformations are commonly used to study purified protein
networks. This technique requires low sample volumes and provides a consistent
set of experimental tools and generic protocols to probe the visco-elastic
properties of soft gels in both the small-strain (linear) and large-strain
(non-linear) regimes, and in addition, normal stresses can be measured. Recent
data collected by Janmey et al. [9] show in particular that sheared
biopolymers exert negative normal forces, a fact that is in contradiction with
the hyperelastic behavior of other well studied elastomers. The broad
availability of experimental data in that geometry has encouraged a large
number of related theoretical and numerical studies [12, 13, 14], focused
primarily on the linear response of the material. However, since simple shear
rheology assumes that the material undergoes purely isochoric deformations in
the limit of small strains, it only allows for partial exploration of material
behavior. In particular, these experiments do not allow one to study
completely the non-linear regime (strain typically larger than 10%) that is
most relevant in many biological situations (single cell or tissue
deformation). And furthermore, it does not allow for a probe of the
dilatational rheology of the networks.

In contrast, at mesoscopic scales, uniaxial and biaxial testing are most
common for tissue mechanical characterization [15, 7, 16, 17] and have been
used to study reconstituted collagen networks [18, 19, 20, 21, 22], the
simplest tissue equivalents [23]. In contrast with simple shear, uniaxial
stretch generically leads to non-isochoric deformations, and hence allows one
to measure quantities such as the material Poisson ratio which can have values
as large as $\nu\approx 3$ for strongly deformed collagen gels. These values
arise in highly anisotropic materials, as reported for instance for solid
foams [24], and it is somewhat surprising to see similar behavior in in vitro
collagen gels which display little or no anisotropy in their undeformed state.

To understand this, we recall that early studies [25] on cell/matrix
interactions show that cells or groups of cells tend to generate tensile
forces on the extracellular environment . When cell colonies were plated on
fibrous materials such as collagenous gels, Harris and Stopak reported the
formation of anisotropic and denser regions connecting these cellular
assemblies, and showed that the matrix structure has a strong influence on
cell motility. Although these observations are well accepted, little is known
about the mechanical response of a fibrous matrix subject to an internal local
strain. Neither of the mechanical characterizations described previously focus
on how deformation changes the microstructure at the fiber scale, an issue of
particular importance in the large strain regime, that is all too easy to
observe (Figure 1).

In this paper, we use collagen type I gels as a model system to address this
question and shed light on the morphological evolution of both the fiber
111Although the expression ”collagen fiber” traditionally refers to large
bundles of collagen fibrils, we will use the words ”fiber” and ”fibril”
interchangeably in this paper to refer to the 0.5 micron diameter bundles. and
the network on an externally imposed stretching strain. Collagen is a
convenient biomaterial for biomechanical studies for a number of reasons: a)
it is readily available in large amounts, which makes it suitable for
milliliter-size gels; b) in vitro reconstituted networks have fibers that are
easily identifiable using confocal microscopy; c) many of its properties have
been extensively studied [26, 27, 28, 29, 30, 18]; d) the large diameter of
the fibers ($\approx 0.5\mu m$ for collagen fibrils [21]) and the stability of
the network [31, 32] make it easy to handle and image over a range of spatial
and temporal scales; e) fibrillogenesis is conveniently controlled in vitro by
pH, temperature and concentration [28, 26].

Figure 1:  Collagen gel morphological changes induced by presence of cells.
(A) Single U87 glioblastoma cell in a collagen network 10 hours after gel
polymerization. bar=50 $\mu$m. (B) Several U87 cells on the surface of a
collagen gel 10 hours after gel polymerization. bar=200 $\mu$m. (C) Two cell
colonies embedded in a collagen matrix 48 hours after gel polymerization.
bar=200 $\mu$m. Fibers (artificial red color) are imaged through confocal
reflectance; cell nuclei (green) are labeled with a GFP-histone heterodimer.

We first verified the presence of cell-induced alignments and densification
with our experimental system. As shown in figure 1A, an isolated human
glioblastoma cell (see Methods) in a collagen network induces stress
variations and modifies the network texture in its vicinity [33, 34]. Several
isolated cells on the surface of a collagen gel produce fiber alignment and
network densification along lines connecting individual cells (figure 1B).
Following Stopak and Harris, we also observed fiber alignment on macroscopic
length scales when we introduce large cell assemblies in the same extra-
cellular environment (figure 1C). Since active matrix remodeling is restricted
to the vicinity of living cells, such an effect can only be accounted for by
the mechanical properties of the network.

With these observations in mind, we employ a specific experimental approach
and develop a set of tools to quantitatively study the coupling between strain
and the morphology of fibrous networks in a range of strain and strain rates
that are typical of many biomechanical situations. Experiments on cell
colonies suggest that such a process can be conveniently studied at the
millimeter scale, and over a time-scale of a few hours. However, instead of
using cells to deform the extra-cellular matrix, we use an external imposed
displacement to stretch collagenous samples and monitor the gel response. In
particular, this experiment allows us to probe a range of dynamical regimes
and independently tune the biochemistry (crosslinking) to study the coupling
of tensile strain to network density and fiber orientation in a controlled
setting and investigate the origin and generality of these mechanical
processes. This also allows us to address the outstanding question of the
mechanical reversibility of these patterns in an extracellular environment.

## Methods

### Network synthesis

In vitro collagen networks are prepared according to a previously described
cell culture-compatible protocol [35], with a final collagen concentration
ranging from 0.5 to 4.0mg/mL. Solutions consist of 10% 10X minimum essential
medium (Invitrogen, Carlsbad, CA), 10% fetal bovine serum (JRH Biosciences,
Lenexa, KS), 1% penicillin-streptomycin (Invitrogen), bovine collagen diluted
to desired concentration (from 3.1mg/mL or 6.4mg/mL batch, Inamed
Biomaterials, Fremont, CA), a few $\mu$L of 1M sodium hydroxide (NaOH, Sigma,
St. Louis, MO) to bring pH to neutral, 50mM sodium bicarbonate (NaHCO3, Sigma)
buffer and deionized water. 800$\mu$L of solution are pipetted onto glass-
bottom Petri dishes (MatTek, Ashland, MA) Samples then polymerize for 30-60
minutes in a cell incubator at 37∘C, 5% CO2. After polymerization, samples are
20mm in diameter and 1-2mm in height.

In addition to untreated in vitro collagen gels, we also prepare polymerized
samples, to which we add glutaraldehyde (GA) - a common cell and tissue
fixative. The effect of this is an increase in the overall stiffness of the
gel by at least an order of magnitude (from a few tens of Pascals to over 1kPa
for a 1mg/mL gel, as evaluted using a cone-plate rheometer), without
noticeably changing the structure of the collagen network (fiber width and
length, mesh size). 2mL of 4% v/v GA (Sigma, St. Louis, MO) in water is
pipetted onto the sample, which is then incubated once more for at least 2
hours. It is subsequently rinsed twice with deionized water. Before use, all
samples are immersed into 2mL of deionized water, to allow the gel to swell,
and to reduce friction.

### Cell experiments

U87 human glioblastoma cells are cultured as described in [35]. After
passaging the cells, tissue equivalents are generated by diluting the cells to
approximately $10^{5}$/mL in an unpolymerized collagen solution at 1.5mg/mL.
After polymerization in a 37∘C temperature- and humidity-controlled incubator,
the spacing between individual cells, as seen in figures 1a and 1b, is on the
order of 100$\mu$m. Cell colonies (or spheroids), with an estimated $10^{3}$
cells, are generated with the hanging droplet method [36] and subsequently
seeded in 500$\mu$L of collagen solution at 1.5mg/mL shortly before it
polymerizes.

### Bulk rheology

We use an AR-G2 (TA Instruments, New Castle, DE) rheometer with a 4∘, 40mm
cone-plate geometry with a 109$\mu$m gap. 1.2mL of collagen solution is
pipetted onto the 37∘C preheated bottom plate of the rheometer and the cone is
lowered onto the sample. We use a solvent trap to prevent the sample from
drying during the measurement. During polymerization, the increase in G′ and
G′′ is probed by continuously oscillating the sample at a fixed 0.5% strain
amplitude and at a frequency of 0.2Hz. The oscillatory strain sweep is
performed at the same frequency and temperature, after the gel has polymerized
for 2-3hrs. The strain amplitude is increased logarithmically until the sample
breaks.

During oscillatory strain sweeps, we simultaneously record the maximum stress
and strain of the sample for each oscillatory cycle. To characterize the onset
of stiffening from the stress-strain data, we define the critical strain as
the value at which the stress $\sigma$ exceeds the product
$G^{\prime}_{0}\gamma$ by more than 10%, where $G^{\prime}_{0}$ is the elastic
modulus in the small-strain linear regime.

### Mechanical setup

We place a polymerized collagen sample onto a glass cover slip and perforate
it with two rough-ended 1mm-diameter glass cylinders (capillaries) (see
figures 2A,B), which are gently pushed all the way to the glass bottom to
prevent the collagen from slipping beneath them. Each cylinder is attached to
two secondary transverse elastic capillary rods, themselves attached to two
linear transducers (Newport, Irvine, CA) controlled by the ESP300 controller
(Newport). The transverse capillaries act as springs that allow to maintain
contact with the bottom cover slip of the dish with constant pressure. The
tips, initially 1cm from each other, can then be moved apart at speeds ranging
from 0.125 to 12.5 $\mu m/s$, corresponding to strain rates of 2.5$\cdot
10^{-5}$ to 2.5$\cdot 10^{-3}$ per second; this range includes measured rates
of cell-induced contraction [37]. The movement of the tips results in the
local stretching of the gel sitting between them. For imaging purposes, the
whole mechanical setup (motors and tips) is clamped to the microscope sample
holder plate.

Figure 2:  Mechanical setup and sample imaging. (A) Side and (B) top views of
the mechanical setup used to deform the network; the collagen gel has a
pancake-like shape, typically 1mm in thickness and 20mm in diameter. As
defined in our experiments, the stretch axis is $x$. Drawn to scale, the two
squares represent the fields-of-view of the wide-field fluorescence images
(5x) and confocal reflectance images (60x). (C) Correlation of multiple slices
over time gives an estimate of the interslice distances, and hence vertical
strain. (D) Collagen network (blue) obtained with confocal reflectance.
Fluorescent tracers (pink) are embedded in the network. Scale bar 20$\mu m$.
(E) Wide-field fluorescence image of the embedded beads. Scale bar 500$\mu m$.

### Confocal imaging

We use a Zeiss LSM 510 Meta (Carl Zeiss Microimaging Inc., Thornwood, NY)
equipped with a 488nm Argon laser line and several photomutliplier tubes. We
set the Meta channel of the microscope (which allows for selection of specific
wavelengths) to detect wavelengths between 474 and 494nm to allow for the fact
that we work in reflectance mode [38], which has the significant advantage of
avoiding the use of fluorescent dyes in the samples. A 60X 1.2-NA Olympus
water immersion objective (Olympus America Inc., Center Valley, PA) is mounted
onto the microscope.

While deforming the sample, we acquire timelapse 2D confocal images at various
heights between 50 and 150$\mu$m from the bottom surface. Tracking in-plane
deformations in multiple slices improves the statistics of our analysis;
moreover, tracking out-of-plane motion via image correlation allows us to
estimate the vertical deformation of the sample (figure 2C). The timelapse
interval is 10s, which corresponds to 0.25% imposed deformation at the typical
strain rate.

### Image segmentation and fiber detection

Fiber orientation is calculated for each 2D confocal slice by several image
processing steps: the raw images are filtered using a 2D Gaussian blur and
subsequently thresholded so that at 0% stretch, 10% of the pixels are above
that threshold; this threshold value is applied to all subsequent images of
the same experiment. A circular window of diameter 15px moves across the
thresholded image, and the 2nd order moment tensor $M$, defined below, is
locally calculated using binary pixel weights: below the threshold level, the
pixel weight $A_{ij}$ is 0, above that level, $A_{ij}$ is 1. Typically, the
moment tensor quantifies the spatial distribution of weight around a center of
mass and its eigenvalues give an indication as to whether weights are
distributed isotropically around the center or not. Similarly here, $M$
quantifies the distribution of pixel intensities around the intensity-weighted
center of mass, with $X_{ij}$ and $Y_{ij}$ the pixel coordinates in the local
circular window:

$M\\!=\\!\\!\\!\\\ \left[\\!\\!\\!\\!\begin{array}[]{cc}\
\\!\\!\\!\scriptstyle\sum A_{ij}(X_{ij}-\bar{X})^{2}&\
\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\scriptstyle\sum
A_{ij}(X_{ij}-\bar{X})(Y_{ij}-\bar{Y})\\\ \ \\!\scriptstyle\sum
A_{ij}(X_{ij}-\bar{X})(Y_{ij}-\bar{Y})&\
\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\scriptstyle\sum
A_{ij}(Y_{ij}-\bar{Y})^{2}\end{array}\\!\\!\\!\\!\right]$ (1)

$r=\frac{\hbox{Max}\left\\{\hbox{Eigenvalues}(M)\right\\}}{\hbox{Min}\left\\{\hbox{Eigenvalues}(M)\right\\}}$
(2)

The ratio $r$ characterizes the aspect ratio of binarized image fragment
enclosed in the sliding window and can be used to detect fibers: a single
fiber passing through the middle of the window, with a diameter less than half
the window size will yield a high aspect ratio; a single fiber in the corner
of the window or multiple fibers in the window will give a low aspect ratio.
We only keep the regions with $r>2$ for which the eigenvector corresponding to
the higher eigenvalue of $M$ provides the local fiber orientation $\phi$.

The choice of window size (15 pixels) is a compromise between increasing
angular resolution at low strains and avoiding multiple fibers in the window.
As shown later on, the fiber density strongly increases at high strain,
resulting in less accuracy in the orientation analysis. In practice, this sets
the limits of the method. See Methods S1, Figure S1 and Figure S2 for
validation of this image processing algorithm.

### Orientational order parameter

To quantify the network anisotropy as a function of the applied deformation,
we calculate an orientation tensor $\Omega$ and an associated nematic order
parameter $\mu$ from the distribution of the fiber orientation $\phi$ at each
time step, defined by:

$\Omega=\left(\begin{array}[]{cc}\left<\cos^{2}(\phi)\right>&\left<\cos\phi\sin\phi\right>\\\
\left<\cos\phi\sin\phi\right>&\left<\sin^{2}\phi\right>\end{array}\right)$ (3)
$\mu=\hbox{Max}\left\\{\hbox{Eigenvalues}(2\;\\!\Omega-\hbox{\bf
Id})\right\\}$ (4)

where $\bf Id$ is the identity matrix, and $<\cdot>$ denotes spatial averaging
over the domain of interest. The order parameter ranges from $\mu=0$ for a
uniform angular (isotropic) distribution to $\mu=1$ for a perfectly aligned
system. Although this parameter, based on a 2D image analysis, only
characterizes the order in the horizontal plane, it does provide a suitable
signature of the microstructure evolution, and in particular its non-linear
behavior.

### Deformation field at the mesoscale

We also characterize, at the mesoscopic millimeter scale, the strain field
$\varepsilon_{xx}$, $\varepsilon_{yy}$ and $\varepsilon_{zz}$ induced by the
imposed displacement of the glass tips. Vertical strain is estimated by
following the individual displacement of confocal z-stacks. For each slice of
a stack taken at a time $t$, this is done by calculating the correlation with
the slices obtained at a neighboring time $t+\Delta t$ (see figure 2C). The
height of the slice which provides the largest correlation value indicates the
new location of the material layer and this information is used to calculate
the vertical component of the strain.

The deformation in the $xy$ plane is measured by a PIV (particle imaging
velocimetry) method. 1$\mu m$ diameter rhodamine Fluospheres (Invitrogen,
Carlsbad, CA), at a volume ratio of 1:100000, are used as tracers to measure
the deformation of the gel at millimeter length-scales. Most of the carboxy-
coated particles stick to the network, as seen in figure 2D. Using a 5x lens
on a Zeiss wide-field microscope and focusing on the middle of the sample, we
image the local density of these particles, which displays heterogeneities as
seen in figure 2E. An image cross-correlation technique is then used to track
these heterogeneities at a scale of $10-50\mu m$ as the network is
progressively deformed. To identify the local displacement of a mesoscopic
region of the gel located at $(x,y)$ from time $t$ to $t+\Delta t$, we extract
a domain of 48x48 pixels surrounding $(x,y)$ at $t$ and look for the best
matching region - i.e. the one that maximizes the cross-correlation function -
in the image obtained at $t+\Delta t$. The cross-correlation function is
defined as:

$\rho_{AB}=\sum{(A_{ij}-A_{avg})(B_{ij}-B_{avg})/(\sigma_{A}\sigma_{B})}$ (5)

where A and B are the pixel intensity values associated with the two regions
of interest, $A_{avg}$ and $B_{avg}$ are the local average pixel intensities,
and $\sigma_{A}$ and $\sigma_{B}$ are the standard deviations of intensity
values of those regions.

This tracking through cross-correlation process is iterated over time to
extract the full trajectory of the material point and corresponds to a
Lagrangian description of the material. For a grid deformation that is fairly
homogeneous in the field of view of the microscope (millimeter scale), we use
a least-squares planar fit of the nodal displacement to quantify the material
deformation at the mesoscale. The deviation from the fit provides a measure of
the error on the deformation field.

These measurements allow us to extract a number of strain and stress
characteristics. In particular, in our geometry, the normal stresses
$\sigma_{yy}$ and $\sigma_{zz}$ are negligible as the gel is not attached on
the lateral sides, and this allows us to estimate the incremental Poisson
ratio, defined as:

$\nu_{xy}=-\partial\varepsilon_{yy}/\partial\varepsilon_{xx}$ (6)

which characterizes the coupling of incremental deformations in orthogonal
directions. In two dimensions and in the small deformation limit, the area of
the grid $A(\epsilon_{xx})$ is related to the Poisson ratio via the following
relationship:

$\frac{1}{A(0)}\frac{dA}{d\epsilon_{xx}}=1-\nu_{xy}$ (7)

Most materials respond to tension with a slight increase in their area (volume
in three dimensions), which for an isotropic material translates into the
condition:

$\nu_{xy}\leq 1$ (8)

The analogous condition in 3 dimensions is $\nu_{xy}\leq 0.5$. At large
deformations for an isotropic material, the criterion is slightly more
complex, but the critical strain, beyond which the change of area (or volume
in 3D) with respect to elongational strain becomes negative, remains of the
same order in practice.

## Results

### Rheological characterization of the samples

The range of collagen concentrations we work with (0.5-4.0 mg/mL) display mesh
sizes from 1-5 $\mu$m (measured through analysis of confocal reflectance
slices as in [35]) and span over two orders of magnitude in shear modulus (see
figure 3A), with values in close agreement with previously reported data [35,
39, 40]. The linear shear modulus $G^{\prime}$ (measured at small deformation)
has a strong dependence on the concentration $c$, with a behavior well
approximated by $G^{\prime}\sim c^{3}$ over the probed range. The onset of
nonlinear strain-stiffening typically occurs at strains of the order of 5%
(see figure 3B and [28]) and has only a weak dependence on collagen
concentration (see figure 3A inset).

Figure 3:  Bulk rheology measurements. (A) Linear elastic modulus as a
function of collagen concentration. A power-law of 3 is shown for comparison.
Inset: the onset of non-linear mechanical behavior, as defined in Methods. (B)
Elastic modulus $G^{\prime}$ during oscillatory strain sweeps and as a
function of collagen concentration. Figure 4:  Typical results of a sample
stretching experiment at micro- and meso-scale. The two montages of 5 images
each show, for two different 1mg/mL samples and at various strains, (A) wide-
field fluorescence images of beads embedded within the network - the super-
imposed grid is the result of the tracking of the deformation field (scale bar
500$\mu$m); (B) direct imaging of the fibers through CRLSM (scale bar
50$\mu$m). In inset, each corresponding $\phi$ histogram, with angle values
going from 0 to $\pi$. For the same samples depicted above: (C) represents the
evolution of the orientation statistics; the color at each point corresponds
to the relative count of fibers oriented along a specific direction at a given
strain $\varepsilon_{xx}$. (D) shows the order parameter $\mu$ resulting from
the data in (C); the curve beyond 15% stretch is grayed out due to the lack of
confidence of the order parameter when the high value of the density prevents
a proper detection of the fibers (see Methods). (E) gives the deformations
$\varepsilon_{yy}$ and $\varepsilon_{zz}$ as a function of the local strain
$\varepsilon_{xx}$. In inset, the incremental Poisson ratio $\nu_{xy}$ as a
function of the imposed deformation.

### Fiber alignment and non-linear Poisson effect

Figures 4A and 4B illustrate the evolution of the network microstucture over
the domain as it is stretched. At low strains ($<5\%$), no particular
alignment can be observed; however, above this threshold both fiber alignment
and network density increase. This is quantified in the figures 4C and 4D
where we show, as a function of the applied strain, the probability
distribution of local in-plane fiber orientation and the resulting order
parameter.

For a prescribed displacement of the capillary tips, we also characterized the
gel microstructure as we move away from the axis connecting the two
capillaries. Both the alignment and fiber density (expressed as fraction of
pixels above a given threshold) decrease (figure 5). This picture is, as
expected, in direct agreement with the alignment pattern induced by cell
colonies pulling on extra-cellular matrix shown in figure 1, where fiber
alignment and density are maximum along the axis joining the colonies, and
decay away from it.

Figure 5:  Characterization of the spatial variations of collagen fiber
alignment and densification. (A-C) Confocal reflectance images of 1mg/mL
collagen network stretched up to 15%; images are located at (A) 0, (B) 0.5 and
(C) 2 millimeters from the stretching axis. bars = 50$\mu m$. (D) Order
parameter and density as a function of distance from the stretching axis, for
the same sample. Letters A, B and C on figures (D) and (E) refer to
corresponding images above. (E) Drawn to scale, locations of images (A), (B)
and (C) with respect to stretching axis.

In figure 4E we show the average induced strain components $\varepsilon_{yy}$
and $\varepsilon_{zz}$ as a function of the externally applied strain
$\varepsilon_{xx}$. When $\varepsilon_{xx}\leq 2\%$ there is little
contraction in the transverse direction so that the Poisson ratio
$\nu_{xy}\sim 0$. As the applied strain increases, the material first thins by
contracting in the $z$ direction when $2\%\leq\varepsilon_{xx}\leq 5\%$, and
only when $\varepsilon_{xx}\geq 5\%$ does it also contract in the transverse
$y$ direction, with observed values of $\nu_{xy}$ as high as 5 (figure 4E
inset). The lag in response between these two directions can be attributed to
the sample geometry as well as a slight initial anisotropy of fiber
orientations in the $yz$ plane [21]; here we consider only the properties of
the planar projection of the network. The large in-plane incremental Poisson
ratio $\nu_{xy}$ quantifies the change in local fiber density and is
consistent with the confocal observations of densification. In order to
quantify and compare these changes with an independent measure of the geometry
of deformation, we define a critical strain $\epsilon_{crit}$ for which
$dA(\epsilon_{crit})/d\varepsilon_{xx}=0$, beyond which the areal strain
(yellow frame on figure 4A) starts to decrease with the applied strain. We
find that $\epsilon_{crit}\sim 5\%$, consistent with the critical strain
observed for fiber alignment.

The critical deformation, as defined above from the kinematic behavior in
local stretching tests, can be directly compared with the strain associated
with the mechanical stiffening measured in rheological experiments (see figure
3). These two quantities, measured independently, show good correlations in
their values and trends. For unfixed collagen samples, critical strain values,
measured either from rheological meansurements or from the kinematics, range
from a few percents at high collagen concentration to 15 % at low
concentration (see figure 6). The critical strain is therefore very weakly
dependant on the concentration, in particular compared with the variation of
the elastic modulus at small deformation that varies over more than two orders
of magnitude in the same concentration range. GA-fixed samples, whose
stiffness is estimated at an order of magnitude higher than their non-fixed
counterparts, also show similar values for $\epsilon_{crit}$. Taken together,
these results show that the strains above which the gel behavior becomes non-
linear as evidenced i) from the elastic modulus for a simple shear geometry
and ii) from the Poisson ratio in the local stretching experiments are related
with each other and only weakly sensitive to physical scales such as the
actual value of the elastic modulus.

Figure 6:  Comparison of critical strains measured from bulk rheology and
2-point stretching. Comparison of critical strain measured from the PIV method
($\epsilon_{crit}$ such that $dA(\epsilon_{crit})/d\varepsilon_{xx}=0$, see
Results) on locally stretched samples with the onset of strain stiffening
($\gamma_{crit}$ such that $\sigma>1.1G^{\prime}_{0}\gamma$, see Results)
obtained from rheological measurements.

### Orientational ordering is an elastic effect

Strain-induced alignment arises a priori from a combination of reversible
elastic effects and irreversible inelastic effects. To disentangle these two
contributions, we apply repeated strain cycles to the sample. All pure type-I
collagen samples display very little reversibility regardless of their
concentration once the imposed strains exceed about 10%; the gel never
recovers its initial configuration once the capillary tips return to their
initial location. In figures 7A,B we show the evolution of the fiber
orientation histograms as the material is cyclically stretched up to strains
of 15%: fiber aligment is permanently imprinted (7C-E).

Figure 7:  Characterization of the reversibility of alignment and
densification in crosslinked and uncrosslinked collagen samples. Response to
cycles of deformation for untreated collagen gels (A-E) and gels treated with
glutaraldehyde (F-K). (A) Histograms of the fiber orientation as a function of
the imposed strain for untreated sample; color at each point corresponds to
the relative count of fibers oriented along a specific direction at a given
strain $\varepsilon_{xx}$. (B) Resulting order parameter for untreated sample.
(C-E) Confocal reflectance images of a 1mg/mL collagen sample cycled 4 times
up to 15% stretch and back to 0%, showing the extent of the reversibility at
the microstructural level: (C) beginning of 1st stretch cycle; (D) middle of
1st stretch cycle; (E) end of 1st stretch cycle; bars = 50$\mu m$. (F)
Histograms of the fiber orientation as a function of imposed strain for
glutaraldehyde-treated sample; color-coding as in (A). (G) Resulting order
parameter and density for glutaraldehyde-treated sample. (H-J) Confocal images
of a 1mg/mL collagen samples with glutaraldehyde, cycled 4 times up to 15%
stretch and back to 0%: (H) beginning of 3rd stretch cycle; (I) middle of 3rd
stretch cycle; (J) end of 3rd stretch cycle; bars = 25$\mu m$. (K) Mesoscopic
response to 15% strain cycles for various strain rates at a constant imposed
strain amplitude. The bold line corresponds to the strain rate used for all
other experiments.

With the addition, post-polymerization, of GA to our samples (see Methods), we
change collagen’s material properties (elastic modulus, plasticity threshold)
without changing the microstructure of the network. For small deformations,
the initial response is similar to that of the unfixed sample, indicating that
fiber alignment and the anomalous Poisson effect are only weakly sensitive to
the fiber material properties. Furthermore, during cycles of applied
deformation - i.e. ramping up the applied strain to 15% and returning back to
0% at the same rate - GA-crosslinked samples exhibit near-reversibility at the
microscale, with fiber images of successive cycles being almost identical
(figures 7H-J). Figures 7F and 7G quantify this reversibility in terms of the
orientation histogram and the order parameter over multiple cycles. Consistent
with this behavior, we find (see figure 7K) that the local deformation field
at the mesoscale evolves along a reversible path for the same strain rate
($2.5\cdot 10^{-4}$/s). Taken together, these results show that fiber
plasticity, though observed for pure collagen samples, is unimportant in
determining alignment at microscopic scales and the large Poisson ratio at
larger scales. That is, strain-induced alignment is primarily an elastic
effect.

The reversible behavior characterized above might then serve as a baseline to
study more subtle effects (e.g. time-dependent and/or irreversible processes)
that can be observed after many cycles or different strain rates. We see, in
particular, a slight decrease in the amplitude of the alignment with the
number of cycles, indicating that fiber plasticity still occurs, although on
much larger timescales. In figure 7K we show the effects of strain rate on the
response of the system for a fixed strain amplitude. At lower strain rates the
system does not recover completely after a full cycle. This offset in the
response can be attributed to a slow creeping process occurring over a time-
scale of a few hours. At larger strain rates, one expects to see dynamic
effects related to viscous dissipation. Previous studies on the poroelasticity
of collagen networks [41, 18, 20, 42, 43], have reported equilibration times
ranging from a few seconds to a few hours, probably due to the diversity of
geometries and setups used for all these measurements. We find that dynamic
effects associated with higher strain rates induce an asymmetry in the
response between loading and unloading, influencing the unloading curve more
than the loading curve. This shows that the material responds with different
time-scales in extension and compression; this, in turn, suggests that
different physical processes are involved during the loading which is
dominated by fiber stretching and unloading which is dominated by fiber
bending.

## Discussion

In this study, we used an experimental and computational approach to quantify
the emergence of fiber alignment as a collagen sample is stretched. This
behavior is consistent with observations of cell-induced morphological changes
in tissue equivalents and sheds light on a biologically relevant material non-
linearity that arises from stress heterogeneities in fiber networks. Fiber
alignment at the microscale results in tissue densification when boundary
conditions allow it; this leads to high values of the Poisson ratio at large
deformations, and is observed for all studied collagen concentrations, with or
without addition of glutaraldehyde.

The reversibility of fiber alignment and gel densification, seen in
crosslinked collagen samples, show that these effects are primarily elastic.
Experiments on a piece of synthetic felt [44] have demonstrated that geometry
alone can account for such a behavior based on the generic non-linearity of
individual fibers, stiff in extension, but soft in compression
(bending/buckling). Under uniaxial loading, a tensile stress is necessarily
balanced in the microstructure by a compressive load on fibers normal to the
stretch direction, leading to a collapse of the material in the normal
direction and a strong enghancement of fiber alignment along the load
direction, as observed here. This argument also explains the correlation
between the moment where stress builds up (onset of non-linearity in
rheological measurements) and the critical strain associated with the Poisson
effect. This is also consistent with the bulk rheology experiments performed
by Janmey et al. [9] who report large negative normal stresses in all tested
fibrous materials. The normal stress in volume-constrained geometries (as in
simple shear flows) is the counterpart of high Poisson ratios in unconstrained
tests such as the local stretch performed in the present study.

Our findings seem to contrast with the experimental results of Tower et al.
[19] on stretched collagen samples, where either irreversible alignment (pure
collagen sample) or early fracture with negligible alignment (GA crosslinked)
occurs, suggesting that fiber plasticity is a key player in the alignment
process. The difference with the present study can be readily explained from a
geometrical standpoint: Tower et al. used a geometry in which the width of the
stretched region is comparable to its length; this might significantly
constrain transverse motion of the material and prevent local volume changes
to their full extent. By contrast, the 2-point stretching device we use
ensures that the material is free to move along the direction normal to the
load.

Our experiments have demonstrated the importance of sample geometry and
boundary conditions on the microstructure and mechanical response of
reconstituted biopolymer gels. For functional tissues, it is known that
mechanical properties are often finely controlled by the texture of the
underlying extracellular matrix as well [45]. Understanding the mechanisms
leading to such organization is an important step in learning how it happens
in the formation of natural tissues and for developing strategies to engineer
suitable tissue equivalents. We have shown here that ECM texturization can be
brought about simply by applying a deformation with a purely mechanical
device, without any intervention due to active modeling by cells. This
external perturbation is applied post-polymerization, in contrast with
previous reports of fiber alignment induced by a flow [46] or a magnetic field
[47, 48] during collagen polymerization. This provides direct support to the
in vivo studies of post-polymerization collagen texturization in developing
tissues [49]. However, the microscopic origin of the permanent texturization
occurring in our uncrosslinked samples remains unclear. Fibril-fibril
junctions are likely to be where plastic deformation occurs, allowing fibrils
to slide with respect to one another and thus inducing irreversible changes in
the network topology. One of the effects of glutaraldehyde crosslinking is to
strengthen these junctions and reduce the amount of plastic reorganization
allowed at the network level.

Whatever their origin, geometrical changes in the network structure are known
to influence, in turn, cell behavior. Our mechanical setup allows for the
dynamical control of network texture in a passive system, but clearly can be
used to study how active cells respond to externally-induced anisotropy. Such
experiments will provide more insights into specific mechanisms of
mechanotransduction and cell behavior, which are crucial to processes such as
morphogenesis, stem cell differentiation, metastasis and wound healing. More
generally, the external control of fiber alignment post polymerization not
only offers a convenient way to design anisotropic tissue equivalents with
collagen, but can also, when applied to other biopolymer systems, shed light
on a range of analogous phenomena, such as actin gel contraction [50] or
platelet-fibrin interactions [51], where microscopic agents interact through
the network and lead to large scale evolution and reorganization of matter.

## References

  * 1. Star A, Tu E, Niemann J, Gabriel JCP, Joiner CS, et al. (2006) Label-free detection of dna hybridization using carbon nanotube network field-effect transistors.  Proc Natl Acad Sci U S A 103: 921–926.
  * 2. Alberts B (2002) Molecular Biology of the Cell.  Garland Science.
  * 3. Gov N, Zilman AG, Safran S (2003) Cytoskeleton confinement and tension of red blood cell membranes.  Phys Rev Lett 90: 228101.
  * 4. Weisel JW (2004) The mechanical properties of fibrin for basic scientists and clinicians.  Biophys Chem 112: 267–276.
  * 5. Palsson B, Bhatia S (2003) Tissue Engineering.  Prentice Hall.
  * 6. Taiz L, Zeiger E (2002) Plant Physiology.  Sinauer Associates.
  * 7. Sacks MS (2000) Biaxial mechanical evaluation of planar biological materials.  Journal of Elasticity 61: 199-246.
  * 8. Storm C, Pastore JJ, MacKintosh FC, Lubensky TC, Janmey PA (2005) Nonlinear elasticity in biological gels.  Nature 435: 191–194.
  * 9. Janmey PA, Georges PC, Hvidt S (2007) Basic rheology for biologists.  Methods Cell Biol 83: 3–27.
  * 10. Heussinger C, Frey E (2006) Floppy modes and nonaffine deformations in random fiber networks.  Phys Rev Lett 97: 105501.
  * 11. Wyart M, Liang H, Kabla A, Mahadevan L (2008) Elasticity of floppy and stiff random networks.  Phys Rev Lett 101: 215501.
  * 12. Head DA, Levine AJ, MacKintosh FC (2003) Deformation of cross-linked semiflexible polymer networks.  Phys Rev Lett 91: 108102.
  * 13. Wilhelm J, Frey E (2003) Elasticity of stiff polymer networks.  Phys Rev Lett 91: 108103.
  * 14. Onck PR, Koeman T, van Dillen T, van der Giessen E (2005) Alternative explanation of stiffening in cross-linked semiflexible networks.  Phys Rev Lett 95: 178102.
  * 15. Billiar KL, Sacks MS (1997) A method to quantify the fiber kinematics of planar tissues under biaxial stretch.  J Biomech 30: 753–756.
  * 16. Ivancic PC, Coe MP, Ndu AB, Tominaga Y, Carlson EJ, et al. (2007) Dynamic mechanical properties of intact human cervical spine ligaments.  Spine J 7: 659–665.
  * 17. Cheng T, Gan RZ (2007) Mechanical properties of stapedial tendon in human middle ear.  J Biomech Eng 129: 913–918.
  * 18. Knapp D, Barocas V, Moon A (1997) Rheology of reconstituted type i collagen gel in confined compression.  Journal of Rheology 41: 971-993.
  * 19. Tower TT, Neidert MR, Tranquillo RT (2002) Fiber alignment imaging during mechanical testing of soft tissues.  Ann Biomed Eng 30: 1221–1233.
  * 20. Girton TS, Barocas VH, Tranquillo RT (2002) Confined compression of a tissue-equivalent: collagen fibril and cell alignment in response to anisotropic strain.  J Biomech Eng 124: 568–575.
  * 21. Roeder BA, Kokini K, Robinson JP, Voytik-Harbin SL (2004) Local, three-dimensional strain measurements within largely deformed extracellular matrix constructs.  J Biomech Eng 126: 699–708.
  * 22. Krishnan L, Weiss JA, Wessman MD, Hoying JB (2004) Design and application of a test system for viscoelastic characterization of collagen gels.  Tissue Eng 10: 241–252.
  * 23. Tranquillo RT (1999) Self-organization of tissue-equivalents: the nature and role of contact guidance.  Biochem Soc Symp 65: 27–42.
  * 24. Lee T, Lakes R (1997) Anisotropic polyurethane foam with poisson s ratio greater than 1.  Journal of Materials Science 32: 2397-2401.
  * 25. Harris AK, Wild P, Stopak D (1980) Silicone rubber substrata: a new wrinkle in the study of cell locomotion.  Science 208: 177–179.
  * 26. Raub CB, Suresh V, Krasieva T, Lyubovitsky J, Mih JD, et al. (2007) Noninvasive assessment of collagen gel microstructure and mechanics using multiphoton microscopy.  Biophys J 92: 2212–2222.
  * 27. Gentleman E, Lay AN, Dickerson DA, Nauman EA, Livesay GA, et al. (2003) Mechanical characterization of collagen fibers and scaffolds for tissue engineering.  Biomaterials 24: 3805–3813.
  * 28. Roeder BA, Kokini K, Sturgis JE, Robinson JP, Voytik-Harbin SL (2002) Tensile mechanical properties of three-dimensional type i collagen extracellular matrices with varied microstructure.  J Biomech Eng 124: 214–222.
  * 29. Brightman AO, Rajwa BP, Sturgis JE, McCallister ME, Robinson JP, et al. (2000) Time-lapse confocal reflection microscopy of collagen fibrillogenesis and extracellular matrix assembly in vitro.  Biopolymers 54: 222–234.
  * 30. Bozec L, Horton M (2005) Topography and mechanical properties of single molecules of type i collagen using atomic force microscopy.  Biophys J 88: 4223–4231.
  * 31. Leikina E, Mertts MV, Kuznetsova N, Leikin S (2002) Type i collagen is thermally unstable at body temperature.  Proc Natl Acad Sci U S A 99: 1314–1318.
  * 32. Persikov AV, Brodsky B (2002) Unstable molecules form stable tissues.  Proc Natl Acad Sci U S A 99: 1101–1103.
  * 33. Friedl P, Z nker KS, Br cker EB (1998) Cell migration strategies in 3-d extracellular matrix: differences in morphology, cell matrix interactions, and integrin function.  Microsc Res Tech 43: 369–378.
  * 34. Kim A, Lakshman N, Petroll WM (2006) Quantitative assessment of local collagen matrix remodeling in 3-d culture: the role of rho kinase.  Exp Cell Res 312: 3683–3692.
  * 35. Kaufman LJ, Brangwynne CP, Kasza KE, Filippidi E, Gordon VD, et al. (2005) Glioma expansion in collagen i matrices: analyzing collagen concentration-dependent growth and motility patterns.  Biophys J 89: 635–650.
  * 36. Kelm JM, Timmins NE, Brown CJ, Fussenegger M, Nielsen LK (2003) Method for generation of homogeneous multicellular tumor spheroids applicable to a wide variety of cell types.  Biotechnol Bioeng 83: 173–180.
  * 37. Pizzo AM, Kokini K, Vaughn LC, Waisner BZ, Voytik-Harbin SL (2005) Extracellular matrix (ecm) microstructural composition regulates local cell-ecm biomechanics and fundamental fibroblast behavior: a multidimensional perspective.  J Appl Physiol 98: 1909–1921.
  * 38. Gunzer M, K mpgen E, Br cker EB, Z nker KS, Friedl P (1997) Migration of dendritic cells in 3d-collagen lattices. visualisation of dynamic interactions with the substratum and the distribution of surface structures via a novel confocal reflection imaging technique.  Adv Exp Med Biol 417: 97–103.
  * 39. Parsons JW, Coger RN (2002) A new device for measuring the viscoelastic properties of hydrated matrix gels.  J Biomech Eng 124: 145–154.
  * 40. Velegol D, Lanni F (2001) Cell traction forces on soft biomaterials. i. microrheology of type i collagen gels.  Biophys J 81: 1786–1792.
  * 41. Wille JJ, Elson EL, Okamoto RJ (2006) Cellular and matrix mechanics of bioartificial tissues during continuous cyclic stretch.  Ann Biomed Eng 34: 1678–1690.
  * 42. Chandran PL, Barocas VH (2004) Microstructural mechanics of collagen gels in confined compression: poroelasticity, viscoelasticity, and collapse.  J Biomech Eng 126: 152–166.
  * 43. Ozerdem B, T zeren A (1995) Physical response of collagen gels to tensile strain.  J Biomech Eng 117: 397–401.
  * 44. Kabla A, Mahadevan L (2007) Nonlinear mechanics of soft fibrous networks.  J R Soc Interface 4: 99–106.
  * 45. Shadwick RE (1999) Mechanical design in arteries.  J Exp Biol 202: 3305–3313.
  * 46. Lee P, Lin R, Moon J, Lee LP (2006) Microfluidic alignment of collagen fibers for in vitro cell culture.  Biomed Microdevices 8: 35–41.
  * 47. Barocas VH, Girton TS, Tranquillo RT (1998) Engineered alignment in media equivalents: magnetic prealignment and mandrel compaction.  J Biomech Eng 120: 660–666.
  * 48. Guo C, Kaufman LJ (2007) Flow and magnetic field induced collagen alignment.  Biomaterials 28: 1105–1114.
  * 49. Stopak D, Wessells NK, Harris AK (1985) Morphogenetic rearrangement of injected collagen in developing chicken limb buds.  Proc Natl Acad Sci U S A 82: 2804–2808.
  * 50. Bendix PM, Koenderink GH, Cuvelier D, Dogic Z, Koeleman BN, et al. (2008) A quantitative analysis of contractility in active cytoskeletal protein networks.  Biophys J 94: 3126–3136.
  * 51. Shah JV, Janmey PA (1997) Strain hardening of fibrin gels and plasma clots.  Rheological Acta 36: 262-268.


# Symmetry-protected non-Abelian geometric phases in optical waveguides with
nonorthogonal modes

Julien Pinske  Stefan Scheel<EMAIL_ADDRESS>Institut für
Physik, Universität Rostock, Albert-Einstein–Straße 23-24, D-18059 Rostock,
Germany

###### Abstract

The generation of non-Abelian geometric phases from a system of evanescently
coupled waveguides is extended towards the framework of nonorthogonal coupled-
mode theory. Here, we study an experimentally feasible tripod arrangement of
waveguides that contain dark states from which a nontrivial
$\mathrm{U}(2)$-mixing can be obtained by means of an adiabatic parameter
variation. We investigate the influence of higher-order contributions beyond
nearest-neighbour coupling as well as self-coupling on the stability of a
$\mathrm{U}(3)$-phase generated from an optical tetrapod setup. Our results
indicate that, despite the mode nonorthogonality, the symmetry of dark states
protects the geometric evolution of light from distortion.

##  I Introduction

Abelian as well as non-Abelian gauge theories feature prominently in modern
theories of fundamental interactions. They also occur frequently in a variety
of geometric settings and are therefore central to much of modern mathematics.
Their formulation in terms of gauge fields is deeply connected to the notion
of a geometric phase. These phase factors arise naturally when considering the
adiabatic evolution of a state vector in Hilbert space Bohm ; Berry ; Wilzeck
. Besides this deeper insight into geometric and topological notions at an
experimentally feasible scale, non-Abelian phases (holonomies) are important
for holonomic and topological quantum computation, where they offer parametric
robustness to increase stability of the computational process HQC ; Pachos .
Consequently, there has been increased interest in the study and
implementation of artificial gauge fields and symmetry groups Cirac .
Successful implementations ranged from single artificial gauge fields in
photonic Segev , superconducting Abdumal , and atomic systems Spielman  to
experimental simulation of lattice gauge theories Zoller .

In recent years, a novel approach was put forward to realise geometric phases
in terms of integrated photonic structures such as laser-written waveguides in
fused silica. This have been proven to be a versatile tool box that combines
the proven capabilities of modern optics with a high degree of interferometric
stability Szameit . For instance, the emergence of a non-Abelian Berry phase
was observed when injecting coherent states of light into topologically guided
modes Chamon ; ChamonExp , whereas other implementations made use of a tripod
arrangement of evanescently coupled waveguides Teuber . In Ref. Ralph , a
proposal based on a photonic bus mode was studied with the view to implement a
controlled-NOT gate on single photons. An extension to an $N$-pod scheme with
arbitrary photon-number states injected was proposed in Ref. Pinske2020 .
However, all the above-mentioned proposals rely on power-orthogonality of
transverse modes. In general, the deviations from mode orthogonality can
significantly distort the dynamics of light in coupled-mode systems.

In our present work, we overcome the limitations of the weak-coupling regime
by investigating the properties of geometric phases within the framework of
nonorthogonal coupled-mode theory (NOCMT). In this setting, the
nonorthogonality of transverse modes comes to light when decreasing the
separation between adjacent waveguides. In this regime, the energy splitting
between an adiabatic subspace and excited states may be large, thus favouring
the generation of an adiabatic quantum holonomy. We present an optical setup
for the generation of $\mathrm{U}(2)$-valued and $\mathrm{U}(3)$-valued
geometric phases arising from tripod and tetrapod arrangements of waveguides,
respectively, into which coherent light is being injected. After an expansion
of the system in terms of normal modes, an analytical computation of the
geometric phase reveals that the symmetry of dark states protects the state of
light from any distortion.

The article is structured as follows. In Sec. II, we review the nonorthogonal
coupled-mode theory (NOCMT) for a network of weak-index contrast waveguides.
Sec. III is dedicated to the study of a tripod arrangement of waveguides in an
experimentally feasible setting. We compute its geometric phase from a normal
mode expansion of the dark states and compare it to a numerical propagation of
the longitudinal fields, thus quantifying the diabatic error. In Sec. IV, the
influence of higher-order coupling and self-coupling is illuminated using the
example of a $\mathrm{U}(3)$-valued geometric phase generated from an optical
tetrapod. Finally, Sec. V contains a summary as well as some concluding
remarks.

##  II Nonorthogonal coupled-mode theory

Since the inception of the concept of coupled modes in electromagnetic systems
Pierce ; Schelkunoff , coupled-mode theory has become a well-established tool
in the description of parametric nonlinear devices, waveguide structures,
optical fibre networks and various other optoelectronic structures Streifer ;
Chuang ; Peall . The starting point of our investigation is the propagation of
coherent light in a linear, lossless, and isotropic medium patterned with an
inhomogeneous lattice of $N$ waveguides with refractive-index profile
$n(\bm{r})=n_{0}+\sum_{j}\Delta n_{j}(\bm{r})$, where $n_{0}$ is the
refractive index of the host material and $\Delta n_{j}$ is the refractive-
index contrast of the $j$th waveguide. The dynamics of an electromagnetic wave
is determined by the inhomogeneous wave equation for the electric field
$\bm{E}(\bm{r},t)$

$\bm{\nabla}^{2}\bm{E}-\frac{n^{2}}{c^{2}}\partial_{t}^{2}\bm{E}=\bm{E}\cdot\bm{\nabla}\mathrm{ln}\,n^{2},$
(1)

where $c$ denotes speed of light in vacuum. Under the assumption of a weak
refractive-index contrast, that is
$\bm{\nabla}\mathrm{ln}\,n^{2}\approx\bm{0}$, Eq. (1) turns into a homogeneous
wave equation. It follows that the Helmholtz equation for the paraxial
propagation (slowly varying envelope) of an electromagnetic plane wave
$\bm{E}(\bm{r},t)=\bm{E}(\bm{r})\mathrm{e}^{\mathrm{i}(\omega t-\beta z)}$
reads

$\left(\bm{\nabla}^{2}_{\bot}+\lambdabar^{-2}n^{2}(\bm{r})-2\mathrm{i}\beta\partial_{z}-\beta^{2}\right)\bm{E}(\bm{r})=\bm{0}$
(2)

where $\lambdabar=c/\omega$, $z$ denotes the propagation coordinate, $\beta$
is the propagation constant, and $\bm{\nabla}^{2}_{\bot}$ denotes the
Laplacian with respect to the transverse coordinates $\bm{r}_{\bot}$. In a
weak index-contrast lattice we have $\Delta n_{j}\ll n_{0}$, hence the vector
character of the Helmholtz equation can be neglected by factoring out a
constant unit vector $\bm{\nu}$, because polarization effects at the interface
of a waveguide become unimportant, and scalar wave theory can be applied Gloge
. For the treatment of coupled-mode theory that includes high-index contrast
waveguides we refer the reader to Ref. Bassi .

In the light of the assumption of weak index contrast, we expand the wave
packet in terms of tight-binding modes as

$\bm{E}(\bm{r})=\sum_{j=1}^{N}a_{j}(z)w_{j}(\bm{r})\bm{\nu},$ (3)

where $w_{j}$ is the transverse field mode that is localised around the
waveguide $\Delta n_{j}$, and $a_{j}$ is the corresponding longitudinal field
amplitude determining the dynamics of the light field. The former constitute a
quasi-complete (complete for divergence free fields) basis in which fields
with compact support over the transverse plane can be expanded. Here we assume
that the shape of the transverse field $w_{j}$ remains constant throughout the
propagation. However, because we consider waveguides that vary their position
in the transverse plane along the propagation direction, the fields $w_{j}$
depend still on $z$. On the other hand, in any typical setup the waveguide
position changes only very slowly over the propagation length (paraxial
approximation), so that we can safely assume $w_{j}(\bm{r})\approx
w_{j}(\bm{r}_{\bot})\mathrm{e}^{\mathrm{i}\beta_{j}z}$ with an individual
propagation constant $\beta_{j}$ of the $j$th waveguide. The latter ones can
be neglected throughout this work by noting that they act merely as offsets
when considering an identical writing process for each waveguide.

We assume the transverse fields to satisfy their own Helmholtz equation in the
$(x,y)$-plane in the absent of the other waveguides, i.e.

$\left(\bm{\nabla}_{\bot}^{2}+\lambdabar^{-2}\left[n_{0}+\Delta
n_{j}(\bm{r})\right]^{2}-\beta^{2}\right)w_{j}(\bm{r}_{\bot})=0.$ (4)

This puts the expansion (3) in the light of a weighted eigenmode expansion. We
can simplify the above expression by noting that for a weak index contrast we
have $\left[n_{0}+\Delta n_{j}\right]^{2}\approx n_{0}^{2}+2n_{0}\Delta n_{j}$
and $\beta^{2}\approx n_{0}^{2}/\lambdabar^{2}$. It thus follows that Eq. (4)
takes the form

$\left(\frac{\lambdabar^{2}}{2n_{0}}\bm{\nabla}_{\bot}^{2}+\Delta
n_{j}(\bm{r})\right)w_{j}(\bm{r}_{\bot})=0.$ (5)

Inserting the ansatz (3) into Eq. (2), making use of Eq. (5) and employing the
same assumptions as before, we obtain

$\sum_{j=1}^{N}\Big{(}\mathrm{i}\lambdabar\partial_{z}a_{j}w_{j}-\sum_{\begin{subarray}{c}m=1\\\
m\neq j\end{subarray}}^{N}\Delta n_{m}a_{j}w_{j}\Big{)}=0.$ (6)

When contracting Eq. (6) with $w_{k}^{*}$, we observe that the transverse
fields overlap with the surrounding waveguides (sites) and their modes gives
rise to transverse interactions (see Fig. 1).

Figure 1:  Schematic front view of evanescently coupled waveguides (black
circles). Each site $j$ supports a fundamental transverse mode $w_{j}$ (shaded
areas). The modes extend towards the nearest neighbouring sites.

Here, we distinguish between two different contributions that occur in Eq. (6)
after contraction, namely

$\begin{split}\kappa_{jk}&=\frac{1}{\lambdabar}\sum_{\begin{subarray}{c}m=1\\\
m\neq j\end{subarray}}^{N}\int_{S_{\infty}}\Delta
n_{m}w_{k}^{*}w_{j}\mathrm{d}^{2}\bm{r}_{\bot},\\\
\sigma_{jk}&=\int_{S_{\infty}}w_{k}^{*}w_{j}\mathrm{d}^{2}\bm{r}_{\bot},\\\
\end{split}$ (7)

where we integrate over the entire $(x,y)$-plane $S_{\infty}$. In Eq. (7),
$\kappa_{jk}\approx\frac{1}{\lambdabar}\int_{S_{\infty}}\Delta
n_{k}w_{k}^{*}w_{j}\mathrm{d}^{2}\bm{r}_{\bot}$ is the evanescent coupling
between waveguides $j$ and $k$, while $\sigma_{jk}$ describes the overlap of
their respective transverse modes. In particular, there is also a self-
coupling $\nu_{j}=\kappa_{jj}$ due to the presence of other waveguides around
the $j$th mode. The latter is usually the smallest contribution to the
dynamical propagation (6). This can be seen from Fig. 2 where these quantities
were computed for two adjacent (identically written) cylindrical waveguides as
a function of their separation. We considered a scenario in which each
waveguide supports only its first transverse mode $w_{j}(\bm{r}_{\bot})$,
given in terms of Bessel functions FOP . These assumptions are not too
restrictive, as monomode operation and mode synchronism are often required as
operating conditions. Furthermore, due to mode overlaps, the transverse modes
form a nonorthogonal basis for describing the light propagation. The
nonorthogonality can be neglected if the distance between adjacent waveguides
is large. Then one works in the regime of OCMT, i.e. $\sigma_{jk}\approx 0$
for $j\neq k$ Szameit ; Haus ; Haus2 .

Figure 2:  Coupling $\kappa$ (in $\mathrm{mm}^{-1}$), overlap $\sigma$, and
self-coupling $\nu$ (in $\mathrm{mm}^{-1}$) between two adjacent waveguides as
a function of their separation $\delta$ (measured from the center). The
results are shown for cylindrical fibres of radius $R=4.8\,\mu\mathrm{m}$ with
bulk index $n_{0}=1.452$ and weak contrast $n_{\mathrm{I}}=6.53\cdot 10^{-4}$
for the inside of each waveguide. The wavelength of the injected light beam
considered, is $\lambda=633\,\mathrm{nm}$. The interpolating functions were
obtained by evaluating Eq. (7) for several values of $\delta$.

With the definitions from Eqs. (7) at hand, the coupled-mode equations become

$\mathrm{i}\sum_{j=1}^{N}\sigma_{jk}\partial_{z}a_{j}=\sum_{j=1}^{N}\kappa_{jk}a_{j}.$
(8)

This form of the paraxial Helmholtz equation resembles a discrete Schrödinger
equation. Equation (8) can be reformulated as a first-order matrix
differential equation in $z$, with $(\bm{\mathrm{\Sigma}})_{jk}=\sigma_{jk}$
and $(\bm{\mathrm{K}})_{jk}=\kappa_{jk}$ being the power and coupling matrix,
respectively. Because $\bm{\mathrm{\Sigma}}$ and $\bm{\mathrm{K}}$ do not
commute in general, the generator $\bm{\mathrm{\Sigma}}^{-1}\bm{\mathrm{K}}$
of the dynamics is not necessarily Hermitian, and
$\bm{\mathrm{a}}^{\dagger}\bm{\mathrm{a}}$ would not be conserved throughout
the propagation. Particularly, we observe that these generators of the NOCMT
lie within the class of real-valued non-symmetric matrices. Nevertheless, the
modified intensity distribution
$\bm{\mathrm{a}}^{\dagger}\bm{\mathrm{\Sigma}}\bm{\mathrm{a}}$ (POVM
measurement) remains constant as long as $\bm{\mathrm{K}}$ is Hermitian, that
is, dissipative effects are negligible as assumed.

The non-Hermitian nature of Eq. (8) can even be lifted completely by
transforming to a set of longitudinal normal modes $\\{b_{j}(z)\\}_{j=1}^{N}$
such that the overall propagation is unitary Haus . Even though we did not
start our derivation from power-orthogonal modes, conservation of energy
demands that such normal modes always exist. Obviously, in contrast to the
waveguide mode $a_{j}$, the corresponding normal mode $b_{j}$ will contain
contributions from adjacent waveguides. To be precise, their relation is given
by the conserved quantity
$\bm{\mathrm{b}}^{\dagger}\bm{\mathrm{b}}=\bm{\mathrm{a}}^{\dagger}\bm{\mathrm{\Sigma}}\bm{\mathrm{a}}$.
Moreover, because the power matrix $\bm{\mathrm{\Sigma}}$ is positive
definite, there exists a (non-unique) matrix $\bm{\mathrm{Q}}$ such that
$\bm{\mathrm{\Sigma}}=\bm{\mathrm{Q}}^{\dagger}\bm{\mathrm{Q}}$. It follows
that the normal modes are given by the transformation
$\bm{\mathrm{b}}=\bm{\mathrm{Q}}\bm{\mathrm{a}}$. With these preparations, Eq.
(8) can be rewritten as

$\mathrm{i}\partial_{z}\bm{\mathrm{b}}=\bm{\mathrm{H}}\bm{\mathrm{b}},$ (9)

where $\bm{\mathrm{H}}$ is given by Haus

$\bm{\mathrm{H}}=\left(\bm{\mathrm{Q}}^{-1}\right)^{\dagger}\bm{\mathrm{K}}\bm{\mathrm{Q}}^{-1},$
(10)

which is Hermitian as long as $\bm{\mathrm{K}}=\bm{\mathrm{K}}^{\dagger}$.
From similarity with Eq. (10) it follows that
$\bm{\mathrm{\Sigma}}^{-1}\bm{\mathrm{K}}$ is always diagonalisable with real
spectrum.

###  II.1 Adiabatic Propagation

While a general state vector $\bm{\mathrm{\Psi}}(z)$ evolves according to Eq.
(9), under the adiabatic assumption any initial preparation
$\bm{\mathrm{\Psi}}(z_{0})$ in a dark subspace (zero-eigenvalue eigenspace)
will be mapped to a state
$\bm{\mathrm{\Psi}}(z_{\mathrm{f}})=\bm{\mathrm{U}}(z_{0},z_{\mathrm{f}})\bm{\mathrm{\Psi}}(z_{0})$
that is also in the dark subspace. Then, we can expand the state
$\bm{\mathrm{\Psi}}$ at every instance $z$ in terms of the dark states, i.e.
$\bm{\mathrm{\Psi}}(z)=\sum_{a}U_{ab}(z)\bm{\mathrm{d}}_{a}(z)$, with initial
condition $\bm{\mathrm{\Psi}}(z_{0})=\bm{\mathrm{d}}_{a}(z_{0})$. Inserting
this ansatz into Eq. (9) and following Ref. Wilzeck , we obtain
$\left(\bm{\mathrm{U}}^{-1}\partial_{z}\bm{\mathrm{U}}\right)_{ba}=(\bm{\mathrm{A}}_{z})_{ab}$
where we used $\bm{\mathrm{H}}\bm{\mathrm{\Psi}}=\bm{\mathrm{0}}$ for states
$\bm{\mathrm{\Psi}}$ in the dark subspace. Here we defined the adiabatic
connection Chern  as
$(\bm{\mathrm{A}}_{z})_{ab}=\bm{\mathrm{d}}_{b}^{\dagger}\partial_{z}\bm{\mathrm{d}}_{a}$.
A formal solution for $\bm{\mathrm{U}}$ is then given in terms of the matrix
exponential

$\bm{\mathrm{U}}(z_{0},z_{\mathrm{f}})=\mathscr{T}\exp\int_{z_{0}}^{z_{\mathrm{f}}}\bm{\mathrm{A}}_{z}\,\mathrm{d}z,$
(11)

where $\mathscr{T}$ denotes $z$-ordering.

For a collection of coupled waveguides, the $z$-dependence of the system is
functionally connected to the distances between the waveguides
$\\{\delta_{\mu}(z)\\}_{\mu}$ that form local coordinates of an abstract
parameter space $\mathscr{M}$ (control manifold). Then, the propagation along
the $z$-direction can equivalently be viewed as a parameter variation along a
curve $\mathcal{C}:[z_{0},z_{\mathrm{f}}]\rightarrow\mathscr{M}$. In the
following sections we are not interested in arbitrary paths, but only those
that (approximately) form loops in $\mathscr{M}$, i.e.
$\mathcal{C}(z_{0})=\mathcal{C}(z_{\mathrm{f}})$. Population transfer between
waveguides associated with such a loop appears to be robust against stochastic
fluctuations of the control parameter, as the transformation
$\bm{\mathrm{U}}(z_{0},z_{\mathrm{f}})$ only depends on the area enclosed by
the loop $\mathcal{C}$ Pachos2000 .

##  III Tripod arrangement of waveguides

In order to examine the properties of geometric phases in relation to
nonorthogonal modes more closely, we shall focus on a benchmark example that
can be implemented using the current state of technology. A schematic
representation of the waveguide network under investigation is shown in Fig.
3. There, the outer waveguides $\mathrm{L}$, $\mathrm{R}$, and $\mathrm{U}$
interact only implicitly via the central waveguide $\mathrm{C}$ (nearest-
neighbour coupling). To be more precise, we consider a situation in which the
transverse fields of the outer waveguides have couplings and overlaps that can
be safely disregarded. Then, the coupling matrix of the system is

$\bm{\mathrm{K}}=\sum_{j\neq\mathrm{C}}\kappa_{j\mathrm{C}}\left(\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{j}\right)+\sum_{j}\nu_{j}\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{j},$

where $j$ enumerates the waveguides
$\\{\mathrm{L},\mathrm{R},\mathrm{U},\mathrm{C}\\}$. Here,
$\\{\bm{\mathrm{w}}_{j}\\}_{j}$ denotes the standard basis in $\mathbb{C}^{4}$
that forms an orthonormal representation of the transverse eigenmodes of the
$j$th waveguide. For the moment, we neglect self-coupling on the matrix
diagonal, as it appears to be the smallest contribution to the overall
propagation. The corresponding power matrix of the system then reads

$\bm{\mathrm{\Sigma}}=\sum_{j\neq\mathrm{C}}\sigma_{j\mathrm{C}}\left(\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{j}\right)+\sum_{j}\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{j},$
(12)

where we assumed normalised transverse modes such that $\sigma_{jj}=1$. The
setup depicted in Fig. 3 is known as a tripod system. The tripod scheme,
arising in different physical settings, is often used as a starting point for
the generation of $\mathrm{U}(2)$-valued phases, see, e.g., Refs. Bergmann ;
Recati ; Duan .

Figure 3:  Front view of the tripod arrangement, arrows indicate interactions
between the outer waveguides and the central waveguide. Within NOCMT, this
includes evanescent coupling as well as the overlap of transverse modes. In
OCMT, this overlap is neglected. The interaction strength depends on the
separation $\delta_{i\mathrm{C}}$ between the corresponding waveguides.

One possible transformation between waveguide and normal modes can be obtained
by a Cholesky factorisation of the matrix (12), that reads

$\bm{\mathrm{Q}}=\sum_{j\neq\mathrm{C}}\left(\sigma_{j\mathrm{C}}\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{j}\right)+\sqrt{1-\bm{\sigma}^{2}}\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}},$
(13)

where we introduced the vector
$\bm{\sigma}=(\sigma_{\mathrm{L}\mathrm{C}},\sigma_{\mathrm{R}\mathrm{C}},\sigma_{\mathrm{U}\mathrm{C}})$.
From the matrix (13) we can compute the Hermitian generator for the evolution
of normal modes according to Eq. (10), and we find

$\bm{\mathrm{H}}=\sum_{j\neq\mathrm{C}}s\kappa_{j\mathrm{C}}\left(\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{j}\right)-2s^{2}\bm{\sigma}\cdot\bm{\kappa}\,\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}},$

with
$\bm{\kappa}=(\kappa_{\mathrm{L}\mathrm{C}},\kappa_{\mathrm{R}\mathrm{C}},\kappa_{\mathrm{U}\mathrm{C}})$
and $s=1/\sqrt{1-\bm{\sigma}^{2}}$. The system’s dark states

$\begin{split}\bm{\mathrm{d}}_{1}&=\sin\theta\bm{\mathrm{w}}_{\mathrm{L}}-\cos\theta\bm{\mathrm{w}}_{\mathrm{R}},\\\
\bm{\mathrm{d}}_{2}&=\cos\theta\sin\varphi\bm{\mathrm{w}}_{\mathrm{L}}+\sin\theta\sin\varphi\bm{\mathrm{w}}_{\mathrm{R}}-\cos\varphi\bm{\mathrm{w}}_{\mathrm{U}},\\\
\end{split}$ (14)

satisfy $\bm{\mathrm{H}}\bm{\mathrm{d}}_{a}=\bm{\mathrm{0}}$, where we defined
$\theta=\arctan(\kappa_{\mathrm{R}\mathrm{C}}/\kappa_{\mathrm{L}\mathrm{C}})$
and
$\varphi=\arctan(\kappa_{\mathrm{U}\mathrm{C}}/\sqrt{\kappa_{\mathrm{L}\mathrm{C}}^{2}+\kappa_{\mathrm{R}\mathrm{C}}^{2}})$.
Quite remarkably, as long as the second-order (next-nearest-neighbour)
contributions to the coupling SzameitSOC  as well as self-coupling are
negligible, the introduction of mode overlaps does not break the degeneracy of
the system, and still gives rise to a two-fold degenerate dark subspace
similar to its counterpart in OCMT Teuber . In particular, the dark states of
$\bm{\mathrm{H}}$ have the same structure as the ones found in OCMT but with
their coefficients now giving the population of the associated normal mode
that differ from the individual modes of each waveguide.

With the dark states known, computing the connection $\bm{\mathrm{A}}_{z}$
becomes a straightforward task and, thus, the matrix exponential (11) can be
evaluated analytically to

$\bm{\mathrm{U}}(z_{0},z_{\mathrm{f}})=\begin{pmatrix}\cos\phi&\sin\phi\\\
-\sin\phi&\cos\phi\\\ \end{pmatrix},$ (15)

written in the basis $\\{\bm{\mathrm{d}}_{a}(z)\\}_{a}$. Here

$\phi(z_{0},z_{\mathrm{f}})=\int_{z_{0}}^{z_{\mathrm{f}}}\frac{\kappa_{\mathrm{U}\mathrm{C}}\left(\kappa_{\mathrm{L}\mathrm{C}}\partial_{z}\kappa_{\mathrm{R}\mathrm{C}}-\kappa_{\mathrm{R}\mathrm{C}}\partial_{z}\kappa_{\mathrm{L}\mathrm{C}}\right)}{\sqrt{\kappa_{\mathrm{L}\mathrm{C}}^{2}+\kappa_{\mathrm{R}\mathrm{C}}^{2}+\kappa_{\mathrm{U}\mathrm{C}}^{2}}\left(\kappa_{\mathrm{L}\mathrm{C}}^{2}+\kappa_{\mathrm{R}\mathrm{C}}^{2}\right)}\mathrm{d}z$
(16)

is a phase factor depending only on the geometry of the associated parameter
variation $\mathcal{C}$. The holonomy matrix (15) is the same unitary
transformation as for the case of orthogonal transverse modes Teuber . This
implies that the results of the NOCMT and the OCMT predict the same output
state as long as the evolution takes place solely in the dark subspace. This
means that the geometric phase generated from the tripod is robust against the
nonorthogonal nature of transverse modes. This is a remarkable result as, in
general, the deviations from mode-orthogonality can significantly distort the
dynamics of light in a coupled-mode system.

###  III.1 Propagation with Gaussian-shaped geometries

Here, we present an experimentally feasible parameter loop $\mathcal{C}$ that
can be realised by current laser-writing techniques Szameit . Suppose one
implements the geometry

$\begin{split}\delta_{\mathrm{L}\mathrm{C}}(z)&=\Delta-\Omega\,\exp\left(-\frac{\left(z-\overline{z}-\tau\right)^{2}}{T^{2}}\right),\\\
\delta_{\mathrm{R}\mathrm{C}}(z)&=\Delta-\Omega\,\exp\left(-\frac{\left(z-\overline{z}+\tau\right)^{2}}{T^{2}}\right),\\\
\delta_{\mathrm{U}\mathrm{C}}(z)&=\Delta_{\mathrm{U}},\\\ \end{split}$ (17)

in which the central waveguide is assumed to remain at the origin of the
transverse plane, while the separation to the $\mathrm{L}$ and $\mathrm{R}$
waveguide changes according to a Gaussian function. In Eq. (17),
$\overline{z}$ denotes half the total propagation length, $\Delta-\Omega$ is
the minimal separation between waveguide $\mathrm{L}$ ($\mathrm{R}$) and the
central waveguide $\mathrm{C}$, with waveguide $\mathrm{U}$ having constant
distance $\Delta_{\mathrm{U}}$. Furthermore, $T$ is the width parameter, and
$\tau$ is the separation of the two Gaussian peaks from the center at
$\overline{z}$. The relevant coupling strengths and overlaps of the respective
waveguides can be determined from their mutual separation, shown in Fig. 2.
Given the variation (17), the respective coupling constants and overlaps are
shown in Fig. 4.

Figure 4:  Couplings $\kappa_{i\mathrm{C}}$ (in $\mathrm{mm}^{-1}$) and
overlaps $\sigma_{i\mathrm{C}}$ between the $i$th outer waveguide and the
central waveguide as a function of the propagation distance $z$ for the case
of the $z$-dependent waveguide position (17). The chosen parameters are
$\Delta=40\,\mu\mathrm{m}$, $\Delta_{\mathrm{U}}=14\,\mu\mathrm{m}$,
$\Omega=26\,\mu\mathrm{m}$, $T=8\,\mathrm{cm}$, and $\tau=2\,\mathrm{cm}$.

A parameter loop starting at the initial point
$\bm{\delta}_{0}=(\Delta,\Delta,\Delta_{\mathrm{U}})$ will induce a nontrivial
mixing of states

$\bm{\mathrm{d}}_{1}(\bm{\delta}_{0})\approx\bm{\mathrm{w}}_{\mathrm{L}},\qquad\bm{\mathrm{d}}_{2}(\bm{\delta}_{0})\approx\bm{\mathrm{w}}_{\mathrm{R}}.$

After integration of the phase factor (16) and subsequent evaluation of the
transformation (15), this mixing can be given explicitly. The adiabatic limit
is applicable if the loop $\mathcal{C}$ given by Eq. (17) is traversed slowly
enough compared to the energy splitting between the dark subspace and the
excited bright subspaces. We found that the separation between dark and bright
states is slightly shifted when compared to the eigenvalue splitting known
from OCMT, due to the contribution from mode overlaps. In order to evaluate to
which extent the adiabatic approximation is justified in NOCMT, we compute the
(expected) gate fidelity $0\leq\overline{F}\leq 1$ in terms of input state
fidelities

$F_{\bm{\Psi}(z_{0})}=|\tilde{\bm{\Psi}}(z_{\mathrm{f}})^{\dagger}\bm{\mathrm{Q}}^{-1}(z_{\mathrm{f}})\bm{\mathrm{U}}(z_{\mathrm{f}})\bm{\Psi}(z_{0})|^{2},$
(18)

with $\bm{\mathrm{U}}(z_{\mathrm{f}})\bm{\Psi}(z_{0})$ being the ideal output
state predicted by the holonomic theory (15) and
$\tilde{\bm{\Psi}}(z_{\mathrm{f}})$ obtained from a numerical propagation
subject to the paraxial Helmholtz equation (8). For the scenario under study,
the gate fidelity of the dark state ensemble
$\\{\\{1/2,\bm{\mathrm{w}}_{\mathrm{L}}\\},\\{1/2,\bm{\mathrm{w}}_{\mathrm{R}}\\}\\}$
amounts to approximately $\overline{F}\approx 98.6\,$%.

We can conclude that robust adiabatic parallel transport of an initial wave
packet can be achieved within the range of experimentally feasible setups,
despite the influence of mode-nonorthogonality. In fact, the expected
distortion from close-coupling dynamics comes to light when one of the bright
states $\bm{\mathrm{w}}_{\mathrm{U}}$ or $\bm{\mathrm{w}}_{\mathrm{C}}$ is
excited. In Fig. 5 we compare the population transfer between the waveguides
for different initial input states. One can clearly see the robustness of the
dark states, while the excited states are exposed to more rapid population
transfer. In fact, when adding the bright states to the input ensemble, a gate
error emerges that is not of diabatic nature, but has its origin in the fact
that the parameter variation (17) does not form a perfect loop. This matching
error of input and output basis is rather small for the dark states but gets
amplified for the bright states due to rapid intensity oscillations, thus
leading to a substantial decrease in gate fidelity, $\overline{F}\approx
88.9\,$%. The relation between mode-nonorthogonality and rapid intensity
oscillations was experimentally verified in Ref. Weimann . In general, the
robustness of dark states can be attributed to their form given by Eq. (14)
which is independent of $\sigma_{i\mathrm{C}}$ and therefore identical to
their counterparts in OCMT Teuber . Hence, despite the overlaps, absence of
light in waveguide $\mathrm{C}$ can still be used as a measure of
adiabaticity.

(a) Intensity in site $\mathrm{L}$ for
$\bm{\mathrm{\Psi}}(z_{0})=\bm{\mathrm{w}}_{\mathrm{L}}$.(b) Intensity in site
$\mathrm{R}$ for $\bm{\mathrm{\Psi}}(z_{0})=\bm{\mathrm{w}}_{\mathrm{R}}$. (c)
Intensity in site $\mathrm{U}$ for
$\bm{\mathrm{\Psi}}(z_{0})=\bm{\mathrm{w}}_{\mathrm{U}}$.(d) Intensity in site
$\mathrm{C}$ for $\bm{\mathrm{\Psi}}(z_{0})=\bm{\mathrm{w}}_{\mathrm{C}}$.
Figure 5:  Simulated propagation of $\bm{\mathrm{\Psi}}(z)$ in term of
waveguide modes [subject to Eq. (8)] for the waveguide position (17) under
investigation. The four plots show the change of intensity throughout
propagation for different initial input states
$\\{\bm{\mathrm{w}}_{j}\\}_{j}$. It can be seen that the bright states
($\mathrm{U}$ and $\mathrm{C}$) are subject to rapid intensity fluctuations
due to deviations from mode orthogonality. The parameter values are
$\Delta=40\,\mu\mathrm{m}$, $\Delta_{\mathrm{U}}=14\,\mu\mathrm{m}$,
$\Omega=26\,\mu\mathrm{m}$, $T=8\,\mathrm{cm}$, and $\tau=2\,\mathrm{cm}$.

##  IV Higher-order coupling in the tetrapod arrangement

We now turn to the influence of self-coupling and higher-order contributions
in a waveguide network which had been neglected in our calculations thus far.
These effects potentially break the degeneracy structure of the system under
investigation. We found numerically that, in the tripod arrangement, these
contributions are completely negligible for the propagation lengths
investigated ($|z_{\mathrm{f}}-z_{0}|\leq 15\,\mathrm{cm}$). However, this
result depends crucially on the distance behaviour shown in Fig. 2, and does
not need to be valid beyond step-index waveguides.

In order to show that coupling beyond nearest-neighbour hopping as well as
self-coupling can become relevant when the distance between waveguides is
substantially decreased, we now turn to a tetrapod arrangement depicted in
Fig. 6. Placing more waveguides around the central site will inevitably result
in a system having no degeneracy at all. The coupling and power matrices of
the (ideal) tetrapod system without self-coupling and higher-order effects are
given by

$\bm{\mathrm{K}}=\sum_{j\neq\mathrm{C}}\kappa_{j\mathrm{C}}\left(\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{j}\right)$
(19)

and

$\bm{\mathrm{\Sigma}}=\sum_{j\neq\mathrm{C}}\sigma_{j\mathrm{C}}\left(\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{\mathrm{C}}+\bm{\mathrm{w}}_{\mathrm{C}}^{\dagger}\bm{\mathrm{w}}_{j}\right)+\sum_{j}\bm{\mathrm{w}}_{j}^{\dagger}\bm{\mathrm{w}}_{j},$
(20)

where $j$ runs through the alphabet
$\\{\mathrm{N},\mathrm{E},\mathrm{S},\mathrm{W},\mathrm{C}\\}$, and hence
$\\{\bm{\mathrm{w}}_{j}\\}_{j}$ might be thought of as the standard basis in
$\mathbb{C}^{5}$. For this system, the adiabatic propagation of normal modes
will generate a $\mathrm{U}(3)$-valued mixing within the three-fold degenerate
dark subspace.

Figure 6:  Front view of the tetrapod arrangement, arrows indicate interaction
between the respective waveguides. Within the NOCMT this includes evanescent
coupling as well as the overlap of transverse modes. In the case of OCMT the
overlap is neglected. The strength of interaction depends on the separation
$\delta_{i\mathrm{C}}$ between the corresponding waveguides. Black arrows were
chosen for nearest-neighbour coupling (first order). Grey solid lines are
reserved for next-nearest neighbour coupling (second order), while grey dotted
lines depict coupling between opposite placed waveguides (third order).

Higher-order coupling becomes relevant as soon as the propagation length of
the evolution is substantially increased. In this regime, we can assume
adiabaticity to apply. Then, the matrices in Eq. (19) and Eq. (20) are no
longer sparse, but have only nonvanishing entries. The diagonal entries of
$\bm{\mathrm{K}}$ now contain the self-couplings $\nu_{j}$. In the following,
we will study the influence of such contributions on a benchmark geometry

$\begin{split}\delta_{\mathrm{N}\mathrm{C}}(z)&=\Delta-\Omega\,\exp\left(-\frac{\left(z-\overline{z}-\tau\right)^{2}}{T^{2}}\right),\\\
\delta_{\mathrm{E}\mathrm{C}}(z)&=\Delta-\Omega\,\exp\left(-\frac{\left(z-\overline{z}\right)^{2}}{T^{2}}\right),\\\
\delta_{\mathrm{S}\mathrm{C}}(z)&=\Delta-\Omega\,\exp\left(-\frac{\left(z-\overline{z}+\tau\right)^{2}}{T^{2}}\right),\\\
\delta_{\mathrm{W}\mathrm{C}}(z)&=\Delta_{\mathrm{W}}.\\\ \end{split}$ (21)

In the chosen waveguide geometry, the evolution of coherent light
approximately starts and ends at the point
$\bm{\delta}_{0}=(\Delta,\Delta,\Delta,\Delta_{\mathrm{W}})$, where the
holonomy will generate a mixing of dark states

$\bm{\mathrm{d}}_{1}(\bm{\delta}_{0})\approx\bm{\mathrm{w}}_{\mathrm{N}},\quad\bm{\mathrm{d}}_{2}(\bm{\delta}_{0})\approx\bm{\mathrm{w}}_{\mathrm{E}},\quad\bm{\mathrm{d}}_{3}(\bm{\delta}_{0})\approx\bm{\mathrm{w}}_{\mathrm{S}}.$

Given the geometry (21), we can derive the $z$-dependent form of the relevant
parameters from Fig. 2. Knowing these parameters allows us, by simple
geometric considerations (cf. Fig. 6), to construct higher-order couplings and
overlaps as well. In Fig. 7, the coupling strengths of these different
contributions is shown for the waveguide $\mathrm{N}$. One observes clearly
how the higher-order couplings have a substantially smaller influence compared
to the nearest-neighbour coupling. Nevertheless, for longer propagation
lengths (e.g. $|z_{\mathrm{f}}-z_{0}|=40\,\mathrm{cm}$) these contributions
will lead to a deviation in the intensity distribution from the evolution
predicted by the propagation through the ideal tetrapod system.

Figure 7:  Couplings (in $\mathrm{mm}^{-1}$) and overlaps between the northern
($\mathrm{N}$) waveguide and $j$th waveguide
($j\in\\{\mathrm{C},\mathrm{E},\mathrm{S}\\}$) of the tetrapod arrangement as
a function of the propagation distance $z$ for the waveguide geometry (21).
Coupling to the central ($\mathrm{C}$) waveguide is of first order, coupling
to the eastern ($\mathrm{E}$) waveguide is of second order, and coupling to
the south ($\mathrm{S}$) waveguide is of third order. The relevant parameters
are given by $\Delta=38\,\mu\mathrm{m}$,
$\Delta_{\mathrm{W}}=12\,\mu\mathrm{m}$, $\Omega=26\,\mu\mathrm{m}$,
$T=15\,\mathrm{cm}$, and $\tau=5\,\mathrm{cm}$.

We also investigated the change of self-coupling $\nu_{j}(z)$ of each
respective waveguide. There we observed that the order of magnitude of the
self-coupling in the outer waveguides lies between those of the second-order
and third-order contributions. However, the self-coupling $\nu_{\mathrm{C}}$
is substantially larger, as the central waveguide mode is closely surrounded
by the other waveguides.

In Fig. 8, the propagation through the ideal tetrapod (solid line) is compared
to a propagation including higher-order couplings and self-coupling (dashed
line). The state was propagated numerically for an initial input
$\bm{\mathrm{w}}_{\mathrm{N}}$. Our results show that these additional effects
(mostly second-order coupling and self-coupling $\nu_{\mathrm{C}}$) lead to
increased population transfer between the waveguides. Due to symmetry-
breaking, this additional population transfer will not be robust against
parametric fluctuations and mode-nonorthogonalities as it is no longer
described by a parallel transport in the degenerate dark subspace of the
tetrapod.

Figure 8:  Intensity distribution for the tetrapod arrangement subject to the
geometry (21). The solid lines belong to the ideal tetrapod [Eq. (19) and Eq.
(20)], while the dashed lines describe propagation including higher-order
couplings as well as self-coupling. The results were computed via a numerical
simulation of Eq. (8) for an initial wave package
$\bm{\mathrm{w}}_{\mathrm{N}}$, with parameters $\Delta=38\,\mu\mathrm{m}$,
$\Delta_{\mathrm{W}}=12\,\mu\mathrm{m}$, $\Omega=26\,\mu\mathrm{m}$,
$T=15\,\mathrm{cm}$, and $\tau=5\,\mathrm{cm}$ for a network of
$40\,\mathrm{cm}$ total length.

##  V Conclusion

In this article, we have shown that the generation of geometric phases in
integrated photonic waveguide structures can be considered within the
framework of nonorthogonal coupled-mode theory. General arguments show that,
including the nonorthogonality of transverse modes leads to a set of tight-
binding equations governing an evolution in which the conventional intensity
distribution is not necessarily conserved. This issue was lifted by means of a
normal mode expansion. An analytical computation of the geometric phase and a
subsequent evaluation of the intensity distribution for the tripod system
showed a robustness of the adiabatic parallel transport against deviations
from mode-orthogonalities as long as adiabaticity holds and higher-order
couplings are negligible. A subsequent study of the tetrapod arrangement
showed that higher-order coupling as well as self-coupling can become relevant
for longer propagation lengths, eventually breaking the degeneracy of the
system, thus perturbing the generation of a $\mathrm{U}(3)$-valued geometric
phase. However, for appropriate propagation lengths our numerical simulations
of the population transfer between the waveguides showed that high fidelity
transformations can be generated when working in the close-coupling regime.

Our work paves the way for the study of adiabatic parallel transport by
photonic structures within the close-coupling regime. The symmetry of the dark
subspace leads to an inherent robustness of geometric phases towards the
distortion of light by mode-nonorthogonalities. Together with their parametric
stability, geometric phases can be considered a powerful tool for the
generation of stable population transfer of light even within the realm of the
nonorthogonal coupled-mode theory.

###### Acknowledgements.

Financial support by the Deutsche Forschungsgemeinschaft (DFG SCHE 612/6-1) is
gratefully acknowledged.

## References

  * (1) Y. Aharonov and D. Bohm, Significance of Electromagnetic Potentials in the Quantum Theory, Phys. Rev. 115, 485 (1959).
  * (2) M. V. Berry, Quantal Phase Factors Accompanying Adiabatic Changes, Proc. R. Soc. London, Ser. A 392, 45 (1984).
  * (3) F. Wilczek and A. Zee, Appearance of Gauge Structure in Simple Dynamical Systems, Phys. Rev. Lett. 52, 2111 (1984).
  * (4) P. Zanardi and M. Rasetti, Holonomic Quantum Computation, Phys. Lett. A 264, 94 (1999).
  * (5) J. K. Pachos, Introduction to Topological Quantum Computation (Cambridge, Cambridge University Press, 2012).
  * (6) M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero et al., Simulating lattice gauge theories within quantum technologies, Eur. Phys. J. D 74, 165 (2020).
  * (7) Y. Lumer, M. A. Bandres, M. Heinrich, L. J. Maczewsky, H. Herzig-Sheinfux, A. Szameit, and M. Segev, Light guiding by artificial gauge fields, Nat. Photon. 13, 339 (2019).
  * (8) A. A. Abdumalikov, Jr., J. M. Fink, K. Juliusson, M. Pechal, S. Berger, A. Wallraff, and S. Filipp, Experimental realization of non-Abelian non-adiabatic geometric gates, Nature (London) 496, 482 (2013).
  * (9) V. Galitski, G. Juzeli$\overline{\text{u}}$nas, and I. B. Spielman, Artificial gauge fields with ultracold atoms, Phys. Today 72, 38 (2019).
  * (10) E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller et al., Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature (London) 534, 516 (2016).
  * (11) A. Szameit and S. Nolte, Discrete optics in femtosecond-laser-written photonic structures, J. Phys. B: At. Mol. Opt. Phys. 43, 163001 (2010).
  * (12) T. Iadecola, T. Schuster, and C. Chamon, Non-Abelian Braiding of Light, Phys. Rev. Lett. 117, 073901 (2016).
  * (13) J. Noh, T. Schuster, T. Iadecola, S. Huang, M. Wang, K. P. Chen, C. Chamon, and M. C. Rechtsman, Nature Physics (London) 16, 989 (2020).
  * (14) M. Kremer, L. Teuber, A. Szameit, and S. Scheel, Optimal design strategy for non-Abelian geometric phases using Abelian gauge fields based on quantum metric, Phys. Rev. Res. 1, 033117 (2019).
  * (15) A. P. Hope, T. G. Nguyen, A. Mitchell, and A. D. Greentree, Adiabatic two-photon quantum gate operations using a long-range photonic bus, J. Phys. B: At. Mol. Opt. Phys. 48, 055503 (2015).
  * (16) J. Pinske, L. Teuber, and S. Scheel, Highly degenerate photonic waveguide structures for holonomic computation, Phys. Rev. A 101, 062314 (2020).
  * (17) J. R. Pierce, J. Appl. Phys 25, 179 (1954).
  * (18) S. A. Schelkunoff, Conversion of Maxwell’s Equations into Generalized Telegraphist’s Equations, Bell Syst. Tech. J. 34, 995 (1955).
  * (19) W. Streifer, M. Osinski, and A. Hardy, Reformulation of the coupled-mode theory of multiwaveguide systems, J. Lightwave Tech. 5, 1 (1987).
  * (20) S.-L. Chuang, A coupled mode formulation by reciprocity and a variational principle, J. Lightwave Tech. 5, 5 (1987).
  * (21) R. G. Peall and R. R. A. Syms, Comparison between strong coupling theory and experiment for three-arm directional couplers in Ti: LiNbO/sub 3, J. Lightwave Tech. 7, 540 (1989).
  * (22) D. Gloge, Weakly Guiding Fibers, Appl. Opt. 10, 2252, (1971).
  * (23) G. Bellanca, P. Orlandi, and P. Bassi, J. Opt. Soc. Am. A 35, 577 (2018).
  * (24) B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (New York, Wiley-VCH, 2007).
  * (25) H. A. Haus, W. P. Huang, S. Kawakami, and N. A. Whitaker, J. Lightwave Tech. 5, 13 (1987).
  * (26) H. A. Haus and W. Huang, Coupled-mode theory, Proc. IEEE, 19, 10 (1991).
  * (27) R. Bott and S. S. Chern, Hermitian vector bundles and the equidistribution of the zeros of their holomorphic sections, Acta Math. 114, 71 (1968).
  * (28) J. K. Pachos and P. Zanardi, Quantum Holonomies for Quantum Computing, Int. J. Mod. Phys. B 15, 1257 (2001).
  * (29) R. G. Unanyan, B. W. Shore, and K. Bergmann, Laser-driven population transfer in four-level atoms: Consequences of non-Abelian geometrical adiabatic phase factors, Phys. Rev. A 59, 2910 (1999).
  * (30) A. Recati, T. Calarco, P. Zanardi, J. I. Cirac, and P. Zoller, Holonomic quantum computation with neutral atoms, Phys. Rev. A 66, 032309 (2002).
  * (31) L. M. Duan, J. Cirac, and P. Zoller, Geometric Manipulation of Trapped Ions for Quantum Computation, Science 292, 1695 (2001).
  * (32) R. Keil, B. Pressl, R. Heilmann, M. Gräfe, G. Weihs, and A. Szameit, Direct measurement of second-order coupling in a waveguide lattice, Appl. Phys. Lett. 107, 241104 (2015).
  * (33) S. Weimann, PhD Thesis, Complex lattice modes in waveguide networks and photonic solids (University of Rostock, Germany, 2018).