# Biomimetic non-ergodic aging by dynamic-to-covalent transitions in physical
hydrogels

Samya Sen  Department of Materials Science & Engineering, Stanford University,
Stanford, CA 94305, USA  Anthony C. Yu  Department of Materials Science &
Engineering, Stanford University, Stanford, CA 94305, USA  Changxin Dong
Department of Materials Science & Engineering, Stanford University, Stanford,
CA 94305, USA  Andrea I. D’Aquino  Department of Materials Science &
Engineering, Stanford University, Stanford, CA 94305, USA  Eric A. Appel
Corresponding author<EMAIL_ADDRESS>Department of Materials Science &
Engineering, Stanford University, Stanford, CA 94305, USA  Department of
Bioengineering, Stanford University, Stanford, CA 94305, USA  Stanford ChEM-H,
Stanford University, Stanford, CA 94305, USA  Institute for Immunity,
Transplantation and Infection, Stanford University School of Medicine,
Stanford, CA 94305, USA  Department of Pediatrics - Endocrinology, Stanford
University School of Medicine, Stanford, CA 94305, USA  Woods Institute for
the Environment, Stanford University, Stanford CA 94305, USA

###### Abstract

Hydrogels are soft materials engineered to suit a multitude of applications
that exploit their tunable mechanochemical properties. Dynamic hydrogels
employing noncovalent, physically crosslinked networks dominated by either
enthalpic or entropic interactions enable unique rheological and stimuli-
responsive characteristics. In contrast to enthalpy-driven interactions that
soften with increasing temperature, entropic interactions result in largely
temperature-independent mechanical properties. By engineering interfacial
polymer-particle interactions, we can induce a dynamic-to-covalent transition
in entropic hydrogels that leads to biomimetic non-ergodic aging in the
microstructure without altering the network mesh size. This transition is
tuned by varying temperature and formulation conditions such as $p$H, which
allows for multivalent tunability in properties. These hydrogels can thus be
designed to exhibit either temperature-independent metastable dynamic
crosslinking or time-dependent stiffening based on formulation and storage
conditions, all while maintaining strucutural features critical for
controlling mass transport, akin to many biological tissues. Such robust
materials with versatile and adaptable properties can be utilized in
applications such as wildfire suppression, surgical adhesives, and depot-
forming injectable drug delivery systems.

##  1 Introduction

Physically crosslinked soft materials exhibit many unique functional and
stimuli-responsive physical and mechanical properties [1, 2, 3, 4]. Over the
years, soft materials employing either enthalpy- or entropy-driven dynamic
noncovalent crosslinking interactions have been designed to have unique
mechanochemical properties tunable to a fine degree based on target behaviors
[5]. These types of dynamic network materials, including many examples of
physical hydrogels, have been employed in a variety of applications that
exploit their robust and versatile properties, including biomaterials [6, 7,
8, 9, 10], drug delivery [11, 12, 13, 2, 3], self-healing materials [14, 15,
16, 17, 18], wearable electronics [16, 19, 20, 18, 21], 3D printing [22, 23,
24], and sprayable coatings for wildland fire suppression [25, 26].

The dynamic physical hydrogel literature is rife with materials that leverage
enthalpy-dominated crosslinking interactions, which are inherently
thermoresponsive _via_ temperature-dependent viscoelastic properties. Such
materials generally soften at elevated temperatures owing to the highly
specific enthalpy-driven binding reactions including ionic, host-guest,
hydrogen bonding, and hydrophobic interactions [2, 6, 15, 16, 19, 20, 17, 3,
27]. More recently, entropy-driven hydrogels with dynamic networks based on
multivalent polymer-particle interactions utilizing carefully tuned
interfacial chemistry have been formulated [25, 5, 28, 8]. Such networks built
with entropy-dominated crosslinking interactions can exhibit temperature-
independent viscoelasticity, and this unique property has been leveraged in
diverse applications from fire suppression [25] to drug delivery [29] where
consistent mechanical properties despite temperature variations are necessary.

In this work, we evaluate an example entropy-driven hydrogel system where
mixing of solutions of cellulose derivatives (HEC and MC) with colloidal
silica particles (CSP) generate dynamic HEC-MC/CSP hydrogels with facile and
scalable formulation (Fig. 1a) This material platform exhibits viscoelastic
rheological properties largely independent of temperature, attributed to the
dominance of the entropic contribution to free energy of the polymer-particle
interactions. We report a curious effect in this system where when aged, the
material gets stiffer over time (i.e., the modulus increases with aging time)
without alteration of the hydrogel microstructure. This aging behavior is
observed consistently over the course of weeks, and is dependent on storage
temperature and formulation pH, suggesting activated processes drive a
transition from dynamic to covalent crosslinking interactions in the system.
Such behaviors are fundamentally akin to aging observed in the extracellular
matrix found in many biological tissues and enable advantageous multifaceted
tunability of hydrogel mechanical properties by exploiting the differential
temperature and $p$H dependence. The biomimetic dynamic-to-covalent
transitions in these synthetic physical hydrogels enable the development of
materials with structure defined by the self-assembly of the gel components
that then evolve into covalent interactions in defined ways, making them
highly distinguishing in numerous application spaces.

Figure 1: Schematic representation of the chemical aging and proliferation of
crosslinks in the hydrogels with age. (a) HEC and MC polymer chains when mixed
with CSP dispersion form a hydrogel, where the CSPs function as transient
crosslink junctions between polymer chains adsorbed onto the surface. The plot
shows the ITC data from Yu _et al._ [25] for calculating the free energy of
binding $\Delta G^{0}$, with the enthalpic, $\Delta H^{0}$, and entropic,
$T\Delta S^{0}$, contributions, for the entropy-dominated hydrogels. (b)
Progression of interactions between polymer chains and particle surface
studied in this work, with dynamic, entropy-driven reversible binding at short
times after mixing, followed by the formation of irreversible enthalpy-driven
covalent bonds between the CSP and the adsorbed HEC-MC chains due to aging.

##  2 Experimental

###  2.1 Materials

The hydrogels studied in this work are prepared following a protocol outlined
in the literature [25]. First, polymer solutions were prepared by dissolving
hydroxyethylcellulose or HEC ($M_{\rm v}\simeq 1,300$ kDa, Sigma) and
methylcellulose or MC ($M_{\rm v}\simeq 90$ kDa, Sigma) in water (or citrate
buffer for $p$H 4 gels) at a concentration of 3% w/w, with stirring and mild
heating. Colloidal silica particles or CSP dispersions were made at 15% w/w by
diluting Ludox TMA 34% w/w dispersions (acid stable, $D_{\rm H}\simeq 22$ nm)
in water (or citrate buffer for $p$H 4 gels). Finally, equal weights of HEC-MC
solutions and CSP dispersions were mixed well by vortexing to obtain 1.5-7.5%
w/w HEC-MC/CSP gels. The properties of the gels are independent of the scale
of mixing as has already been shown in earlier studies [25, 26].

###  2.2 Rheometry

Rheological characterization of the gels was performed on a separated motor-
transducer rheometer (ARES-G2, TA Instruments) at Stanford Soft and Hybrid
Materials Facility (SMF). Stress relaxation tests were done by applying a step
strain $\gamma_{0}=2\%$ to materials and the evolution of the shear
viscoelastic modulus $G(t;t_{\rm age},T)$ was monitored over 1000 s, where
$t_{\rm age}$ is the age of the sample and $T$ is the absolute temperature at
which the gels were stored. For rheometry, the samples were loaded between 25
mm hatched parallel plates to mitigate wall slip artifacts, and the
temperature was controlled by an Advanced Peltier System. Five sample ages
were used for the aging study: 0 ($\sim 1~{}{\rm h}$), 1, 3, 7, and 15 days.
Three temperatures of storage were used to study the effect of temperature on
the aging: 4, 25, 37∘C (at $p$H 7). Effect of $p$H was tested by making gels
in sodium citrate buffer with $p$H 4 compared with gels made in water, stored
at 25∘C. All tests were duplicated.

###  2.3 FTIR spectroscopy

The extent of reaction was inferred using Fourier Transform Infrared
Spectroscopy (FTIR) on each sample. Attenuated Transmission Reflectance (ATR)
setup was used to probe each sample on a Nicolet-iS50 spectrometer at Stanford
Soft and Hybrid Materials Facility (SMF). Each sample was aged for the
requisite duration and lyophilized at the end of the aging process; the
lyophilized powder was loaded on the diamond window of the ATR setup and the
reflectance was observed using a KBr beamsplitter for wavenumbers
$\bar{\nu}=400$-5000 cm-1. Sample data was collected in duplicates. The
absorbance spectrum of the material in this range of wavenumbers was used to
further quantify the decay or proliferation of bonds due to aging.

##  3 Results

In this work, we first discuss the rheological signatures of aging in our HEC-
MC/CSP hydrogel system and comment on the thermodynamic nature of the aging
process, along with the temperature and $p$H dependence. We then provide
structural insights into this aging process obtained from diffusivity
measurements, which suggests negligible alterations to the network structure
during aging. Based on this, we describe the physics of the aging process and
hypothesize the possible interfacial interactions that may underlie this
unique behavior. We then provide corroborating evidence for the hypothesized
chemical changes using infrared spectroscopy data that sheds more light on the
nature of interfacial chemistry between the polymer-particle components of the
hydrogel system, and hence rationalize the evolution of chemical bonds in the
material with aging time. Following this, we model the chemical reaction
underlying the aging process. We use second-order reaction kinetics to model
the dynamic-to-covalent transition and infer microstructural changes in the
material by connecting these kinetics to bulk mechanical properties using
network theory. Finally, we summarize the findings of this paper and comment
on the outlook for engineering soft materials with such novel and highly
tunable mechanochemical properties.

It is common practice to study the stiffness of aging systems using mechanical
perturbation, such as rheometry. Often, percolating systems, aggregating
dispersions, soft glasses, and crowded suspensions are probed using stress
relaxation, which provides long-time information about material properties
while keeping the sample deformation within the linear regime [30, 31, 32, 33,
34, 35, 36, 37, 38, 39]. In this work, we first determined the stiffness of
our HEC-MC/CSP gels from the short-time plateau modulus for each gel
formulation, $G_{0}$, as a function of aging time to quantify the aging
process. We then studied the effect of temperature and $p$H on the aging
process by monitoring the evolution of plateau modulus values over time. The
plateau modulus can be connected to the microstructure of the gels using
network theory, thus quantifying the various chemical reactions and structural
evolution networks undergo when aged (vide infra).

###  3.1 Temperature-dependent aging

Stress relaxation data for gels at different timepoints following aging at
three different temperatures demonstrates that, for a given temperature, the
gels stiffen with age (Fig. 2). As $t_{\rm age}$ increases, the relaxation
curve shifts towards higher modulus, resulting in higher values of $G_{0}$.
Additionally, the relaxation behavior of the gels also changes. The stress
relaxation transients become broader and the characteristic timescale of
relaxation becomes longer, which is evident from the slower decay of
relaxation modulus with age. This observation suggests that the material
possesses multiple timescales of relaxation, or a relaxation spectrum, and not
just one single characteristic timescale [40, 41]. Polymeric materials,
especially those with transient crosslinks, often exhibit such behavior. These
signatures of the stress relaxation behavior indicate a broadening in the
dynamics of the gels with age, linked with the stiffening behavior and
proliferation of stronger crosslinks through the material.

(a)

(b)

(c)

Figure 2: Stress relaxation for gels at different aging times at (a) 4∘C, (b)
25∘C, and (c) 37∘C. For a given temperature, gels get stiffer with age, and
its relaxation timescale becomes longer. The same features are observed as
temperature increases for a given age of the gel.

As the temperature increases, the aging becomes more pronounced. Indeed, the
modulus increases faster with age, and the relaxation behavior becomes even
broader. Clearly, temperature plays a key role in the structural change of the
gels, and the reaction responsible for such a change proceeds faster at higher
temperatures. These behaviors are characteristic of endothermic covalent
reaction kinetics, which we hypothesized may occur over time between the
alcohol groups on the cellulose chains and the silanol groups on the CSP
surface, which are in close proximity to each other due to the strong yet
dynamic non-covalent polymer-particle interactions. As such a condensation
reaction would be expected to be $p$H-dependent, we then evaluated the effect
of $p$H on aging.

###  3.2  _p_ H-dependent aging

Stress relaxation data for gels prepared at two $p$H values ($p{\rm H}\simeq
4$ and 7) and stored at 25∘C were obtained for different aging times (Fig. 3).
Similar to our previous observations, the modulus increases and the relaxation
spectrum gets broader with age. However, both the relative magnitudes of
$G_{0}$ and the broadening of the relaxation spectrum is decreased at lower
$p$H.

(a)

(b)

Figure 3: Stress relaxation for gels at different aging times at (a) $p{\rm
H}=4$, and (b) $p{\rm H}=7$, all gels at $25^{\circ}$C. For a given $p$H, gels
get stiffer with age, and its relaxation timescale becomes longer. The same
features are observed as $p$H increases for a given age of the gel.

These data indicate that the reaction responsible for aging in the gels
proceeds slower in more acidic conditions, while increasing temperature speeds
up the reaction. These observations suggest that a potential driver of the
aging process is an interfacial covalent reaction between the silanol groups
on the CSPs and the alcohol groups on the cellulose polymers. Such a reaction
would be expected to proceed faster at higher temperatures (whereby the
enthalpic contribution to the free energy of change dominates), and faster in
a more alkaline environment as deprotonated silanol groups are more reactive,
favoring a condensation reaction with the cellulose alcohol groups. In the
following section, we discuss the thermodynamic nature of these aging
processes.

###  3.3 Non-ergodicity in aging gels

Typically, for aging gel systems, the underlying process of aging is
classified as non-ergodic (i.e., in violation of time-invariant physics) [33].
Briefly, non-ergodic physical processes proceed such that upon changing the
starting point (i.e., equivalent to shifting the location of $t=0$), one
changes the physics of the problem, thus breaking time-translation symmetry.
The evolution of the relaxation spectrum of our HEC-MC/CSP gels with age
corroborates this observation. The non-ergodicity is captured using a
nonlinear time-axis shift that, in theory, should collapse the data. The
shifting is based on the time of aging itself, $t_{\rm age}$, and a shifting
exponent, $\mu$, which generates master curves for the aging gels [32, 33,
34]. The time-axis is scaled, in a sort of coordinate transformation wherein
the passage of time is either “stretched” or “compressed”, using a shift
factor $\theta$, where $\theta\propto t_{\rm age}^{\mu}$. This form of
shifting works well for non-ergodic systems. It is important to note that
$\mu$ and therefore $\theta$ are functions of temperature, and successful
shifting can be conducted for data at different temperatures (Fig. 4).

(a)

(b)

(c)

Figure 4: Master curves of stress relaxation data shifted to account for
aging, for gels at different aging times at (a) 4∘C, (b) 25∘C, and (c) 37∘C.
Such shifting captures non-ergodic effects in the material.

For stress relaxation data at each temperature, applying the transformation
$t_{\rm age}:\rightarrow t_{\rm age}/\theta$ using the appropriate value of
$\mu$ in $\theta\propto t_{\rm age}^{\mu}$ collapses the data onto master
curves (n.b., the modulus data has also been normalized by its plateau value
to achieve this). The same has been repeated with success for the data at two
different $p$H values (Fig. 5).

(a)

(b)

Figure 5: Master curves for stress relaxation for gels at different aging
times at (a) $p{\rm H}=4$, and (b) $p{\rm H}=7$. Such shifting captures non-
ergodic effects in the material.

(a)

(b)

Figure 6: Power law correlation between shift factor and aging time,
$\theta(t_{\rm age},T)\propto t_{\rm age}^{\mu(T)}$. (a) Extent of correlation
between gels at different temperatures, with higher power law slopes at higher
temperatures. (b) Extent of correlation between gels at different $p$H, with
higher power law slope at higher $p$H.

The successful shifting is characteristic of non-ergodic systems, and is
displayed by systems that age through the growth of interconnected structure
elements, in this case being crosslinks between the cellulose chains
facilitated at the surface of the CSPs. Such aging systems also tend to be
strongly enthalpic (and therefore strongly dependent on temperature), with
entropy contributing much less to the free energy.

The correlation between $\theta$ and $t_{\rm age}$ is a power law (Fig.6), and
we see that $\mu$ is higher at higher $T$, which suggests a stronger
correlation and therefore stronger aging. These observations are in full
agreement with the rheological data discussed earlier. The same trend is
maintained for $p$H dependence, where $\mu$ is higher at higher $p$H, again
suggesting a stronger correlation and therefore stronger aging. To further
support these observations, we sought to gain structural insight into the
effect of aging on the network characteristics by studying the diffusion of
polymeric cargo through the gel, as discussed next.

###  3.4 Structural insights into hydrogel aging from diffusivity

The observed trends in rheology under the effect of temperature and $p$H is
suggestive of the nature of interfacial reactions underlying the aging of the
hydrogel. Accordingly, we hypothesized that the number density of crosslinks
between polymer chains at particle surfaces would remain largely the same upon
mixing the components to make the hydrogel. Aging would not form new crosslink
junctions, but only convert a fraction of dynamic noncovalently bound sites
into covalently bound anchors. Since the average mesh size of the polymer
matrix is determined by the number density of crosslinks, we hypothesized this
system would maintain a rather consistent mesh size, and the stiffening only
stems from an aging of each bond rather than an increase in the number of
bonds. This process is illustrated schematically in Fig. 7. Note that
immediately upon mixing, transient multivalent crosslinks form within the
polymer network _via_ chain adsorption at the particle surfaces. The particle
surfaces, therefore, are already populated with polymer chains and the number
density of polymer-particle interactions has already reached its maximum
extent. Yet, while these polymer-particle interactions are dynamic and non-
covalent, they can act as anchoring sites. In this paper, when we refer to
crosslinks, we shall mean network anchor points which irreversibly connect
chains _via_ covalent bonding with particle surfaces. It must be kept in mind
that, on average, the total number of network anchor points remains the same,
yet the fraction of covalent crosslinks increases while the fraction of
dynamic bonds decreases as the material ages as polymer-particle interactions
undergo a dynamic-to-covalent transition.

Figure 7: (a) Diffusivity data for three cargo sizes obtained from FRAP
experiments as functions of aging time. (b) Schematic for the hypothesized
dynamic-to-covalent transition scheme for the hydrogel mesh, leading to aging.

To evaluate our hypothesis, we used diffusivity measurements of polymeric
solutes through the hydrogel using fluorescence recovery after photobleaching
(FRAP) tests (Fig. 7a). For FRAP measurements, we used dextran probes
conjugated with a fluorescin-isothiocyanate dye (FITC-dextran) as the solute
which diffuses through the hydrogel mesh. Three average molecular weights of
the polymer were used (4, 40, and 150 kDa), which correspond to three solute
“sizes” pertaining to their radii of gyration $R_{\rm g}$ in solution (1.4,
4.5, and 8.8 nm respectively). The diffusion coefficient $D$ as a function of
aging time $t_{\rm age}$ of each solute obtained from FRAP is compared against
the ideal Stokes-Einstein diffusivity, $D_{0}=k_{\rm B}T/6\pi\eta_{\rm
s}R_{\rm g}$, where $\eta_{\rm s}$ is the solvent viscosity, taken as 1
mPa$\cdot$s for water, $k_{\rm B}$ is the Boltzmann constant, and $T$ is the
absolute temperature of the system (25∘C). We see that $D/D_{0}$ negligibly
changes over time (i.e. $<15\%$ change in diffusivity over 15 days) for all
three solute sizes. In contrast, rheological measurements showed that the
shear modulus increased by more than 250% in the same duration. Since
$D/D_{0}$ of passively diffusing solutes like polymer globules is a function
of mesh size [42], these observations suggest that the average mesh size of
the hydrogel matrix does not meaningfully change during aging while the
modulus increases significantly by virtue of strengthening of the crosslinking
interactions. Below we propose a reaction scheme underlying the slow
transition from non-covalent reversible polymer-particle binding in the
hydrogels to irreversibly covalent bonded gels while preserving the number
density of crosslinks, which is synergistic with the rheology and diffusion
data.

###  3.5 Spectroscopic characterization of aging within HEC-MC/CSP gels

We hypothesized that each dynamic, noncovalent bond between the cellulosic
hydroxyl groups and the silanol groups on the surface of the CSPs is slowly
converted to a (particle)-Si-O-(polymer) bond with age, and consequently the
topology of the crosslinked chains does not change (Fig. 8). In this scheme,
the total number of polymer-particle interactions does not change as the
hydrogels age, but these interactions only strengthen with time as they
undergo the proposed dynamic-to-covalent transition. Such aging has been
observed previously in the literature for physical gels of gelatin in water,
where chains form a noncovalent network upon dissolution and these junctions
undergo a perpetual aging process that results in gel stiffening [34, 43].
Consistent with our hypothesis, the number density of crosslinks within these
gelatin-based systems does not change significantly with time, thus preserving
the effective mesh size of the network. We then sought to gather direction
spectroscopic evidence of the chemical changes occurring during aging and
discuss how it supports the proposed reaction scheme and mechanism
hypothesized from rheological data described above.

Figure 8: Hypothesized interfacial interactions between HEC-MC chains and CSP
surface. Shortly after mixing, the polymer chains interact reversibly with the
particle surface, which includes hydrogen bonding effects. This slowly
transitions to irreversible conjugation as the material ages, the rate of
which is dependent on temperature and $p$H.

In this section, we monitor the aging processes in our gels with FTIR
spectroscopy to generate further insight into the dynamic-to-covalent
crosslinking process and quantify the rate of structural stiffening. We focus
on the wavenumber range 3200-3600 cm-1, since the absorption peak $(a_{\rm
m})$ of the cellulosic hydroxyl groups occurs at approximately
$\bar{\nu}\approx 3410~{}{\rm cm}^{-1}$. Per our hypothesis, the condensation
reaction between the alcohol groups on the cellulose polymers and the silanol
group on the CSP surface would lead to the disappearance of a fraction of
these -OH bonds and commensurate formation of (particle)-Si-O-(polymer) bonds.
Consequently, we would expect to observe a decay in this $3410~{}{\rm
cm}^{-1}$ absorption peak as the sample ages. We looked for this signature in
aging gels to study both the effects of temperature and $p$H. We only show
absorbance data for the relevant range of wavenumbers here, but transmittance
data over the full spectral range accessed during our FTIR tests are shown in
the SI (Figs. S3 and S4).

Figure 9: FTIR data for gels through aging at different temperatures. (a)
Absorbance data for gels at three different temperatures, where aging leads to
a decay in the absorbance peak at 3410 cm-1, $a_{\rm m}$, which is accelerated
as the temperature increases. (b) Decay behavior of $a_{\rm m}$ as function of
$t_{\rm age}$ and temperature, showing faster disappearance of the cellulosic
hydroxyl peak at higher temperatures. Fit lines are for an exponential
function. (c) Peak decay timescale $\tau_{\rm m}$ from fitting of the FTIR
data to an exponential decay function, showing that the timescale of decay
shortens at elevated temperatures. Figure 10: FTIR data for gels through aging
at different $p$H values. (a) Absorbance data for gels at $p$H of 4 and 7,
where aging leads to a decay in the absorbance peak at 3410 cm-1, $a_{\rm m}$,
which is accelerated as the $p$H of the gel increases. (b) Decay behavior of
$a_{\rm m}$ as function of $t_{\rm age}$ and $p$H, showing slower
disappearance of the cellulosic hydroxyl peak in more acidic conditions. Fit
lines are for an exponential function. (c) Peak decay timescale $\tau_{\rm m}$
from fitting FTIR data to an exponential decay function, showing that the
timescale of decay slows in more acidic conditions.

For a given temperature, we see that the absorption spectrum flattens with
aging time (Fig. 9a), and the absorbance value of peak at $3410~{}{\rm
cm}^{-1}$, $a_{\rm m}$, diminishes with $t_{\rm age}$. As the temperature at
which the gels are stored is increased from 4 to 25 to 37∘C, the peak
diminishes more rapidly (Fig. 9b). This decrease in the peak absorbance can be
captured by a simple exponential decay. The decay timescale, $\tau_{\rm m}$,
clearly decreases as the temperature of gel aging increases, indicating faster
aging at elevated temperatures (Fig. 9c). These observations are in direct
agreement with the rheology and kinetics modeling data described above, which
showed an increase in stiffness with age resulting from a proliferation of
covalent crosslinks over time.

Similarly, as the $p$H of the gels is increased from 4 to 7, the absorbance
spectrum flattens out more rapidly (Fig. 10a), and the peak of the spectrum,
$a_{\rm m}$, decays more quickly (Fig. 10b). This decrease in the peak
absorbance has been captured by a simple exponential decay, and the decay
timescale also decreases with increased $p$H, thus showing an accelerated
decay and therefore faster aging in more alkaline conditions (Fig. 10c). These
FTIR spectroscopy studies support the proposed reaction scheme and enable
quantification of the rate and extent of dynamic-to-covalent transitions
occurring over time for various aging conditions. Finally, for further
quantification, we model the kinetics of this proposed reaction, and infer the
progression of crosslink formation during aging through the rheological data
shown here. We relate the chemical crosslinking strength to plateau modulus of
the gels, giving insight into the energetics of the process.

###  3.6 Modeling of aging processes

We propose a simple modeling scheme for the chemical reaction underlying the
observed aging behaviors to quantify the process. The reaction kinetics
associated with the crosslink proliferation is connected back to rheological
signatures using network theory, which allows us to estimate energetics of the
covalent crosslinking process.

The evolution of the hydrogel network can be modeled using a second-order
reaction kinetics model. For macromolecular crosslinking systems involving two
reactants, the gelation kinetics have been compared to the predictions of a
simple second-order reaction kinetics scheme in the literature [43, 44, 45].
In our HEC-MC/CSP gels, where the cellulosic polymers and CSP form the two
components forming crosslinks by reacting according to the reaction given by

$\displaystyle\mathcal{A}+\mathcal{B}$
$\displaystyle\xrightleftharpoons[k_{\rm d}]{k_{\rm a}}\mathcal{C},$ (1)

with $\mathcal{A}$ and $\mathcal{B}$ representing the two reactants (cellulose
polymers and colloidal silica particles), and $\mathcal{C}$ represents the
dynamically crosslinked gel product, with $k_{\rm a}$ and $k_{\rm d}$ the rate
constants for crosslink association and dissociation respectively. While the
gel would be expected to have a relatively small dissociation rate constant
compared to the association rate constant, we assume it to be non-zero for the
most general case. The differential equation for the rate kinetics appropriate
to this scheme is given by

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\xi$ $\displaystyle=k_{\rm
a}(a_{0}-\xi)(b_{0}-\xi)-k_{\rm d}\xi,$ (2)

where $a_{0}$ and $b_{0}$ are the reactive site concentrations of
$\mathcal{A}$ and $\mathcal{B}$ respectively, and $\xi$ is some measure of the
degree of crosslinking (e.g., the volume density of crosslinks, often called
the crosslink concentration). Eq. 2 can be rewritten as

$\displaystyle\dfrac{{\rm d}}{{\rm d}t}\xi$ $\displaystyle=k_{\rm
a}(a-\xi)(b-\xi),$ (3)

and the solution to this differential equation is

$\displaystyle\xi$ $\displaystyle=a\dfrac{1-\exp\left[k_{\rm
a}(a-b)t\right]}{1-\dfrac{a}{b}\exp\left[k_{\rm a}(a-b)t\right]},$ (4)

which can be rewritten as

$\displaystyle\xi$ $\displaystyle=a\dfrac{1-\exp\left(\kappa
t\right)}{1-\rho\exp\left(\kappa t\right)},$ (5)

where $\kappa=k_{\rm a}(a-b)$ and $\rho=a/b$. Here, $\kappa$ is a rescaled
rate of reaction, and $\rho$ is the molar ratio of the two reactants. Using
$\kappa$, which is expected to be a function of temperature, we can capture
the rate of crosslink formation and therefore the speed and extent of aging in
the gels, and thus connect it to the shear modulus of the material using
network theory (vide infra).

###  3.7 Relating kinetics to rheology using network theory

From the stress relaxation data, the plateau modulus $G_{0}$ is defined as

$\displaystyle G_{0}(t_{\rm age},T)\equiv\lim_{t\rightarrow 0}G(t;t_{\rm
age},T),$ (6)

and it is a function of aging time $t_{\rm age}$ and dependent on temperature
$T$. From plots of $G_{0}$ at different temperatures (Fig. 11a) and $p$H (Fig.
11b), it is clear that the gels stiffen with aging time, and the stiffening is
stronger at higher temperatures and higher $p$H. These observations are in
full agreement with the experimental results described above.

(a)

(b)

Figure 11: (a) Effect of temperature on aging. Plateau modulus as a function
of aging time for both temperatures, showing that ages progresses slower at a
lower temperature. (b) Effect of $p$H on aging. Plateau modulus as a function
of aging time for both temperatures, showing that ages progresses slower at a
lower $p$H. Also shown are fits to the second order reaction model, Eq. 8.

This modulus can also be related to the degree of crosslinking using network
theory, where the network stiffness is proportional to the degree of
crosslinking [43, 41, 40, 46, 47, 48, 49, 44, 45]. In short, we have

$\displaystyle G_{0}(t_{\rm age},T)\propto\xi(t_{\rm age},T).$ (7)

Leveraging the second-order model solution for $\xi$ described above, we can
rewrite Eq. 5 into our model for the aging modulus as

$\displaystyle G_{0}(t_{\rm age},T)$ $\displaystyle=G_{\rm
b}(T)+G_{\xi}(T)\dfrac{1-\exp\left[\kappa(T)t_{\rm
age}\right]}{1-\rho\exp\left[\kappa(T)t_{\rm age}\right]},$ (8)

where $G_{\rm b}$ is the background viscoelastic modulus independent of the
crosslinks (from nonreacting and purely viscoelastic components) which does
not change with $t_{\rm age}$ but is dependent on temperature, and $G_{\xi}$
is the crosslink-dependent modulus, which changes with time/age, and is also
dependent on temperature. When fitting Eq. 8 to the modulus data for our gels
(Fig. 11), we see that the simple second order reaction kinetics model
captures the aging modulus very well. This approach is useful since the fit
parameter $\kappa$ is common for the modulus and the crosslink density
equations, which means we can use the fit to rheology data to infer
information about crosslinking reaction energetics.

###  3.8 Activation energy of aging reaction

From the solution of the second order reaction (Eq. 5), we had defined

$\displaystyle\kappa=k_{\rm a}(a-b),$ (9)

and also, by definition, the equilibrium rate constant of the reaction is

$\displaystyle K_{\rm eq}\equiv\frac{k_{\rm a}}{k_{\rm d}},$ (10)

which gives the relation of the equilibrium rate constant and the rate of
reaction to be

$\displaystyle\kappa=K_{\rm eq}k_{\rm d}(a-b).$ (11)

Assuming an Arrhenius rate activated process, we have

$\displaystyle K_{\rm eq}(T)=A\exp\left(-\frac{E_{\rm a}}{RT}\right),$ (12)

which finally gives

$\displaystyle\kappa(T)=Ak_{\rm d}(a-b)\exp\left(-\frac{E_{\rm
a}}{RT}\right)\equiv B\exp\left(-\frac{E_{\rm a}}{RT}\right).$ (13)

From this, we can write

$\displaystyle\ln\kappa(T)=\ln B-\frac{E_{\rm a}}{RT}.$ (14)

With this equation we can now calculate the activation energy of the covalent
crosslinking reaction driving aging in these systems from the reaction rate
parameter $\kappa$ obtained from fitting rheology data to our second order
model. We can fit the $\ln\kappa(T)$ vs. $1/T$ data to a linear function, and
calculate the activation energy $E_{\rm a}$ from the slope of this line (Fig.
12).

Figure 12: Activation energy of the irreversibly conjugating HEC-MC/CSP
hydrogel reaction, obtained by connecting the rheological fit parameter
$\kappa$ using the second order model to the Arrhenius dependence of
equilibrium rate constant.

From the fit, the activation energy is calculated to be $E_{\rm a}=52.9~{}{\rm
kJ}~{}{\rm mol}^{-1}$. Compared to the free energy of dynamic crosslink
formation in these HEC-MC/CSP gels [25], which was determined to be $\Delta
G_{0}=-18.9~{}{\rm kJ}~{}{\rm mol}^{-1}$, the transition to covalent
crosslinking has a much higher activation energy, further supporting that this
is a temperature activated process.

##  4 Discussion

In this work, we developed an entropy-driven physical hydrogel system that
undergoes a unique dynamic-to-covalent transition in the crosslinking that
mimics behaviors observed in many naturally occurring extracellular matrix
proteins. Immediately upon mixing, cellulose derivatives (HEC and MC) and CSPs
engage in multivalent, dynamic interactions as the polymers adsorb to the
interface of the particles, leading to the formation of HEC-MC/CSP hydrogels.
Yet, these polymer-particle interactions undergo condensation reactions
leading to covalent bond formation between the cellulose polymer chains and
the CSP surface. This transitions drives aging of the materials, whereby they
get stiffer over time i.e., the modulus increases with aging time), yet
without any observable alteration of the hydrogel structure and mesh size.
These chemical reactions occurring only where the dynamic polymer-particle
interactions are already formed are accelerated rates at elevated temperatures
and higher $p$H values. These underlying chemical changes responsible for
aging in the HEC-MC/CSP gels is supported with both rheology and spectroscopy
data, suggests a slow proliferation of covalent crosslinking _via_ interfacial
interactions over time.

Such unique behavior suggests that the extent of crosslinking _via_
irreversibly conjugated bonds and consequently the mechanical strength of this
system can thus be tuned in multiple ways. Firstly, different temperatures age
the gels differently, and a suitable duty cycle of temperatures can be
employed to achieve various rates of aging. Secondly, diminishing the density
of hydroxyl groups on the cellulose polymers by using alternative cellulosic
variants (e.g., methoxy rather than hydroxyl functionality) should directly
inhibit the dynamic-to-covalent transitions and associated aging. One such
example is using hydroxyethylmethylcellulose (HEMC) instead of HEC or MC as it
has fewer hydroxyl groups along the polymer chain as they are replaced by
methoxy groups, which dramatically attenuates aging (Fig. S1). Thirdly,
varying the $p$H can also modify the aging behavior. Therefore, depending on
the needs in any given application area, slower aging gels can be subjected to
elevated temperatures for speeding up the aging process which is otherwise
stunted due to substitution or $p$H.

Using time-aging and temperature-stiffening hydrogels based on covalent
crosslinking energetics through engineered interfacial interactions, a
complementary paradigm to temperature-softening enthalpy-dominated dynamic
hydrogels is thus established. We have engineered a transition from an
entropy-dominated dynamic crosslinking system to a covalently crosslinked
hydrogel system without altering the structure of the matrix, the strength of
which can further be tuned based on chemistry (choice of cellulose
modification), formulation ($p$H), and temperature. The biomimetic non-ergodic
aging through dynamic-to-covalent transitions in these synthetic physical
hydrogels enable the development of materials with structures defined by the
self-assembly of the gel components that then evolve into versatile covalent
hydrogels with tunable mechanical properties. The properties of these systems
can be modulated through minor changes in formulation and storage conditions,
making them amenable to numerous applications.

## Acknowledgement

The rheology and spectroscopy data in this paper was collected at the Stanford
Soft and Hybrid Materials Facility (SMF) (supported by the NSF grant National
Science Foundation grant ECCS-1542152). We thank the Gordon & Betty Moore
Foundation for their financial support of this work as part of our efforts to
develop wildland fire retardants to improve the execution of prescribed burns
and enhance forest management.

## References

  * Appel et al. [2012] Appel, E. A.; del Barrio, J.; Loh, X. J.; Scherman, O. A. Supramolecular polymeric hydrogels. _Chem. Soc. Rev._ 2012, _41_ , 6195–6214.
  * Mann et al. [2017] Mann, J. L.; Yu, A. C.; Agmon, G.; Appel, E. A. Supramolecular polymeric biomaterials. _Biomater. Sci._ 2017, _6_ , 10–37.
  * Webber et al. [2016] Webber, M. J.; Appel, E. A.; Meijer, E. W.; Langer, R. Supramolecular biomaterials. _Nat. Mater._ 2016, _15_ , 13–26.
  * Seiffert [2015] Seiffert, S. _Supramolecular Polymer Networks and Gels_ , 1st ed.; Springer, 2015.
  * Yu et al. [2021] Yu, A. C.; Lian, H.; Kong, X.; Hernandez, H. L.; Qin, J.; Appel, E. A. Physical networks from entropy-driven non-covalent interactions. _Nat. Comm._ 2021, _12_ , 746.
  * Rodell et al. [2013] Rodell, C. B.; Kaminski, A. L.; Burdick, J. A. Rational design of network properties in guest-host assembled and shear-thinning hyaluronic acid hydrogels. _Biomacromolecules_ 2013, _14_ , 4125–4134.
  * Rodell et al. [2016] Rodell, C. B.; Dusaj, N. N.; Highley, C. B.; Burdick, J. A. Injectable and cytocompatible tough double-network hydrogels through tandem supramolecular and covalent crosslinking. _Adv. Mater._ 2016, _28_ , 8419–8424.
  * Hernandez et al. [2018] Hernandez, H. L.; Grosskopf, A. K.; Stapleton, L. M.; Agmon, G.; Appel, E. A. Non-Newtonian polymer-nanoparticle hydrogels enhance cell viability during injection. _Macromol. Biosci._ 2018, _19_ , 1800275.
  * Grosskopf et al. [2019] Grosskopf, A. K.; Roth, G. A.; Smith, A. A. A.; Gale, E. C.; Hernandez, H. L.; Appel, E. A. Injectable supramolecular polymer–nanoparticle hydrogels enhance human mesenchymal stem cell delivery. _Bioeng. Transl. Med._ 2019, e70147.
  * Wong Po Foo et al. [2009] Wong Po Foo, C. T.; Lee, J. S.; Mulyasasmita, W.; Parisi-Amon, A.; Heilshorn, S. C. Two-component protein-engineered physical hydrogels for cell encapsulation. _Proc. Natl. Acad. Sci. U. S. A._ 2009, _106_ , 22067–22072.
  * Buwalda et al. [2017] Buwalda, S. J.; Vermonden, T.; Hennink, W. E. Hydrogels for therapeutic delivery: current developments and future directions. _Biomacromolecules_ 2017, _18_ , 316–330.
  * Li and Mooney [2016] Li, J.; Mooney, D. J. Designing hydrogels for controlled drug delivery. _Nat. Rev. Mater._ 2016, _1_ , 16071.
  * Mitragotri et al. [2014] Mitragotri, S.; Burke, P. A.; Langer, R. Overcoming the challenges in administering biopharmaceuticals: formulation and delivery strategies. _Nat. Rev. Drug Discov._ 2014, _13_ , 655–672.
  * Ma et al. [2021] Ma, J.; Porath, L. E.; Haque, M. F.; Sett, S.; Rabbi, K. F.; Nam, S.; Miljkovic, N.; Evans, C. M. Ultra-thin, self-healing vitrimer coatings for durable hydrophobicity. _Nat. Comm._ 2021, _12_ , 5210\.
  * Clarke et al. [2017] Clarke, D. E.; Pashuck, E. T.; Bertazzo, S.; Weaver, J. V. M.; Stevens, M. M. Self-healing, self-assembled beta-sheet peptide poly(gamma-glutamic acid) hybrid hydrogels. _J. Am. Chem. Soc._ 2017, _139_ , 7250–7255.
  * Rao et al. [2016] Rao, Y.-L.; Chortos, A.; Pfattner, R.; Lissel, F.; Chiu, Y.-C.; Feig, V.; Xu, J.; Kurosawa, T.; Gu, X.; Wang, C.; He, M.; Chung, J. W.; Bao, Z. Stretchable Self-Healing Polymeric Dielectrics Cross-Linked Through Metal–Ligand Coordination. _J. Am. Chem. Soc._ 2016, _139_ , 6020–6027.
  * Cordier et al. [2008] Cordier, P.; Tournilhac, F.; Soulie-Ziakovic, C.; Leibler, L. Self-healing and thermoreversible rubber from supramolecular assembly. _Nature_ 2008, _451_ , 977–980.
  * Li et al. [2016] Li, C.-H.; Wang, C.; Keplinger, C.; Zuo, J.-L.; Jin, L.; Sun, Y.; Zheng, P.; Cao, Y.; Lissel, F.; Linder, C.; You, X.-Z.; Bao, Z. A highly stretchable autonomous self-healing elastomer. _Nat. Chem._ 2016, _8_ , 618–624.
  * Oh et al. [2016] Oh, J. Y. et al. Intrinsically stretchable and healable semiconducting polymer for organic transistors. _Nature_ 2016, _539_ , 411–415.
  * Yan et al. [2018] Yan, X. et al. Quadruple H-Bonding Cross-Linked Supramolecular Polymeric Materials as Substrates for tretchable, Antitearing, and Self-Healable Thin Film Electrodes. _J. Am. Chem. Soc._ 2018, _140_ , 5280–5289.
  * Someya et al. [2016] Someya, T.; Bao, Z.; Malliaras, G. G. The rise of plastic bioelectronics. _Nature_ 2016, _540_ , 379–385.
  * Ouyang et al. [2016] Ouyang, L.; Highley, C. C.; Rodell, C. B.; Sun, W.; Burdick, J. A. 3D printing of shear-thinning hyaluronic acid hydrogels with secondary cross-linking. _A. C. S. Biomater. Sci. Eng._ 2016, _2_ , 1743–1751.
  * Dubbin et al. [2016] Dubbin, K.; Hori, Y.; Lewis, K. K.; Heilshorn, S. C. Dual-stage crosslinking of a gel-phase bioink improves cell viability and homogeneity for 3D bioprinting. _Adv. Healthc. Mater._ 2016, _5_ , 2488–2492.
  * Murphy and Atala [2014] Murphy, S. V.; Atala, A. 3D bioprinting of tissues and organs. _Nat. Biotechnol._ 2014, _32_ , 773–785.
  * Yu et al. [2016] Yu, A. C.; Chen, H.; Chan, D.; Agmon, G.; Stapleton, L. M.; Sevit, A. M.; Tibbitt, M. W.; Acosta, J. D.; Zhang, T.; Franzia, P. W.; Langer, R.; Appel, E. A. Scalable manufacturing of biomimetic moldable hydrogels for industrial applications. _Proc. Natl. Acad. Sci. U. S. A._ 2016, _113_ , 14255–14260.
  * Yu et al. [2019] Yu, A. C.; Lopez Hernandez, H.; Kim, A. H.; Stapleton, L. M.; Brand, R. J.; mellor, E. T.; Bauer, C. P.; McCurdy, G. D.; Wolff, A. J.; Chan, D.; Criddle, C. S.; Acosta, J. D.; Appel, E. A. Wildfire prevention through prophylactic treatment of high-risk landscapes using viscoelastic retardant fluids. _Proc. Natl. Acad. Sci. U. S. A._ 2019, _116_ , 20820–20827.
  * Appel et al. [2010] Appel, E. A.; Biedermann, F.; Rauwald, U.; Jones, S. T.; Zayed, J. M.; Scherman, O. A. Supramolecular cross-linked networks via host-guest complexation with cucurbit[8]uril. _J. Am. Chem. Soc._ 2010, _132_ , 14251–14260.
  * Yu et al. [2019] Yu, A. C.; Smith, A. A. A.; Appel, E. A. Structural considerations for physical hydrogels based on polymer–nanoparticle interactions. _Mol. Syst. Des. Eng._ 2019, _5_ , 401–407.
  * Correa et al. [2021] Correa, S.; Grosskopf, A. K.; Hernandez, H. L.; Chan, D.; Yu, A. C.; Stapleton, L. M.; Appel, E. A. Translational applications of hydrogels. _Chem. Rev._ 2021, _121_ , 11385–11457.
  * Sollich et al. [1997] Sollich, P.; Lequeux, F.; Hébraud, P.; Cates, M. E. Rheology of Soft Glassy Materials. _Phys. Rev. Lett._ 1997, _78_ , 2020–2023.
  * Sollich [1998] Sollich, P. Rheological constitutive equation for a model of soft glassy materials. _Phys. Rev. E_ 1998, _58_ , 738–759.
  * Cloitre et al. [2000] Cloitre, M.; Borrega, R.; Leibler, L. Rheological Aging and Rejuvenation in Microgel Pastes. _Phys. Rev. Lett._ 2000, _85_ , 4819–4822.
  * Fielding et al. [2000] Fielding, S. M.; Sollich, P.; Cates, M. E. Aging and rheology in soft materials. _J. Rheol._ 2000, _44_ , 323–369.
  * Ronsin et al. [2009] Ronsin, O.; Caroli, C.; Baumberger, T. Interplay between Shear Loading and Structural Aging in a Physical Gelatin Gel. _Phys. Rev. Lett._ 2009, _103_ , 138302.
  * Joshi and Petekidis [2018] Joshi, Y. M.; Petekidis, G. Yield stress fluids and ageing. _Rheol. Acta_ 2018, _57_ , 521–549.
  * Fielding [2020] Fielding, S. M. Elastoviscoplastic rheology and aging in a simplified soft glassy constitutive model. _J. Rheol._ 2020, _64_ , 723–738.
  * Sen and Ewoldt [2022] Sen, S.; Ewoldt, R. H. Thixotropic spectra and Ashby-style charts for thixotropy. _J. Rheol._ 2022, _66_ , 1041–1053.
  * Sen [2022] Sen, S. Meaningful descriptions of thixotropy and extensibility for yield-stress fluid drop impact on thin films. Ph.D. thesis, University of Illinois Urbana-Champaign, 2022.
  * Wang and Ewoldt [2022] Wang, Y.; Ewoldt, R. H. New insights on carbon black suspension rheology – anisotropic thixotropy and anti-thixotropy. _J. Rheol._ 2022, _66_ , 937–953.
  * Bird et al. [1987] Bird, R. B.; Armstrong, R. C.; Hassager, O. _Dynamics of Polymeric Liquids_ ; Wiley, 1987; Vol. 1.
  * Ferry [1980] Ferry, J. D. _Viscoelastic Properties of Polymers_ , 3rd ed.; Wiley, 1980\.
  * Axpe et al. [2019] Axpe, E.; Chan, D.; Offeddu, G. S.; Chang, Y.; Merida, D.; Lopez Hernandez, H.; Appel, E. A. A multiscale model for solute diffusion in hydrogels. _Macromol._ 2019, _52_ , 6889–6897.
  * Normand et al. [2000] Normand, V.; Muller, S.; Ravey, J.-C.; Parker, A. Gelation Kinetics of Gelatin: A Master Curve and Network Modeling. _Macromol._ 2000, _33_ , 1063–1071.
  * Miller et al. [2019] Miller, K. A.; Morado, E. G.; Samanta, S. R.; Walker, B. A.; Nelson, A. Z.; Sen, S.; Tran, D. T.; Whitaker, D. J.; Ewoldt, R. H.; Braun, P. V.; Zimmerman, S. C. Acid-triggered, acid-generating, and self-amplifying degradable polymers. _J. Am. Chem. Soc._ 2019, _141_ , 2838–2842.
  * Xu et al. [2020] Xu, Y.; Sen, S.; Wu, Q.; Zhong, X.; Ewoldt, R. H.; Zimmerman, S. C. Base-triggered self-amplifying degradable polyurethanes with the ability to translate local stimulation to continuous long-range degradation. _Chem. Sci._ 2020, _11_ , 3326–3331.
  * Rubinstein and Colby [2003] Rubinstein, M.; Colby, R. H. _Polymer Physics_ ; Oxford University Press, 2003\.
  * Graessly [2003] Graessly, W. W. _Polymeric Liquids & Networks: Structure and Properties_, 1st ed.; Garland Science, 2003.
  * Graessly [1975] Graessly, W. W. Statistical mechanics of random coil networks. _Macromol._ 1975, _8_ , 186–190.
  * Chompff and Duiser [1966] Chompff, A. J.; Duiser, J. A. Vicoelasticity of networks consisting of crosslinked or entangled macromolecules. I. Normal modes and mechanical spectra. _J. Chem. Phys._ 1966, _45_ , 1505–1514.


# Efficient Preparation of Nonabelian Topological Orders in the Doubled
Hilbert Space

Shang Liu  sliu.physgmail.com Kavli Institute for Theoretical Physics,
University of California, Santa Barbara, California 93106, USA

###### Abstract

Realizing nonabelian topological orders and their anyon excitations is an
esteemed objective. In this work, we propose a novel approach towards this
goal: quantum simulating topological orders in the doubled Hilbert space – the
space of density matrices. We show that ground states of all quantum double
models (toric code being the simplest example) can be _efficiently_ prepared
in the doubled Hilbert space. In contrast, this is not the case in the
conventional Hilbert space: Ground states of only some of these models are
known to be efficiently preparable. Additionally, we find that nontrivial
anyon braiding effects, both abelian and nonabelian, can be realized in the
doubled Hilbert space, although the intrinsic nature of density matrices
restricts possible excitations.

##  I Introduction

In this Noisy Intermediate-Scale Quantum (NISQ) era, quantum simulation has
attracted significant amount of research efforts [1, 2]. One major question
being actively explored is how to realize certain useful or novel quantum
states on a realistic platform.

In the pursuit of this endeavor, it has been realized that _nonunitary
operations_ are extremely useful tools in many quantum simulation tasks.
Indeed, as shown in Refs. 3, 4, 5, 6, certain long-range entangled many-body
quantum states, such as the ferromagnetic cat state and some topologically
ordered states, can be prepared with constant layers of local unitary gates
and local measurements. In contrast, if only local unitaries are allowed, one
would need a unitary circuit whose depth increases with the system size.
However, although quantum measurements are efficient for generating long-range
entanglement, they have one obvious shortcoming – the randomness of outcomes.
In a limited set of cases discussed in the references above, it is possible to
efficiently cure the randomness from measurements. In other situations, it
seems that one has to do inefficient post-selections or seek alternative
protocols, unless when an intrinsically random state is wanted [7].

In this work, we propose a novel route of quantum simulation where
_deterministic_ nonunitary operations are possible: quantum simulation in the
_doubled Hilbert space_ – the space of density matrices. More explicitly, each
density matrix $\rho$ of a quantum state can be mapped to a pure state
$\ket{\rho}$ in the doubled space. The wave function of $\ket{\rho}$ is
nothing but the density matrix elements $\rho_{ij}$. Applying a quantum
channel on $\rho$ is equivalent to acting a generically nonunitary operator on
$\ket{\rho}$. Therefore, it is possible that interesting quantum states hard
to prepare in the conventional Hilbert space can be efficiently prepared as
states $\ket{\rho}$ in the doubled Hilbert space.

We explicitly demonstrate the advantage of the doubled Hilbert space by
realizing topologically ordered states, which is perhaps the most exciting
subfield of quantum simulation due to its direct relevance to fault-tolerant
quantum computation [8, 9, 10, 11, 12][13, 14, 3, 4, 5, 6, 15, 16]. More
specifically, we show that in the doubled Hilbert space, one can _efficiently_
prepare ground states of Kitaev’s quantum double models [17], or $G$-gauge
theories, for all finite groups $G$. In the conventional Hilbert space,
however, the same task is known to be possible only for a subset of groups
(solvable groups) [3, 4, 6].

This is not a free lunch though. The advantage of doubled Hilbert space comes
with a cost: Intrinsic properties of density matrices restrict the allowed
states and not all excited states of these topological orders are accessible.
Nonetheless, we show that nontrivial anyon braiding, a hallmark of topological
order, can be realized in the doubled Hilbert space. This is true for both
abelian and nonabelian braiding.

We shall mention that doubled-space topological order also appears in Refs.
15, 18. The difference is that these works assume a topologically ordered
state already exists in the conventional Hilbert space, while in this work, we
begin with a product state.

In the following, we will start from an abelian example: building a toric code
ground state in the doubled Hilbert space and realizing abelian anyon
braiding. We then discuss more general quantum doubled models and the
phenomenon of nonabelian anyon braiding.

##  II Example of Toric Code

Given a many-body density matrix $\rho$ and a real-space tensor product basis
$\\{\ket{i}\\}$ of the Hilbert space, we can write the density matrix as
$\rho=\sum_{i,j}\rho_{ij}\ket{i}\bra{j}$, and then map it to a pure state in
the doubled Hilbert space:

$\displaystyle\ket{\rho}$
$\displaystyle:=\frac{1}{\sqrt{\operatorname{Tr}(\rho^{2})}}\sum_{i,j}\rho_{ij}\ket{i}\ket{j}.$
(1)

We will refer to these two copies of Hilbert spaces as the ket and bra spaces,
respectively.

In this paper, we will introduce a strategy for building up topologically
ordered states in the doubled Hilbert space. In this section, we will first
consider the simplest example – that of the toric code model [17]. More
complicated topological orders, especially nonabelian ones, will be discussed
in later sections. Given a square lattice with qubits living on the
edges/links, the toric code model is defined by the following Hamiltonian

$\displaystyle H_{\rm TC}$
$\displaystyle=-\sum_{v}\raisebox{-19.0pt}{\includegraphics{dia_Z2TCVertex}}-\sum_{f}\raisebox{-19.0pt}{\includegraphics{dia_Z2TCFace}},$
$\displaystyle\equiv-\sum_{v}\mathcal{X}_{v}-\sum_{f}\mathcal{Z}_{f}.$ (2)

Throughout this paper, we use $v$ ($e$, $f$) to label vertices (edges, faces).
$\mathcal{X}_{v}$ ($\mathcal{Z}_{f}$) is the product of all Pauli-$X$
(Pauli-$Z$) operators on the neighboring edges of vertex $v$ (face $f$).
Generalizations of the Hamiltonian to non-square lattices are straightforward.
For simplicity, we can assume the lattice to be periodic, although open
boundary conditions are also not hard to analyze. $H_{\rm TC}$ is a
frustration free Hamiltonian: Ground states of it satisfies each single term,
i.e. $\mathcal{X}_{v}=\mathcal{Z}_{f}=1$ for all $v,f$.

###  II.1 Ground State

Let us now describe our protocol for building a toric code ground state in the
doubled Hilbert space. We again consider a periodic square lattice with qubits
living on edges/links, and take the initial state $\ket{0}$ to be the product
state $\prod_{e}\ket{\uparrow}_{e}$. The initial density matrix is denoted as
$\rho_{0}=\ket{0}\bra{0}$, and the corresponding state in the doubled space is
$\ket{\rho_{0}}=\prod_{e,e^{\prime}}\ket{\uparrow}_{e}\ket{\uparrow}_{e^{\prime}}$.
For each vertex $v$ on the lattice, we define a local quantum channel
$\mathcal{E}_{v}[\cdot]$ which acts on an arbitrary density matrix $\sigma$ as

$\displaystyle\mathcal{E}_{v}[\sigma]=\frac{1}{2}(\sigma+\mathcal{X}_{v}\sigma\mathcal{X}_{v}).$
(3)

Let $\mathcal{E}:=\prod_{v}\mathcal{E}_{v}$ be the composition of all
$\mathcal{E}_{v}$; the order does not matter because the $\mathcal{X}_{v}$
operators all commute. The next step of our protocol is to take $\rho_{0}$
through the quantum channel $\mathcal{E}$. Due to the locality of
$\mathcal{X}_{v}$, this quantum channel can be realized by applying a finite-
depth local unitary circuit on $\ket{0}$ together with some auxiliary spins.

Let $\rho$ be the resulting density matrix, and let $\ket{\rho}$ be the
corresponding pure state in the doubled space. We have

$\displaystyle\ket{\rho}\propto\prod_{v}\left[\frac{1}{2}(1+\mathcal{X}_{v}\otimes\mathcal{X}_{v})\right]\ket{0}.$
(4)

We claim and will now explain that $\ket{\rho}$ has the toric code topological
order. We may regard $\ket{\rho}$ as living on two copies of the original
lattice, as illustrated in Fig. 1. For each edge $e$ (vertex $v$, face $f$) in
the ket space lattice, we denote its counterpart in the bra space lattice by
$\bar{e}$ ($\bar{v}$, $\bar{f}$). With this notation, we observe that
$\ket{\rho}$ has $+1$ eigenvalue under the following mutually commuting
operators/stabilizers:

  * •

$Z_{e}Z_{\bar{e}}$ for all $e$,

  * •

$\mathcal{Z}_{f}$ and $\mathcal{Z}_{\bar{f}}$ for all $f$.

  * •

$\mathcal{X}_{v}\mathcal{X}_{\bar{v}}$ for all $v$.

The first set of stabilizers implies that each pair of edges $e$ and $\bar{e}$
can be identified in the $Z$-basis. Under this identification, the bilayer
lattice reduces to a single layer, and the third set of stabilizers is
equivalent to $\mathcal{X}_{v}$. Combining with the second set of stabilizers,
we see that $\ket{\rho}$ indeed minimizes the toric code Hamiltonian $H_{\rm
TC}$. There is also an alternative perspective: Just imagine gluing each pair
of vertices $v$ and $\bar{v}$ as shown by the dashed line in Fig. 1. We obtain
a lattice where each pair of (new) vertices are connected by two edges. The
stabilizers listed above imply that $\ket{\rho}$ minimizes the toric code
Hamiltonian defined on this new lattice.

Figure 1: Illustraction of the doubled lattice.

###  II.2 Restricted Excitations

The instrinsic properties of density matrices, including Hermiticity,
positivity, and trace unity, constraint possible operations in the doubled
Hilbert space. As a consequence, the allowed excitations above the
aforementioned ground state in the doubled space are restricted. For example,
$\operatorname{Tr}(\rho)=1$ implies that $\ket{\rho}$ has a nontrivial overlap
with the state $\ket{I}:=D^{-1/2}\sum_{i}\ket{i}\ket{i}$, where
$\\{\ket{i}|i=1,2,\cdots,D\\}$ is an orthonormal basis of the $D$-dimensional
conventional Hilbert space. Since $\ket{I}$ satisfies (has $+1$ eigenvalue
under) the first and third sets of stabilizers mentioned previously, i.e.
$Z_{e}Z_{\bar{e}}$ and $\mathcal{X}_{v}\mathcal{X}_{\bar{v}}$, no matter what
operations we do on $\ket{\rho}$, the resulting state has no hope to
completely violate (have $-1$ eigenvalue under) these stabilizers.

###  II.3 Abelian Anyon Braiding

Although the accessible excitations are restricted as explained above, we will
show that nontrivial anyon braiding, a hallmark of topological order, can
still be implemented in the doubled Hilbert space. Here, in the case of toric
code, we describe an example of _abelian_ braiding. Later, using general
quantum doubled models, we will also discuss the more exciting phenomenon of
_nonabelian_ braiding.

Figure 2: (a-c) Definition of the $\mathcal{O}_{X}$, $\mathcal{O}_{Z}$, and
$\mathcal{W}_{X}$ operators used for demonstrating abelian braiding. (d) One
view of the doubled Hilbert space.

Let $\rho$ be the density matrix that has toric code topological order in the
doubled Hilbert space. To describe our protocol, define $\mathcal{O}_{X}$,
$\mathcal{O}_{Z}$ and $\mathcal{W}_{X}$ operators as shown in Fig. 2a-c. Also
define $U=(\mathcal{O}_{X}+\mathcal{O}_{Z})/\sqrt{2}$. We propose comparing
the following two states in the doubled space.

  1. 1.

$(\mathcal{W}_{X}\otimes\mathcal{W}_{X})\ket{\rho}=\ket{\rho}$.

  2. 2.

$(U\otimes U)(\mathcal{W}_{X}\otimes\mathcal{W}_{X})(U\otimes
U)\ket{\rho}=:\ket{\rho^{\prime}}$.

Here we have used the fact that $U$ and $\mathcal{W}_{X}$ are both real
unitary operators. The $U\otimes U$ operator does not respect the
$Z_{e}Z_{\bar{e}}=1$ constraints, hence we can no longer identify the edges
from the two lattice layers. We can still use the alternative picture
mentioned previously: Only identify the vertices from the two layers, and
obtain a lattice of the form in Fig. 2d. $U\otimes U$ creates excitations on
this lattice, and in particular, there is a chance that a charge (a vertex
with $\mathcal{X}=-1$) is generated at $v_{2}$. $\mathcal{W}_{X}\otimes W_{X}$
creates a pair of fluxes (faces with $\mathcal{Z}=-1$) from the vacuum, winds
one flux around $v_{2}$, and annihilates the flux pair back to the vacuum.

To see how $\ket{\rho^{\prime}}$ is related to anyon braiding, we write

$\displaystyle(U\otimes U)\ket{\rho}=\ket{\rho_{+}}+\ket{\rho_{-}},$ (5)
$\displaystyle\ket{\rho_{+}}=\frac{1}{2}(1+\mathcal{O}_{X}\otimes\mathcal{O}_{X})\ket{\rho},$
(6)
$\displaystyle\ket{\rho_{-}}=\frac{1}{2}(\mathcal{O}_{X}\otimes\mathcal{O}_{Z}+\mathcal{O}_{Z}\otimes\mathcal{O}_{X})\ket{\rho},$
(7)

where $(\mathcal{O}_{Z}\otimes\mathcal{O}_{Z})\ket{\rho}=\ket{\rho}$ has been
used. $\ket{\rho_{+}}$ contains no excitation near $v_{2}$, but
$\ket{\rho_{-}}$ contains a charge at $v_{2}$. After braiding a flux, we have

$\displaystyle(\mathcal{W}_{X}\otimes\mathcal{W}_{X})(U\otimes
U)\ket{\rho}=\ket{\rho_{+}}-\ket{\rho_{-}},$ (8)

where the minus sign is the result of nontrivial (semionic) mutual statistics
between a charge and a flux. Finally applying $(U\otimes U)$ again, we find

$\displaystyle\ket{\rho^{\prime}}$ $\displaystyle=(U\otimes
U)(\mathcal{W}_{X}\otimes\mathcal{W}_{X})(U\otimes U)\ket{\rho}$
$\displaystyle=(\mathcal{O}_{X}\otimes\mathcal{O}_{X})\ket{\rho},$ (9)

which contains a pair of separated magnetic fluxes and thus differs from
$\ket{\rho}$.

The difference between $\rho$ and $\rho^{\prime}$ can be detected using the
following method: Prepare an ancillary qubit in the $\ket{\uparrow}$ state,
and then send the whole system
$(\ket{\uparrow}\bra{\uparrow})\otimes\rho^{\prime}$ through the following
quantum channel.

$\displaystyle\mathcal{E}_{f}[\sigma]$
$\displaystyle:=[Z\otimes(1+\mathcal{Z}_{f})/2]\sigma[Z\otimes(1+\mathcal{Z}_{f})/2]$
$\displaystyle+[X\otimes(1-\mathcal{Z}_{f})/2]\sigma[X\otimes(1-\mathcal{Z}_{f})/2],$
(10)

where $f$ is one of the two shaded faces in Fig. 2a. The resulting ancillary
qubit is in the $\ket{\downarrow}$ state. In contrast, if
$(\ket{\uparrow}\bra{\uparrow})\otimes\rho$ is sent through $\mathcal{E}_{f}$,
the ancillary qubit will still be in the $\ket{\uparrow}$ state. This
concludes our example of abelian braiding.

We note that although the operator $U$ has a nonlocal nature, it can be step-
by-step built up using local unitaries. Let $U_{n}:=(X_{1}+Z_{1}Z_{2}\cdots
Z_{n})/\sqrt{2}$ be the $n$-qubit version of this operator. $U_{1}$ is the
single-qubit Hadamard gate. General $U_{n}$ can be constructed by the
recursive relation $U_{n}=\mathrm{CNOT}_{n,n-1}U_{n-1}\mathrm{CNOT}_{n,n-1}$.
Here, $\mathrm{CNOT}_{i,j}:=(-1)^{(1-Z_{i})(1-X_{j})/4}$ is a two-qubit gate
and one can check that under the conjugate of this gate,
$(Z_{j},X_{j})\mapsto(Z_{i}Z_{j},X_{j})$.

##  III General Quantum Double Models

Toric code is the simplest instance of the more general _quantum doubled
models_ [17]. In this section, we will show that ground states of all these
models admit efficient preparation in the doubled Hilbert space. This is in
sharp contrast with the case of conventional Hilbert space, where ground
states of only a subset of these models are known to be efficiently preparable
[3, 4, 6].

There is a quantum double model for each finite group $G$ which is generally
nonabelian. The model can live on an arbitrary lattice on an arbitrary
orientable two-dimensional surface. Physical degrees of freedom, called spins,
live on the edges, and the local Hilbert space of each spin is spanned by the
orthonormal group element basis $\\{\ket{g}|g\in G\\}$. We choose a direction
for each edge. Reversing the direction of a particular edge is equivalent to
the basis change $\ket{z}\mapsto\ket{z^{-1}}$ for the corresponding spin. Let
$v$ be a vertex, and $f$ be an adjacent face, we define the local gauge
transformations $A_{v}(g)$ and magnetic flux operators $B_{(v,f)}(h)$ as
follows.

$\displaystyle
A_{v}(g)\Bigg{|}{\raisebox{-19.0pt}{\includegraphics{dia_VertexState}}}\Bigg{\rangle}=\Bigg{|}{\raisebox{-19.0pt}{\includegraphics{dia_VertexGaugeTrans}}}\Bigg{\rangle},$
(11) $\displaystyle
B_{(v,f)}(h)\Bigg{|}{\raisebox{-19.0pt}{\includegraphics{dia_FaceState}}}\Bigg{\rangle}=\delta_{z_{1}z_{2}z_{3}z_{4},h}\Bigg{|}{\raisebox{-19.0pt}{\includegraphics{dia_FaceState}}}\Bigg{\rangle}.$
(12)

Here we use a tetravalent vertex and a square face as examples, and the
generalizations should be straightforward. We further define two projectors:

$\displaystyle A_{v}:=|G|^{-1}\sum_{g\in G}A_{v}(g),\quad
B_{f}:=B_{(v,f)}(1).$ (13)

Note that $B_{f}$ does not depend on the choice of the adjacent vertex $v$.
The quantum double Hamiltonian is then given by [17]

$\displaystyle H_{\rm QD}=-\sum_{v}A_{v}-\sum_{f}B_{f}.$ (14)

The projectors in $H_{\rm QD}$ all commute with each other, and ground states
of $H_{\rm QD}$ satisfy $A_{v}=B_{f}=1$. In the special case $G={Z}_{2}$, if
we identify $\ket{\uparrow}$ ($\ket{\downarrow}$) with the group element
states $\ket{1}$ ($\ket{-1}$), then $H_{\rm QD}$ is equivalent to $H_{\rm TC}$
in Eq. 2 up to a positive factor and an additive constant.

###  III.1 Ground State

The procedure of preparing a ground state of $H_{\rm QD}$ in the doubled space
is similar to the previous toric code example. Consider an arbitrary
orientable two-dimensional lattice with group element states living on the
edges/links. We take the initial state $\ket{0}$ to be the product state
$\prod_{e}\ket{1}_{e}$, where $1\in G$ denotes the identity group element. For
each vertex $v$ on the lattice, we define a local quantum channel
$\mathcal{E}_{v}[\cdot]$ by

$\displaystyle\mathcal{E}_{v}[\sigma]=\frac{1}{|G|}\sum_{g\in G}A_{v}(g)\sigma
A_{v}^{\dagger}(g).$ (15)

Notice that in the group element basis, each $A_{v}(g)$ is a real unitary
matrix, thus $|G|^{-1}\sum_{g\in G}A_{v}(g)^{\dagger}A_{v}(g)=1$, and the
above map is indeed a valid trace-preserving quantum channel. Let
$\mathcal{E}:=\prod_{v}\mathcal{E}_{v}$; the order does not matter because the
$A_{v}(g)$ operators at different vertices commute with each other. The next
step of our protocol is to take the initial density matrix through
$\mathcal{E}$. Let $\rho$ be the resulting density matrix, and let
$\ket{\rho}$ be the corresponding pure state in the doubled space. We have

$\displaystyle\ket{\rho}\propto\prod_{v}\left[|G|^{-1}\sum_{g}A_{v}(g)\otimes
A_{v}(g)\right]\ket{0}.$ (16)

As before, we may regard $\ket{\rho}$ as living on two copies of the original
lattice; see Fig. 1. We observe that $\ket{\rho}$ has $+1$ eigenvalue under
the following mutually commuting projectors:

  * •

$\sum_{g\in G}\ket{g}_{e}\ket{g}_{\bar{e}}\bra{g}_{e}\bra{g}_{\bar{e}}$ for
all $e$,

  * •

$B_{f}$ (ket space) and $B_{\bar{f}}$ (bra space) for all $f$.

  * •

$|G|^{-1}\sum_{g}A_{v}(g)A_{\bar{v}}(g)$ for all $v$.

As in the previous example, there are two ways of seeing the topological order
of $\ket{\rho}$: (1) Using the first set of projectors, we can identify each
pair of edges $e$ and $\bar{e}$ in the group element basis, thus reducing the
bilayer lattice to a single layer. $\ket{\rho}$ is a ground state of $H_{\rm
QD}$ on this reduced lattice. (2) Just identify each pair of vertices $v$ and
$\bar{v}$ as shown by the dashed line in Fig. 1. $\ket{\rho}$ is also a ground
state of $H_{\rm QD}$ on this lattice with doubled edges.

One may wonder what is the purification of $\rho$. In other words, how to
actually implement the quantum channel $\mathcal{E}$? Let us start from the
local quantum channel $\mathcal{E}_{v}$. One way to implement this channel is
to add an auxiliary group-element spin to the vertex $v$, then prepare this
spin in the state $\ket{+}:=|G|^{-1/2}\sum_{g\in G}\ket{g}$, and finally apply
the following local unitary operator.

$\displaystyle U_{v}:=\sum_{g\in G}\ket{g}_{v}\bra{g}_{v}\otimes A_{v}(g),$
(17)

which acts on both the vertex and the surrounding edges. One may check that
for an arbitrary density matrix $\sigma$ living on the edges,

$\displaystyle\operatorname{Tr}_{v}\left[U(\ket{+}_{v}\bra{+}_{v}\otimes\sigma)U^{\dagger}\right]=\mathcal{E}_{v}[\sigma],$
(18)

therefore realizing the desired local quantum channel. Now, to implement the
total quantum channel $\mathcal{E}$, one just need to add auxiliary spins to
all vertices, all prepared in the $\ket{+}$ state, and then apply the _finite-
depth_ local unitary circuit $U=\prod_{v}U_{v}$. Going back to our state
preparation protocol, this means that $\rho$ can be obtained from the pure
state $\ket{\psi}:=U[(\prod_{v}\ket{+}_{v})\otimes\ket{0}]$ by tracing out all
vertex spins. The state $\ket{\psi}$ is actually an interesting state by
itself. When $G=\mathbb{Z}_{2}$, after an $X$-$Z$ basis change on all edges,
$\ket{\psi}$ is the so-called cluster state, a symmetry protected topological
state with a 0-form and a 1-form $\mathbb{Z}_{2}$ symmetries [19].
Interestingly, the close relation between this cluster state and the toric
code topological order has been showing up repeatedly [3, 20, 21].

###  III.2 Nonabelian Anyon Braiding

Similar to the toric code example, intrinsic properties of density matrices
restrict possible excitations above the aforementioned quantum double ground
states. Nonetheless, we will show that _nonabelian_ anyon braiding can be
implemented in the doubled Hilbert space.

Let us begin with a nontechnical sketch of the nonabelian braiding phenomenon.
In a general quantum double model, we consider a type of unitary string
operators that can create anyon pairs from the vacuum, as illustrated in Fig.
3a. If we move around one anyon from such a pair (Fig. 3b), and finally let it
approach the other anyon, then the two anyons can perfectly annihilate (Fig.
3c), as long as the anyon trajectory does not enclose any other anyon or a
hole of the space. On the other hand, if the trajectory of this anyon pair
encloses another anyon, then it is possible that the anyon pair does not
annihilate, but instead fuses into a new nontrivial anyon. This is illustrated
in Fig. 3d-f. In the following, we elaborate the full details of this
phenomenon and show that it is realizable in the doubled Hilbert space.

Figure 3: Schematic illustration of nonabelian braiding. Circles represent
anyon excitations, and lines represent string operators used for creating
anyons. The lines can be smoothly deformed without changing the underlying
state as long as they do not touch other anyons. Comparing panels a-c and d-f,
we see the effect of nontrivial braiding.

Before going into the doubled space, let us first try to understand an example
of nonabelian braiding in the normal setting. Point-like excitations, or
anyons, of a $2d$ topological order are generated by string-like operators.
For our purpose here, let us consider a class of string operators
$A_{\mathcal{L}}(g)$ whose action is as follows.

$\displaystyle
A_{\mathcal{L}}(g)\bigg{|}\raisebox{-10.0pt}{\includegraphics{dia_CombState}}\bigg{\rangle}$
$\displaystyle=$
$\displaystyle\bigg{|}\raisebox{-10.0pt}{\includegraphics{dia_CombStateTrans}}\bigg{\rangle}.$
(19)

A number of notations need explanations here: As the figure suggests, the
operator $A_{\mathcal{L}}(g)$ acts within a comb-like region. $\mathcal{L}$ is
an oriented string coinciding with the base of the comb (horizontal line
above) and ends at the vertex $v$ (hence its orientation agrees with the
horizontal arrows above). $x_{1,k-1}$ is a short-hand notation for
$x_{1}x_{2}\cdots x_{k-1}$. The number of comb teeth $l$ may not equal to the
string length $k$, depending on the shape of $\mathcal{L}$ and the lattice
structure. We note that in the group element basis, $A_{\mathcal{L}}(g)$ is a
real unitary operator.

A remarkable property of $A_{\mathcal{L}}(g)$ is that it commutes with all
$A_{v}$ and $B_{f}$ terms away from the two ends of $\mathcal{L}$. Hence, when
$\mathcal{L}$ is an open string, the operator creates a pair of point-like
excitations. In the special case where $\mathcal{L}$ forms a closed
contractible loop, the operator $A_{\mathcal{L}}(g)$ acts trivially on all
ground states of the quantum double model – creating no excitation.
Physically, this corresponds to first creating a pair of anyons and then
annihilating them back to the vacuum.

Figure 4: (a) Illustration of nonabelian braiding. The closed loop
$\mathcal{C}$ is assumed to be contractible. (b) A spin configuration used for
the computation of nonabelian braiding.

To demonstrate the phenomenon of nonabelian braiding, consider two string
operators $A_{\mathcal{C}}(g)$ and $A_{\mathcal{L}}(h)$ illustrated in Fig.
4a. where $\mathcal{C}$ is a contractible closed loop and $\mathcal{L}$ is
open. Let $\ket{\Omega}$ be an arbitrary ground state of the quantum double
model, we will show that in general
$A_{\mathcal{C}}(g)A_{\mathcal{L}}(h)\ket{\Omega}\neq
A_{\mathcal{L}}(h)A_{\mathcal{C}}(g)\ket{\Omega}$. Hence there is a nontrivial
braiding. Moreover, in contrast to
$A_{\mathcal{L}}(h)A_{\mathcal{C}}(g)\ket{\Omega}$ that has two anyons at the
ends of $\mathcal{L}$, $A_{\mathcal{C}}(g)A_{\mathcal{L}}(h)\ket{\Omega}$ has
one more anyon located at $v_{1}$, the base point of $\mathcal{C}$. This
phenomenon is not possible for abelian braiding.

Consider a particular spin configuration as shown in Fig. 4b, and assume the
absence of flux, i.e. $B_{p}=1$ for all faces $p$. Let us first analyze the
action of string operators on such a state: Acting $A_{\mathcal{L}}(h)$, the
$y_{l}$ spin is mapped to $\tilde{y}_{l}:=y_{l}x_{1,k}^{-1}h^{-1}x_{1,k}$.
Next, acting $A_{\mathcal{C}}(g)$, we have

$\displaystyle z_{1}$ $\displaystyle\mapsto gz_{1},$ $\displaystyle z_{m}$
$\displaystyle\mapsto(y_{1}\cdots\tilde{y}_{l}\cdots
y_{n-1})^{-1}g(y_{1}\cdots\tilde{y}_{l}\cdots y_{n-1})z_{m}.$ (20)

Before the action of these two operators, the flux around the face $f$
vanishes: $z_{1}wz_{m}^{-1}y_{n}=1$. After the action, the product of group
elements around the same face becomes

$\displaystyle gz_{1}wz_{m}^{-1}(y_{1}\cdots\tilde{y}_{l}\cdots
y_{n-1})^{-1}g^{-1}(y_{1}\cdots\tilde{y}_{l}\cdots y_{n-1})y_{n}$
$\displaystyle=$ $\displaystyle g(y_{1}\cdots\tilde{y}_{l}\cdots
y_{n})^{-1}g^{-1}(y_{1}\cdots\tilde{y}_{l}\cdots y_{n}).$ (21)

Using $y_{1}y_{2}\cdots y_{n}=1$, we observe that
$y_{1}\cdots\tilde{y}_{l}\cdots y_{n}$ is in the conjugacy class of $h^{-1}$,
i.e. it equals to $a^{-1}h^{-1}a$ for some $a\in G$. Define
$\tilde{h}^{-1}:=a^{-1}h^{-1}a$, then $\tilde{h}$ is in the conjugacy class of
$h$, and the result in the previous equation simplifies to
$g\tilde{h}g^{-1}\tilde{h}^{-1}$ – the commutator between $g$ and $\tilde{h}$.
In general, this commutator is nontrivial, leading to a violation of the
$B_{f}$ term in the Hamiltonian. Importantly, for any other nearby face
$f^{\prime}\neq f$, $B_{f^{\prime}}$ is preserved. This implies that the
excitation at the face $f$ is a true anyon that can not be annihilated by
local unitary operators, because it can be detected remotely by multiplying
group elements along a large loop. The ground state $\ket{\Omega}$ is a
superposition of many flux-free spin configurations, where $\tilde{h}$ can
take all possible values among the conjugacy class of $h$. Generically, the
$B_{f}=1$ condition will be partially or fully violated in
$A_{\mathcal{C}}(g)A_{\mathcal{L}}(h)\ket{\Omega}$. As a example, consider the
simplest nonabelian group $S_{3}$. It contains three conjugacy classes:
$C_{1}=\\{1\\}$, $C_{2}$ consisting of three transpositions, and $C_{3}$
consisting of two permutations of order $3$. If we take $g\in C_{2}$ and $h\in
C_{3}$, then $g$ does not commute with any possible $\tilde{h}$, and the
resulting state $A_{\mathcal{C}}(g)A_{\mathcal{L}}(h)\ket{\Omega}$ has
$B_{f}=0$.

Finally, let us explain how to demonstrate nonabelian braiding in the doubled
space. Recall that the string operators $A_{\mathcal{L}}(h)$ are real unitary
operators, so we can directly apply them to the density matrix: $\rho\mapsto
A_{\mathcal{L}}(h)\rho A_{\mathcal{L}}^{\dagger}(h)$, or equivalently,
$\ket{\rho}\mapsto A_{\mathcal{L}}(h)\otimes A_{\mathcal{L}}(h)\ket{\rho}$. We
previously mentioned that in the doubled Hilbert space, there is a picture of
reducing a bilayer lattice to a single layer. Under this reduction,
$A_{\mathcal{L}}(h)\otimes A_{\mathcal{L}}(h)$ is mapped to a single copy of
$A_{\mathcal{L}}(h)$. We have thus understood how to apply the string
operators on the doubled-space topological order. In order to observe
nonabelian braiding, we need to compare the following two states (in the
doubled space)

$\displaystyle\ket{\rho_{1}}$ $\displaystyle=[A_{\mathcal{C}}(g)\otimes
A_{\mathcal{C}}(g)][A_{\mathcal{L}}(h)\otimes A_{\mathcal{L}}(h)]\ket{\rho},$
(22) $\displaystyle\ket{\rho_{2}}$ $\displaystyle=[A_{\mathcal{L}}(h)\otimes
A_{\mathcal{L}}(h)][A_{\mathcal{C}}(g)\otimes A_{\mathcal{C}}(g)]\ket{\rho},$
(23)

and show that $\ket{\rho_{1}}$ contains one more anyon excitaton. A possible
strategy is as follows: Prepare an ancillary qubit in the $\ket{\uparrow}$
state, and then send the full density matrix
$(\ket{\uparrow}\bra{\uparrow})\otimes\rho_{i}$ ($i=1,2$) through the
following quantum channel.

$\displaystyle\mathcal{E}_{f}[\sigma]$ $\displaystyle:=(Z\otimes
B_{f})\sigma(Z\otimes B_{f})$
$\displaystyle+[X\otimes(1-B_{f})]\sigma[X\otimes(1-B_{f})].$ (24)

In the case of $\rho_{2}$, the resulting ancillary qubit will still be in the
$\ket{\uparrow}$ state. On the contrary, in the case of $\rho_{1}$, the
ancillary qubit will have a nontrivial probability of in the
$\ket{\downarrow}$ state.

We note that the string operator $A_{\mathcal{L}}(g)$ defined in Eq. 19 can be
step-by-step built up using local unitaries. To see this, define a unitary
operator $E_{u_{1},u_{2}}$ as follows.

$\displaystyle
E_{u_{1},u_{2}}\bigg{|}\raisebox{-10.0pt}{\includegraphics{dia_ESupport}}\bigg{\rangle}=\bigg{|}\raisebox{-10.0pt}{\includegraphics{dia_ETrans}}\bigg{\rangle}.$
(25)

As the figure already suggests, we choose the two vertices $u_{1}$ and $u_{2}$
such that they sandwich the edge labeled by $x_{k}$ in Eq. 19. Then one can
check that

$\displaystyle
E_{u_{1},u_{2}}^{-1}A_{\mathcal{L}}(g)E_{u_{1},u_{2}}=A_{\mathcal{L}^{\prime}}(g),$
(26)

where $\mathcal{L}^{\prime}$ is longer than $\mathcal{L}$ by one lattice unit.

##  IV Discussion

We have shown that the preparation of many topologically ordered states are
simpler in the doubled Hilbert space than the conventional Hilbert space. We
have also found that nontrivial anyon braiding effects, both abelian and
nonabelian, can be realized in the doubled space. What we are unable to figure
out is whether such doubled-space topological orders are useful for fault-
tolerant quantum computation. The intrinsic restrictions on the accessible
states is a big problem and future works are needed to elucidate whether this
problem sets an essential obstacle for practical applications. Nonetheless, we
hope our work illustrates the potential of quantum simulation in the doubled
Hilbert space and can inspire more interesting results in this direction.

## Acknowledgments

I am grateful to Yanting Cheng, Chengshu Li, and Pengfei Zhang for detailed
suggestions on the writing of this manuscript. I would like to also thank Yimu
Bao and Nat Tantivasadakarn for fruitful discussions. I am supported by the
Gordon and Betty Moore Foundation under Grant No. GBMF8690 and the National
Science Foundation under Grant No. NSF PHY-1748958. This work was initiated at
Aspen Center for Physics, which is supported by National Science Foundation
grant PHY-2210452. Part of this work was conducted at the Institute for
Advanced Study at Tsinghua University. I would like to express my gratitude
for the hospitality.

## References

  * Altman  _et al._ [2021] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, M. Greiner, K. R. Hazzard, R. G. Hulet, A. J. Kollár, B. L. Lev, M. D. Lukin, R. Ma, X. Mi, S. Misra, C. Monroe, K. Murch, Z. Nazario, K.-K. Ni, A. C. Potter, P. Roushan, M. Saffman, M. Schleier-Smith, I. Siddiqi, R. Simmonds, M. Singh, I. Spielman, K. Temme, D. S. Weiss, J. Vučković, V. Vuletić, J. Ye, and M. Zwierlein, Quantum simulators: Architectures and opportunities, PRX Quantum 2, 017003 (2021).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature (London) 607, 667 (2022).
  * Tantivasadakarn  _et al._ [2021] N. Tantivasadakarn, R. Thorngren, A. Vishwanath, and R. Verresen, Long-range entanglement from measuring symmetry-protected topological phases, arXiv e-prints , arXiv:2112.01519 (2021), arXiv:2112.01519 [cond-mat.str-el] .
  * Verresen  _et al._ [2021a] R. Verresen, N. Tantivasadakarn, and A. Vishwanath, Efficiently preparing Schrödinger’s cat, fractons and non-Abelian topological order in quantum devices, arXiv e-prints , arXiv:2112.03061 (2021a), arXiv:2112.03061 [quant-ph] .
  * Tantivasadakarn  _et al._ [2023a] N. Tantivasadakarn, R. Verresen, and A. Vishwanath, Shortest Route to Non-Abelian Topological Order on a Quantum Processor, Phys. Rev. Lett. 131, 060405 (2023a), arXiv:2209.03964 [quant-ph] .
  * Tantivasadakarn  _et al._ [2023b] N. Tantivasadakarn, A. Vishwanath, and R. Verresen, Hierarchy of Topological Order From Finite-Depth Unitaries, Measurement, and Feedforward, PRX Quantum 4, 020339 (2023b), arXiv:2209.06202 [quant-ph] .
  * Lee  _et al._ [2022] J. Y. Lee, W. Ji, Z. Bi, and M. P. A. Fisher, Decoding Measurement-Prepared Quantum Phases and Transitions: from Ising model to gauge theory, and beyond, arXiv e-prints , arXiv:2208.11699 (2022), arXiv:2208.11699 [cond-mat.str-el] .
  * Satzinger  _et al._ [2021] K. J. Satzinger, Y. J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, C. Gidney, I. Aleiner, F. Arute, K. Arya, J. Atalaya, _et al._ , Realizing topologically ordered states on a quantum processor, Science 374, 1237 (2021), arXiv:2104.01180 [quant-ph] .
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021), arXiv:2104.04119 [quant-ph] .
  * Google Quantum AI and Collaborators [2023] Google Quantum AI and Collaborators, Non-Abelian braiding of graph vertices in a superconducting processor, Nature (London) 618, 264 (2023).
  * Iqbal  _et al._ [2023a] M. Iqbal, N. Tantivasadakarn, T. M. Gatterman, J. A. Gerber, K. Gilmore, D. Gresh, A. Hankin, N. Hewitt, C. V. Horst, M. Matheny, T. Mengle, B. Neyenhuis, A. Vishwanath, M. Foss-Feig, R. Verresen, and H. Dreyer, Topological Order from Measurements and Feed-Forward on a Trapped Ion Quantum Computer, arXiv e-prints , arXiv:2302.01917 (2023a), arXiv:2302.01917 [quant-ph] .
  * Iqbal  _et al._ [2023b] M. Iqbal, N. Tantivasadakarn, R. Verresen, S. L. Campbell, J. M. Dreiling, C. Figgatt, J. P. Gaebler, J. Johansen, M. Mills, S. A. Moses, J. M. Pino, A. Ransford, M. Rowe, P. Siegfried, R. P. Stutz, M. Foss-Feig, A. Vishwanath, and H. Dreyer, Creation of Non-Abelian Topological Order and Anyons on a Trapped-Ion Processor, arXiv e-prints , arXiv:2305.03766 (2023b), arXiv:2305.03766 [quant-ph] .
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of Rydberg atoms on a kagome lattice, Proceedings of the National Academy of Science 118, e2015785118 (2021), arXiv:2011.12295 [cond-mat.quant-gas] .
  * Verresen  _et al._ [2021b] R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of Toric Code Topological Order from Rydberg Blockade, Physical Review X 11, 031005 (2021b), arXiv:2011.12310 [cond-mat.str-el] .
  * Bao  _et al._ [2023] Y. Bao, R. Fan, A. Vishwanath, and E. Altman, Mixed-state topological order and the errorfield double formulation of decoherence-induced transitions, arXiv e-prints , arXiv:2301.05687 (2023), arXiv:2301.05687 [quant-ph] .
  * Fan  _et al._ [2023] R. Fan, Y. Bao, E. Altman, and A. Vishwanath, Diagnostics of mixed-state topological order and breakdown of quantum memory, arXiv e-prints , arXiv:2301.05689 (2023), arXiv:2301.05689 [quant-ph] .
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Annals of Physics 303, 2 (2003), arXiv:quant-ph/9707021 [quant-ph] .
  * Liu [2023] S. Liu, Anyon Quantum Dimensions from an Arbitrary Ground State Wave Function, arXiv e-prints , arXiv:2304.13235 (2023), arXiv:2304.13235 [cond-mat.str-el] .
  * Yoshida [2016] B. Yoshida, Topological phases with generalized global symmetries, Phys. Rev. B 93, 155131 (2016), arXiv:1508.03468 [cond-mat.str-el] .
  * Ashkenazi and Zohar [2022] S. Ashkenazi and E. Zohar, Duality as a feasible physical transformation for quantum simulation, Phys. Rev. A 105, 022431 (2022).
  * Liu and Ji [2022] S. Liu and W. Ji, Towards Non-Invertible Anomalies from Generalized Ising Models, arXiv e-prints , arXiv:2208.09101 (2022), arXiv:2208.09101 [cond-mat.str-el] .