# Scaling Properties of Liquid Dynamics Predicted from a Single Configuration:
Pseudoisomorphs for Harmonic-Bonded Molecules

Zahraa Sheydaafar  Jeppe C. Dyre  Thomas B. Schrøder<EMAIL_ADDRESS>Glass and
Time, IMFUFA, Department of Science and Environment, Roskilde University, P.O.
Box 260, DK-4000 Roskilde, Denmark

(August 31, 2024)

###### Abstract

Isomorphs are curves in the thermodynamic phase diagram of invariant excess
entropy, structure, and dynamics, while pseudoisomorphs are curves of
invariant structure and dynamics, but not of the excess entropy. The latter
curves have been shown to exist in molecular models with flexible bonds [Olsen
et al., J. Chem. Phys. 145, 241103 (2016)]. We here present three force-based
methods to trace out pseudoisomorphs based on a single configuration and test
them on the asymmetric dumbbell and 10-bead Lennard-Jones chain models with
bonds modeled as harmonic springs. The three methods are based on requiring
that particle forces, center-of-mass forces, and torques, respectively, are
invariant in reduced units. For each of the two investigated models we
identify a method that works well for tracing out pseudoisomorphs.

##  I Introduction

Defining the excess entropy ${S}_{\rm ex}$ as the entropy minus that of an
ideal gas at same density and temperature, isomorphs are curves in the
thermodynamic phase diagram along which structure, dynamics, and ${S}_{\rm
ex}$ are invariant to a good approximation Gnan  _et al._ (2009). Systems with
isomorphs, termed R-simple Malins  _et al._ (2013); Flenner  _et al._ (2014);
Prasad and Chakravarty (2014); Schrøder and Dyre (2014); Heyes  _et al._
(2015); Khrapak  _et al._ (2016); Kaskosz  _et al._ (2023), are characterized
by a strong correlation between the canonical-ensemble equilibrium
fluctuations of potential energy $U$ and virial $W$ as quantified by the
Pearson correlation coefficient (sharp brackets denote canonical averages and
$\Delta$ the deviation from the mean),

$R=\dfrac{\langle\Delta W\Delta U\rangle}{\sqrt{\langle(\Delta
W)^{2}\rangle\langle(\Delta U)^{2}\rangle}}\,.$ (1)

The criterion for being R-simple is $R>0.9$ Gnan  _et al._ (2009), but even
systems with somewhat lower R values may have good isomorphs Ingebrigtsen  _et
al._ (2012).

Isomorph invariance of structure and dynamics refers to the use of units where
the energy unit is $e_{0}\equiv k_{B}T$, the length unit is
$l_{0}\equiv\rho^{-1/3}$ in which $\rho$ is the particle number density, and
the time unit is $t_{0}\equiv\rho^{-1/3}\sqrt{m/k_{B}T}$ in which $m$ is a
particle mass. Quantities made dimensionless by reference to this unit system
are termed “reduced” and marked by a tilde; for instance the reduced particle
position $\mathbb{r}$ is given by
${\tilde{\mathbb{r}}}\equiv\mathbb{r}/l_{0}=\rho^{1/3}\mathbb{r}$.

Because ${S}_{\rm ex}$ is the part of entropy that refers to particle
positions, isomorphs are configurational adibats. These can be traced out in
the thermodynamic phase diagram by using the generally valid statistical-
mechanical relation Gnan  _et al._ (2009):

$\gamma\equiv\left(\frac{\partial\ln
T}{\partial\ln\rho}\right)_{\mathrm{S_{ex}}}=\frac{\langle\Delta U\Delta
W\rangle}{\langle(\Delta U)^{2}\rangle}\,.$ (2)

Evaluating the right-hand side by equilibrium $NVT$ simulations, an isomorph
is traced out by solving the differential equation Eq. (2) numerically using,
e.g., the Euler or Runge-Kutta methods Attia  _et al._ (2021).

Figure 1: Scatter plots of equilibrium fluctuations of the virial and the
potential energy. (a) Results for the asymmetric dumbbell model with harmonic
bonds and (b) for the 10-bead flexible Lennard-Jones chain model with harmonic
bonds. For both models the correlation coefficient $R$ is much smaller than
the criterion for an R-simple system, $R>0.9$ Gnan  _et al._ (2009).

Isomorphs have been identified and validated for both atomic Gnan  _et al._
(2009); Bøhling  _et al._ (2012); Bacher  _et al._ (2018); Heyes  _et al._
(2019); Attia  _et al._ (2021); Castello  _et al._ (2021) and molecular
systems Ingebrigtsen  _et al._ (2012); Veldhorst  _et al._ (2014), and
isomorph-theory predictions have also been verified in experiments on glass-
forming van der Waals molecular liquids Xiao  _et al._ (2015); Hansen  _et
al._ (2018). In molecular systems, isomorphs are found when bonds are modeled
as constraints Ingebrigtsen  _et al._ (2012); Veldhorst  _et al._ (2014), but
not when the bonds are flexible Veldhorst  _et al._ (2015); Olsen  _et al._
(2016); Kaskosz  _et al._ (2023). As an example, Fig. 1 shows scatter plots of
virial versus potential energy for the asymmetric dumbbell (ASD) and 10-bead
Lennard-Jones chain (LJC) models, both with harmonic springs as commonly used
in simulations Koperwas and Paluch (2022); Kaskosz  _et al._ (2023). Neither
model is R-simple; the virial potential-energy correlation coefficients are
$0.579$ and $0.284$, respectively. As expected, these models do not have
isomorphs (data not shown), i.e., structure and dynamics are not invariant
along the curves of constant ${S}_{\rm ex}$. Nevertheless, it has been found
that both models have curves in the phase diagram of invariant structure and
dynamics; termed “pseudoisomorphs” Veldhorst  _et al._ (2014); Kaskosz  _et
al._ (2023).

Figure 2: Invariance of the dynamics along state points determined by Olsen et
al. Olsen  _et al._ (2016) to define pseudoisomorphs. (a), (b), (c) show
results for the asymmetric dumbbell model with harmonic bonds; (d), (e), (f)
show results for the 10-bead Lennard-Jones chain model with harmonic bonds.
(a) and (d) report the reduced mean-square displacement of the center of mass
as a function of reduced time. (b) and (e) report the center-of-mass
incoherent intermediate scattering function as a function of reduced time at
the wave-vector corresponding to the first peak of the static structure
factor. (c) and (f) report the normalized end-to-end vector time-
autocorrelation function, probing the decay of molecular orientation.

Since pseudoisomorphs are not configurational adiabats, Eq. (2) cannot be used
to identify such curves in the phase diagram. In 2016 Olsen et al. presented a
method for tracing out pseudoisomorphs involving the following steps Olsen
_et al._ (2016): 1) An equilibrium configuration is quenched to the nearest
local minimum in the high-dimensional potential energy landscape, the so-
called inherent state Stillinger and Weber (1983); 2) the Hessian matrix is
diagonalized to find the vibrational spectrum of the inherent state; 3) the
high-frequency part of the spectrum (related to the springs) is identified; 4)
the scaling properties upon a density change of the remaining part of the
spectrum is used to identify the pseudoisomorph. This method is very
computationally demanding, but works well. Olsen et al. demonstrated that the
it traces out pseudoisomorphs for the ASD and the flexible LJC models
demonstrating good, though not perfect invariance of the dynamics probed via
the incoherent intermediate scattering function. For comparison with later
results, Fig. 2 shows the dynamics probed via the mean-square displacement,
incoherent intermediate scattering function, and orientational time-
autocorrelation function, for the state points identified in Olsen  _et al._
(2016).

The present paper proposes simpler methods for generating pseudoisomorphs
based on the scaling properties of the forces of a single configuration
Schrøder (2022). This idea has been shown to work well not only for atomic
systems like the Kob-Andersen binary Lennard-Jones mixture, but also for
molecular models like the ASD model and the Lewis-Wahnstrom OTP model with
rigid bonds Ingebrigtsen  _et al._ (2012); Koperwas  _et al._ (2020); Schrøder
(2022); Sheydaafar  _et al._ (2023).

##  II Models and Simulation details

The ASD model is a toy model of toluene Vrabec  _et al._ (2001); Galbraith and
Hall (2007); Schrøder  _et al._ (2009a); Chopra  _et al._ (2010a, b);
Ingebrigtsen  _et al._ (2012); Fragiadakis and Roland (2017); Santiago (2018);
Dombrowski and Klotsa (2020). We simulated a system of 5000 ASD molecules
defined as two different-sized spheres, a large (A) and a small (B) one Vrabec
_et al._ (2001); Schrøder  _et al._ (2009b); Milinkovic  _et al._ (2013). The
spheres interact via Lennard-Jones (LJ) potentials, and in the units defined
by the large sphere ($\sigma_{AA}\equiv 1$, $\epsilon_{AA}\equiv 1$, and
$m_{A}\equiv 1$) the other LJ-parameters are $\sigma_{AB}=0.894$,
$\sigma_{BB}=0.788$, $\epsilon_{AB}=0.342$, $\epsilon_{BB}=0.117$,
$m_{B}=0.195$. Bonds are modeled as harmonic springs with equilibrium length
$0.584$ and spring constant $k=3000$.

The LJC model is a generic coarse-grained polymer model Bennemann  _et al._
(1998); Aichele  _et al._ (2003); Puosi and Leporini (2011); Shavit  _et al._
(2013). We simulated 1000 10-bead LJC molecules. Non-bonded particles interact
via a standard LJ potential, cutting and shifting the forces at $2.5\sigma$.
All particles are of same type and the potential parameters are set to unity:
$\sigma=1,\epsilon=1$. Bonds are modeled as harmonic springs with equilibrium
length $1$ and spring constant $k=3000$.

All Molecular Dynamics simulations were performed in the $NVT$ ensemble with a
Nose-Hoover thermostat in periodic boundary conditions, using RUMD that is an
open-source molecular dynamics package optimized for GPU computing Bailey  _et
al._ (2017) (http://rumd.org).

##  III Identifying pseudoisomorphs via force-based methods

Single-configuration force-based methods for generating isomorphs were
introduced recently Schrøder (2022); Sheydaafar  _et al._ (2023). The idea is
the following. Given a configuration $\mathbf{R}_{1}$ at state point
$(\rho_{1},T_{1})$, a uniform scaling to density $\rho_{2}$ is performed
leading to
$\mathbf{R}_{2}=\left(\rho_{1}/\rho_{2}\right)^{1/3}\mathbf{R}_{1}$. For
molecules, two variants of this scaling can be applied: “center-of-mass
scaling” where the molecular center-of-masses are scaled while all
orientations and internal degrees of freedom are kept fixed, or “atomic
scaling” where a uniform scaling is applied to all atoms thus modifying also
the intramolecular bond lengths. After scaling a configuration, the forces
associated with the two configurations, $\mathbf{F}(\mathbf{R}_{1})$ and
$\mathbf{F}(\mathbf{R}_{2})$, are compared ($\mathbf{F}$ is the long vector of
all forces). The temperature $T_{2}$ at the density $\rho_{2}$ is identified
from the condition of invariant reduced forces. Specifically,
$|\mathbf{\tilde{F}}(\mathbf{R}_{1})|=|\mathbf{\tilde{F}}(\mathbf{R}_{2})|$
implies Schrøder (2022)

$T_{2}=\frac{|\mathbf{F}(\mathbf{R}_{2})|}{|\mathbf{F}(\mathbf{R}_{1})|}\left(\frac{\rho_{1}}{\rho_{2}}\right)^{1/3}T_{1}\,.$
(3)

If the two force vectors are parallel, this leads to the forces being
identical in reduced units,
$\mathbf{\tilde{F}}(\mathbf{R}_{1})=\mathbf{\tilde{F}}(\mathbf{R}_{2})\,$.
Assuming this is representative of all relevant configurations, it follows
that structure and dynamics are invariant in reduced units because the same
reduced-unit equation of motion applies at the two state points Schrøder and
Dyre (2014); Dyre (2018). We do not expect the force vectors before and
scaling to be perfectly parallel, however, but Eq. (3) can still be used.
Reference Schrøder (2022) proposed that for force-based method to work well,
the standard Pearson as well as Spearman (ordering) correlation coefficients
of the force components should both be larger than 0.95.

Different variants of the method are arrived at by different interpretations
of what exactly $\mathbf{F}(\mathbf{R})$ represents, e.g., the forces on all
the atoms or just on the center-of-mass forces. We consider below also a
variant based on invariance of torques in reduced units,
$\boldsymbol{\tilde{\tau}_{1}}=\boldsymbol{\tilde{\tau}_{2}}$, which leads to

$T_{2}=\frac{\boldsymbol{|\tau_{2}|}}{\boldsymbol{|\tau_{1}|}}T_{1}\,.$ (4)

Figure 3:  Force and torque correlations of a single configuration of the
harmonic ASD model using center-of-mass scaling. Each point represents scaled
and unscaled forces/torques. (a) Atomic-force method. The figure shows the
x-coordinates of the reduced particle forces plotted against the same
quantities of the uniformly scaled configuration. The temperature
$T_{2}=0.197$ is identified by means of Eq. (3). (b) Molecular-force method.
The figure shows the center-of-mass “molecular” forces before and after
scaling, which have no contributions from the springs. In this case
$T_{2}=0.299$. (c) Torque method. The figure shows the correlation between the
torques of the molecules before and after scaling ($T_{2}=0.310$).

The three methods are applied to the ASD model in Fig. 3, using center-of-mass
scaling, i.e., $\mathbb{R}_{2,cm}=(\rho_{1}/\rho_{2})^{1/3}\mathbb{R}_{1,cm}$.
The state point $(\rho_{1},T_{1})=(0.785,0.174)$ is used as reference and
$\rho_{2}=0.856$, i.e., a 9% density increase is considered. Figure 3(a) shows
a scatter plot of the atomic-force components before and after scaling for a
single configuration; here $T_{2}=0.197$ is found by applying atomic forces in
Eq. (3). Figure 3(b) shows a similar plot based on the center-of-mass
“molecular” forces. In this case the quite different $T_{2}=0.299$ is arrived
at and a significantly better correlation is obtained. Finally, Fig. 3(c)
shows the torque correlations of molecules of unscaled and scaled
configurations. The correlation is here not far from that of the molecular
forces; using Eq. (4) gives $T_{2}=0.310$.

Figure 4:  Testing for invariance of the reduced dynamics of the harmonic-
spring ASD model using the three methods generating pseudoisomorphs by scaling
a single configuration of the reference state point
$(\rho_{1},T_{1})=(0.785,0.174)$. Each method investigates the reduced center-
of-mass mean-square displacement (upper figures), center-of-mass incoherent
intermediate scattering function at the wave vectors corresponding to the
first peak of the static structure factor (middle figures), and directional
time-autocorrelation function probed via the autocorrelation function of the
normalized bond vectors (bottom figures). (a), (b), (c) show results of the
atomic-force method requiring invariant reduced forces between all atoms,
i.e., including the harmonic bond contributions into Eq. (3). (d), (e), (f)
show results of the molecular-force method requiring invariant reduced center-
of-mass forces between the molecules. (g), (h), (i) show results of the torque
method requiring invariant reduced torques on the molecules (Eq. (4)).

For each of the three methods, Fig. 4 compares results for the dynamics. While
the methods in principle require only a single configuration, for carefully
comparing the methods, $T_{2}$-values were obtained by averaging over 195
independent configurations. We find that the molecular-force method gives
almost invariant dynamics, except at the lowest density (black curve) where
negative pressure and phase separation is observed. The torque method gives
similar, though slightly inferior results, while the atomic force method does
not work well.

Next we turn to the 10-bead harmonic-spring LJC model for which one can use
not only the three above methods, but also a fourth one based on invariant
reduced segmental forces. These are defined as

$\mathbf{F}_{Seg,j}\equiv\dfrac{1}{d_{j}}\mathbf{F}_{j}+\dfrac{1}{d_{j+1}}\mathbf{F}_{j+1},$
(5)

in which $d_{j}$ and $d_{j+1}$ are the number of bonds that particles $j$ and
$j+1$ are involved in. Clearly

$\mathbf{F}_{Mol}=\sum_{j=1}^{9}\mathbf{F}_{Seg,j}\,.$ (6)

The results are quite different from the ASD results. Thus the dynamics are to
a good approximation invariant for the atomic- and segmental-force methods
(Fig. 5 (g), (h), (i)), except at the lowest density from reference point
$(\rho_{1},T_{1})=(1.00,0.700)$ at which the virial becomes negative. On the
other hand, neither the molecular force method (d, e, f) nor the torque method
(j, k, l) produce good pseudoisomorphs.

Figure 5: Dynamics of the harmonic-spring LJC model along proposed
pseudoisomorphs. The reference state point is $(\rho_{1},T_{1})=(1.00,0.700)$
and the density increase is $17\%$. The atomic-force and molecular-force
methods both generate good pseudoisomorphs when increasing density as compared
to the reference state point. For both methods, decreasing the density leads
to a higher temperature, which obviously leads to too fast dynamics. The
molecular-force method does not achieve a good pseudoisomorph. The torque
method results in better invariance than the molecular-force method for a
range of densities, but worse than the atomic-force and segmental-force
methods.  Figure 6: Dynamics of the ASD model along proposed pseudoisomorphs.
The densities are here higher than in Fig. 4. The density increase is $19\%$.
For all three methods center-of-mass scaling was employed. The reference state
point $(\rho_{1},T_{1})=(0.932,0.465)$ is taken from Ref. Olsen  _et al._
(2016). The molecular-force methods works best, but the invariance is not as
good as at lower densities (Fig. 4).

##  IV Pseudoisomorphs in the ASD Model at Higher Densities

In Fig. 4 we applied the three force-based single-configuration methods to the
ASD model and found best invariance of the dynamics using the molecular-force
and, to a slightly less degree, torque methods. Figure 6 shows the results of
applying the three methods at higher densities. The invariance of the dynamics
is not as good as at lower densities (Fig. 4). The molecular-force method
works best, but no method works as well as the method of Ref. Olsen  _et al._
, 2016, compare Fig. 2.

###  IV.1 What is the problem?

[]fig7a []fig7b
---

Figure 7: (a) Equilibrium distribution of ASD bond lengths at state points
identified by the molecular-force method (see Figure 6(d,e,f)). Bonds are
compressed when density is increased. This means that when a configuration is
center-of-mass scaled from the reference state point to a (say) higher
density, the resulting configuration is not representative of the equilibrium
configurations at this higher density. (b) Quenching the system according to
the constraint conditions (fixed centre-of-mass and orientational direction of
all molecules). Black: bond length distribution of original unscaled
configuration at the reference state point. Orange: applying the quenching
procedure to the unscaled configuration leads to virtually the same
distribution. Red: atomic scaling to the new higher density leads to bond
lengths that are considerably smaller than the ones in equilibrium at the new
state point (green curve). Blue: center-of-mass scaling to the new higher
density followed by the quenching procedure leads to a bond length
distribution that is close to the one in equilibrium at the new state point
(green curve). Note: the procedure adds the same length ’l’ to all bonds,
which lead to shift in the distribution leaving the shape invariant.

Figure 7(a) shows the distributions of bond lengths for simulations at the
state points generated by the molecular-force method. The bonds are compressed
when the density is increased, an effect that is not seen at the lower
densities of Fig. 4. This means that when a configuration from the reference
state point is scaled to a higher density, the scaled configuration is not
representative of equilibrium configurations at the new state point because
the bonds are too long, an issue that becomes more serious at higher
densities.

To eliminate the effects of the harmonic bonds, we introduce the following
procedure: For a given configuration fix the center-of-mass and orientation of
each molecule. For this “constrained” system add a length $l$ to all bond
lengths and minimize the potential energy with respect to $l$. We refer to
this procedure as “quenching”, but note that the minimization only involves a
single variable, the length $l$ added to all bonds. When applied to an
unscaled configuration, the bond length distribution remains largely unchanged
(black and orange curves in Fig. 7 (b)). Now, consider the scaling of the
original configuration to higher density, with the aim to achieve a bond
length distribution close to that of the equilibrium at the new state point
(green curve). Using center-of-mass scaling of the molecules leaves the bond
lengths unchanged (black curve), whereas atomic scaling lead to too small bond
lengths (red curve). Applying the quenching procedure after center-of-mass
scaling shifts the distribution to shorter bond lengths (from the black line
to the blue line), and the average bond length approaches that of the
equilibrium distribution at the higher density (green line). This effectively
eliminates the non-scaling degrees of freedom thought to cause the poor
invariance observed in Fig. 6.

[]fig8a []fig8b
---
[]fig8c []fig8d

Figure 8: Correlation of atomic-force components, (a) and (b), and molecular-
force components, (c) and (d), for the ASD model. A configuration is taken
from a simulation at $(\rho_{1},T_{1})=(0.932,0.465)$ and $T_{2}$ as
previously calculated by means of Eq. (3). For both the atomic-force and
molecular-force methods we find significantly better correlations between
scaled and non-scaled forces components when applying a quench. Note that only
the molecular forces fulfill the criterion that the Pearson and Spearman
correlation coefficients are both above 0.95 Schrøder (2022).

Based on pairs of quenched configurations, one may again apply the atomic-
force, molecular-force, and torque methods to generate state points with,
hopefully, same dynamics. Specifically, we proceeded as follows. A single
configuration is selected from an equilibrium simulation at the reference
state point (density $\rho_{1}$). This configuration is scaled uniformly to
the density of interest, $\rho_{2}$. Both scaled and unscaled configurations
were then quenched as described above in order to eliminate the bond
vibrational degrees of freedom. After this the relevant forces / torques were
evaluated and $T_{2}$ determined from Eq. (3) and Eq. (4), respectively.
Figure 8(a) shows the forces of a single scaled configuration versus those of
the unscaled configuration before quenching, while (b) demonstrates better
correlation after quenching; (c) and (d) show the correlations between the
center-of-mass forces before and after quenching. The quench method leads to
significantly better correlation, with correlation coefficient increasing from
0.850 to 0.934 in the atomic-force case and from 0.975 to 0.990 in the
molecular (center-of-mass) force case.

Figure 9:  Testing the ASD model for invariance of the reduced dynamics using
the three quench methods for generating pseudoisomorphs based on a single
configuration from $(\rho_{1},T_{1})=(0.932,0.465)$. In contrast to Fig. 6,
the scaled and non-scaled configurations were here quenched to a potential-
energy minimum in order to eliminate the harmonic bond degrees of freedom;
after which $T_{2}$ was determined as in Fig. 6.  Figure 10:  Testing the ASD
model for invariance of the reduced-unit radial distribution function for the
three different methods, without and with quenching, using
$(\rho_{1},T_{1})=(0.932,0.465)$ as reference state point. (a), (b), (c) show
results for state points generated by the atomic-force, molecular-force, and
torque methods; (d), (e), (f) show results for state points generated by the
same methods after minimization. In all cases, structure is invariant to a
good approximation.

The effect of the quenching procedure is tested in Fig. 9, which is analogous
to Fig. 6 except that the temperatures $T_{2}$ are calculated based on
quenched configurations. The best results are obtained with the molecular-
force method, which was also the one that worked best in Fig. 4. For this
method we find excellent collapse of the reduced center-of-mass mean-square
displacement as a function of time, as well as of the center-of-mass
incoherent intermediate scattering function. The directional time-
autocorrelation function shows a slightly worse collapse, but is nevertheless
significantly better than without quenching. Comparing the results of the
torque method with and without quenching shows that quenching also here
significantly improves the invariances. – Using atomic scaling gives the same
results for the molecular force and torque methods after the system has been
quenched (Table 1), while the atomic-force method gives results different from
Fig. 9.

Table 1: Predicted temperature using atomic scaling after quenching (ASD model). Note that the molecular-force and torque methods give the same results as when center-of-mass scaling is used (Fig. 9). $\rho=0.932$ is the reference density in all cases.  density | $T(F_{Atomic})$ | $T(F_{Mol})$ | Torque
---|---|---|---
0.886 | 0.444 | 0.352 | 0.345
0.932 | 0.465 | 0.465 | 0.465
0.969 | 0.494 | 0.573 | 0.581
1.009 | 0.540 | 0.710 | 0.730
1.060 | 0.624 | 0.917 | 0.957

##  V Summary

Genuine isomorphs do not exist in systems without strong virial potential-
energy correlations, e.g., systems of molecules with harmonic bonds. For the
ASD model and the flexible LJC model with harmonic bonds curves nevertheless
exist along which the structure and dynamics are invariant to a good
approximation, so-called pseudoisomorphs Olsen  _et al._ (2016). These curves
do not have invariant excess entropy, but otherwise behave much like isomorphs
by having approximately invariant structure and dynamics in reduced units. A
previously discussed method for tracing out pseudoisomorphs in the
thermodynamic phase diagram works well, but is quite complicated to apply in
practiceOlsen  _et al._ (2016). The present paper explored the possibility of
tracing out pseudoisomorphs based on the simple requirement of invariant
length of the reduced force vector of a single configuration after scaling.
Although the focus above was on the dynamics, we note that for all three
methods discussed the structure is also invariant to a good approximation,
both with and without quenching (Fig. 10).

For the harmonic-bond ASD model the best method for tracing out
pseudoisomorphs is the molecular-force method with center-of-mass scaling and
quenching of the harmonic bonds (Fig. 9). At low densities, the quenching can
be skipped (Fig. 4). We conjecture that this method will work on all small
molecules with harmonic bonds, assuming of course that the systems in question
have pseudoisomorphs. The molecular-force method does not work well for the
10-bead harmonic-bond LJ-chain model, however, for which the best method is
the atomic-force method with atomic scaling (Fig. 5).

From a practical point of view one would like to have a single method that
works for all molecular models with pseudoisomorphs, including the two models
investigated here. The differences found can be interpreted as reflecting
differences regarding which forces are important for the motion of the
molecules. Despite the harmonic bonds, the ASD is relatively stiff and its
motion is largely governed by the center-of-mass “molecular” force and the
torque. If the scaling were perfect, both the molecular-force and the torque
methods would lead to invariant dynamics. The fact that we find the molecular-
force method to work better than the torque method can be interpreted as the
molecular forces being more important than the torques for the motion of the
molecules. It would be interesting to investigate whether molecular models
exist for which the torque method works better than the molecular-force
method. The molecular-force method does not work well for the LJC model. We
interpret this as reflecting the fact that the center-of-mass motion does not
determine the motion of this complex, flexible molecule, which should rather
be be thought of in terms of the motion of its beads.

In summary, simple single-configuration force-based methods exist for tracing
out pseudoisomorphs, which work as well as the complicated method of
diagonalizing the Hessian before and after scaling.Olsen  _et al._ (2016) No
single-configuration method of general validity has been identified, however.
More work is needed to investigate these and other methods. In this regard, we
find it encouraging that the introduction of quenches after scaling improves
all the methods studied.

###### Acknowledgements.

This work was supported by the VILLUM Foundation’s Matter grant (VIL16515).

## References

  * Gnan  _et al._ (2009) N. Gnan, T. B. Schrøder, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, Pressure-energy correlations in liquids. IV. “Isomorphs” in liquid phase diagrams, J. Chem. Phys. 131, 234504 (2009).
  * Malins  _et al._ (2013) A. Malins, J. Eggers, and C. P. Royall, Investigating isomorphs with the topological cluster classification, J. Chem. Phys. 139, 234505 (2013).
  * Flenner  _et al._ (2014) E. Flenner, H. Staley, and G. Szamel, Universal Features of Dynamic Heterogeneity in Supercooled Liquids, Phys. Rev. Lett. 112, 097801 (2014).
  * Prasad and Chakravarty (2014) S. Prasad and C. Chakravarty, Onset of simple liquid behaviour in modified water models, J. Chem. Phys. 140, 164501 (2014).
  * Schrøder and Dyre (2014) T. B. Schrøder and J. C. Dyre, Simplicity of condensed matter at its core: Generic definition of a Roskilde-simple system, J. Chem. Phys. 141, 204502 (2014).
  * Heyes  _et al._ (2015) D. M. Heyes, D. Dini, and A. C. Branka, Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and other properties along fluid-solid coexistence, Phys. Status Solidi (b) 252, 1514 (2015).
  * Khrapak  _et al._ (2016) S. A. Khrapak, B. Klumov, L. Couedel, and H. M. Thomas, On the long-waves dispersion in Yukawa systems, Phys. Plasmas 23, 023702 (2016).
  * Kaskosz  _et al._ (2023) F. Kaskosz, K. Koperwas, A. Grzybowski, and M. Paluch, The origin of the density scaling exponent for polyatomic molecules and the estimation of its value from the liquid structure, J. Chem. Phys. 158, 144503 (2023).
  * Ingebrigtsen  _et al._ (2012) T. S. Ingebrigtsen, T. B. Schrøder, and J. C. Dyre, Isomorphs in model molecular liquids, J. Phys. Chem. B 116, 1018 (2012).
  * Attia  _et al._ (2021) E. Attia, J. C. Dyre, and U. R. Pedersen, Extreme case of density scaling: The Weeks-Chandler-Andersen system at low temperatures, Phys. Rev. E 103, 062140 (2021).
  * Bøhling  _et al._ (2012) L. Bøhling, T. S. Ingebrigtsen, A. Grzybowski, M. Paluch, J. C. Dyre, and T. B. Schrøder, Scaling of viscous dynamics in simple liquids: Theory, simulation and experiment, New J. Phys. 14, 113035 (2012).
  * Bacher  _et al._ (2018) A. K. Bacher, T. B. Schrøder, and J. C. Dyre, The EXP pair-potential system. II. Fluid phase isomorphs, J. Chem. Phys. 149, 114502 (2018).
  * Heyes  _et al._ (2019) D. M. Heyes, D. Dini, L. Costigliola, and J. C. Dyre, Transport coefficients of the Lennard-Jones fluid close to the freezing line, J. Chem. Phys. 151, 204502 (2019).
  * Castello  _et al._ (2021) F. L. Castello, P. Tolias, and J. C. Dyre, Testing the isomorph invariance of the bridge functions of Yukawa one-component plasmas, J. Chem. Phys. 154, 034501 (2021).
  * Veldhorst  _et al._ (2014) A. A. Veldhorst, J. C. Dyre, and T. B. Schrøder, Scaling of the dynamics of flexible Lennard-Jones chains, J. Chem. Phys. 141, 054904 (2014).
  * Xiao  _et al._ (2015) W. Xiao, J. Tofteskov, T. V. Christensen, J. C. Dyre, and K. Niss, Isomorph theory prediction for the dielectric loss variation along an isochrone, J. Non-Cryst. Solids 407, 190 (2015).
  * Hansen  _et al._ (2018) H. W. Hansen, A. Sanz, K. Adrjanowicz, B. Frick, and K. Niss, Evidence of a one-dimensional thermodynamic phase diagram for simple glass-formers, Nat. Commun. 9, 518 (2018).
  * Veldhorst  _et al._ (2015) A. A. Veldhorst, J. C. Dyre, and T. B. Schrøder, Scaling of the dynamics of flexible Lennard-Jones chains: Effects of harmonic bonds, J. Chem. Phys. 143, 194503 (2015).
  * Olsen  _et al._ (2016) A. E. Olsen, J. C. Dyre, and T. B. Schrøder, Communication: Pseudoisomorphs in liquids with intramolecular degrees of freedom, J. Chem. Phys. 145, 241103 (2016).
  * Koperwas and Paluch (2022) K. Koperwas and M. Paluch, Computational evidence for the crucial role of dipole cross-correlations in polar glass-forming liquids, Phys. Rev. Lett. 129, 025501 (2022).
  * Stillinger and Weber (1983) F. H. Stillinger and T. A. Weber, Dynamics of structural transitions in liquids, Phys. Rev. A 28, 2408 (1983).
  * Schrøder (2022) T. B. Schrøder, Predicting scaling properties from a single fluid configuration, Phys. Rev. Lett. 129, 245501 (2022).
  * Koperwas  _et al._ (2020) K. Koperwas, A. Grzybowski, and M. Paluch, Virial–potential-energy correlation and its relation to density scaling for quasireal model systems, Phys. Rev. E 102, 062140 (2020).
  * Sheydaafar  _et al._ (2023) Z. Sheydaafar, J. C. Dyre, and T. B. Schrøder, Scaling properties of liquid dynamics predicted from a single configuration: Small rigid molecules, J. Phys. Chem. B 127, 3478 (2023).
  * Vrabec  _et al._ (2001) J. Vrabec, J. Stoll, and H. Hasse, A set of molecular models for symmetric quadrupolar flu, J. Phys. Chem. B 105, 12126 (2001).
  * Galbraith and Hall (2007) A. L. Galbraith and C. K. Hall, Solid–liquid phase equilibria for mixtures containing diatomic Lennard–Jones molecules, Fluid Phase Eq. 262, 1 (2007).
  * Schrøder  _et al._ (2009a) T. B. Schrøder, N. P. Bailey, U. R. Pedersen, N. Gnan, and J. C. Dyre, Pressure-energy correlations in liquids. III. Statistical mechanics and thermodynamics of liquids with hidden scale invariance, J. Chem. Phys. 131, 234503 (2009a).
  * Chopra  _et al._ (2010a) R. Chopra, T. M. Truskett, and J. R. Errington, Excess-entropy scaling of dynamics for a confined fluid of dumbbell-shaped particles, Phys. Rev. E 82, 041201 (2010a).
  * Chopra  _et al._ (2010b) R. Chopra, T. M. Truskett, and J. R. Errington, Excess entropy scaling of dynamic quantities for fluids of dumbbell-shaped particles, J. Chem. Phys. 133, 104506 (2010b).
  * Fragiadakis and Roland (2017) D. Fragiadakis and C. M. Roland, A test for the existence of isomorphs in glass-forming materials, J. Chem. Phys. 147, 084508 (2017).
  * Santiago (2018) I. Santiago, Nanoscale active matter matters: Challenges and opportunities for self-propelled nanomotors, Nano Today 19, 11 (2018).
  * Dombrowski and Klotsa (2020) T. Dombrowski and D. Klotsa, Kinematics of a simple reciprocal model swimmer at intermediate Reynolds numbers, Phys. Rev. Fluids 5, 063103 (2020).
  * Schrøder  _et al._ (2009b) T. B. Schrøder, U. R. Pedersen, N. P. Bailey, S. Toxvaerd, and J. C. Dyre, Hidden Scale Invariance in Molecular van der Waals Liquids: A Simulation Study, Phys. Rev. E 80, 041502 (2009b).
  * Milinkovic  _et al._ (2013) K. Milinkovic, M. Dennison, and M. Dijkstra, Phase diagram of hard asymmetric dumbbell particles, Phys. Rev. E 87, 032128 (2013).
  * Bennemann  _et al._ (1998) C. Bennemann, W. Paul, K. Binder, and B. Dünweg, Molecular-dynamics simulations of the thermal glass transition in polymer melts: $\alpha$-relaxation behavior, Phys. Rev. E 57, 843 (1998).
  * Aichele  _et al._ (2003) M. Aichele, Y. Gebremichael, F. W. Starr, J. Baschnagel, and S. C. Glotzer, Polymer-specific effects of bulk relaxation and stringlike correlated motion in the dynamics of a supercooled polymer melt, J. Chem. Phys. 119, 5290 (2003).
  * Puosi and Leporini (2011) F. Puosi and D. Leporini, Scaling between relaxation, transport, and caged dynamics in polymers: From cage restructuring to diffusion, J. Phys. Chem. B 115, 14046 (2011).
  * Shavit  _et al._ (2013) A. Shavit, J. F. Douglas, and R. A. Riggleman, Evolution of collective motion in a model glass-forming liquid during physical aging, J. Chem. Phys. 138, 12A528 (2013).
  * Bailey  _et al._ (2017) N. P. Bailey, T. S. Ingebrigtsen, J. S. Hansen, A. A. Veldhorst, L. Bøhling, C. A. Lemarchand, A. E. Olsen, A. K. Bacher, L. Costigliola, U. R. Pedersen, H. Larsen, J. C. Dyre, and T. B. Schrøder, RUMD: A general purpose molecular dynamics package optimized to utilize GPU hardware down to a few thousand particles, Scipost Phys. 3, 038 (2017).
  * Dyre (2018) J. C. Dyre, Perspective: Excess-entropy scaling, J. Chem. Phys. 149, 210901 (2018).


# Clifford Dressed Time-Dependent Variational Principle

Antonio Francesco Mello  International School for Advanced Studies (SISSA),
via Bonomea 265, 34136 Trieste, Italy  Alessandro Santini  International
School for Advanced Studies (SISSA), via Bonomea 265, 34136 Trieste, Italy
Guglielmo Lami  Laboratoire de Physique Théorique et Modélisation, CNRS UMR
8089, CY Cergy Paris Université, 95302 Cergy-Pontoise Cedex, France  Jacopo De
Nardis  Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089, CY
Cergy Paris Université, 95302 Cergy-Pontoise Cedex, France  Mario Collura
International School for Advanced Studies (SISSA), via Bonomea 265, 34136
Trieste, Italy  INFN Sezione di Trieste, 34136 Trieste, Italy

###### Abstract

We propose an enhanced Time-Dependent Variational Principle (TDVP) algorithm
for Matrix Product States (MPS) that integrates Clifford disentangling
techniques to efficiently manage entanglement growth. By leveraging the
Clifford group, which maps Pauli strings to other Pauli strings while
maintaining low computational complexity, we introduce a Clifford dressed
single-site 1-TDVP scheme. During the TDVP integration, we apply a global
Clifford transformation as needed to reduce entanglement by iteratively
sweeping over two-qubit Clifford unitaries that connect neighboring sites in a
checkerboard pattern. We validate the new algorithm numerically using various
quantum many-body models, including both integrable and non-integrable
systems. Our results demonstrate that the Clifford dressed TDVP significantly
improves entanglement management and computational efficiency, achieving
higher accuracy, extended simulation times, and enhanced precision in computed
observables compared to standard TDVP approaches. Additionally, we propose
incorporating Clifford gates directly within the two-site 2-TDVP scheme.

## Introduction. —

Simulating non-equilibrium dynamics of closed quantum many-body systems is a
fundamental challenge in modern physics, whose difficulty primarily arises
from the exponential size of the underlying Hilbert space and from the complex
quantum correlations involved. The dynamics of these systems, where
entanglement is generally expected to grow linearly in time (following a time-
dependent volume law scaling of the entanglement entropy), demand powerful
numerical methods to be addressed [1, 2, 3, 4, 5]. Matrix Product States (MPS)
provide a robust numerical framework for simulating one-dimensional (1D)
quantum many-body systems, offering a compact representation of 1D quantum
states [6, 7, 8, 9].

A standard approach for simulating the non-equilibrium dynamics of an MPS wave
function involves the Trotterization of the evolution operator, thus
transforming the Hamiltonian dynamics into a discrete quantum circuit of
unitary gates acting locally between neighboring lattice sites (qubits) in the
system. This method results in the well-known Time-Evolving Block Decimation
(TEBD) scheme [6], which efficiently handles real-time evolution by breaking
down the global evolution into a sequence of local unitary gates.

Unlike TEBD, the Time-Dependent Variational Principle (TDVP) method evolves
the MPS parameters directly within the MPS variational manifold [10, 11]. TDVP
ensures that the quantum state remains an optimal approximation within the MPS
ansatz throughout the evolution, providing a more accurate representation of
the system’s dynamics, especially for long-time simulations and systems with
strong entanglement growth. This approach is crucial for studying phenomena
such as quantum quenches [12], transport properties [13], and the growth of
entanglement [14, 15], where traditional methods would struggle with the
exponentially growing Hilbert space.

Unfortunately, regardless of the chosen approach, a quantum system in a non-
equilibrium state experiences a scrambling of quantum information, manifested
as entanglement spreading throughout the system [16, 17]. This process
inevitably increases the complexity of the state. In the MPS representation,
this results in the unbounded growth of the bond dimension. Consequently,
beyond a certain threshold time, this exponential growth leads to a loss of
accuracy in our approximate representation of the wave function, as the MPS
can no longer efficiently capture the system’s entanglement.

However, if we can manipulate the wave function to strategically reduce some
of the quantum correlations and keep the entanglement level under control,
this would be a promising strategy to extend the capabilities of current
classical algorithms for simulating the unitary dynamics of quantum systems
[18]. Indeed, a remarkable insight comes from the field of quantum computing
and information: circuits composed exclusively of Clifford gates (Hadamard, S
and Controlled-NOT gates) can be efficiently simulated on a classical
computer.

Quantum states generated by Cliffords, known as stabilizer states, can exhibit
volume law entanglement yet remain classically tractable (Gottesman Knill
theorem) [19, 20, 21]. The algorithm allowing this is known as tableau
algorithm, as it consists in storing and evolving efficiently a compact table
representing the Pauli strings that stabilize the state [22, 23]. This facts
highlight that while entanglement is a fundamental quantum resource, its
presence alone does not necessarily make a computational problem intractable
for classical devices.

Beyond entanglement, another essential quantum resource is ‘non-
stabilizerness’ or ‘magic’, which significantly influences the complexity of
quantum problems [24]. Similar to entanglement, non-stabilizerness has been
rigorously defined within the paradigm of quantum resource theories, with
various measures proposed to quantify it [25, 26, 27, 28, 29, 30, 31, 32, 33].
These measures have become instrumental in understanding the computational
power of quantum systems, emphasizing that non-stabilizerness is a key factor
in determining the classical hardness of simulating quantum systems [34, 35,
36, 37, 35, 38, 39, 40, 41, 42, 43].

Successfully integrating MPS with stabilizer formalism [23, 44, 45] poses
significant challenges. Recently, initial progress in this direction has been
made with the introduction of an efficient algorithm for finding the
stabilizer Pauli strings of an MPS [40, 46], the stabilizer tensor network
ansatz [47], and subsequently the Clifford enhanced MPS ($\mathcal{C}$MPS)
[48]. The latter are states obtained by applying a Clifford unitary to an MPS.
An optimization algorithm for these states has been proposed in Ref. [49] in
case of ground state search, effectively augmenting the Density Matrix
Renormalization Group (DMRG) with Clifford circuits. Hybrid Clifford tensor
networks algorithms have been introduced also in Refs. [50] and [51].

In this letter, we advance the research in this direction by exploring a novel
approach that integrates a Clifford-based disentanglement scheme with the TDVP
algorithm to reduce the complexity of simulating out-of-equilibrium quantum
systems.

## Preliminaries. —

We examine a system of $N$ qubits, although our findings can be easily
extended to qudits. We designate the local computational basis
$\\{\ket{0},\ket{1}\\}$ as the eigenstates of Pauli matrix $\hat{\sigma}^{3}$,
such that $\hat{\sigma}^{3}\ket{s}=(-1)^{s}\ket{s}$. A general Pauli string
can be formed as
$\hat{\Sigma}^{\boldsymbol{\mu}}=\hat{\sigma}_{1}^{\mu_{1}}\hat{\sigma}_{2}^{\mu_{2}}\cdots\hat{\sigma}_{N}^{\mu_{N}}\in\mathcal{P}_{N}$,
where the boldface superscript $\boldsymbol{\mu}$ represents the set
$\\{\mu_{1},\dots,\mu_{N}\\}$, with $\mu_{j}\in\\{0,1,2,3\\}$ and zero stays
for the local identity operator. In the following, for simplicity, subscripts
indicating the qubits will be omitted whenever they are not required. Pauli
strings form a complete basis for any operator $\hat{O}$ acting on the many-
body Hilbert space $\mathcal{H}=\\{\ket{0},\ket{1}\\}^{\otimes N}$.
Specifically, we can express $\hat{O}$ as
$\hat{O}=\sum_{\boldsymbol{\mu}}O_{\boldsymbol{\mu}}\hat{\Sigma}^{\boldsymbol{\mu}}$,
where the coefficients are determined by $O_{\boldsymbol{\mu}}={\rm
Tr}(\hat{O}\hat{\Sigma}^{\boldsymbol{\mu}})/2^{N}$. This uses the
orthogonality condition ${\rm
Tr}(\hat{\Sigma}^{\boldsymbol{\mu}}\hat{\Sigma}^{\boldsymbol{\nu}})=2^{N}\delta_{\boldsymbol{\mu}\boldsymbol{\nu}}$.

When an operator undergoes a generic unitary transformation, i.e.
$\hat{U}\hat{O}\hat{U}^{{\dagger}}$, its complexity in the Pauli basis
representation can grow arbitrarily. However, when this transformation is
chosen from the set of Clifford unitaries $\mathcal{C}_{N}$, the complexity of
that operator remains bounded. Indeed, the Clifford group $\mathcal{C}_{N}$ on
$N$ qubits, consists of unitary operators that map Pauli strings into Pauli
strings under conjugation. Specifically, for any Clifford unitary
$\hat{C}\in\mathcal{C}_{N}$ and any Pauli operator
$\hat{\Sigma}^{\boldsymbol{\mu}}\in\mathcal{P}_{N}$, we have
$\hat{C}\hat{\Sigma}^{\boldsymbol{\mu}}\hat{C}^{\dagger}=\pm\hat{\Sigma}^{\boldsymbol{\nu}}$.
Clifford transformations, which are generated by the Hadamard gate, phase
gate, and CNOT gate, are fundamental to fault-tolerant quantum computation and
quantum error correction [52, 19, 20, 21, 53].

Let us stress that the CNOT gate, which is the only responsible for the
entanglement generation, is a member of this group of operators. In practice,
any unitary matrix $\hat{C}\in\mathcal{C}_{N}$ can produce significant
entanglement when applied to a many-body state, yet it preserves the
intrinsically low complexity associated with Pauli string operators. In turn,
the operator entanglement of a typical local Hamiltonian will be bounded by
$\sim\log(N)$.

Figure 1: Clifford dressed evolution. (a) The evolution of the state
$\ket{\phi_{0}}$, driven by the successive application of $U_{1}$ and $U_{2}$,
is transformed according to a modified evolution path
$\tilde{U}_{2}\tilde{U}_{1}$. In this new trajectory, each intermediate state
is strategically adjusted within its respective Clifford orbit to
significantly reduce the entanglement in its MPS representation. This approach
ensures that the MPS wave-function maintains a lower entanglement throughout
the evolution, enhancing the efficiency and accuracy of the simulation. (b)
Clifford dressed TDVP steps, the disentangling routine optimally selects local
Clifford two-qubit gates by minimizing the bipartite entanglement entropy,
then sweeps along the chain up to a number of Clifford bi-layers of depth $D$.
(c) At the end of the disentangling procedure, we obtain the Clifford dressed
state $|\tilde{\psi}(t_{m})\rangle$. The optimal Clifford transformation
$\hat{C}_{m}$ is encoded using the stabilizer tableau and applied to the
diagonal MPO Hamiltonian $\hat{H}_{m-1}$, preserving its diagonal structure.
(d) The effective operators for the next TDVP step are iteratively constructed
from the MPS of the Clifford dressed state in mixed canonical form, and the
MPO of the Clifford dressed Hamiltonian.

## Clifford disentangled 1-TDVP. —

The typical scenario we consider involves the dynamics generated by the
Hamiltonian $H_{0}$ after preparing the system in a short-range correlated
initial state (e.g., a product state) at time $t=0$. The subscript zero
indicates that the Hamiltonian is initially in its “bare” form, meaning it has
not yet been modified by Clifford disentanglers.

We are generally interested in evolving the state and monitoring local
observables which are experimentally relevant. This includes computing the
expectation value of Pauli strings
$\bra{\psi(0)}e^{i\hat{H}_{0}t}\hat{\Sigma}^{\boldsymbol{\mu}}e^{-i\hat{H}_{0}t}\ket{\psi(0)}$.
Typically, we track the time evolution of the state by discretizing time into
small intervals $dt$. The state, approximated as an MPS with a fixed bond
dimension $\chi$, is evolved using the single-site 1-TDVP scheme [11].

However, the entanglement entropy typically grows unboundedly, rendering the
MPS representation an inadequate wave-function ansatz beyond a relatively
short time threshold. This is where Clifford-inspired disentangling strategies
come into play.

Basically, at each step $m\in\\{1,\dots\\}$, corresponding to time
$t_{m}=m\,dt$, we can apply a suitable Clifford transformation $\hat{C}_{m}$
to the state $\ket{\psi(t_{m})}$ to reduce entanglement of the new state
$|\tilde{\psi}(t_{m})\rangle=\hat{C}_{m}\ket{\psi(t_{m})}$.

This disentangling routine, or “entanglement cooling” [54, 48, 55, 56, 57],
applied to the states induces a corresponding transformation of the
Hamiltonian: $\hat{H}_{m}=\hat{C}_{m}\hat{H}_{m-1}\hat{C}_{m}^{\dagger}$. The
latter can be efficiently performed using the stabilizer tableau formalism
[58].

In practice, we can rewrite the original discrete evolution
$\prod_{m=0}^{\lfloor t/dt\rfloor-1}e^{-i\hat{H}_{0}dt}\ket{\psi(0)}$ as a
Clifford dressed evolution (see Fig. 1(a)):

$\prod_{m=0}^{\lfloor
t/dt\rfloor-1}e^{-i\hat{H}_{m}dt}\hat{C}_{m}\ket{\psi(0)},$ (1)

where $\hat{C}_{0}=\hat{I}$ is the identity operator, and the product is
understood to be in reverse order. Eq. (1) iteratively constructs the final
$\mathcal{C}$MPS which results in $\hat{C}_{\lfloor
t/dt\rfloor-1}\cdots\hat{C}_{1}\,e^{-i\hat{H}_{0}t}\ket{\psi(0)}$ [48].
Additionally, at each time step $t_{m}$, any Pauli string we are interested in
measuring must be transformed as
$\hat{C}_{m}\cdots\hat{C}_{1}\hat{\Sigma}^{\boldsymbol{\mu}}\hat{C}_{1}^{\dagger}\cdots\hat{C}_{m}^{\dagger}$.

Note that we are assuming the application of a Clifford disentangler at each
time step of the TDVP integrator. However, this is not strictly necessary. We
can choose to apply the entanglement cooling routine less frequently, for
instance, every $k$ time steps, thereby reducing the computational overhead
while still mitigating the growth of entanglement.

The disentangling routine constructs the optimal Clifford operator iteratively
by sweeping over two-qubit Clifford unitaries that connect neighboring sites
of the chain in a checkerboard pattern [54, 48]. The number of bi-layers is
denoted by $D$ (see Fig. 1(b)). The optimal Clifford two-qubit gate is
selected by searching through the subset of 720 Clifford tableaux with
positive sign [58] and randomly choosing one of those that minimizes the von
Neumann entanglement entropy. The sign of the tableau does not affect the
local Singular Value Decomposition (SVD) of the MPS tensors. However, since
the sign choice may influence subsequent minimizations, we reintroduce a
random sign configuration from the $2^{4}$ possible options. The procedure is
repeated sequentially, sweeping back and forth over the entire chain.

We implement the 1-TDVP scheme because it is a symplectic integrator.
Specifically, to evolve the state from $t_{m}$ to $t_{m+1}$, we project the
dressed Hamiltonian $\hat{H}_{m}$ using the MPS tensors of the Clifford
enhanced state
$|\tilde{\psi}(t_{m})\rangle=\mathbb{A}^{s_{1}}_{L}\cdots\mathbb{A}^{s_{n-1}}_{L}\mathbb{A}^{s_{n}}_{C}\mathbb{A}^{s_{n+1}}_{R}\cdots\mathbb{A}^{s_{N}}_{R}\ket{s_{1},\dots,s_{N}},$
which is in the mixed canonical form with respect to the central site $n$.
This allows us to define the effective Hamiltonian (see Fig. 1(c))

$H^{\text{eff}}_{m}(n)=\sum_{k}^{O(N)}J_{k}\,\mathbb{L}^{m}_{k}(n-1)\otimes\hat{\sigma}_{n}^{\mu^{m}_{k}}\otimes\mathbb{R}^{m}_{k}(n+1),$
(2)

which remains diagonal in the operator auxiliary dimension, as we keep
$\hat{H}_{m}$ diagonal in the Pauli basis. Here, $J_{k}$ are the couplings
associated with each Pauli string. A similar transformation is applied for
$\hat{K}^{\text{eff}}_{m}(n)$ (we refer to Ref. [11] for more details on the
TDVP algorithm in the context of MPS).

Figure 2: (a) Average magnetization and (b) half-chain entanglement entropy
for the critical transverse field Ising model. (c) Half-chain magnetization
and (d) half-chain entanglement entropy for the XX model. (e) Average
magnetization and (f) half-chain entanglement entropy for the next-to-nearest-
neighbors Ising chain. Black lines represent free-fermions (upper and central
panels) or ED results (lower panel). Solid green lines are TDVP data. Markers
indicate the Clifford disentangler, applied every $k$ time steps (see legend
in the upper panel). In all panels we set $N=20,\ \chi=128$. Figure 3: Error
on expectation value of observables (same observables as in Fig.2). (a)(c)(e)
Integrated error $\epsilon(t)$, and (b)(d)(f) integrated error at final time
$\epsilon_{T}$ ($T=8$). (a)(b) Critical transverse field Ising model. (c)(d)
XX model. (e)(f) Next-to-nearest-neighbors Ising chain. Solid lines are TDVP
data. Markers indicate the Clifford disentangler, applied every $k=10$ time
steps. We set $N=20$ and explore different bond dimensions $\chi$ (see legend
in the upper panel).

## Numerical experiments. —

We benchmark our algorithm on the following bare Hamiltonian

$\displaystyle\hat{H}_{0}=$ $\displaystyle
J^{x}_{1}\sum_{j=0}^{N-2}\hat{\sigma}^{1}_{j}\hat{\sigma}^{1}_{j+1}+J^{y}_{1}\sum_{j=0}^{N-2}\sigma^{2}_{j}\sigma^{2}_{j+1}+$
$\displaystyle+$ $\displaystyle
J^{x}_{2}\sum_{j=0}^{N-3}\sigma^{1}_{j}\sigma^{1}_{j+2}+h\sum_{j=0}^{N-1}\hat{\sigma}^{3}_{j}.$
(3)

Specifically, after fixing the MPS bond dimension $\chi$, that is the amount
of employable resources, we compare the standard 1-TDVP with our novel
strategy where the Clifford disentangler routine is invoked every $k$ time
steps (for various values of $k$).

Firstly, we consider two integrable cases by setting $J_{2}^{x}=0$ in Eq. (3).
Specifically, the critical quantum Ising chain $J^{y}_{1}=0$, $J^{x}_{1}=-h=1$
prepared in the fully polarized state $\ket{00\dots 0}$ along the $\hat{z}$
direction, and, the XX model with $J^{x}_{1}=J^{y}_{1}=1$ and $h=0$ prepared
in the Néel state $\ket{0101\dots 01}$. Finally, we consider the non-
integrable next-to-nearest-neighbors Ising model by setting
$J_{2}^{x}=J_{1}^{x}=-h=1$ and $J^{y}_{1}=0$ prepared in $\ket{00\dots 0}$.

During the time evolution, we observe the dynamics of the bipartite
entanglement entropy at the midpoint of the chain, which shows discontinuities
each time the Clifford disentangler is applied.

To evaluate the algorithm’s effectiveness in reproducing the expectation
values of local observables (which may become non-local after the Clifford
dressing), we measure the half-chain magnetization $\hat{\sigma}^{3}_{N/2}$
for the XX model. Notice that this model admits a $U(1)$ symmetry, therefore,
the total magnetization is conserved. Instead, for the Ising model, in both
cases we measure the average magnetization
$\hat{M}_{z}/N=\sum_{j=0}^{N-1}\hat{\sigma}^{3}_{j}/N$. These quantities,
which are not conserved, experience a nontrivial evolution. The XX and
critical nearest-neighbors Ising models are exactly solvable, allowing us to
compare our numerical findings with the exact solutions obtained via free-
fermion techniques. For the next-to-nearest-neighbors Ising model, instead, we
resort to comparison with exact diagonalization (ED) calculations.

We define a measure of the error up to time $t$ as the integrated distance
between the exact expectation value of an observable $O(t)$ and the
corresponding value obtained with either TDVP or Clifford enhanced TDVP
denoted as $\tilde{O}_{\chi}(t)$ (for fixed $\chi$), i.e.

$\epsilon(t)=\frac{1}{t}\int_{0}^{t}\absolutevalue{\tilde{O}_{\chi}(t^{\prime})-O(t^{\prime})}\,dt^{\prime}\,.$
(4)

We also define the integrated error at final time as
$\epsilon_{T}=\epsilon(T)$. As mentioned, in our case $O$ corresponds to
either the average magnetization or the on site half-chain magnetization.

Fig. 2 shows the results for the observable evolution and the half-chain
entanglement entropy in all analyzed models for fixed bond dimension
$\chi=128$ and system size $N=20$. The corresponding panels in Fig. 3 show the
evolution of the error $\epsilon(t)$ for $k=10$.

For the critical Ising chain in Figs. 2(a) and 2(b), Our approach limits the
growth of entanglement entropy and enables us to achieve higher precision in
magnetization dynamics compared to the standard 1-TDVP, which becomes
ineffective around $t\approx 5$ due to the entanglement exceeding the bond
dimension $\chi$. The most remarkable results are attained for the XX model,
Figs. 2(c) and 2(d), where at around $t\approx 2$ the 1-TDVP simulation breaks
down. Clifford dressed evolution closely matches the exact solution for
significantly longer times, demonstrating remarkable accuracy. Notice how the
improved evolution obtains several order of magnitudes in precision for the
same amount of resources. Finally, Figs. 2(e) and 2(f) refer to the non-
integrable setup, where disentangling the evolved state becomes more
challenging than in the previous scenarios. Nevertheless, the evolution of the
Clifford dressed observable seems to match the one provided by ED for
$\chi\leq 128$. The integrated error, contrarily to the other case studies,
does not allow to claim a significant advantage.

## Conclusion and Outlook. —

In this work, we introduce a novel algorithm to enhance the simulation of
quantum time-evolution. We improve the standard single-site TDVP algorithm for
MPS by incorporating optimized Clifford disentanglers in order to mitigate the
growth of entanglement. The algorithm can be considered as an effective way to
perform time evolution for the recently introduced states known as Clifford
enhanced MPS [48]. Our numerical experiments demonstrate the effectiveness of
the method in managing controlled entanglement growth while accurately
reproducing the expectation values of physical observables.

Several possible future research are apparent. First, improving the way
Clifford disentanglers are optimized could yield further improvements. Various
possibilities can be explored, such as introducing finite temperature in a
Metropolis-like algorithm, or simultaneously optimizing multiple local
Clifford gates. These approaches could help in avoiding sub-optimal local
minima. Secondly, akin to the approach taken in Ref. [49] for DMRG, a
potential extension could leverage the two-site TDVP (2-TDVP). In this
scenario, the two-qubit Clifford disentangler would be applied to disentangle
the state just before performing the SVD. Once the optimal disentangler is
identified, the effective Hamiltonian is mapped under conjugation, affecting
the Pauli strings only locally. This approach would maintain locality while
efficiently managing entanglement, albeit at the cost of sacrificing the
symplectic nature of the integrator.

A fundamental theoretical question that remains unresolved is whether it is
feasible to quantify the effectiveness of the ansatz used, thereby
distinguishing between quantum states that can be successfully simulated using
our method (i.e. those that can be significantly Clifford disentangled) and
those for which this is not possible. This would be crucial to gain a deeper
understanding of the benefits and limitations method. A possibility is a
connection with the entanglement spectrum statistics, which should exhibit
signatures of entanglement irreversibility [56, 57].

## Acknowledgments. –

We acknowledge the use of Stim [58] for the stabilizer formalism operations,
ITensor for the tensor networks simulations [59] and QuTiP for the ED
simulations [60]. We are particularly grateful to Martina Frau, Alessio
Paviglianiti, Poetri Tarabunga, Emanuele Tirrito, Tobias Haug, Gerald Fux,
Marcello Dalmonte and Rosario Fazio for inspiring discussions on the topic.
This work was supported by the ANR-22-CPJ1-0021-01, the ERC Starting Grant
101042293 (HEPIQ), the PNRR MUR project PE0000023-NQSTI, and the PRIN 2022
(2022R35ZBF) - PE2 - “ManyQLowD”.

## References

  * Feynman [1982] R. P. Feynman, Simulating physics with computers, International Journal of Theoretical Physics 21, 467 (1982).
  * 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 607, 667 (2022).
  * Xu  _et al._ [2023] X. Xu, S. Benjamin, J. Sun, X. Yuan, and P. Zhang, A herculean task: Classical simulation of quantum computers (2023), arXiv:2302.08880 [quant-ph] .
  * Pashayan  _et al._ [2020] H. Pashayan, S. D. Bartlett, and D. Gross, From estimation of quantum probabilities to simulation of quantum circuits, Quantum 4, 223 (2020).
  * Anshu and Arunachalam [2023] A. Anshu and S. Arunachalam, A survey on the complexity of learning quantum states (2023), arXiv:2305.20069 [quant-ph] .
  * Vidal [2004] G. Vidal, Efficient Simulation of One-Dimensional Quantum Many-Body Systems, Phys. Rev. Lett. 93, 040502 (2004).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Paeckel  _et al._ [2019a] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019a).
  * Silvi  _et al._ [2019] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Physics Lecture Notes 10.21468/scipostphyslectnotes.8 (2019).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Ranabhat and Collura [2022] N. Ranabhat and M. Collura, Thermalization of long range Ising model in different dynamical regimes: a full counting statistics approach (2022).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium $xxz$ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Alba and Calabrese [2018] V. Alba and P. Calabrese, Entanglement dynamics after quantum quenches in generic integrable systems, SciPost Phys. 4, 017 (2018).
  * Alba and Calabrese [2016] V. Alba and P. Calabrese, Entanglement and thermodynamics after a quantum quench in integrable systems, Proceedings of the National Academy of Sciences 114, 7947 (2016).
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, Journal of Statistical Mechanics: Theory and Experiment 2005, P04010 (2005).
  * Amico  _et al._ [2008] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Entanglement in many-body systems, Rev. Mod. Phys. 80, 517 (2008).
  * Paeckel  _et al._ [2019b] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Time-evolution methods for matrix-product states, Annals of Physics 411, 167998 (2019b).
  * Gottesman [1997] D. Gottesman, Stabilizer Codes and Quantum Error Correction (1997).
  * Gottesman [1998a] D. Gottesman, Theory of fault-tolerant quantum computation, Phys. Rev. A 57, 127 (1998a).
  * Gottesman [1998b] D. Gottesman, The Heisenberg Representation of Quantum Computers, arXiv preprint arXiv:9807006 10.48550/ARXIV.QUANT-PH/9807006 (1998b).
  * Dehaene and De Moor [2003] J. Dehaene and B. De Moor, Clifford group, stabilizer states, and linear and quadratic operations over GF(2), Physical Review A 68, 10.1103/physreva.68.042318 (2003).
  * Aaronson and Gottesman [2004] S. Aaronson and D. Gottesman, Improved simulation of stabilizer circuits, Physical Review A 70, 10.1103/physreva.70.052328 (2004).
  * Gu  _et al._ [2024] A. Gu, S. F. E. Oliviero, and L. Leone, Magic-induced computational separation in entanglement theory (2024), arXiv:2403.19610 [quant-ph] .
  * Gross  _et al._ [2021] D. Gross, S. Nezami, and M. Walter, Schur–weyl duality for the clifford group with applications: Property testing, a robust hudson theorem, and de finetti representations, Communications in Mathematical Physics 385, 1325–1393 (2021).
  * Bu  _et al._ [2023] K. Bu, W. Gu, and A. Jaffe, Stabilizer testing and magic entropy (2023), arXiv:2306.09292 [quant-ph] .
  * Chen  _et al._ [2024] L. Chen, R. J. Garcia, K. Bu, and A. Jaffe, Magic of random matrix product states, Phys. Rev. B 109, 174207 (2024).
  * Rall  _et al._ [2019] P. Rall, D. Liang, J. Cook, and W. Kretschmer, Simulation of qubit quantum circuits via pauli propagation, Phys. Rev. A 99, 062337 (2019).
  * Grewal  _et al._ [2024] S. Grewal, V. Iyer, W. Kretschmer, and D. Liang, Improved stabilizer estimation via bell difference sampling (2024), arXiv:2304.13915 [quant-ph] .
  * Haug and Piroli [2023a] T. Haug and L. Piroli, Quantifying nonstabilizerness of matrix product states, Physical Review B 107, 10.1103/physrevb.107.035148 (2023a).
  * Haug and Kim [2023] T. Haug and M. Kim, Scalable Measures of Magic Resource for Quantum Computers, PRX Quantum 4, 010301 (2023).
  * Leone  _et al._ [2022] L. Leone, S. F. E. Oliviero, and A. Hamma, Stabilizer Renyi Entropy, Phys. Rev. Lett. 128, 050402 (2022).
  * Leone and Bittel [2024] L. Leone and L. Bittel, Stabilizer entropies are monotones for magic-state resource theory (2024), arXiv:2404.11652 [quant-ph] .
  * Niroula  _et al._ [2024] P. Niroula, C. D. White, Q. Wang, S. Johri, D. Zhu, C. Monroe, C. Noel, and M. J. Gullans, Phase transition in magic with random quantum circuits (2024), arXiv:2304.10481 [quant-ph] .
  * Lami and Collura [2023] G. Lami and M. Collura, Nonstabilizerness via Perfect Pauli Sampling of Matrix Product States, Phys. Rev. Lett. 131, 180401 (2023).
  * Mello  _et al._ [2024a] A. F. Mello, G. Lami, and M. Collura, Retrieving non-stabilizerness with Neural Networks (2024a), arXiv:2403.00919 [quant-ph] .
  * Tarabunga  _et al._ [2023] P. S. Tarabunga, E. Tirrito, T. Chanda, and M. Dalmonte, Many-body magic via pauli-markov chains—from criticality to gauge theories, PRX Quantum 4, 040317 (2023).
  * Heinrich and Gross [2019] M. Heinrich and D. Gross, Robustness of Magic and Symmetries of the Stabiliser Polytope, Quantum 3, 132 (2019).
  * Fux  _et al._ [2023] G. E. Fux, E. Tirrito, M. Dalmonte, and R. Fazio, Entanglement-magic separation in hybrid quantum circuits (2023), arXiv:2312.02039 [quant-ph] .
  * Tarabunga  _et al._ [2024] P. S. Tarabunga, E. Tirrito, M. C. Bañuls, and M. Dalmonte, Nonstabilizerness via matrix product states in the pauli basis (2024), arXiv:2401.16498 [quant-ph] .
  * Tirrito  _et al._ [2024] E. Tirrito, P. S. Tarabunga, G. Lami, T. Chanda, L. Leone, S. F. E. Oliviero, M. Dalmonte, M. Collura, and A. Hamma, Quantifying nonstabilizerness through entanglement spectrum flatness, Phys. Rev. A 109, L040401 (2024).
  * Haug and Piroli [2023b] T. Haug and L. Piroli, Stabilizer entropies and nonstabilizerness monotones, Quantum 7, 1092 (2023b).
  * Frau  _et al._ [2024] M. Frau, P. S. Tarabunga, M. Collura, M. Dalmonte, and E. Tirrito, Non-stabilizerness versus entanglement in matrix product states (2024), arXiv:2404.18768 [quant-ph] .
  * García  _et al._ [2017] H. J. García, I. L. Markov, and A. W. Cross, On the Geometry of Stabilizer States, arXiv preprint arXiv:1711.07848 10.48550/ARXIV.1711.07848 (2017).
  * Elliott [2008] M. B. Elliott, Stabilizer states and local realism (2008).
  * Lami and Collura [2024] G. Lami and M. Collura, Learning the stabilizer group of a Matrix Product State (2024), arXiv:2401.16481 [quant-ph] .
  * Masot-Llima and Garcia-Saez [2024] S. Masot-Llima and A. Garcia-Saez, Stabilizer tensor networks: universal quantum simulator on a basis of stabilizer states (2024), arXiv:2403.08724 [quant-ph] .
  * Lami  _et al._ [2024] G. Lami, T. Haug, and J. D. Nardis, Quantum State Designs with Clifford Enhanced Matrix Product States (2024), arXiv:2404.18751 [quant-ph] .
  * Qian  _et al._ [2024] X. Qian, J. Huang, and M. Qin, Augmenting density matrix renormalization group with clifford circuits (2024), arXiv:2405.09217 [cond-mat.str-el] .
  * Mello  _et al._ [2024b] A. F. Mello, A. Santini, and M. Collura, Hybrid stabilizer matrix product operator (2024b), arXiv:2405.06045 [quant-ph] .
  * Paviglianiti  _et al._ [2024] A. Paviglianiti, G. Lami, M. Collura, and A. Silva, Estimating Non-Stabilizerness Dynamics Without Simulating It (2024), arXiv:2405.06054 [quant-ph] .
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Bravyi and Kitaev [2005] S. Bravyi and A. Kitaev, Universal quantum computation with ideal Clifford gates and noisy ancillas, Phys. Rev. A 71, 022316 (2005).
  * True and Hamma [2022] S. True and A. Hamma, Transitions in Entanglement Complexity in Random Circuits, Quantum 6, 818 (2022).
  * Odavic  _et al._ [2023] J. Odavic, G. Torre, N. Mijic, D. Davidovic, F. Franchini, and S. M. Giampaolo, Random unitaries, Robustness, and Complexity of Entanglement, Quantum 7, 1115 (2023).
  * Chamon  _et al._ [2014] C. Chamon, A. Hamma, and E. R. Mucciolo, Emergent Irreversibility and Entanglement Spectrum Statistics, Physical Review Letters 112, 10.1103/physrevlett.112.240501 (2014).
  * Shaffer  _et al._ [2014] D. Shaffer, C. Chamon, A. Hamma, and E. R. Mucciolo, Irreversibility and entanglement spectrum statistics in quantum circuits, Journal of Statistical Mechanics: Theory and Experiment 2014, P12007 (2014).
  * Gidney [2021] C. Gidney, Stim: a fast stabilizer circuit simulator, Quantum 5, 497 (2021).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013), arXiv:1211.6518 [quant-ph] .