††thanks: This manuscript has been authored by UT-Battelle, LLC, under
contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US
government retains and the publisher, by accepting the article for
publication, acknowledges that the US government retains a nonexclusive, paid-
up, irrevocable, worldwide license to publish or reproduce the published form
of this manuscript, or allow others to do so, for US government purposes. DOE
will provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-
public-access-plan).

# String order melting of spin-1 particle chains in superconducting transmons
using optimal control

Paul Kairys<EMAIL_ADDRESS>Bredesen Center for Interdisciplinary
Research and Graduate Education, University of Tennessee, Knoxville, Tennessee
Quantum Science Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee
Travis S. Humble<EMAIL_ADDRESS>Bredesen Center for Interdisciplinary
Research and Graduate Education, University of Tennessee, Knoxville, Tennessee
Quantum Science Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee

(March 2021; March 2021)

###### Abstract

Utilizing optimal control to simulate a model Hamiltonian is an emerging
strategy that leverages the intrinsic physics of a device with digital quantum
simulation methods. Here we evaluate optimal control for probing the non-
equilibrium properties of symmetry-protected topological (SPT) states
simulated with superconducting hardware. Assuming a tunable transmon
architecture, we cast evolution of these SPT states as a series of one- and
two-site pulse optimization problems that are solved in the presence of
leakage constraints. From the generated pulses, we numerical simulate time-
dependent melting of the perturbed SPT string order across a six-site model
with an average state infidelity of $10^{-3}$. The feasibility of these pulses
as well as their efficient application indicate that high-fidelity simulations
of string-order melting are within reach of current quantum computing systems.

##  I Introduction

Understanding the static and dynamical properties of quantum states is of
paramount interest to the physical sciences. One route to enable these studies
uses quantum devices and quantum information processing in protocols known as
quantum simulation georgescu_quantum_2014 . There are several unique
approaches to quantum simulation; from purely digital approaches on universal
quantum computers to purely analog approaches using tailor-made quantum
devices. Recent work suggests that there are also a number of intermediate
paradigms capable of realizing quantum simulation salathe_digital_2015 ;
parra-rodriguez_digital-analog_2020 ; lamata_digital-analog_2018 ;
martin_digital-analog_2020 ; celeri_digital-analog_2021 ; lamata_digital-
analog_2017 ; arrazola_digital-analog_2016 ; galicia_enhanced_2020 ;
babukhin_hybrid_2020 ; kairys_parametrized_2021 .

Quantum simulation based on quantum optimal control (QOC) permits one to take
advantage of the natural device dynamics and Hilbert space while also enabling
the use of digital decomposition methods kairys_parametrized_2021 . In
principle, this permits efficient use of coherent resources within quantum
hardware. One of the most common applications of QOC is to identify the device
controls which realize a desired unitary evolution glaser_training_2015 .
Unitary evolution can then be used to preform state preparation, dynamical
evolution, or even mitigate errors to improve information extraction
glaser_training_2015 .

As a leading case study for quantum simulation, topological physics provide a
route for both scientific discovery and engineering as well as the validation
of quantum devices and simulation protocols. One of the primary drivers of
interest in these exotic phases of matter are their potential use within
quantum information processing as robust quantum memories fowler_surface_2012
; terhal2015quantum . Some of these topological phases are well understood via
both analytical and numerical methods that may be used for validating quantum
simulation while others represent key open problems for the field
terhal2015quantum .

Consider the example of one-dimensional symmetry-protected topological (SPT)
phases, which exhibit interesting static and dynamical properties that are
still under much study senthil_symmetry-protected_2015 . One of these phases,
known as the Haldane phase, is realized in chains of interacting spin-1
particles that can be characterized by phenomena such as edge modes,
degenerate entanglement spectra, and dilute (or hidden) antiferromagnetic
order affleck_rigorous_1987 ; pollmann_entanglement_2010 ;
pollmann_detection_2012 ; calvanese_strinati_destruction_2016 . A number of
studies have shown that the latter phenomena, also called string order, can
undergo a dynamic process called melting that leads to infinitesimally fast
vanishing of string order under symmetry-breaking quantum quenches
calvanese_strinati_destruction_2016 .

Here we test the feasibility of studying string order melting using quantum
simulation based on QOC. Our approach considers the simulation performed on
superconducting quantum devices composed of coupled transmons operated as
three-level (qutrit) systems. Transmon technology has matured significant over
the past decade and such devices are currently available in configurations of
up to $100$ transmons. Notably, these hardware systems are driven by analog
control pulses that make them well suited for quantum simulation
arute_quantum_2019 ; blok2021quantum ; mckay2018qiskit . In modeling the
physics of these superconducting transmon devices, we demonstrate that control
optimization can determine the local unitary evolutions that generate a
symmetry-breaking quantum quench. These results are validated using exact
numerical simulation which confirm that feasibility of simulating string order
melting in superconducting quantum devices.

The remainder is outlined as follows: In Section II, the phenomena of string
order melting is defined and strategies to observe the phenomena are
discussed, based on Ref. calvanese_strinati_destruction_2016 . In Section III,
we outline our quantum simulation protocol of the quench dynamics and string
order measurement. In Section IV, the device architecture is introduced and
aspects of controlling this architecture are discussed. We report a set of
optimal controls for this device architecture in Section V and we validate
that these controls enable the study of string order melting in Section VI.
Finally, we conclude and discuss avenues for future experimental
demonstrations in Section VII.

##  II String order melting

Studying dynamical many-body quantum systems is a classically challenging task
because it requires integration of the Schrödinger equation on a Hilbert space
which grows exponentially with increasing particle number. This difficulty
motivates quantum simulation protocols to be used when studying such phenomena
georgescu_quantum_2014 . In this work we propose a quantum simulation protocol
to study string order melting analogous to the classical simulations performed
by Calvanese et al. in Ref. calvanese_strinati_destruction_2016  via classical
numerical methods.

The simulation of string-order melting requires evolving a quantum state with
string order under symmetry-breaking time evolutions and observing the
dynamics of string order as a function of string length, time, and direction.
One of the best known models with string order is given by the AKLT model,
named for Affleck, Kennedy, Lieb, and Tasaki. The AKLT Hamiltonian is defined
as affleck_rigorous_1987

$\displaystyle\hat{H}_{AKLT}$
$\displaystyle=\sum_{i=1}^{N}P(s_{i}+s_{i+1}=2)=\sum_{i=1}^{N}\mathcal{P}_{i,i+1},$
(1)

where $N$ is the number of spin-1 particles, $s_{i}$ is the total spin of the
particle on site $i$, and $P(s_{i}+s_{i+1}=2)=\mathcal{P}_{i,i+1}$ is the
projector onto the subspace between two particles $i$ and $i+1$ with total
spin equal to two: $s_{i}+s_{i+1}=2$.

It was shown that this Hamiltonian has a unique set of ground states
corresponding to the the mutual eigenvectors with eigenvalue $0$ for all the
projectors $\mathcal{P}_{i,i+1}$. These states are known as valence bond
solids and are naturally defined by projecting pairs of spin-1/2 particles in
a singlet state into the spin-1 triplet subspace formed between singlets
affleck_rigorous_1987 ; affleck1988valence . The ground states of the AKLT
Hamiltonian are contained within a phase called the Haldane phase which is a
symmetry-protected topological phase with string order preserved by certain
symmetries affleck_rigorous_1987 ; calvanese_strinati_destruction_2016 ;
pollmann_entanglement_2010 .

To define string order one first defines an operator

$\displaystyle\hat{\mathcal{O}}^{\alpha}_{k,l}:=\hat{S}_{k}^{\alpha}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi\hat{S}_{n}^{\alpha}}\bigg{]}\hat{S}_{l}^{\alpha}$
(2)

where indices $k,n,l$ are lattice sites and $\hat{S}^{\alpha}_{i}$ is the
spin-1 operator acting on lattice site $i$ which generates rotations around
the $\alpha=x,y,z$ axes. The expectation value of
$\hat{\mathcal{O}}^{\alpha}_{k,l}$ can be used to define an order parameter
known as a string order parameter:

$\displaystyle\mathcal{O}^{\alpha}_{\text{string}}(\psi)=\lim_{|l-k|\rightarrow\infty}\bra{\psi}\hat{\mathcal{O}}^{\alpha}_{k,l}\ket{\psi}.$
(3)

When $\mathcal{O}^{\alpha}_{\text{string}}(\psi)\neq 0$ the system is said to
have string order and this order parameter is used to quantify string order
melting. As defined in Ref. calvanese_strinati_destruction_2016 , string order
melting refers to the the decay of string order in the long-range limit of
$|l-k|\rightarrow\infty$ at infinitesimal times. This phenomena is in stark
contrast to typical Landau theory order parameters which cannot vanish
instantly because of the continuity of time evolution
calvanese_strinati_destruction_2016 . Thus string order melting represents a
fascinating class of dynamical quantum physics in which quantum simulators may
be particularly useful in studying.

Following Ref. calvanese_strinati_destruction_2016 , we consider a Hamiltonian
under which the AKLT model is quenched to be a modified version of the spin-1
XXZ Hamiltonian:

$\displaystyle\hat{H}(\lambda,b)$
$\displaystyle=\sum_{i=1}^{N}\bigg{[}\hat{S}^{x}_{i}\hat{S}^{x}_{i+1}+\hat{S}^{y}_{i}\hat{S}^{y}_{i+1}+\lambda\hat{S}^{z}_{i}\hat{S}^{z}_{i+1}\bigg{]}$
$\displaystyle+b\sum_{i=1}^{N}\hat{S}^{x}_{i}$ (4)

where $\lambda,b$ are two competing energy scales to be parameterized. It is
known that for $b=0$ and $\lambda\leq\lambda_{c}\approx 1.186$ the ground
state of the Hamiltonian is within the Haldane phase, the same phase as the
AKLT ground states ueda2008finite . The string order found within the Haldane
phase is known to be preserved under perturbations invariant to the action of
a symmetry group pollmann_entanglement_2010

$\mathcal{G}=\\{\hat{1},e^{i\pi\sum_{i}\hat{S}^{x}_{i}},e^{i\pi\sum_{i}\hat{S}^{y}_{i}},e^{i\pi\sum_{i}\hat{S}^{z}_{i}}\\}.$
(5)

This symmetry group is a specific representation of the dihedral group $D_{2}$
and any evolution which is invariant under the action of all elements of
$\mathcal{G}$ will preserve string order. The transverse field term
$b\sum_{i=1}^{N}\hat{S}^{x}_{i}$ in Eq. (II) is not invariant under the action
of all elements of $\mathcal{G}$ and therefore will lead to the loss of string
order. However, in this case, the perturbation is invariant under the action
of subsets of $\mathcal{G}$, which therefore preserves string order in the $x$
direction calvanese_strinati_destruction_2016 . The perturbation thus leads to
melting of string order in the $z$ and $y$ directions only. We consider the
implementation of quantum evolution under Eq. (II).

The simulation proceeds by initially preparing a ground state of the AKLT
model on an open chain, evolving the state under Eq. (II) with $\lambda=0.2$
for a maximum time $T=2.5$ (here $\hbar=1$). At discrete time steps spaced by
$\delta t=0.1$, the string order observables
$\hat{\mathcal{O}}^{\alpha}_{k,l}$ are calculated and the dependence on the
length of the operator $l-k$, the direction $\alpha$, and the transverse field
perturbation, $b$, are probed. These parameters are precisely the ones we will
use to demonstrate that string order melting can be observed by quantum
simulation enabled with optimal control however, other types of symmetry-
breaking quenches could be explored in the future pollmann_entanglement_2010 ;
calvanese_strinati_destruction_2016 .

We compare the observed simulation results against a numerical simulation
using exact diagonalization for $N=6$ spin-1 sites. From a physical point of
view, this small size does not allow one to observe the large-$N$ effects
which define string order melting. However, it does provide validation that
the optimal controls found in this work enable quantitative simulation of the
phenomena associated with string order melting.

##  III Quantum simulation of string order melting

We consider generating dynamics of the Hamiltonian in Eq. (II) with the
transverse field strength parameter given as $b=\Delta bn_{x}$, where $\Delta
b$ is the interval of $b$ on which we wish to explore quench dynamics and
$n_{x}$ is an integer that determines the total magnitude of $b$. This re-
parametrization allows us to decompose the global time evolution of the
quantum quench into a product of local evolutions via Trotterization
kairys_parametrized_2021 . In this work, the time evolution operator can be
defined as

$\displaystyle U(\lambda,b,T_{s})$
$\displaystyle=\exp\bigg(-\frac{iT_{s}}{\hbar}H(\lambda,b)\bigg{missing})$
$\displaystyle=\lim_{q\rightarrow\infty}\bigg{(}\prod_{i=1}^{N}U_{XYZ}^{i,i+1}\bigg{(}\frac{T_{s}}{q}\bigg{)}U_{X}\bigg{(}\frac{T_{s}}{q}\bigg{)}^{n_{x}}\bigg{)}^{q},$
(6)

where the $XYZ$ unitary determined by $\lambda$ is given by

$U_{XYZ}^{i,i+1}\bigg{(}\frac{T_{s}}{q}\bigg{)}=\exp\bigg(-\frac{iT_{s}}{q\hbar}\bigg{[}\hat{S}^{x}_{i}\hat{S}^{x}_{i+1}+\hat{S}^{y}_{i}\hat{S}^{y}_{i+1}+\lambda\hat{S}^{z}_{i}\hat{S}^{z}_{i+1}\bigg{]}\bigg{missing})$
(7)

and the $X$ field unitary is given by

$\displaystyle U_{X}\bigg{(}\frac{T_{s}}{q}\bigg{)}$
$\displaystyle=\exp\bigg(-\frac{iT_{s}\Delta
b}{q\hbar}\sum_{i=0}^{N}\hat{S}^{x}_{i}\bigg{missing})$ (8)
$\displaystyle=\prod_{i=0}^{N}\exp\bigg(-\frac{iT_{s}\Delta
b}{q\hbar}\hat{S}^{x}_{i}\bigg{missing}).$ (9)

For the numerical demonstrations below, we consider $\Delta b=0.2$ and
$\lambda=0.2$.

By truncating the limit in Eq. (6) one obtains a $q$th-order approximation to
the global evolution operator. Alternatively, we choose to use another common
definition of the Trotter decomposition order given by a step size
$\tau=T_{s}/q$. In this work we use optimal control methods to determine a set
of device controls which generate the individual, local, Trotter evolutions.
Then, by composing the local optimal controls in sequence, one is able
approximate to the desired global quantum dynamics.

We now show how to measure the string observables. We first introduce our
notation for a system of $N$ spin-1 particles and demonstrate how measurements
of the expectation value of the string-order operator can be evaluated.
Consider the spin-1 operators $\hat{S}^{\alpha}$ with directional component
$\alpha=x,y,z$. These operators obey the angular momentum commutation
relations
$[\hat{S}^{\alpha},\hat{S}^{\beta}]=i\hbar\varepsilon_{\alpha\beta\gamma}\hat{S}^{\gamma}$,
where $\varepsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol. The
eigenvectors of each local operator are given by
$\hat{S}^{\alpha}\ket{s^{\alpha}}=s\ket{s^{\alpha}}$ where $s=0,$ or $\pm$
labels the vectors with eigenvalues $0$ and $\pm 1$, respectively. One can
convert from eigenstates of $\hat{S}^{\beta}$ to eigenstates of
$\hat{S}^{\alpha}$ via the unitary operator $\hat{u}^{\alpha,\beta}$ defined
as:

$\hat{u}^{\alpha,\beta}=\ket{-^{\alpha}}\bra{-^{\beta}}+\ket{0^{\alpha}}\bra{0^{\beta}}+\ket{+^{\alpha}}\bra{+^{\beta}}.$
(10)

A tensor product of eigenstates with the same component $\alpha$ can be
labeled via a string $\mathbf{s}$ as
$\ket{\mathbf{s}^{\alpha}}=\otimes_{i=1}^{N}\ket{s_{i}^{\alpha}}$, where
$s_{i}^{\alpha}$ is the $i$th element of the string $\mathbf{s}^{\alpha}$.
Here $\ket{\mathbf{s}^{\alpha}}$ is a state in the composite Hilbert space of
$N$ spin-1 particles and the set of all states formed by all possible strings
$\mathbf{s}$ form a complete orthonormal basis for the composite Hilbert
space, i.e.,
$\hat{I}=\sum_{\mathbf{s}}\ket{\mathbf{s}^{\alpha}}\bra{\mathbf{s}^{\alpha}}$.
Moreover, one can transform between tensor product basis states using
$\hat{U}^{\alpha,\beta}=\bigotimes_{i=1}^{N}\hat{u}^{\alpha,\beta}$:

$\ket{\mathbf{s}^{\alpha}}=\hat{U}^{\alpha,\beta}\ket{\mathbf{s}^{\beta}}=\bigotimes_{i=1}^{N}\hat{u}^{\alpha,\beta}\ket{\mathbf{s}^{\beta}}$
(11)

We will now consider the expansion of an expectation value of a string order
operator upon a particular direction $\alpha$ for an arbitrary state
$\ket{\psi}$ (See Appendix A for full derivation)

$\displaystyle\bra{\psi}\hat{\mathcal{O}}^{\alpha}_{k,l}\ket{\psi}$
$\displaystyle=\bra{\psi}\bigg{(}\hat{S}_{k}^{\alpha}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi\hat{S}_{n}^{\alpha}}\bigg{]}\hat{S}_{l}^{\alpha}\bigg{)}\ket{\psi}$
(12)
$\displaystyle=\sum_{\mathbf{s}}\bigg{(}s_{k}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi
s_{n}}\bigg{]}s_{l}\bigg{)}\bigg{|}\bra{\mathbf{s}^{z}}\hat{U}^{\alpha
z}\ket{\psi}\bigg{|}^{2}$ (13)

where we have $\hat{U}^{z,\alpha}$ to convert from the $\hat{S}^{\alpha}$
basis to the $\hat{S}^{z}$ basis. This decomposition shows that by applying
$\hat{U}^{z,\alpha}$ after preparing the state $\ket{\psi}$ permits
measurements in the standard $\hat{S}^{z}$ basis that are the same as
measuring in the $\hat{S}^{\alpha}$ basis. Then, when we observe a measurement
outcome string $\mathbf{s}$ we can calculate the weights in the sums of Eq.
(13). Once we have the weights and the probabilities of measuring a particular
string, we can easily estimate
$\bra{\psi}\hat{\mathcal{O}}^{\alpha}_{k,l}\ket{\psi}$ on a classical
computer. Therefore, to estimate the string order operator expectation value
on a device requires implementing $\hat{U}^{xz}$ and $\hat{U}^{yz}$. We define
$\hat{U}^{xz}$ and $\hat{U}^{yz}$ in the standard $\hat{S}_{z}$ basis in Eqs.
(23) and (24), respectively.

##  IV Device architecture and problem mapping

There are a variety of routes to realize quantum simulation of spin-1 systems.
One route is to construct tailored quantum devices with spin-1 degrees of
freedom, such as cold atoms, trapped ions, or strongly correlated
superconducting circuits yip2003dimer ; garcia2004implementation ;
imambekov2003spin ; hilker2017revealing ; cohen2014proposal ;
cohen2015simulating ; senko2015realization ; albarran2018spin . However such
devices have either not been constructed or their construction and control may
be infeasible with current technology. Development and engineering of such
novel quantum devices would require significant, and potentially unknown, time
and cost.

An alternative route is to adapt currently developed quantum hardware to
perform the desired simulation, leveraging a tremendous body of work in the
understanding and engineering of such devices. There are a number of existing
device paradigms which can enable quantum simulation of spin-1 systems via
qubit-based digital quantum computing. Unfortunately, there is a distinct
overhead in using qubits to model spin-1 systems that will be particularly
inconvenient for near- and mid-term realizations. This overhead arises because
the local Hilbert space of each spin-1 particle is of dimension three,
requiring at least two qubits to represent each spin-1 particle.

To eliminate this overhead it is therefore natural to consider quantum devices
in which the basic physical element has local Hilbert space of at least three.
One of the leading device paradigms which can satisfy this requirement are
superconducting devices based on transmons krantz_quantum_2019 . Each transmon
is a nonlinear oscillator and computations can be performed in the low energy
subspace of these systems, allowing one to create, in principle, systems of
qudits krantz_quantum_2019 ; wu_high-fidelity_2020 . Recently, systems of
multiple interacting transmons operated as qutrits have been demonstrated,
suggesting a feasible route for enabling quantum simulations of interacting
spin-1 systems blok2021quantum .

We choose a mapping of the device-model Hilbert spaces that identifies the
eigenstates of the local spin-1 $z$ operator $\hat{S}^{z}$
($\ket{-},\ket{0},\ket{+}$) with the eigenstates of the local excitation
number operator $\hat{n}$ ($\ket{0},\ket{1},\ket{2}$) of each transmon:

$\displaystyle\ket{-}$ $\displaystyle\rightarrow\ket{0}$ (14)
$\displaystyle\ket{0}$ $\displaystyle\rightarrow\ket{1}$
$\displaystyle\ket{+}$ $\displaystyle\rightarrow\ket{2}.$

This composite Hilbert space formed by a system of $N$ spin-1 particles can be
realized by a system of $N$ transmons. For all discussions that follow this
mapping will be used.

Having selected the basic quantum information element, we now state our
assumptions about the device architecture. A variety of superconducting device
architectures based on transmons have been developed krantz_quantum_2019
including multiple types of transmons and modalities of transmon interactions.
We consider an architecture of tunable-frequency transmons mediated by tunable
couplers, as explored previously s li2022realization ; yan_tunable_2018 ;
sung_realization_2020 . Similar architectures have been generalized to large
devices with more than 50 transmons operating at high fidelity
arute_quantum_2019 .

Within this architecture, transmons are modeled as coupled anharmonic
(Duffing) oscillators with interactions between transmons mediated by a
tunable coupler krantz_quantum_2019 ; li2022realization ; yan_tunable_2018 ;
sung_realization_2020 . The effective device Hamiltonian then becomes

$\displaystyle H$
$\displaystyle=\sum_{i=1}^{N}\bigg{[}\omega_{i}\hat{n}_{i}+\frac{\delta_{i}}{2}\hat{n}_{i}(\hat{n}_{i}-1)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\delta\omega_{i}(t)\hat{n}_{i}+\varepsilon_{i}(t)(\hat{a}^{\dagger}_{i}+\hat{a}_{i})\bigg{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\sum_{\langle
i,j\rangle}g_{i,j}(t)(\hat{a}^{\dagger}_{i}\hat{a}_{j}+\hat{a}_{i}\hat{a}^{\dagger}_{j})$
(15)

where the operators $\hat{a}_{i}^{\dagger},\hat{a}_{i}$ are Bosonic creation
and annihilation operators, respectively, and
$\hat{n}_{i}=\hat{a}^{\dagger}_{i}\hat{a}_{i}$ is the Bosonic number operator.
The parameters $\omega_{i},\delta_{i}$ are the idling frequency and
anharmonicity of transmon $i$. The time-dependent functions in Eq. (IV), are
the frequency detunings $\delta\omega_{i}(t)$ of each transmon, the local
microwave controls are given by $\varepsilon_{i}(t)$, and the tunable coupling
is given by $g_{i,j}(t)$.

Figure 1: Optimal microwave controls for a single transmon that generates the
single-site Trotter evolution for the quench dynamics for two different
Trotter-step sizes $\tau=0.1$ and $\tau=0.01$. (Top panels) The optimal pulse
envelopes found via optimal control modulated at the transmon $0\rightarrow 1$
frequency $\omega_{01}$. The frequency spectrum of these pulses in the lab
frame are shown to the right of the time-domain pulses, with key transition
frequencies labeled. (Bottom panels) The optimal pulse envelopes modulated at
the transmon $1\rightarrow 2$ transition frequency $\omega_{12}$. The
frequency spectrum of these pulses in the lab frame are shown to the right of
the time-domain pulses. The solid lines and dashed-dot lines refer to the
optimal pulses that generate the Trotter-step unitary Eq. (8) at $\tau=0.1$
and $\tau=0.01$, respectively. Figure 2: Optimal microwave controls for a
single transmon that generates the rotation operators $\hat{U}^{zx}$ and
$\hat{U}^{zy}$, enabling measurement of the string order parameter in the
assumed device architecture. (Top panels) The optimal pulse envelopes found
via optimal control modulated at the transmon $0\rightarrow 1$ frequency
$\omega_{01}$. The frequency spectrum of these pulses in the lab frame are
shown to the right of the time-domain pulses, with key transition frequencies
labeled. (Bottom panels) The optimal pulse envelopes modulated at the transmon
$1\rightarrow 2$ transition frequency $\omega_{12}$. The frequency spectrum of
these pulses in the lab frame are shown to the right of the time-domain
pulses. The solid lines and dashed-dot lines refer to the optimal pulses that
generate the unitaries Eq. (23) and Eq. (24).

The density and layout of the transmons in a real system depends on a number
of practical limitations such as calibration complexity and noise
arute_quantum_2019 . since we are only considering a one-dimensional chain of
spin-1 particles, it suffices to assume that each transmon is connected to at
most two neighboring transmons in order to form the needed one dimensional
topology. This ensures that the unitaries implemented via optimal controls are
local within the device. Additional connections to ancilla transmons may be
useful for state preparation or observable measurements but we do not consider
these possibilities here.

We consider a general task of embedding a spin-$s$ particle into the first
$2s+1$ levels of an oscillator. This can be done compactly by defining a map
between the eigenvalues $n$ of the excitation number operator $\hat{n}$ and
the eigenvalues of the spin-$s$ $z$ operator $\hat{S}_{z}$:
$\hat{S}_{z}\ket{n}=(n+s)\ket{n}$ for $n\leq(2s+1)$, where
$\hat{n}\ket{n}=n\ket{n}$ are the eigenvectors of the local excitation number
operator.

We define the $z$ component of the total-spin operator for $M$ spins as
$\hat{J}_{z}=\sum_{i=1}^{M}S_{z}^{(i)}$, which satisfies the eigenvalue
equation $\hat{J}_{z}\ket{m}=m\ket{m}$. Now, we consider the action of
$\hat{J}_{z}$ on a state with a fixed excitation number
$\ket{\psi}=\otimes_{i=1}^{M}\ket{n_{i}}$ with $n_{i}\leq(2s+1)$:

$\displaystyle\hat{J}_{z}\ket{\psi}$
$\displaystyle=\otimes_{i=1}^{M}\hat{S}_{z}\ket{n_{i}}$ (16)
$\displaystyle=\otimes_{i=1}^{M}(n_{i}+s)\ket{n_{i}}$ (17)
$\displaystyle=\bigg{[}\prod_{i=1}^{M}(n_{i}+s)\bigg{]}\ket{\psi}.$ (18)

In other words, the $z$ component of total spin of state $\ket{\psi}$ is
$\prod_{i=1}^{M}(n_{i}+s)$. Most importantly, we note that the eigenvalue of
$\hat{J}_{z}$ will be the same for every state with the same number of
excitations because the eigenvalue $m=\prod_{i=1}^{M}(n_{i}+s)$ is invariant
to permutations of the number of excitations between sites.

Intuitively, this means that all states with a fixed $z$ component of total
spin $m$ lie within the subspace spanned by states with a fixed particle
number. This is important when considering implementing evolutions which will
preserve the $z$ component of total spin, which are common. For example,
typical two-local interactions of the form
$\hat{S}_{+}^{(i)}\hat{S}_{-}^{(j)}+\hat{S}_{+}^{(i)}\hat{S}_{-}^{(j)}\propto
S_{x}^{(i)}S_{x}^{(j)}+S_{y}^{(i)}S_{y}^{(j)}$ preserve the $z$ component of
total spin (as do Heisenberg-type interactions), and, therefore, any action
generated by such terms would be constrained to operate within a fixed
particle-number block of the oscillator system.

These considerations limit the dynamics for the two-site Trotter step defined
in Eq. (7) to a block-diagonal representation of the total particle number
basis for the transmon system. Therefore, analog device controls which
preserve total particle number are the most natural controls to use to attempt
to generate the unitary defined in Eq. (7).

Examining the assumed device Hamiltonian reveals that all the terms in the
Hamiltonian except for the microwave control will preserve particle number.
Thus we generate the desired two-site Trotter step unitary, Eq. (7), using
only transmon frequency and coupling controls. This choice of device controls
reduces simulation complexity and potentially reduces the complexity of pulse
calibration/characterization because the dynamics during gate operation are
constrained to blocks of total particle number.

##  V Optimal control results

Figure 3: Optimal controls for a two-transmon system that generates the two-
site Trotter evolution for the quench dynamics (Eq. (7)) for two different
Trotter-step sizes $\tau=0.1$ and $\tau=0.01$. (Top panels) The optimal
coupling controls found via optimal. The frequency spectrum of these pulses in
the lab frame are shown to the right of the time-domain pulses. (Middle
panels) The optimal detuning controls for transmon $1$ found via optimal. The
frequency spectrum of these pulses in the lab frame are shown to the right of
the time-domain pulses. (Bottom panels) The optimal detuning controls for
transmon $2$ found via optimal. The frequency spectrum of these pulses in the
lab frame are shown to the right of the time-domain pulses. The solid lines
and dashed-dot lines refer to the optimal pulses that generate the Trotter-
step unitary Eq. (7) at $\tau=0.1$ and $\tau=0.01$, respectively. Figure 4:
Time dynamics of the string observable in different directions
$\alpha=x,y,~{}z$ (left to right) for a 6-site AKLT model as a function of
string length when $b=0.2$. The solid lines are exact dynamics generated by
numerical diagonalization of the global 6-site quench Hamiltonian, Eq. (II).
The markers represent the approximate Trotterized dynamics generated via
repeated application of the final-time unitaries identified via optimal
control for two different Trotter step sizes $\tau=0.1$ and $\tau=0.01$.

In Section III, we outlined a set of four unitary operations to simulate and
observe the phenomena of string order melting within a superconducting
transmon device. The first two sets of unitaries, Eqs. (7) and (8), simulate
the time dynamics of string order melting. The second two, $\hat{U}^{zx}$(Eq.
(23)) and $\hat{U}^{zy}$(Eq. (24)), rotate into the correct basis so that
string order information can be extracted from the device.

In this section, we analyze a set of optimized device controls capable of
generating these unitaries. We begin by analyzing the unitaries defined on
single transmons: Eqs. (8), (23), and (24) and then proceed to discuss optimal
controls for a two-transmon system which generate Eq. (7). The numerical
methods used to generate these controls are detailed in Appendix B.

The Trotter-step unitary for a single-site transverse field perturbation is
defined in Eq. (8) and can be generated by modulating microwave-frequency
controls at the transition frequency of the transmon. In Fig. 1, we show a set
of optimal pulses at the two main transmon frequencies, $\omega_{01}$ and
$\omega_{12}$. These plots compare optimal controls for Trotter-step sizes
$\tau=0.1$ and $\tau=0.01$ with infidelities (as defined in Eq. (32)) of
$\approx 5\times 10^{-10}$ and $\approx 1\times 10^{-10}$, respectively.

From the results shown in Fig. 1, we find that the optimal control pulses for
a $50$ ns control time have high fidelity for two key reasons: 1) the pulses
are well localized in the frequency domain which prevents transitioning
population to higher levels and 2) the amplitude of the control pulses are
very low, which suppresses the off-resonant excitations.

The next results are for the basis changes needed to measure string order in
the transmon system: $\hat{U}^{zx}$(Eq. (23)) and $\hat{U}^{zy}$(Eq. (24)).
The optimal microwave controls which generate both unitaries were defined to
be $50$ ns and are shown in Fig. 2. In both cases, we observe similar pulse
characteristics such as amplitude limits and shifts in the frequency domain
around the primary drive frequencies. The primary difference between the
pulses is the small phase offset of the $\omega_{01}$ drive and large phase
offset of the $\omega_{12}$ drive. The phase of the pulses that generate
$\hat{U}^{zy}$ are offset from those that generate $\hat{U}^{zx}$ because the
phase of the microwave control pulses selects the axis of rotation in the
$x-y$ plane. Another slight difference is the infidelity of the pulses, which
is $\approx 3\times 10^{-6}$ and $\approx 2\times 10^{-5}$ for $\hat{U}^{zx}$
and $\hat{U}^{zy}$, respectively.

The difference in infidelity is due to different amounts of residual leakage
outside of the computational subspace. It can be seen in the power spectrum of
the optimal $\omega_{12}$ drive for $\hat{U}^{zy}$ has about an order of
magnitude more amplitude on the $\omega_{23}$ transition frequency,
contributing more to leakage. Why the optimizer was unable to find a better
solution is most likely due to the presence of a local minima in the
optimization landscape induced by our pulse constraints
riviello_searching_2015  and would probably be resolved by more exhaustive
numerical searches. However, we do not perform such exhaustive search because
these pulses need only be applied once after the state evolution and therefore
achieving significantly higher fidelity only negligibly affects the overall
simulation accuracy.

For the two-transmon optimal controls, in Fig. 3 we show a set of optimal
controls for Trotter-step sizes $\tau=0.1$ and $\tau=0.01$ with infidelities
of $\approx 4\times 10^{-6}$ and $\approx 4\times 10^{-8}$, respectively. The
similarity between the optimal pulses is primarily in the transmon detuning
(lower two rows) because the optimal controls for $\tau=0.1$ were used as the
starting point for the optimization of $\tau=0.01$. However, there is a
striking difference between the optimal coupling controls in the average
amplitude of the coupling, which corresponds to $0$ GHz in the frequency
domain. This decrease can be empirically understood as leading to less
generated rotation around the desired axis during the control time.

Overall, the optimal controls we identified meet a number of criteria
necessary for implementation in real devices. All of the pulses we found are
well within experimental bandwidth limitations and the control amplitudes are
also within feasible limits krantz_quantum_2019 ; sung_realization_2020 ;
li2022realization . Moreover, each control pulse requires about $60$
parameters to describe via the control ansatz we used (defined in Appendix B);
perhaps fewer if a more optimal basis is chosen. This suggests that
implementing and calibrating these pulses should be feasible in modern
superconducting devices egger_adaptive_2014 ; arute_quantum_2019 ;
sung_realization_2020 ; li2022realization .

##  VI Validation of melting dynamics with optimal control

Figure 5: Time dynamics of the $l=6$ string observable in different directions
$\alpha=x,y,~{}z$ (left to right) for a 6-site AKLT model as a function of
perturbation strength $b$. The solid lines are exact dynamics generated by
numerical diagonalization of the global 6-site quench Hamiltonian, Eq. (II).
The markers represent the approximate Trotterized dynamics generated via
repeated application of the final-time unitaries identified via optimal
control for two different Trotter step sizes $\tau=0.1$ and $\tau=0.01$. The
single-transmon optimal controls were intended to generate a perturbation of
$b=0.2$, thus to realize $b=0.2n$ the single-transmon optimal controls are
applied $n$ times per Trotter layer. Figure 6: Time dynamics of state
infidelity for a 6-site AKLT state generated by the Trotterized quench
dynamics found via optimal control for increasing strengths of the symmetry-
breaking perturbation. The two different marker types represent dynamics
approximating the true dynamics at two different Trotter step sizes $\tau=0.1$
and $\tau=0.01$.

We now present results from numerical simulations to validate that the optimal
controls identified enable the simulation of string order melting. We first
extract the unitaries generated by the optimal controls in the previous
section. For computational convenience, we choose to make a Markovian
assumption about the device dynamics, where we neglect population generated
outside of the computational subspace. This leads to a non-unitary simulation
but we find that the errors introduced by our optimal controls are less than
those induced by Trotter error.

First, in Fig. 4 we plot the time dynamics of the string observables for each
spatial direction as a function of time and string length for a weak
transverse field perturbation $b=0.2$. We observe that there is an initial
decrease in the string order parameter in every direction between $t=0$ and
$t=1$ which is in agreement with Ref. calvanese_strinati_destruction_2016  for
short strings. Importantly, we observe that the Trotter step of $\tau=0.1$ is
capable of tracking the qualitative dynamics for each observable with only
limited quantitative error.

Next, in Fig. 5, we plot the dynamics of the $l=6$ string order operator in
each spatial direction. We observe that for increasing strength of the
perturbation increases the rate of string order destruction for
$\langle\hat{\mathcal{O}}^{y}(t)\rangle$ and
$\langle\hat{\mathcal{O}}^{z}(t)\rangle$ but leads to less destruction of
string order for $\langle\hat{\mathcal{O}}^{x}(t)\rangle$. This is understood
because the perturbation being applied only partially breaks the protecting
symmetry and no string order melting occurs in the $x$ direction. This is,
again in qualitative agreement with the observations made in Ref.
calvanese_strinati_destruction_2016  for small string lengths.

Finally, in Fig. 6 we show the state fidelity as a function of perturbation
and for the Trotterized dynamics generated by the optimal controls. As
expected, we see that the smaller Trotter-step size yields a lower state
fidelity and we observe that the infidelity is much less dependent on the
perturbation strength for $\tau=0.01$. Importantly, we observe a non-monotonic
change in state fidelity for the larger Trotter step size and a monotonic one
for $\tau=0.01$. This indicates that choosing an optimal Trotter step size for
this problem may be difficult, depending on the timescale of dynamics one
wishes to observe. Because string order melting is a fast dynamical phenomena
a smaller $\tau$ may be better because the accuracy of the simulation as
determined by state fidelity seems to change logarithmically at short times.

##  VII Conclusion

In this work we have presented a path towards quantum simulation of string
order dynamics in a superconducting transmon architecture. We have used
quantum optimal control to verify that the necessary quantum processes can be
implemented and our numerical simulations based on these optimal controls
indicate that experimental observation of string order melting should be
feasible in the a near generation of superconducting devices. Our approach to
designing these simulations uses superconducting transmons to simulate spin-1
particles, microwave control lines to generate excitations, and tunable
couplers and tunable frequency transmons to drive particle-conserving
evolution that generates entanglement between transmons.

We employed numerical optimization methods to generate device controls that
are capable of driving the designed unitary evolutions with high fidelity. The
characteristics of these controls indicate that optimal pulses could be
implemented within existing devices, even when accounting for known hardware
limitations. We have used these controls to simulate the time dynamics of the
AKLT state in the presence of a quantum quench, and we have verified that
these simulations predict a number of qualitative and quantitative markers of
string order melting for SPT states.

The analyses of these simulations have further quantified the dependence of
state fidelity on evolution time, perturbation strength, and Trotter
decomposition order. We found that a choosing smaller time step $\tau$ is
significantly less sensitive to the perturbation strength and generally leads
to a monotonic decrease in state fidelity. However, it is important to note
that a smaller choice of $\tau$ requires a proportionally longer evolution
time because the number of repeated applications of the Trotter step. In this
case, $\tau=0.01$ requires a factor of 10 more gates to realize the same
simulation time of $\tau=0.1$. In the presence of decoherence or control
errors, the well known trade off in errors generated by Trotterization versus
gate depth are expected to limit the total time simulated salathe_digital_2015
; knee2015optimal .

Our results support simulating the dynamics of string order melting in a
transmon architecture as experimentally feasible, but there are still a number
of challenges for implementation that require further consideration. The most
immediate route of future work is through additional numerical simulations.

A clearer understanding of the relationships between the duration of the
controls, their infidelity, and ultimately their robustness to noise is
needed. We have fixed the pulse duration for single (50 ns) and two-transmon
(100 ns) optimal control simulations but optimal controls with comparable
fidelities may exist at smaller control times caneva_optimal_2009 ;
kirchhoff_optimized_2018 . Optimal pulses with smaller control times would
enable the observation of longer time dynamics within the same device
coherence times. Furthermore, quantifying or improving the robustness of
controls against noise within the applied field will aide in their
implementations egger_adaptive_2014 ; lysne2020small .

Another challenge is state preparation of the entangled spin-1 states within
the Haldane phase, such as the AKLT state. There are a variety of routes to
achieve this that may be feasible for near-term devices, for example,
variational state preparation and adiabatic or digitized adiabatic evolution
barends_digitized_2016 ; kandala_hardware-efficient_2017 . In particular,
because the Haldane phase is symmetry-protected topological phase, it can be
adiabatically connected to a product state via an evolution that breaks all of
the symmetries that protect the phase and therefore can be efficiently
prepared via adiabatic methods pollmann_entanglement_2010 .

Another approach to state preparation that is less efficient, but perhaps
complimentary, is to prepare the AKLT state by projecting into the ground-
state subspace via a sequence of parity measurements. This can be achieved by
noting that the AKLT Hamiltonian is a sum of Hermitian projectors that can be
converted to involutory Hermitian operators. These projections can be
implemented by controlled evolution using an ancilla qubit to perform a type
of parity measurement that projects the state into a subspace where the AKLT
state lies. This non-unitary route to preparing the AKLT state would also
enable direct measurement of the terms within the Hamiltonian or even a route
for dissipative state preparation kraus2008preparation ; verstraete2009quantum
.

In conclusion, we have shown how quantum simulation using quantum optimal
control offers a unique approach to study the static and dynamical properties
of a model Hamiltonian. Using numerical simulations, we have confirmed that
the quantum simulation of dynamical topological phenomena such as string order
melting are feasible in current superconducting devices. This affords many new
opportunities for the implementation and optimization of quantum simulations.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers,
Quantum Science Center and the U.S. Department of Energy, Office of Science,
Early Career Research Award. This research used resources of the Compute and
Data Environment for Science (CADES) at the Oak Ridge National Laboratory,
which is supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC05-00OR22725.

This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy (DOE). The U.S. Government
retains and the publisher, by accepting the article for publication,
acknowledges that the U.S. Government retains a nonexclusive, paid-up,
irrevocable, worldwide license to publish or reproduce the published form of
this manuscript, or allow others to do so, for U.S. Government purposes.

##  Appendix A String order parameter measurement

In this section we present the derivation of Eq. (13) for the expectation
value of a string order operator for an arbitrary state $\ket{\psi}$ in the
direction $\alpha$:

$\displaystyle\bra{\psi}\hat{\mathcal{O}}^{\alpha}_{k,l}\ket{\psi}$
$\displaystyle=\bra{\psi}\bigg{(}\hat{S}_{k}^{\alpha}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi\hat{S}_{n}^{\alpha}}\bigg{]}\hat{S}_{l}^{\alpha}\bigg{)}\ket{\psi}=\sum_{\mathbf{s}}\bra{\psi}\bigg{(}\hat{S}_{k}^{\alpha}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi\hat{S}_{n}^{\alpha}}\bigg{]}\hat{S}_{l}^{\alpha}\bigg{)}\ket{\mathbf{s}^{\alpha}}\bra{\mathbf{s}^{\alpha}}\ket{\psi}$
(19)
$\displaystyle=\sum_{\mathbf{s}}\bra{\psi}\bigg{(}s_{k}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi
s_{n}}\bigg{]}s_{l}\bigg{)}\ket{\mathbf{s}^{\alpha}}\bra{\mathbf{s}^{\alpha}}\ket{\psi}=\sum_{\mathbf{s}}\bigg{(}s_{k}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi
s_{n}}\bigg{]}s_{l}\bigg{)}\bra{\psi}\ket{\mathbf{s}^{\alpha}}\bra{\mathbf{s}^{\alpha}}\ket{\psi}$
(20)
$\displaystyle=\sum_{\mathbf{s}}\bigg{(}s_{k}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi
s_{n}}\bigg{]}s_{l}\bigg{)}\bra{\psi}U^{z\alpha\dagger}U^{z\alpha}\ket{\mathbf{s}^{\alpha}}\bra{\mathbf{s}^{\alpha}}U^{z\alpha\dagger}U^{z\alpha}\ket{\psi}$
(21)
$\displaystyle=\sum_{\mathbf{s}}\bigg{(}s_{k}\bigg{[}\prod_{n=k+1}^{l-1}e^{i\pi
s_{n}}\bigg{]}s_{l}\bigg{)}\bra{\psi}U^{\alpha
z\dagger}\ket{\mathbf{s}^{z}}\bra{\mathbf{s}^{z}}U^{\alpha z}\ket{\psi}.$ (22)

Using the definitions of the spin-1 operators and the commutation relations,
we express the unitary operators $\hat{U}^{z,\alpha}$ in the basis of
eigenvectors of $\hat{S}^{z}$:

$\displaystyle\hat{U}^{zx}$
$\displaystyle=\frac{1}{2}\ket{-}\bigg{(}\bra{+}-\sqrt{2}\bra{0}+\bra{-}\bigg{)}$
(23) $\displaystyle+\frac{1}{\sqrt{2}}\ket{0}\bigg{(}-\bra{+}+\bra{-}\bigg{)}$
$\displaystyle+\frac{1}{2}\ket{+}\bigg{(}\bra{+}+\sqrt{2}\bra{0}+\bra{-}\bigg{)}$

and

$\displaystyle\hat{U}^{zy}$
$\displaystyle=\frac{1}{2}\ket{-}\bigg{(}-\bra{+}+i\sqrt{2}\bra{0}+\bra{-}\bigg{)}$
(24) $\displaystyle+\frac{1}{\sqrt{2}}\ket{0}\bigg{(}\bra{+}+\bra{-}\bigg{)}$
$\displaystyle+\frac{1}{2}\ket{+}\bigg{(}-\bra{+}-i\sqrt{2}\bra{0}+\bra{-}\bigg{)}.$

##  Appendix B Numerical methods

We utilize numerical optimal control techniques to simulate the dynamics
generated by Eq. (IV) and iteratively optimize the time-dependent device
controls to minimize the infidelity between the evolved system dynamics and a
target unitary operator. These methods and implementations follow what has
been discussed in previous work kairys_parametrized_2021 ; kairys2021efficient
and we review the details for completeness and convenience.

We identify optimal controls that generate four unitary operators as defined
in Eqs. (6), (23), (24), and (7). The first three unitaries are each defined
on a single-transmon Hilbert space and therefore we use a single transmon
Hamiltonian derived from Eq. (IV) that omits the frequency detuning control
and uses only microwave control parameters:

$\displaystyle H_{D}$
$\displaystyle=\omega_{1}\hat{n}_{1}+\frac{\delta_{1}}{2}\hat{n}_{1}(\hat{n}_{1}-1)+\varepsilon_{1}(t)(\hat{a}^{\dagger}_{1}+\hat{a}_{1}).$
(25)

For the two-site Trotter step unitary, we use a two-transmon Hamiltonian
derived from Eq. (IV) omitting the microwave control lines because the target
unitary, Eq. (7), is block-diagonal in the total particle number basis, such
that

$\displaystyle H_{D}$
$\displaystyle=\omega_{1}\hat{n}_{1}+\frac{\delta_{1}}{2}\hat{n}_{1}(\hat{n}_{1}-1)+\delta\omega_{1}(t)\hat{n}_{1}$
(26)
$\displaystyle+\omega_{2}\hat{n}_{2}+\frac{\delta_{2}}{2}\hat{n}_{2}(\hat{n}_{2}-1)+\delta\omega_{2}(t)\hat{n}_{2}$
(27)
$\displaystyle+~{}g_{1,2}(t)(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}_{1}\hat{a}^{\dagger}_{2}).$
(28)

In both cases, quantum information processing is performed in reference to a
pre-calibrated rotating frame associated with the idling frequency of each
transmon. This transformation is given by the unitary transform

$\displaystyle R(t)$ $\displaystyle=\bigotimes_{i}R_{i}(t)$ (29)
$\displaystyle=\bigotimes_{i}\exp[\frac{it}{\hbar}\omega_{i}\hat{n}_{i}\bigg{]}$
(30)

which yields an effective Hamiltonian in this rotating frame,
$H^{R}=i\hbar(\partial_{t}R(t))R^{\dagger}(t)+R(t)H_{D}R^{\dagger}(t)$, and
the corresponding time-ordered evolution operator for control time $T_{c}$ of
a transmon system is

$U_{D}(T_{c})=\mathcal{T}\exp[-\frac{i}{\hbar}\int_{0}^{T_{c}}d\tau
H^{R}(\tau)\bigg{]}.$ (31)

We then use the Gradient Optimization of Analytic Controls (GOAT) algorithm to
evaluate gradients of the objective function with respect to control
parameters machnes_tunable_2018 . Our primary objective function to be
minimized is the unitary infidelity of a quantum processes using a projective
$SU$ measure derived from the Hilbert-Schmidt inner product palao_optimal_2003
:

$g(\vec{\alpha})=1-\frac{1}{d^{2}}\bigg{|}Tr(U^{\dagger}_{T}P_{c}U_{D}(\vec{\alpha},T_{c})P_{c})\bigg{|}^{2}$
(32)

where $U_{T}$ is the target unitary operator we wish to prepare,
$U_{D}(\vec{\alpha},T_{c})$ is the unitary evolution operator for the device
with control parameters $\vec{\alpha}$ and control time $T_{c}$, $P_{c}$ is a
projection onto the desired computational subspace and $d$ is the dimension of
the computational subspace.

Within the GOAT algorithm, each control field $f(t)$ (e.g.,
$\delta\omega_{i}(t)$, $\varepsilon_{i}(t)$, $g_{i,j}(t)$) is decomposed into
a (not-necessarily orthonormal) function basis which is parameterized by a set
of real numbers – enabling traditional numerical optimization techniques to be
used to optimize the device controls.

Here we describe each control field via a functional form which includes both
the control parameters to be optimized as well as additional functions which
constrain the optimization to a class of practical pulses. For all of our
optimizations we describe each control field $f(t)$ (e.g.,
$\delta\omega_{i}(t)$, $\varepsilon_{i}(t)$, $g_{i,j}(t)$) as

$\displaystyle
f(\vec{\alpha},t)=\Omega(t)\cos(\overline{\omega}t)S(h(\vec{\alpha},t))$ (33)

Where we have defined a carrier frequency $\overline{\omega}$, a window
function $\Omega(t)$ to ensure that the pulse turns on and off smoothly, a
saturation function $S(x)$ to ensure that the optimal pulses stay within a
pre-specified amplitude range, and the parameterized function
$h(\vec{\alpha},t)$. The window function is a flat-top cosine defined as

$\displaystyle\Omega(t)=\begin{cases}\frac{1-\cos(\pi
t/\tau_{r})}{2}\Omega_{m}&0\leq t\leq\tau_{r}\\\ \Omega_{m}&\tau_{r}\leq
t\leq(\tau_{c}-\tau_{r})\\\
\frac{1-\cos(\pi(\tau_{c}-t)/\tau_{r})}{2}\Omega_{m}&(\tau_{c}-\tau_{r})\leq
t\leq\tau_{c},\end{cases}$ (34)

where $\tau_{c}=T_{c}$ is the total control time, $\Omega_{m}$ scales the
magnitude of the control pulse, and $\tau_{r}$ is the ramp time, which was
chosen to be $0.3\tau_{c}$ to reduce spectral leakage tripathi_operation_2019
. The saturation function is a generalized logistic function defined as

$S(x)=-B-\frac{2B}{1-3\exp(-\frac{4x}{B})}$ (35)

where $B=0.08\text{ GHz}/2\pi$ for microwave controls, $B=0.5\text{ GHz}/2\pi$
for frequency detuning controls, $B=0.01\text{ GHz}/2\pi$ for coupling
controls. These parameters were chosen to agree roughly with control
limitations currently observed in superconducting transmon devices of the
assumed architecture sung_realization_2020 ; li2022realization . Finally, we
expand each control field as a linear combination of $N$ analytic functions,
which we choose to be sinusoidal functions with varying amplitude, frequency,
and phase:

$h(\vec{\alpha},t)=\bigg{[}\sum_{n}^{N}\alpha_{n,1}\sin(\alpha_{n,2}t+\alpha_{n,3})\bigg{]}.$
(36)

For the microwave controls operating on a single transmon we optimize two
drive channels as in Eq. (33) each at a different carrier frequency
$\overline{\omega}=\omega_{01},\omega_{12}$, which are the transition
frequencies between the transmon levels $\ket{0}\rightarrow\ket{1}$ and
$\ket{1}\rightarrow\ket{2}$, respectively. We also assume the control time for
a single qutrit operation is $T_{c}=50$ ns. For each channel we set $N=10$ in
Eq. (36). This means that each of the single-transmon control pulses require
optimization of $60$ parameters.

For the coupling control and frequency detuning controls we choose a carrier
frequency of $\overline{\omega}=0$, which places the dynamics within those
typically used to generate resonant two-qubit gates on superconducting
hardware sung_realization_2020 ; krantz_quantum_2019 . We assume a control
time of $100$ ns for the two-transmon evolution because this coincides with
experimental timescales realizable in current devices sung_realization_2020 ;
li2022realization . Moreover, we set $N=7$ in Eq. (36) for each control
channel which gives a total of $63$ total parameters for optimization.

For the microwave controls we penalize controls with intermediate-time leakage
to higher subspaces by computing a functional and adding it to the infidelity
in Eq. (32):

$\displaystyle\mathcal{L}(\vec{\alpha})$
$\displaystyle=\frac{1}{T_{c}}\int_{0}^{T_{c}}d\tau\Tr(P_{c}U^{\dagger}(\vec{\alpha},\tau)P_{d}U(\vec{\alpha},\tau)P_{c})$
(37)

where we have defined $P_{c}$ as the projector onto the computational subspace
of the transmon (In our simulations of spin-1 particles, this means
span$(\\{\ket{0},\ket{1},\ket{2}\\})$ and the operator $P_{d}$ assigns weights
to specific leakage levels depending on the importance that little
intermediate population lie in that state. In our simulations we set
$P_{d}=0.1\ket{3}\bra{3}+1\ket{4}\bra{4}$ to weakly penalize leakage to
$\ket{3}$ and strongly on $\ket{4}$.

The precise choice of weights for the leakage subspaces was not significant as
we constrained the optimizer to look for pulses well localized within the
frequency domain at the allowed transition frequencies. However the leakage
penalty does help ensure that the optimal pulses identified will generalize to
a true anharmonic quantum oscillator with an infinite number of energy levels.

In order to perform gradient based optimization we require the gradient of
this function with respect to a parameter $\alpha$ we derive the gradient of
the leakage penalty as

$\displaystyle\partial_{\alpha}\mathcal{L}$
$\displaystyle=\partial_{\alpha}\frac{1}{T_{c}}\int_{0}^{T_{c}}d\tau\Tr(P_{c}U^{\dagger}(\vec{\alpha},\tau)P_{d}U(\vec{\alpha},\tau)P_{c})$
(38)
$\displaystyle=\frac{1}{T_{c}}\int_{0}^{T_{c}}d\tau\partial_{\alpha}\Tr(P_{c}U^{\dagger}(\vec{\alpha},\tau)P_{d}U(\vec{\alpha},\tau)P_{c})$
(39) $\displaystyle=\frac{1}{T_{c}}\int_{0}^{T_{c}}d\tau$

and we see that calculating this gradient requires only knowledge of
$\partial_{\alpha}U(\vec{\alpha},t)$, which is obtained through the GOAT
method machnes_tunable_2018 . A leakage penalty was not added to the two-site
operator optimization because non-computational states are actively being used
for quantum information processing in that case kairys_parametrized_2021 .

In all simulations, we do not implement the rotating wave approximation to
ensure a more accurate estimate of gate fidelity. Moreover, we model each
transmon as a 5-level system to fully account for leakage in the two-transmon
evolutions and to accurately quantify leakage in the single transmon
simulations.

We perform optimizations to identify controls that generate the Trotter-step
operators (Eqs. (7) and (8)) at two Trotter step sizes: $\tau=0.1$ and
$\tau=0.01$. We seed the optimization for Trotter step size $\tau=0.1$ with a
random initial guess and $\tau=0.01$ with the optimal controls obtained for
$\tau=0.1$.

We use experimental hardware parameters for the frequency and anharmonicity of
each transmon derived from Ref. blok2021quantum  in which a system of
transmons are used as qutrits. Specifically we define $\omega_{1}=5.634\text{
GHz}/2\pi$, $\delta_{1}=-0.266\text{ GHz}/2\pi$, $\omega_{2}=5.447\text{
GHz}/2\pi$, and $\delta_{2}=-0.270\text{ GHz}/2\pi$.

We implement the GOAT algorithm using the programming language Julia and
various open-source packages Bezanson_Julia_A_fresh_2017 . Our implementation
uses the Julia package DifferentialEquations.jl to numerically solve the
coupled GOAT equations of motion using a order 5/4 Runge-Kutta method with
adaptive time stepping rackauckas2017differentialequations . For the gradient-
based control optimization of $\vec{\alpha}$, we use a limited-memory Broyden-
Fletcher-Goldfarb-Shanno (L-BFGS) algorithm with a backtracking line-search
method which are implemented in the Optim.jl package and LineSearches.jl
package, respectively mogensen2018optim . We limit each optimization to 2000
iterations of L-BFGS and define a stopping criteria when the infinity-norm of
the gradient falls below 1e-9 or the relative change in the objective function
is below 1e-8. For further details on the derivations of gradients via the
GOAT algorithm we refer the reader to our previous work
kairys_parametrized_2021 .

Finally, to validate that our resulting optimal controls should enable the
observation of string order melting in real superconducting hardware we
perform a numerical simulation to observe the dynamics of string order under
the optimized controls. Specifically, we compute the dynamics of a 6-site
spin-1 system using the resulting optimal controls and compare the results
obtained via exact numerical integration of the model. We initialize the
system in the AKLT state obtained via diagonalization of the AKLT Hamiltonian
Eq. (1). Then, we evolve the state under the exact quench dynamics given by
Eq. (II) and compare these with the dynamics generated by two sets of
Trotterized dynamics generated via optimal controls for two different Trotter
step sizes: $\tau=0.1$ and $\tau=0.01$. We then compute the string order
operator expectation value and state fidelity to draw our final conclusions.

## References

  * [1] I. M. Georgescu, S. Ashhab, and F. Nori.  Quantum simulation.  Reviews of Modern Physics, 86(1):153–185, March 2014.
  * [2] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff.  Digital Quantum Simulation of Spin Models with Circuit Quantum Electrodynamics.  Physical Review X, 5(2):021027, June 2015.
  * [3] A. Parra-Rodriguez, P. Lougovski, L. Lamata, E. Solano, and M. Sanz.  Digital-analog quantum computation.  Physical Review A, 101(2):022305, February 2020.
  * [4] Lucas Lamata, Adrian Parra-Rodriguez, Mikel Sanz, and Enrique Solano.  Digital-analog quantum simulations with superconducting circuits.  Advances in Physics: X, 3(1):1457981, January 2018.
  * [5] A. Martin, L. Lamata, E. Solano, and M. Sanz.  Digital-analog quantum algorithm for the quantum Fourier transform.  Physical Review Research, 2(1):013012, January 2020.
  * [6] L. C. Celeri, D. Huerga, F. Albarrán-Arriagada, E. Solano, and M. Sanz.  Digital-analog quantum simulation of fermionic models.  arXiv:2103.15689 [quant-ph], March 2021.  arXiv: 2103.15689.
  * [7] L. Lamata.  Digital-analog quantum simulation of generalized Dicke models with superconducting circuits.  Scientific Reports, 7(1):43768, March 2017.
  * [8] I. Arrazola, J. S. Pedernales, L. Lamata, and E. Solano.  Digital-Analog Quantum Simulation of Spin Models in Trapped Ions.  Scientific Reports, 6(1):30534, July 2016.
  * [9] A. Galicia, B. Ramon, E. Solano, and M. Sanz.  Enhanced connectivity of quantum hardware with digital-analog control.  Physical Review Research, 2(3):033103, July 2020.
  * [10] D. V. Babukhin, A. A. Zhukov, and W. V. Pogosov.  Hybrid digital-analog simulation of many-body dynamics with superconducting qubits.  Physical Review A, 101(5):052337, May 2020.
  * [11] Paul Kairys and Travis S. Humble.  Parametrized hamiltonian simulation using quantum optimal control.  Phys. Rev. A, 104:042602, Oct 2021.
  * [12] Steffen J. Glaser, Ugo Boscain, Tommaso Calarco, Christiane P. Koch, Walter Köckenberger, Ronnie Kosloff, Ilya Kuprov, Burkhard Luy, Sophie Schirmer, Thomas Schulte-Herbrüggen, Dominique Sugny, and Frank K. Wilhelm.  Training Schrödinger’s cat: quantum optimal control: Strategic report on current status, visions and goals for research in Europe.  The European Physical Journal D, 69(12):279, December 2015.
  * [13] Austin G. Fowler, Matteo Mariantoni, John M. Martinis, and Andrew N. Cleland.  Surface codes: Towards practical large-scale quantum computation.  Phys. Rev. A, 86:032324, Sep 2012.
  * [14] Barbara M Terhal.  Quantum error correction for quantum memories.  Reviews of Modern Physics, 87(2):307, 2015.
  * [15] T. Senthil.  Symmetry-Protected Topological Phases of Quantum Matter.  Annual Review of Condensed Matter Physics, 6(1):299–324, 2015.
  * [16] Ian Affleck, Tom Kennedy, Elliott H. Lieb, and Hal Tasaki.  Rigorous results on valence-bond ground states in antiferromagnets.  Physical Review Letters, 59(7):799–802, August 1987.
  * [17] Frank Pollmann, Ari M. Turner, Erez Berg, and Masaki Oshikawa.  Entanglement spectrum of a topological phase in one dimension.  Physical Review B, 81(6):064439, February 2010.
  * [18] Frank Pollmann and Ari M. Turner.  Detection of symmetry-protected topological phases in one dimension.  Physical Review B, 86(12):125441, September 2012.
  * [19] Marcello Calvanese Strinati, Leonardo Mazza, Manuel Endres, Davide Rossini, and Rosario Fazio.  Destruction of string order after a quantum quench.  Physical Review B, 94(2):024302, July 2016.
  * [20] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. GSL. Brandao, D. A. Buell, et al.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, October 2019.
  * [21] MS Blok, VV Ramasesh, T Schuster, K O’Brien, JM Kreikebaum, D Dahlen, A Morvan, Beni Yoshida, NY Yao, and I Siddiqi.  Quantum information scrambling on a superconducting qutrit processor.  Physical Review X, 11(2):021010, 2021.
  * [22] David C McKay, Thomas Alexander, Luciano Bello, Michael J Biercuk, Lev Bishop, Jiayin Chen, Jerry M Chow, Antonio D Córcoles, Daniel Egger, Stefan Filipp, et al.  Qiskit backend specifications for openqasm and openpulse experiments.  arXiv preprint arXiv:1809.03452, 2018.
  * [23] Ian Affleck, Tom Kennedy, Elliott H Lieb, and Hal Tasaki.  Valence bond ground states in isotropic quantum antiferromagnets.  In Condensed matter physics and exactly soluble models, pages 253–304. Springer, 1988.
  * [24] Hiroshi Ueda, Hiroki Nakano, and Koichi Kusakabe.  Finite-size scaling of string order parameters characterizing the haldane phase.  Physical Review B, 78(22):224402, 2008.
  * [25] SK Yip.  Dimer state of spin-1 bosons in an optical lattice.  Physical review letters, 90(25):250402, 2003.
  * [26] Juan J Garcia-Ripoll, Miguel A Martin-Delgado, and J Ignacio Cirac.  Implementation of spin hamiltonians in optical lattices.  Physical review letters, 93(25):250405, 2004.
  * [27] Adilet Imambekov, Mikhail Lukin, and Eugene Demler.  Spin-exchange interactions of spin-one bosons in optical lattices: Singlet, nematic, and dimerized phases.  Physical Review A, 68(6):063602, 2003.
  * [28] Timon A Hilker, Guillaume Salomon, Fabian Grusdt, Ahmed Omran, Martin Boll, Eugene Demler, Immanuel Bloch, and Christian Gross.  Revealing hidden antiferromagnetic correlations in doped hubbard chains via string correlators.  Science, 357(6350):484–487, 2017.
  * [29] Itsik Cohen and Alex Retzker.  Proposal for verification of the haldane phase using trapped ions.  Physical review letters, 112(4):040503, 2014.
  * [30] I Cohen, P Richerme, Z-X Gong, C Monroe, and A Retzker.  Simulating the haldane phase in trapped-ion spins using optical fields.  Physical Review A, 92(1):012334, 2015.
  * [31] C Senko, P Richerme, J Smith, A Lee, I Cohen, A Retzker, and C Monroe.  Realization of a quantum integer-spin chain with controllable interactions.  Physical Review X, 5(2):021026, 2015.
  * [32] F Albarrán-Arriagada, L Lamata, E Solano, G Romero, and JC Retamal.  Spin-1 models in the ultrastrong-coupling regime of circuit qed.  Physical Review A, 97(2):022306, 2018.
  * [33] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver.  A quantum engineer’s guide to superconducting qubits.  Applied Physics Reviews, 6(2):021318, June 2019.
  * [34] Xian Wu, S. L. Tomarken, N. Anders Petersson, L. A. Martinez, Yaniv J. Rosen, and Jonathan L. DuBois.  High-Fidelity Software-Defined Quantum Logic on a Superconducting Qudit.  Physical Review Letters, 125(17):170502, October 2020.
  * [35] Shaowei Li, Daojin Fan, Ming Gong, Yangsen Ye, Xiawei Chen, Yulin Wu, Huijie Guan, Hui Deng, Hao Rong, He-Liang Huang, et al.  Realization of fast all-microwave controlled-z gates with a tunable coupler.  Chinese Physics Letters, 39(3):030302, 2022.
  * [36] F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver.  Tunable Coupling Scheme for Implementing High-Fidelity Two-Qubit Gates.  Physical Review Applied, 10(5):054062, November 2018.
  * [37] Y. Sung, L. Ding, J. Braumüller, A. Vepsäläinen, B. Kannan, M. Kjaergaard, A. Greene, G. O. Samach, C. McNally, D. Kim, et al.  Realization of high-fidelity cz and $zz$-free iswap gates with a tunable coupler.  Phys. Rev. X, 11:021058, Jun 2021.
  * [38] Gregory Riviello, Katharine Moore Tibbetts, Constantin Brif, Ruixing Long, Re-Bing Wu, Tak-San Ho, and Herschel Rabitz.  Searching for quantum optimal controls under severe constraints.  Physical Review A, 91(4):043401, April 2015.
  * [39] D. J. Egger and F. K. Wilhelm.  Adaptive Hybrid Optimal Quantum Control for Imprecisely Characterized Systems.  Physical Review Letters, 112(24):240503, June 2014.
  * [40] George C Knee and William J Munro.  Optimal trotterization in universal quantum simulators under faulty control.  Physical Review A, 91(5):052327, 2015.
  * [41] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro.  Optimal Control at the Quantum Speed Limit.  Physical Review Letters, 103(24):240501, December 2009.
  * [42] S. Kirchhoff, T. Kebler, P. J. Liebermann, E. Assémat, S. Machnes, F. Motzoi, and F. K. Wilhelm.  Optimized cross-resonance gate for coupled transmon systems.  Physical Review A, 97(4):042348, April 2018.
  * [43] Nathan K Lysne, Kevin W Kuper, Pablo M Poggi, Ivan H Deutsch, and Poul S Jessen.  Small, highly accurate quantum processor for intermediate-depth quantum simulations.  Physical review letters, 124(23):230501, 2020.
  * [44] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. Las Heras, R. Babbush, A. G. Fowler, B. Campbell, Yu Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, E. Solano, H. Neven, and John M. Martinis.  Digitized adiabatic quantum computing with a superconducting circuit.  Nature, 534(7606):222–226, June 2016.  arXiv: 1511.03316.
  * [45] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, September 2017.
  * [46] Barbara Kraus, Hans P Büchler, Sebastian Diehl, Adrian Kantian, Andrea Micheli, and Peter Zoller.  Preparation of entangled states by quantum markov processes.  Physical Review A, 78(4):042307, 2008.
  * [47] Frank Verstraete, Michael M Wolf, and J Ignacio Cirac.  Quantum computation and quantum-state engineering driven by dissipation.  Nature physics, 5(9):633–636, 2009.
  * [48] Paul Kairys and Travis S Humble.  Efficient quantum gate discovery with optimal control.  In 2021 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 413–418. IEEE, 2021.
  * [49] S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm.  Tunable, Flexible, and Efficient Optimization of Control Pulses for Practical Qubits.  Physical Review Letters, 120(15):150401, April 2018.
  * [50] J. P. Palao and R. Kosloff.  Optimal control theory for unitary transformations.  Physical Review A, 68(6):062308, December 2003.
  * [51] Vinay Tripathi, Mostafa Khezri, and Alexander N Korotkov.  Operation and intrinsic error budget of a two-qubit cross-resonance gate.  Physical Review A, 100(1):012301, 2019.
  * [52] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah.  Julia: A fresh approach to numerical computing.  SIAM Review, 59(1):65–98, 9 2017.
  * [53] C. Rackauckas and Q. Nie.  Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia.  Journal of Open Research Software, 5(1), 2017.
  * [54] P. K. Mogensen and A.N. Riseth.  Optim: A mathematical optimization package for Julia.  Journal of Open Source Software, 3(24):615, 2018.