# Accelerated Newton-Raphson GRAPE methods for optimal control

David L. Goodwin<EMAIL_ADDRESS>Chemistry Research Laboratory,
University of Oxford, Mansfield Road, Oxford OX1 3TA, UK  Institute for
Biological Interfaces 4 – magnetic resonance, Karlsruhe Institute for
Technology (KIT), Karlsruhe, Germany  Mads Sloth Vinding<EMAIL_ADDRESS>Center of Functionally Integrative Neuroscience (CFIN), Department of Clinical
Medicine, Faculty of Health, Aarhus University, Denmark

###### Abstract

A Hessian based optimal control method is presented in Liouville space to
mitigate previously undesirable polynomial scaling of computation time. This
new method, an improvement to the state-of-the-art Newton-Raphson GRAPE
method, is derived with respect to two exact time-propagator derivative
techniques: auxiliary matrix and ESCALADE methods. We observed that compared
to the best current implementation of Newton-Raphson GRAPE method, for an
ensemble of 2-level systems, with realistic conditions, the new auxiliary
matrix and ESCALADE Hessians can be 4-200 and 70-600 times faster,
respectively.

Optimal Control, Hessian, Newton-Raphson, GRAPE

###### pacs:

02.30.Yy, 02.60.Pn, 03.67.-a, 87.80.Lg

##  I Introduction

The problem of transferring the state of a dynamical system to a desired
target state while minimising the remaining distance and costs is often solved
with optimal control theory [1, 2]. Applications include quantum sensing [3,
4, 5, 6], quantum computing [7, 8, 9], and nuclear magnetic resonance (nmr)
spectroscopy [10, 11, 12, 13] and imaging (mri) [14, 15].

A number of different approaches to optimal control (OC) has lead to the
development of different methods: Lagrangian methods [16, 17, 18, 19];
minimal-time OC [10]; gradient ascent pulse engineering (grape) [12];
sophisticated gradient-free searches [20, 21]; Krylov-Newton methods [22]; OC
with a basis of analytic controls [23, 24]; a tensor product approach for
large quantum systems [25]. The method outlined in this work is based on a
piecewise-constant control pulse approximation [26, 27, 28, 29] of grape [12,
30, 31] using a gradient following numerical optimisation.

Although finding an optimal solution to the problem of controlling a single
2-level system from a defined initial system state to a desired target state,
a _state-to-state_ problem, is considered straight-forward and computationally
inexpensive with modern methods and computing power; an OC problem can become
numerically and computationally arduous [32], particularly for applications
that account for practical hardware configurations and limitations [11, 33,
34]. Additionally, the computational expense can increase dramatically when
optimising over an ensemble of systems, such as the case in solid-state
nuclear magnetic resonance [13, 35, 36], where the ensemble also includes
crystalline orientations of a powder average [37]. The particular application
of OC of interest to the authors is that of a neural network based method for
mri [38, 39, 40] and a method of morphic OC [41, 42], requiring hundreds of
thousands of optimised pulse shapes to form their optimal libraries.

Keeping the example of mri, the utility of OC is highlighted when considering
the legal constraints of power deposition safeguards [43] and the obvious
financial rewards of reducing the time a patient stays inside the mri machine.

Modern techniques can mitigate protracted numerical convergence with a
quadratically convergent optimisation method, requiring the calculation of a
Hessian matrix, giving a large saving in the number of required serial
optimisation iterations [31, 35, 44]. However, it is known within the
community that calculation of the Hessian matrix does not scale well to a
control problem with a large number of controllable amplitudes [45] (also
shown in Fig. 1).

This communication presents a jump in computational efficiency of multiple-
orders of magnitude with a rework of the original method, devised to calculate
the exact Hessian matrix [46, 31, 35], and the recently published exact,
matrix-free, method of efficient spin control using analytical Lie algebraic
derivatives (escalade) [44]. The manuscript will present the mathematical
formulation of these exact Newton-Raphson grape methods in the irreducible
spherical-tensor basis of a Liouville space and show a newly devised method to
calculate the Hessian matrix with $\mathcal{O}(N)$ scaling, compared to the
previous $\mathcal{O}(N^{2})$ scaling (Fig. 1). Results show the comparative
speedup of this method to the original in the context of state-to-state mri
problems, requiring a generalisation of the escalade method with calculations
of optimal $\mathrm{z}$-controls.

##  II Exact Newton-Raphson optimal control

###  II.1 Optimal control in Liouville space

A quantum system can be described by a density operator, $\hat{\rho}(t)$; a
time-dependent system state. The evolution of this state is dictated by the
Liouville-von Neumann equation,

$\frac{\partial}{\partial
t}\hat{\rho}(t)=-i\big{[}\hat{H}(t),\hat{\rho}(t)\big{]}$ (1)

where $\hat{H}(t)$ is a time-dependent Hamiltonian, an operator in a Hilbert
space with a spectrum of allowed energy levels. The usual factor of $\hbar$ is
dropped here, resulting in the eigenspectrum of $\hat{H}$ expressed in angular
frequency units. The methods presented in this manuscript are particular to a
Liouville space, also named the adjoint representation[30, 31, 35]. A system
state in a Liouville space is represented by a vector, $|\hat{\rho}\rangle$,
obtained by stacking columns of the density operator $\hat{\rho}$ with Eq. (1)
becoming

$\displaystyle\frac{\partial}{\partial
t}\big{|}\hat{\rho}(t)\big{\rangle}=-i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}(t)\big{|}\hat{\rho}(t)\big{\rangle},$
$\displaystyle\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}(t)\triangleq\openone\\!\otimes\\!\hat{H}(t)-\hat{H}(t)^{\intercal}\\!\otimes\\!\openone$
(2)

where the identity matrix, $\openone$, and Hamiltonian have the same
dimension.

The form of a bilinear control problem is to split that which is controllable,
_the control_ , from that which is not, _the drift_. The Liouvillian for a
control problem with $\mathrm{x}$-, $\mathrm{y}$-, and $\mathrm{z}$-controls
on a 2-level system can be written as

$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}(t)=\underbrace{\omega\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{z}}}_{\text{drift}}+\underbrace{c_{\mathrm{x}}(t)\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{x}}+c_{\mathrm{y}}(t)\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{y}}+c_{\mathrm{z}}(t)\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{z}}}_{\text{control}}$
(3)

where the angular frequency $\omega$ is the time-independent resonant
frequency offset, $c_{\mathrm{x,y,z}}(t)$ are time-dependent control
amplitudes, and
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{x,y,z}}$
are Pauli matrices of a Liouville space.

The OC method of grape [12] uses piecewise constant control pulses, where
control pulses are constant over a small time interval, $\Delta t$ [28, 29]:

$\displaystyle
c_{k}(t)\to\begin{bmatrix}c_{k,1}&c_{k,2}&\cdots&c_{k,N}\end{bmatrix},$
$\displaystyle k\in\\{\mathrm{x},\mathrm{y},\mathrm{z}\\}$ (4)

using the notation $c_{k,n}\equiv c_{k}(t_{n})$ for convenience, and
$t_{N}=N\Delta t$. This discrete formulation allows numerical solution to Eq.
(2), given an initial state of the system $|\hat{\rho}_{0}\rangle$, through
time-ordered propagation

$|\hat{\rho}_{N}\rangle=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!N}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!N-1}\cdots\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!2}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!1}|\hat{\rho}_{0}\rangle$
(5)

where
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}$
are time-propagators of an isolated time-slice and are defined through the
exponential map

$\displaystyle|\hat{\rho}_{n}\rangle=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}|\hat{\rho}_{n-1}\rangle,$
$\displaystyle\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}=\mathrm{e}^{-i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.61948pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.61948pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{n}\Delta
t}$ (6)

The matrix exponential of Eq. (6) usually calculated with the Padé
approximant, Taylor series, or Krylov propagation [35].

An additional method to calculate the time-propagators of a 2-level system is
by explicitly calculating the elements of the matrix. In a spherical-tensor
basis, this matrix is the Wigner-matrix [47]:

$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}=\begin{bmatrix}\alpha^{2}&\sqrt{2}\alpha\beta&\beta^{2}\\\
-\sqrt{2}\alpha\beta^{\ast}&\alpha\alpha^{\ast}-\beta\beta^{\ast}&\sqrt{2}\alpha^{\ast}\beta\\\
{\beta^{\ast}}^{2}&-\sqrt{2}\alpha^{\ast}\beta^{\ast}&{\alpha^{\ast}}^{2}\end{bmatrix}$
(7)

which is formulated in terms of the complex elements

$\displaystyle\alpha=\cos{\phi}+i\frac{z}{r}\sin{\phi},$
$\displaystyle\beta=\frac{y}{r}\sin{\phi}+i\frac{x}{r}\sin{\phi},$ (8)

where the shorthands $x=c_{\mathrm{x},n}$, $y=c_{\mathrm{y},n}$, and
$z=c_{\mathrm{z},n}+\omega$ are used, $\phi=\frac{1}{2}r\Delta t$ is a polar
angle of rotation, and $r=\sqrt{x^{2}+y^{2}+z^{2}}$ is the polar radius.

As a notational convenience for what follows, the following _effective
propagators_ are defined as the effect of the pulse between time slices $m$
and $n$:

$\displaystyle\mathbf{U}_{m}^{n}\triangleq\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n-1}\cdots\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!m+1}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!m},$
$\displaystyle\forall(1\leqslant m\leqslant n\leqslant N)$ (9)

where a backward (time-reversed) propagation can be denoted by
$\mathbf{U}_{n}^{m}={\mathbf{U}_{m}^{n}}^{\dagger}$.

Optimal control requires a metric to optimise; the _fidelity_ [48, 12],
$\mathcal{F}$, a measure how well the pulses perform a desired control task.
The task of the OC problem is to find a set of control amplitudes,
$c_{\mathrm{x,y,z}}(t)$, that maximise the fidelity e.g. the real part of an
inner product:

$\max_{c_{\mathrm{x,y,z}}(t)}\\!\big{(}\mathcal{F}\big{)}=\max_{c_{\mathrm{x,y,z}}(t)}\\!\Big{(}\mathrm{Re}\big{\langle}\hat{\sigma}\big{|}\mathbf{U}_{1}^{N}|\hat{\rho}_{0}\big{\rangle}\Big{)}$
(10)

This form of the fidelity metric is defined in terms of state-to-state
problems, where a system is in a defined initial state,
$|\hat{\rho}_{0}\rangle$, and the control task is to take this state to a
desired one, $|\hat{\sigma}\rangle$. With the notation introduced in Eq. (9),
$\mathbf{U}_{1}^{N}$ is interpreted as the effective propagator over the
shaped pulse. In addition to Eq. (5), the system is propagated backwards from
the desired target state,

$|\hat{\chi}_{n}\rangle=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}^{\dagger}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n+1}^{\dagger}\cdots\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!N-1}^{\dagger}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!N}^{\dagger}|\hat{\sigma}\rangle$
(11)

which is termed the _adjoint state_ of the control problem [35].

Grape is a gradient following numerical optimisation method and requires
derivatives of the fidelity with respect to the controls. In turn, this
requires the directional propagator derivatives, $D_{\\!k,n}^{a}$; with the
subscripts $k$ denoting the derivative in the direction of
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{k}$
and $n$ denoting the derivative operating on the time-propagator
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}$.
For each time-slice, $n$, and for each control direction
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{k}\in\\{\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{x}},\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{y}},\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{\mathrm{z}}\\}$,
the $a^{\text{th}}$ order derivative takes the form

$\nabla^{a}\\!\mathcal{F}(c_{k,n})=\mathrm{Re}\underbrace{\big{\langle}\hat{\chi}_{n+1}}_{\mathclap{N\text{
times}}}\\!\big{|}\overbrace{\\!D_{\\!k,n}^{a}\hat{\rho}_{n-1}\big{\rangle}}^{\mathclap{N\text{
times}}},$ (12)

where bra-ket notation explicitly shows vector structures i.e.
$|D_{\\!k,n}^{a}\hat{\rho}_{n-1}\rangle=D_{\\!k,n}^{a}|\hat{\rho}_{n-1}\rangle$.
The number of forward and backward propagations is indicated for each control
channel to produce a gradient vector, $\nabla\\!\mathcal{F}$, or the diagonal
elements of a Hessian matrix, $\nabla^{2}\\!\mathcal{F}$.

This is sufficient for a fidelity gradient, scaling linearly with $N$, but a
fidelity Hessian also requires mixed second order derivatives [31], where the
off-diagonal Hessian elements are

$\nabla^{2}\\!\mathcal{F}(c_{k,n},c_{j,m})=\mathrm{Re}\underbrace{\\!\big{\langle}\hat{\chi}_{n+1}{D_{\\!k,n}}}_{\mathclap{N\text{
times}}}\\!\big{|}\overbrace{\\!\mathbf{U}_{m+1}^{n-1}\big{|}D_{\\!j,m}\hat{\rho}_{m-1}\big{\rangle}}^{\mathclap{\frac{1}{2}N(N-1)\text{
times}}}$ (13)

Clearly, the form of Eq. (13) has a central propagator that cannot be absorbed
into the bra or ket because it depends on both $t_{n}$ and $t_{m}$, and
therefore the computation scales polynomially with $\frac{1}{2}N(N-1)$ (the
factor $\frac{1}{2}$ comes from the symmetric property of a Hessian [35]).
This is known within the OC community [45] and is highlighted in Fig. 1. The
linear plots of the fidelity and gradient, on log-log axes, show these
calculations are efficient with increasing $N$, whereas the Hessian
calculation time is not linear on these log-log axes. The subject of this
manuscript is to mitigate this undesirable scaling, resulting in a linearly
scaling Hessian calculation (dashed lines in Fig. 1), after the following
section outlines calculation of directional propagator derivatives.

Figure 1: Average wall-clock time of fidelity, gradient, and Hessian for
increasing $N$. Solid lines show the Newton-Raphson grape method [31] and
dashed lines show the proposed _accelerated_ Newton-Raphson grape method. Both
methods use the auxiliary matrix method [46] to calculate propagator
derivatives.

###  II.2 Directional propagator derivatives with auxiliary matrices

As has been published previously [46, 31, 35], exact propagator derivatives
required by Eq. (12) and (13) can be calculated by exponentiating an auxiliary
matrix: resulting in an upper triangular block matrix [49] with a time-
propagator,
$\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}$,
on the block diagonal, and with the directional derivative of that propagator,
calculated in the direction of a control operator, in the upper triangular
block. The first and second order propagator derivatives can be extracted from

$\exp{\begin{bmatrix}\mathbf{A}_{n}&\mathbf{C}_{j}&\bm{0}\\\
\bm{0}&\mathbf{A}_{n}&\mathbf{C}_{k}\\\
\bm{0}&\bm{0}&\mathbf{A}_{n}\end{bmatrix}}=\begin{bmatrix}\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}&D_{\\!j,n}&\frac{1}{2}D_{\\!jk,n}^{2}\\\
\bm{0}&\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}&D_{\\!k,n}\\\
\bm{0}&\bm{0}&\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}\end{bmatrix}$
(14)

where the block matrix is formed from
$\mathbf{A}_{n}=-i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{n}\Delta
t$, as in Eq. (6), and
$\mathbf{C}_{k}=-i\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{L}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{L}}}}_{k}\Delta
t$ are functions of the control operators with
$j,k\in\\{\mathrm{x},\mathrm{y},\mathrm{z}\\}$.

###  II.3 Directional propagator derivatives with escalade

Whereas Eq. (14) is exact, the matrix exponential is expensive. Foroozandeh
and Singh recently devised a method that is free from this expensive matrix
exponential calculation, in the OC method of Efficient Spin Control using
Analytical Lie Algebraic Derivatives (escalade) [44]. This method can also be
extended with interaction propagator splitting [50].

The efficient method of escalade to calculate the directional derivatives in
Eqs. (12) and (13), in this single spin model, is to construct a matrix with
rows containing all elements needed to construct the propagator derivatives

$\begin{bmatrix}\mathbf{\Theta}_{\\!\mathrm{x}}\\\
\mathbf{\Theta}_{\\!\mathrm{y}}\\\
\mathbf{\Theta}_{\\!\mathrm{z}}\end{bmatrix}=\mathrm{vec}[\openone]+\frac{\sin^{2}{\phi}}{\phi
r}\mathrm{vec}[\mathbf{S}]+\frac{2\phi-\sin{2\phi}}{2\phi
r^{2}}\mathrm{vec}[\mathbf{S}^{2}]$ (15)

where $\mathrm{vec}[\openone]$ is a vectorised identity matrix, with a
vectorisation operation on a matrix $\mathbf{A}$ such that
$\mathbf{A}=\mathrm{vec}^{-1}[\mathrm{vec}[\mathbf{A}]]$. The skew-symmetric
matrix $\mathbf{S}$ is

$\mathbf{S}=\begin{bmatrix}0&z&-y\\\ -z&0&x\\\ y&-x&0\end{bmatrix}$ (16)

and the symmetric matrix $\mathbf{S}^{2}$ is calculated algebraically. In
turn, the directional propagator derivatives can be written as

$\displaystyle
D_{\\!k,n}=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}\Big{[}\mathrm{vec}^{-1}\big{[}\mathbf{\Sigma}\mathbf{\Theta}_{\\!k}\big{]}\Big{]}$
(17) $\displaystyle\begin{aligned}
D_{\\!jk,n}^{2}=&\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}\Big{[}\mathrm{vec}^{-1}\big{[}\mathbf{\Sigma}\mathbf{\Theta}_{\\!jk}\big{]}\Big{]}+\\\
&\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n}\Big{[}\mathrm{vec}^{-1}\big{[}\mathbf{\Sigma}\mathbf{\Theta}_{\\!j}\big{]}\Big{]}\Big{[}\mathrm{vec}^{-1}\big{[}\mathbf{\Sigma}\mathbf{\Theta}_{\\!k}\big{]}\Big{]}\end{aligned}$
(18)

where
$\mathbf{\Sigma}=\begin{bmatrix}\mathrm{vec}[\mathbf{C}_{\mathrm{x}}]&\mathrm{vec}[\mathbf{C}_{\mathrm{y}}]&\mathrm{vec}[\mathbf{C}_{\mathrm{z}}]\end{bmatrix}$
is a 3-column matrix to be multiplied with the 3-row matrix of
$\mathbf{\Theta}_{\\!k}$. Using a similar notation to Eq. (15), the
$\mathbf{\Theta}_{\\!jk}$ matrices required by the second order derivatives in
Eq. (18) are

$\displaystyle\begin{bmatrix}\mathbf{\Theta}_{\\!j\mathrm{x}}\\\
\mathbf{\Theta}_{\\!j\mathrm{y}}\\\
\mathbf{\Theta}_{\\!j\mathrm{z}}\end{bmatrix}=$
$\displaystyle\frac{\cos(2\phi)-1}{r^{2}}\mathrm{vec}\bigg{[}\mathbf{S}\frac{\mathrm{d}\mathbf{S}}{\mathrm{d}j}+\frac{\mathrm{d}\mathbf{S}}{\mathrm{d}j}\mathbf{S}\bigg{]}$
$\displaystyle+\frac{2\phi-\sin(2\phi)}{r^{3}}\mathrm{vec}\bigg{[}\frac{\mathrm{d}\mathbf{S}}{\mathrm{d}j}\bigg{]}$
$\displaystyle+\frac{2k(1-\cos(2\phi)-\phi\sin(2\phi))}{r^{3}}\mathrm{vec}[\mathbf{S}]$
$\displaystyle+\frac{3k\sin(2\phi)-2k\phi(2+\cos(2\phi))}{r^{4}}\mathrm{vec}[\mathbf{S}^{2}]$
(19)

where three such equations are required for
$j\in\\{\mathrm{x},\mathrm{y},\mathrm{z}\\}$ and the derivatives of
$\mathbf{S}$ are derived algebraically from Eq. (16).

Since the directional derivatives of Eq. (17) and Eq. (18), together with
time-propagators of Eq. (7), do not involve any matrix operations, other than
a few trivial multiplications, escalade offers substantial computational gains
relative to the auxmat method in Eq. (14).

##  III Accelerated Newton-Raphson optimal control

Moving away from a chosen calculation method of directional propagator
derivatives, the remaining bottleneck of the Newton-Raphson methods presented
above is the off-diagonal Hessian elements; the mixed derivatives of Eq. (13).

As a starting point, escalade includes an additional efficiency, which also
applies to the auxmat method, where the central effective propagator in Eq.
(13) can be split to
$\mathbf{U}_{m+1}^{n-1}=\mathbf{U}_{m+1}^{n-1}\mathbf{U}_{1}^{m}\mathbf{U}_{m}^{1}=\mathbf{U}_{1}^{n-1}\mathbf{U}_{m}^{1}$,
Eq. (13) becomes

$\nabla^{2}\\!\mathcal{F}(c_{k,n},c_{j,m})=\mathrm{Re}\big{\langle}\hat{\chi}_{n+1}D_{\\!k,n}\big{|}\mathbf{U}_{1}^{n-1}\mathbf{U}_{m}^{1}\big{|}D_{\\!j,m}\hat{\rho}_{m-1}\big{\rangle}$
(20)

An interpretation of these two central effective propagators is: the right
hand side, $\langle\hat{\chi}_{n+1}D_{\\!k,n}|$, is multiplied by the
effective propagator $\mathbf{U}_{1}^{n-1}$, presenting the directional
derivative of $t_{n}$ to be evaluated at $t_{0}$; the left hand side,
$|D_{\\!j,m}\hat{\rho}_{m-1}\rangle$, is multiplied by the time-reversed
effective propagator $\mathbf{U}_{m}^{1}$, presenting the directional
derivative of $t_{m}$ to also be evaluated at $t_{0}$.

To outline how this can be overcome, a representation of a _trajectory_ is
introduced:

$\big{[}\bm{\rho}\big{]}_{m}^{n}\triangleq\begin{bmatrix}|\hat{\rho}_{n}\rangle&|\hat{\rho}_{n-1}\rangle&\cdots&|\hat{\rho}_{m+1}\rangle&|\hat{\rho}_{m}\rangle\end{bmatrix}$
(21)

which is an array of column vectors and the whole trajectory from the dynamics
in Eq. (3) is contained in $[\bm{\rho}]_{0}^{N}$. Trajectory analysis is
useful for visualising pulse dynamics [51]. Taking the concept of a single
matrix containing all trajectory information, a _directional derivative
trajectory_ , evaluated at $t_{0}$ as with the right side of Eq. (20), can be
defined as

$\big{[}\partial_{j}^{[0]}\\!\bm{\rho}\big{]}_{0}^{n}\triangleq\begin{bmatrix}|\mathbf{U}_{n}^{1}D_{\\!j,n}\hat{\rho}_{n-1}\rangle&\cdots&|\mathbf{U}_{2}^{1}D_{\\!j,2}\hat{\rho}_{1}\rangle&|\mathbf{U}_{1}^{\dagger}D_{\\!j,1}\hat{\rho}_{0}\rangle\end{bmatrix}$
(22)

where the superscript [0] is used to indicate evaluation at $t_{0}$.

With the realisation that a _directional derivative trajectory_ in Eq. (22) is
a matrix in itself, $n$ Hessian elements can be calculated with one matrix-
vector product with

$\begin{bmatrix}\nabla^{2}\\!\mathcal{F}(c_{k,n},c_{j,n})\\\ \vdots\\\
\nabla^{2}\\!\mathcal{F}(c_{k,n},c_{j,2})\\\
\nabla^{2}\\!\mathcal{F}(c_{k,n},c_{j,1})\end{bmatrix}^{\\!\mathrm{T}}\\!\\!\\!\\!=\mathrm{Re}\underbrace{\big{\langle}\hat{\chi}_{n+1}D_{\\!k,n}\mathbf{U}_{1}^{n-1}}_{\mathclap{N\text{
times
backward}}}\\!\big{|}\overbrace{\\!\big{[}\partial_{j}^{[0]}\\!\bm{\rho}\big{]}_{0}^{n}}^{\mathclap{N\text{
times forward}}}$ (23)

Given that the total effective propagator $\mathbf{U}_{1}^{N}$ can be
calculated with forward propagation, $\mathbf{U}_{1}^{n-1}$ can be updated
during subsequent backward propagation with
$\mathbf{U}_{1}^{n-1}=\mathchoice{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{6.77017pt}}\smash{\hat{P}}}}{\hat{\vphantom{\rule{1.0pt}{7.06079pt}}\smash{\hat{P}}}}_{\\!\\!n-1}^{\dagger}\mathbf{U}_{1}^{n}$.
The single vector-matrix product of Eq. (23), per time-slice, is expected to
be much more efficient than the $n-1$ vector-vector products, per time-slice,
in Eq. (13).

##  IV Acceleration of a broadband MRI example

Figure 2: The average speedup of the accelerated escalade method (_acc_
-escalade) and the accelerated auxiliary matrix method (_acc_ -auxmat),
compared to the auxiliary matrix method. The number of time-steps, $N$, the
time-slice width, $\Delta t$, and the pulse duration, $T$, are parameters
affecting computation time. In each of the three plots one parameter is kept
constant: (A) $N=128$, (B) $\Delta t=7.8125$ $\mathrm{\SIUnitSymbolMicro s}$;
and (C) $T=1$ $\mathrm{ms}$. The dotted lines in each plot show the speedup of
_acc_ -escalade relative to _acc_ -auxmat.

The authors choose to set the test of the method presented above in the
context of mri because there is a very real need to have computationally fast
OC methods to run _on-the-fly_ with a need for patient specific solutions.
Often, OC-facilitated mri exploits simultaneous/parallel control systems e.g.
eight radiofrequency channels ($\mathrm{x}$\- and $\mathrm{y}$-controls) [14,
52] and three magnetic field gradients and/or numerous local shims
($\mathrm{z}$-controls) [52, 53, 40].

The optimisation is set as robust for an ensemble of $101$ spin-$\frac{1}{2}$
systems, with offset bandwidth
$\frac{\omega}{2\pi}\in[-\frac{1}{2},+\frac{1}{2}]~{}$\mathrm{kHz}$$,
simulated using a single block-diagonal Liouvillian with each block operating
on a single ensemble member [37, 41].

The time-slice, $\Delta t$, is often fixed to the hardware digitalisation
dwell-time: typically a few $\mathrm{\SIUnitSymbolMicro s}$. Modern mri
systems handle $N>1000$, but the pulse duration, $T=N\Delta t$, is usually
kept around a few $\mathrm{ms}$ to avoid degradation of pulse performance due
to transverse relaxation. As a set of speed tests, these three variables are
incremented and set as a range of OC problems and are shown in Fig. 2. The
vertical dashed lines indicate a set of variables physically relevant to mri.

Figure 2 shows the average speedup of the accelerated auxiliary matrix method
(_acc_ -auxmat) and the accelerated escalade method (_acc_ -escalade), both
relative to the standard auxiliary matrix method. A separate fidelity,
gradient, and Hessian calculation is performed for each method from a random
control waveform, with $c_{\mathrm{x}}$, $c_{\mathrm{y}}$, and
$c_{\mathrm{z}}$ controls in the order of $2\pi\times
10^{3}~{}$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$. The average is over 84
different random control pulses, run in parallel on 28 cpu cores.

Figure 2A shows that the effect of increasing $T$ while also increasing
$\Delta t$ is, approximately, a constant speedup of _acc_ -auxmat and _acc_
-escalade Hessian. This can be attributed to the effective propagator
splitting of Eq. (20). The small speedup of _acc_ -auxmat fidelity and
gradient calculations is due to an algorithmic efficiency from coding in a way
that lends itself to the accelerated Hessian methods.

The main result of this manuscript is the speedup of Hessian calculations of
the _acc_ -auxmat and _acc_ -escalade methods, with speedup increasing as $N$
increases, in Fig. 2B and Fig. 2C, to over $\times 100$ when a large $N$ is
used. Furthermore, this trend does not appear to dampen, indicating further
speedup when $N\gg 1000$.

From the evidence in Fig. 2B and Fig. 2C, the accelerated Hessian calculations
of Eq. (23) do indeed remove the $\mathcal{O}(N^{2})$ scaling of the original
Newton-Raphson method [31], reducing to a linear scaling $\mathcal{O}(N)$.

##  V Conclusion

A new mathematical formulation of the Newton-Raphson grape method [31] has
been presented in Liouville space and applies to optimal control problems with
unitary evolution. The escalade method [44] recast the cumbersome problem of
finding derivatives – using trigonometric evaluations of vectorised arrays
rather than matrix-matrix products and computationally expensive matrix
exponentials – and factorised the central propagator to avoid
$\mathcal{O}(N^{2})$ scaling, now reduced to $\mathcal{O}(N)$, in the
computation of the Hessian. The escalade method is shown in a new light of
Liouville space, with the additional derivation of $\mathrm{z}$-controls,
which are important to mri.

The main result is a new formulation of the Newton-Raphson grape method in
Liouville space, being an optimisation method with quadratic convergence to an
optimal solution, reducing the expensive polynomial scaling of the control
problem to a linear scaling when increasing the number of piecewise-constant
pulses in optimal solution. Speedup increases as the number of time-slices,
$N$, increases: $\times 4$ to $\times 200$, for $N=100$ and $N=1000$
respectively. Furthermore, employing escalade within this new Hessian
calculation method shows further speedup of $\times 70$ to $\times 600$, for
$N=100$ and $N=1000$ respectively, when compared to the original Newton-
Raphson grape method.

###### Acknowledgements.

DLG thanks Pranav Singh and Mohammadali Foroozandeh for useful discussions on
the escalade method, Burkhard Luy for his continued support, and Martin
Goodwin for proofreading. MSV would like to thank Villum Fonden, Eva og Henry
Fraenkels Mindefond, Harboefonden, and Kong Christian Den Tiendes Fond.

## References

  * Glaser  _et al._ [2015] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Training schrödinger’s cat: quantum optimal control, Eur. Phys. J. D 69, 279 (2015).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: a european community view, New J. Phys. 20, 080201 (2018).
  * Saywell  _et al._ [2018] J. C. Saywell, I. Kuprov, D. Goodwin, M. Carey, and T. Freegarde, Optimal control of mirror pulses for cold-atom interferometry, Phys. Rev. A 98, 023625 (2018).
  * Saywell [2020] J. C. Saywell,  _Optimal control of cold atoms for ultra-precise quantum sensors_ , Ph.D. thesis, University of Southampton, UK (2020).
  * Saywell  _et al._ [2020] J. Saywell, M. Carey, M. Belal, I. Kuprov, and T. Freegarde, Optimal control of raman pulse sequences for atom interferometry, J. Phys. B: At. Mol. Opt. Phys. 53, 085006 (2020).
  * Saywell  _et al._ [2021] J. Saywell, M. Carey, N. Dedes, I. Kuprov, and T. Freegarde, Can optimised pulses improve the sensitivity of atom interferometers?, in  _Quantum Technology: Driving Commercialisation of an Enabling Science II_, Vol. 11881, edited by M. J. Padgett, K. Bongs, A. Fedrizzi, and A. Politi, International Society for Optics and Photonics (SPIE, 2021) pp. 83–92.
  * Zhang  _et al._ [2012] J. Zhang, R. Laflamme, and D. Suter, Experimental implementation of encoded logical qubit operations in a perfect quantum error correcting code, Phys. Rev. Lett. 109, 100503 (2012).
  * Waldherr  _et al._ [2014] G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, and J. Wrachtrup, Quantum error correction in a solid-state hybrid spin register, Nature 506, 204 (2014).
  * Dolde  _et al._ [2014] F. Dolde, V. Bergholm, Y. Wang, I. Jakobi, B. Naydenov, S. Pezzagna, J. Meijer, F. Jelezko, P. Neumann, T. Schulte-Herbrüggen, J. Biamonte, and J. Wrachtrup, High-fidelity spin entanglement using optimal control, Nat. Commun. 5, 3371 (2014).
  * Khaneja  _et al._ [2001] N. Khaneja, R. Brockett, and S. J. Glaser, Time optimal control in spin systems, Phys. Rev. A 63, 032308 (2001).
  * Skinner  _et al._ [2003] T. E. Skinner, T. O. Reiss, B. Luy, N. Khaneja, and S. J. Glaser, Application of optimal control theory to the design of broadband excitation pulses for high-resolution NMR, J. Magn. Reson. 163, 8 (2003).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, J. Magn. Reson. 172, 296 (2005).
  * Tošner  _et al._ [2006] Z. Tošner, S. J. Glaser, N. Khaneja, and N. C. Nielsen, Effective Hamiltonians by optimal control: Solid-state NMR double-quantum planar and isotropic dipolar recoupling, J. Chem. Phys. 125, 184502 (2006).
  * Vinding  _et al._ [2017a] M. S. Vinding, D. Brenner, D. H. Y. Tse, S. Vellmer, T. Vosegaard, D. Suter, T. Stöcker, and I. I. Maximov, Application of the limited-memory quasi-Newton algorithm for multi-dimensional, large flip-angle RF pulses at 7T, Magn. Reson. Mater. Phy. 30, 29 (2017a).
  * Reeth  _et al._ [2019] E. V. Reeth, H. Ratiney, K. Tse Ve Koon, M. Tesch, D. Grenier, O. Beuf, S. J. Glaser, and D. Sugny, A simplified framework to optimize MRI contrast preparation, Magn. Reson. Med. 81, 424 (2019).
  * Somlói  _et al._ [1993] J. Somlói, V. A. Kazakov, and D. J. Tannor, Controlled dissociation of I2 via optical transitions between the X and B electronic states, Chem. Phys. 172, 85 (1993).
  * Zhu  _et al._ [1998] W. Zhu, J. Botina, and H. Rabitz, Rapidly convergent iteration methods for quantum optimal control of population, J. Chem. Phys. 108, 1953 (1998).
  * Maday and Turinici [2003] Y. Maday and G. Turinici, New formulations of monotonically convergent quantum control algorithms, J. Chem. Phys. 118, 8191– (2003).
  * Eitan  _et al._ [2011] R. Eitan, M. Mundt, and D. J. Tannor, Optimal control with accelerated convergence: Combining the Krotov and quasi-Newton methods, Phys. Rev. A 83, 053426 (2011).
  * Doria  _et al._ [2011] P. Doria, T. Calarco, and S. Montangero, Optimal control technique for many-body quantum dynamics, Phys. Rev. Lett. 106, 190501 (2011).
  * Rach  _et al._ [2015] N. Rach, M. M. Müller, T. Calarco, and S. Montangero, Dressing the chopped-random-basis optimization: A bandwidth-limited access to the trap-free landscape, Phys. Rev. A 92, 062343 (2015).
  * Ciaramella  _et al._ [2015] G. Ciaramella, A. Borzì, G. Dirr, and D. Wachsmuth, Newton methods for the optimal control of closed quantum spin systems, SIAM J. Sci. Comput. 37, A319 (2015).
  * Machnes  _et al._ [2018] S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm, Tunable, flexible, and efficient optimization of control pulses for practical qubits, Phys. Rev. Lett. 120, 150401 (2018).
  * Lucarelli [2018] D. Lucarelli, Quantum optimal control via gradient ascent in function space and the time-bandwidth quantum speed limit, Phys. Rev. A 97, 062346 (2018).
  * Quiñones Valles  _et al._ [2019] D. Quiñones Valles, S. Dolgov, and D. Savostyanov, Tensor product approach to quantum control, in  _Integral Methods in Science and Engineering: Analytic Treatment and Numerical Approximations_, edited by C. Constanda and P. Harris (Springer, 2019) pp. 367–379.
  * Conolly  _et al._ [1986] S. Conolly, D. Nishimura, and M. A., Optimal control solutions to the magnetic resonance selective excitation problem, IEEE Trans. Med. Imaging 5, 106 (1986).
  * Mao  _et al._ [1986] J. Mao, T. H. Mareci, K. N. Scott, and E. R. Andrew, Selective inversion radiofrequency pulses by optimal control, J. Magn. Reson. 70, 310 (1986).
  * Viola and Lloyd [1998] L. Viola and S. Lloyd, Dynamical suppression of decoherence in two-state quantum systems, Phys. Rev. A 58, 2733 (1998).
  * Viola  _et al._ [1999] L. Viola, E. Knill, and S. Lloyd, Dynamical decoupling of open quantum systems, Phys. Rev. Lett. 82, 2417 (1999).
  * de Fouquieres  _et al._ [2011] P. de Fouquieres, S. G. Schirmer, S. J. Glaser, and I. Kuprov, Second order gradient ascent pulse engineering, J. Magn. Reson. 212, 412 (2011).
  * Goodwin and Kuprov [2016] D. L. Goodwin and I. Kuprov, Modified Newton-Raphson GRAPE methods for optimal control of spin systems, J. Chem. Phys. 144, 204107 (2016).
  * Pechen and Tannor [2011] A. N. Pechen and D. J. Tannor, Are there traps in quantum control landscapes?, Phys. Rev. Lett. 106, 120402 (2011).
  * Kobzar  _et al._ [2004] K. Kobzar, T. E. Skinner, N. Khaneja, S. J. Glaser, and B. Luy, Exploring the limits of broadband excitation and inversion pulses, J. Magn. Reson. 170, 236 (2004).
  * Kobzar  _et al._ [2008] K. Kobzar, T. E. Skinner, N. Khaneja, S. J. Glaser, and B. Luy, Exploring the limits of broadband excitation and inversion: II. Rf-power optimized pulses, J. Magn. Reson. 194, 58 (2008).
  * Goodwin [2017] D. L. Goodwin,  _Advanced optimal control methods for spin systems_ , Ph.D. thesis, University of Southampton, UK (2017).
  * Tos̆ner  _et al._ [2021] Z. Tos̆ner, M. J. Brandl, J. Blahut, S. J. Glaser, and B. Reif, Maximizing efficiency of dipolar recoupling in solid-state NMR using optimal control sequences, Sci. Adv. 7, eabj5913 (2021).
  * Kuprov [2016] I. Kuprov, Fokker-Planck formalism in magnetic resonance simulations, J. Magn. Reson. 270, 124 (2016).
  * Vinding  _et al._ [2019] M. S. Vinding, B. Skyum, R. Sangill, and T. E. Lund, Ultrafast (milliseconds), multidimensional RF pulse design with deep learning, Magn. Reson. Med. 82, 586 (2019).
  * Vinding  _et al._ [2021a] M. S. Vinding, C. S. Aigner, S. Schmitter, and T. E. Lund, DeepControl: 2DRF pulses facilitating inhomogeneity and B0 off-resonance compensation in vivo at 7T, Magn. Reson. Med. 85, 3308 (2021a).
  * Vinding  _et al._ [2021b] M. S. Vinding, D. L. Goodwin, I. Kuprov, and T. E. Lund, Optimal control gradient precision trade-offs: Application to fast generation of DeepControl libraries for MRI, J. Magn. Reson. 333, 107094 (2021b).
  * Goodwin  _et al._ [2020] D. L. Goodwin, M. R. M. Koos, and B. Luy, Second order phase dispersion by optimised rotation pulses, Phys. Rev. Research 2, 033157 (2020).
  * Haller  _et al._ [2022] J. D. Haller, D. L. Goodwin, and B. Luy, SORDOR pulses: expansion of the Böhlen-Bodenhausen scheme for low-power broadband magnetic resonance, Magn. Reson. 3, 53 (2022).
  * Dempsey  _et al._ [2002] M. F. Dempsey, B. Condon, and D. M. Hadley, MRI safety review, in  _Seminars in Ultrasound, CT and MRI_, Vol. 23 (2002) pp. 392–401.
  * Foroozandeh and Singh [2021] M. Foroozandeh and P. Singh, Optimal control of spins by analytical Lie algebraic derivatives, Automatica 129, 109611 (2021).
  * Kallies [2018] W. Kallies,  _Concurrent optimization of robust refocused pulse sequences for magnetic resonance spectroscopy_ , Ph.D. thesis, Technische Universität München, DE (2018).
  * Goodwin and Kuprov [2015] D. L. Goodwin and I. Kuprov, Auxiliary matrix formalism for interaction representation transformations, optimal control, and spin relaxation theories, J. Chem. Phys. 143, 084113 (2015).
  * Siminovitch [1997] D. J. Siminovitch, Rotations in NMR: Part I. Euler-Rodrigues parameters and quaternions, Concepts Magn. Reson. 9, 149 (1997).
  * Glaser  _et al._ [1998] S. J. Glaser, T. Schulte-Herbrüggen, M. Sieveking, O. Schedletzky, N. C. Nielsen, O. W. Sørensen, and C. Griesinger, Unitary control in quantum ensembles: Maximizing signal intensity in coherent spectroscopy, Science 280, 421 (1998).
  * Van Loan [1978] C. F. Van Loan, Computing integrals involving the matrix exponential, IEEE Trans. Automat. Contr. 23, 395 (1978).
  * Goodwin  _et al._ [2022] D. L. Goodwin, P. Singh, and M. Foroozandeh, Adaptive optimal control of entangled qubits (2022), accepted for publication in Sci. Adv.
  * Kuprov [2013] I. Kuprov, Spin system trajectory analysis under optimal control pulses, J. Mag. Reson. 233, 107 (2013).
  * Vinding  _et al._ [2017b] M. Vinding, B. Guérin, T. Vosegaard, and N. Nielsen, Local SAR, global SAR, and power-constrained large-flip-angle pulses with optimal control and virtual observation points, Magnetic Resonance in Medicine 77, 374 (2017b).
  * Stockmann and Wald [2018] J. P. Stockmann and L. L. Wald, In vivo B0 field shimming methods for MRI at 7 T, NeuroImage 168, 71 (2018).


# Unveiling the early stage evolution of local atomic structures in the
crystallization process of a metallic glass

Lin Ma (马琳)1 These two authors contribute equally to this work. Xiao-Dong Yang
(杨晓东)111footnotemark: 1 Feng Yang (杨锋)1
Xin-Jia Zhou (周鑫嘉)1 Zhen-Wei Wu (武振伟)1
1Institute of Nonequilibrium Systems, School of Systems Science, Beijing
Normal University,
100875 Beijing, China Corresponding author Email<EMAIL_ADDRESS>
###### Abstract

The early stage evolution of local atomic structures in a multicomponent
metallic glass during its crystallization process has been investigated via
molecular dynamics simulation. It is found that the initial thermal stability
and earliest stage evolution of the local atomic clusters show no strong
correlation with their initial short-range orders, and this leads to an
observation of a novel symmetry convergence phenomenon, which can be
understood as an atomic structure manifestation of the ergodicity.
Furthermore, in our system we have quantitatively proved that the crucial
factor for the thermal stability against crystallization exhibited by the
metallic glass, is not the total amount of icosahedral clusters, but the
degree of global connectivity among them.

Keywords: metallic glass, crystallization, molecular dynamics simulation,
local atomic clusters

PACS: 64.70.pe, 68.55.A-, 64.60.aq

##  1 Introduction

Since the advents of multicomponent metallic glasses (MGs) in the early 1990s,
the initial atomic-scale structure evolution and nanocrystallization under
annealing or mechanical loading treatments of their glass-forming liquids,
have been extensively investigated for both fundamental and practical reasons
Wang2019Dynamic,
Sheng2006Atomic,Cheng2009Atomic,Pauly2010Transformation,Inoue2000Stabilization,Wang2011Atomic,Liu2008Atomistic,Wang2003In,Luo2004Icosahedral,Chen2006Molecular,Martin2004Nanocrystallization,DellaValle1994Microstructural,Kelton2003First,Cheng2011Atomic,Li2009Structural,
Wu2013Effect, Jiang2016Effect, Mokshin2008Shear,
Mokshin2010Crystal,Wu2021Quantitative,Cui2022Spatial. Various unique internal
structures of MGs underlie their interesting properties, which render MGs
potentially useful for distinct
applicationsLiu2008Atomic,Gao2021Hydrogen,Ye2016Generalized,Chen2014Atomic.
For the excellent mechanical properties of glasses with nanocrystalline
inclusions and as well fundamental meanings to solid state physics, numerous
studies have been conducted to comprehend how amorphous alloys crystallize at
nanoscales. However, due to the lack of directly time-resolved data of atomic
structure evolutions at the early stage of crystallization, clear pictures of
how localized atomic structures rearranging themselves, developing original
translational symmetry, and eventually forming nanocrystals with three-
dimensional (3-D) periodicity, remain less well-understood at present.
Especially the questions that which local parts of the MGs will the firstly
get into the structure with embryonic translational symmetry, what kind of
features these local parts have, and how they will evolve in the future, are
still mysterious at present.

Through high-resolution transmission electron microscopy (HRTEM) and advanced
structural analysis techniques, researchers have accumulated massive precious
experimental data to study and reveal the MGs’ structure evolution under
certain
treatmentsPauly2010Transformation,Inoue2000Stabilization,Wang2011Atomic,Chen2006Molecular,Martin2004Nanocrystallization.
A series of meaningful achievements have been obtained. But HRTEM cannot give
us the answers of these problems mentioned above. A clearly physical picture
that can describe the details of the evolutions associated with different
atomic-scale structures and the features of these localized structures has not
been obtained now. The quantitative description of this process is still in
its early stages.

The icosahedral short-range order (ISRO) in MGs is generally adopted as a
structural indicator for the fundamental process underlying structural
relaxationCheng2008Indicators,Wu2013Correlation,Li2017Five,Qiao2023A. However,
in some MGs the ISRO is absent. As shown in Ref. [30], atomic symmetry is a
more general concept in the glass-forming alloys in which both local fivefold
and translational symmetries are present. Both experimental and simulated
researches have shown that local clusters symmetry plays an important role in
determining many properties of MGs. Therefore, it is worth to check that
whether the degrees of the local fivefold symmetry (LFFS) and the local
translational symmetry (LTS) connected with the analysis of local icosahedral
clusters, especially the spatial topological order of them, are more crucial
for the understanding of the structural evolution and thermal stability of
MGs.

In this study, to address above challenges we performed molecular dynamics
(MD) simulation with the LAMMPS codePli95, and the degrees of LFFS and LTS
were used as the structural indicators to identify and characterize the
irreversible structural evolutions in the crystallization process. The spatial
correlation of the icosahedral clusters in the glass structures and its
influence on the thermal stability of MGs against crystallization has also
been investigated. We find that the initial stability and the very beginning
evolution of the local atomic structure are basically unaffected by the degree
of LFFS of these regions. Moreover, it is found that not the quantity of the
icosahedral clusters but the connectivity among them, will be the crucial
factor for the thermal stability against crystallization in amorphous alloys.

##  2 Methodology

In our MD simulations, a model system of Zr65Al30Cu5 containing $10000$ atoms
in a cubic box with periodic boundary conditions was adopted. The realistic
embedded-atom method potential used is the one developed based on the data in
Ref. [32]. The integration time step of all simulations was set to $2$ fs. In
the process of sample preparation, it was first melted and equilibrated at
$T=2500$K for $1$ ns, and then cooled down to $290$K with a cooling rate of
$10^{12}$K/s ($10^{13}$K/s), during which the cell size was adjusted to
maintain zero pressure in the NPT ensemble. Finally, the sample was relaxed at
$290$K for another $1$ ns. To ensure the reliability of the MD-generated
atomic configurations, the pair correlation functions (PCFs) of the model were
calculated to declare the amorphous states of the samples. For investigating
the crystallization process, annealing treatment was applied on the sample at
$T=650$K (above the kinetic glass transition temperature $T_{g}=615$K) in the
NPT ensemble for doing the following analysis.

##  3 Results and Discussion

Figure 1:  Fig. 1. (color online) The calculated PCFs of the MD-generated
atomic configurations obtained at a cooling rate of $10^{13}$K/s. For clarity,
the noisy and scattered Cu-Cu partial PCF is not included.

The PCFs of the amorphous Zr65Al30Cu5 alloy obtained through MD simulation are
presented in Fig. 1. The overall PCF exhibits a shoulder on the second peak,
and the Zr-Zr partial PCF displays the largest distance for the first peak,
consistent with the previous experimental findings on these
alloysYang2006Design. Additionally, the Cu-Zr and Zr-Zr partial PCFs show a
characteristic splitting of the second peak into two subpeaks, which is a
common feature observed for metallic glassesMendelev2007Using. These results,
to some extent, affirm the reliability of the MD-generated atomic
configurations for the subsequent analysis. Figure 2(a) presents the potential
energy (PE) evolution of the entire system as a function of scaled-time during
the isothermal annealing process. The onset of crystallization is evident by
the PEs exhibiting a drop at a scaled-time of $0.5\sim 0.6$. This decrease
signifies the initiation of crystallization. Notably, Figure 2(a) also
illustrates distinct variation trends in the PEs of two samples generated with
different cooling rates. It is worth to note that we do not emphasise the
relationship between the rate of cooling and the rate of crystallization here.
The two different cooling rates were used to produce two different
configurations for the subsequent analysis. In other words, the only purpose
here choosing two slightly different cooling rates is to generate two samples
that are similar in short-range orders but with a bit subtle difference in
medium-range orders (as it can be seen below). And we also performed the iso-
configurational runs to take the thermal fluctuation effect into account. It
was found that only $\sim$8 cases out of 40 that the fast-quenched sample took
much longer time to crystalize than the slow-quenched one, and the simply
isoconfigurational averaged crystallization-time (in practical defined as the
middle-term of the PE’s drops) of the slow-quenched sample is $3.6$ns which is
indeed statistically longer than the one of fast-quenched sample’s $2.8$ns.
More impressively, as it can be seen in Fig. 2(b), the distribution of the
crystallization-time for the slow-quenched sample is bimodal, in which there
is a peak (This peak remains after increasing the statistical sample size to
80, ruling out the possibility that this is a statistical rise and fall.)
located at long-time intervals, what is significantly different from the one
of the fast-quenched samples. Apparently, all above results indicate that the
initial atomic structure indeed plays an important role on the thermal
stability of our metallic liquid. The underlying reasons for this intriguing
phenomenon will be further explored in subsequent sections of this work.

Figure 2:  Fig. 2. (color online) (a) The annealing time dependence of the
system potential energy during the isothermal annealing treatment at $650$K.
For clarity, the time is scaled by $t_{a}$ which is the total annealing time
for the two samples respectively. (b) The distribution of isoconfigurational
averaged crystallization-time.

The atomic structure changes throughout the annealing process are indicated by
the temporal evolution of the PEs, as demonstrated by the evolving PCFs in
Fig. 3, which illustrate a transition from an amorphous state to a crystalline
solid. To comprehensively investigate the local structural rearrangement
preceding crystallization, we employed the Voronoi tessellation
analysisFinney1970Random,Borodin1999Local. This method involves dividing the
three-dimensional space into cells centered around each atom. A plane is drawn
to bisect each line connecting the central atom and its nearest-neighbors,
resulting in the formation of a Voronoi cell enclosed by these inner planes.
The Voronoi index $<n_{3},n_{4},n_{5},n_{6}>$ is utilized to characterize the
atomic clusters surrounding each atom, where $n_{i}$ represents the number of
$i$-edged faces of the corresponding Voronoi cell. The presence of an
$i$-edged face reflects the local symmetry of the central atom with its
nearest neighbors in a specific direction. Notably, faces with three, four,
and six edges exhibit LTS features, while pentagon faces indicate the presence
of LFFS. The degree of LFFS in local structures is quantified using a specific
metric defined as $d_{5}=n_{5}/\sum_{i}n_{i}$ Peng2011Structural.
Additionally, by analyzing the coordination information of atoms obtained from
the Voronoi analysis, the system’s global bond-angle distributions can be
measured. We selected a number of different annealing time periods and
analyzed the local structures of the corresponding atomic configurations. The
results shown in Fig. 4 reveal the changes in bond angle distributions, which
demonstrate a gradual increase in LTS features within the system throughout
annealing. The initial bond angle distribution clearly indicates the amorphous
state of the sample. As the annealing process progresses, characteristic
angles of $90^{\circ}$, $120^{\circ}$, and $180^{\circ}$ gradually emerge and
become more prominent, suggesting an enhancement of ordering in the system.
The variation in the fraction of $i$-edged faces in the Voronoi polyhedrons
(Fig. 4) also supports the aforementioned trend. The evolution of these local
symmetries behaves similarly in the samples with different cooling rates,
despite the fact that the sample generated with a lower cooling rate evolves
more gradually while the one generated with a higher cooling rate crystallizes
more suddenly after a bit gentle embryo state processing. Below we will
discuss this behavior in more details.

Figure 3:  Fig. 3. (color online) Changes in the total PCFs during the
crystallization of two MG samples obtained at different cooling rates. (a)
Cooling rate = $10^{13}$ K/s. (b) Cooling rate = $10^{12}$ K/s. Figure 4:
Fig. 4. (color online) The variations of the bond-angle distributions and
degree of LFFS and LTS of the MGs system during the annealing process.
(Cooling rate=$10^{13}$K/s, in (a) and (c); cooling rate=$10^{12}$K/s, in (b)
and (d).)

To provide a detailed description of the evolution of local structures with
different symmetry characteristics, we used the degree of LFFS $d_{5}$ in
Voronoi cells associated with local atomic clusters to label corresponding
regions and then examine the subsequent variation of $d_{5}$ of these regions
in the annealing process. As shown in Fig. 5, in contrast to some previous
observations where regions with a higher degree of LFFS exhibited greater
stability, and regions with lower degree of LFFS underwent more frequent
irreversible rearrangements, our results indicate that the thermal stability
and evolution speed of local structures are basically unaffected by the degree
of $d_{5}$ of these regions during annealing. Notably, we observed a
phenomenon of symmetry convergence (Fig. 5), where local structures with
different initial LFFS characteristics undergo transformations in their
localized environments, leading the entire system towards a state where the
degree of LFFS is approximately 0.4. After this convergence phenomenon,
$d_{5}$ begins to fall fast, indicating that the system is crystallising, the
growth and development of nanocrystals occur randomly within the samples.
Based on the aforementioned discussion, our aim is to elucidate the key
mechanism behind the nucleation of high-density nanosized crystals during the
initial stages of crystallization in MGs under supercooled liquid
conditionsWang2011Atomic,Liu2008Atomistic,Schroers1999Continuous,Schroers2000Crystallization,Pauly2010Criteria.
We propose that the convergence of local structure symmetry results in the
system transitioning into an embryo stage, a process that can be understood as
the atomic-structure manifestation of ergodicity, which makes sure that every
local region in the sample has the potential to nucleate and grow into
crystals from a liquid state. To further substantiate our viewpoints, we
provide a more detailed description of the evolution of local atomic
structures with dissimilar symmetries. Utilizing a sample with a cooling rate
of $10^{13}$ K/s, we employ the degree of $d_{5}$ to label different local
atomic regions. In addition to examining the variation in $d_{5}$ for the
corresponding atomic clusters, we also investigate the changes in triangular
($d_{3}=n_{3}/\sum_{i}n_{i}$), quadrilateral ($d_{4}=n_{4}/\sum_{i}n_{i}$),
and hexagonal faces ($d_{6}=n_{6}/\sum_{i}n_{i}$) with LTS features during the
evolution process (Fig. 6). The results further confirm that the thermal
stability and earliest evolution of atomic clusters are not correlated with
their original local symmetries.

Figure 5:  Fig. 5. (color online) The evolution of $d_{5}$ of different local
atomic regions classified by their corresponding original degree of $d_{5}$
(from $0.0$ to $1.0$) for the samples with two different cooling rate during
the annealing process. (a) Cooling rate=$10^{13}$K/s. (b)Cooling
rate=$10^{12}$K/s. $\Delta t$ is the configuration sampling time interval.
Figure 6:  Fig. 6. (color online) The change of $d_{3}$, $d_{4}$, $d_{6}$, and
$d_{5}$ of different local atomic regions classified by their corresponding
original degree of $d_{5}$ ((a)$d_{5}=1.0$, (b) $d_{5}=0.8$, (c)$d_{5}=0.6$,
(d)$d_{5}=0.4$, (e)$d_{5}=0.2$, (f)$d_{5}=0.0$.) for the sample with a cooling
rate of $10^{13}$K/s during the annealing process. The time interval $\Delta
t$ on the x-axis corresponds to $40$ ps.

After elucidating the initial-stage evolution of the local atomic clusters,
our focus now shifts to understanding the factors behind the disparate changes
observed in the PE curves which correspond to the differences in thermal
stability of samples produced at different cooling rates. Figure 7(a) presents
the average fraction of pentagonal faces, reflecting the LFFS of the entire
system, along with the presence of triangular, quadrilateral, and hexagonal
faces exhibiting the LTS features. The result suggests that the different
cooling rates we considered here have a negligible impact on the prevalence of
both LFFS and LTS in the corresponding samples. As well, the proportions of
icosahedral clusters $<0,0,12,0>$ in the two as-quenched samples, as well as
their variations during the annealing process were shown in Fig. 7(b). The
results unequivocally demonstrate that, despite the similar initial fractions
of $<0,0,12,0>$ in these two samples generated at different cooling rates,
their subsequent evolution over time exhibits substantial disparity.
Consequently, we naturally infer that the underlying mechanism accounting for
the aforementioned diversity lies in the distinctive packing and correlation
of icosahedral clusters within the glass structures. This investigation
addresses a crucial question regarding the connectivity and arrangement of
local icosahedral clusters, which profoundly influence the thermal stability
and crystallization process in these two samples.

Figure 7:  Fig. 7. (color online) (a) The fractions of i-edged faces in
polyhedra of two samples associated with different cooling rates. (b) The
variation of the fractions of the icosahedral clusters during the isothermal
annealing treatment.

To investigate the network and extended clusters (ECs) formed by the
interconnection of the $<0,0,12,0>$ clusters, here we introduced the method of
graph theory approach. The central atoms of all icosahedral clusters will be
treated as nodes in a network, and if two different central atoms have the
same one atom as their common nearest-neighbors, the two represented nodes are
considered to be connected. The number of the sharing atoms will be used as
the weight of this link. Now the question is transformed into a mathematical
problem, obtaining and analyzing the connected components (i.e. the ECs in the
glass structures) in a graph. The connectivity number ($N_{con}$) which
characterizes the space length scales of ECs is defined as the number of nodes
in a connected component in the corresponding graph. Distributions of
$N_{con}$ at several time points during the isothermal annealing process have
been shown in Fig. 8. It can be seen that the icosahedral clusters in the
samples generated under a lower cooling rate exhibit stronger spatial
correlations, and they tend to be connected with each other. The results
indicate that the low cooling rate will enhance the spatial correlations of
the icosahedral clusters and consequently increase the stability of the ECs
network in the sample, which would be the key factor of dynamic slowing down
to prevent the crystallization. In order to make different dimensional ECs
have comparability, we define strength factor for an EC as $\sum
w_{i}/N_{con}$, where $N_{con}$ is the number of nodes in an EC, and $w_{i}$
is the sum of weight of all links in this connected component abstracted for
the corresponding EC. The global strength factor of the network representing
the atomic structure will be defined as the sum of the strength factors over
all ECs. The outcome implies that the two adopted cooling rates lead to small
difference (about $7.3\%$, see Fig. 7) in the quantity of icosahedral clusters
in two samples, but large influence (about $27.4\%$, see Fig. 8) on the global
strength factor of the samples. The results demonstrate once more that the
thermal stability against crystallization do not depend on the total quantity
of the icosahedral clusters, but on the degree of connectivity among them in
our cases. As well, the variation of the system strength factor also has been
given in Fig. 8. It can be seen that systems subjected to different cooling
rates exhibit distinct evolving behavior during the crystallization process,
which indicates that the system strength factor is a good indicator for the
thermal stability of the MGs and strongly correlated with the evolutionary
time-scales in the crystallization process of them. And this also can be used
to rationalize what we found in Fig. 2(b), the bimodal feature, by the fact
that in the slow-quenched samples the stochastic nucleus’s development and the
subsequent crystal growth will be greatly influenced by how close it is to a
compact extended icosahedral-clusters, where the connectivity strength is
relatively high.

The crystallization is a complex physical process consisting of two major
events, nucleation and crystal growth which are driven by thermodynamic,
chemical, dynamic, as well as local structure properties. It has been well
documented in the literatures that the presence of icosahedral clusters and of
connected icosahedra plays an important role on the dynamic slowing-down of
the metallic liquids during supercoolingCheng2011Atomic,Wu2013Correlation.
Therefore, what we found here can be attributed to the fact that a more
compact icosahedral network will stabilize the dynamics of the liquids a bit
and then hinders the crystal growth process. At the same time, even more
interesting and complicated, the present of the icosahedra network possessing
lower formation energyWu2016critical could also increase the free energy
barrier of the nucleation process in the view of it makes the interface
between the crystal nuclei and the surrounding medium a bit rough (effectively
large interface area), which leads to a potentially more interface free energy
cost and then against the crystallizationSun2013effects,Desgranges2018unusual.
In other words, the increase of crystalline embryos and the decrease of
icosahedral connectivity during annealing seem to imply that there is a mutual
constraint between the two factors. On the one hand, with the growth of
crystalline embryos, the connectivity of icosahedra is damaged, which drives
the whole crystallization process; on the other hand, the connectivity of
icosahedral clusters inhibits the formation of crystalline embryos. Especially
when a much slower cooling rate was adopted, the crystal-like short-range
orders may start to organize to form some extended structures, which would
definitely influence the following crystallization process when the sample was
annealed at some temperatures. The interplay between the embryos of
crystallization and the connectivity of icosahedra on the thermal stability of
the liquid will become even more complicated at that
timeDesgranges2018unusual,Desgranges2019can, and the mechanism of the
interaction between these two factors in the metallic glass-formimng systems
is a thought-provoking issue that deserves an in-depth investigation in the
future.

Figure 8:  Fig. 8. (color online) The distributions of $N_{con}$ at several
time points ((a)$0.0$ ns, (b)$1.0$ ns, (c)$2.0$ ns.) during the isothermal
annealing process. (d) The change in system strength factors of two samples
with dissimilar cooling rate during the isothermal annealing process.

##  4 Conclusion

In this study, we provided a refined description of the evolutionary details
of the atomic structures with different local symmetries in the
crystallization process of a metallic glass. At the very beginning of the
devitrifaction, the results suggest that the thermal stability and the
earliest atomic structure evolution have no correlation with their short-range
order (SRO) characters, which leads to a novel phenomenon of local symmetry
convergence being found. We propose that this could be the key mechanism for
the nucleation of high-density nanocrystals at the beginning of the
crystallization process of MGs under the deep supercooled
conditionsWang2011Atomic,Liu2008Atomistic,Schroers1999Continuous,Schroers2000Crystallization,Pauly2010Criteria.
Furthermore, it has been clarified that in our system the thermal stability
against crystallization does not depend on the total amount of the icosahedral
SRO, but on the degree of global connectivity among them, which is a result
might promote future studies of the theory for nucleation and growth of
crystals in the disordered systems. At last, it is worth to note that the
global strength factor which quantifies the degree of global connectivity of
the icosahedral networks would also be an interesting alternative to
investigate other relevant physical processes of the
MGsLi2017Local,Wu2018Streched.

## Acknowledgment

Z.W.W. thanks the members of the “Beijing Metallic Glass Club” for the long-
term fruitful discussions. This work was supported by the National Natural
Science Foundation of China (Grant Nos. 52031016 and 11804027). The work was
carried out at National Supercomputer Center in Tianjin, and the calculations
were performed on Tianhe new generation supercomputer.

## References

  * [1] Wang W H 2019 Prog. Mater. Sci. 106 100561 ISSN 0079-6425 URL https://www.sciencedirect.com/science/article/pii/S0079642519300295
  * [2] Sheng H W, Luo W K, Alamgir F M, Bai J M and Ma E 2006 Nature 439(7075) 419–425 URL https://doi.org/10.1038/nature04421
  * [3] Cheng Y Q, Ma E and Sheng H W 2009 Phys. Rev. Lett. 102(24) 245501 URL https://link.aps.org/doi/10.1103/PhysRevLett.102.245501
  * [4] Pauly S, Gorantla S, Wang G, Kühn U and Eckert J 2010 Nat. Mater 9(6) 473–477 URL https://doi.org/10.1038/nmat2767
  * [5] Inoue A 2000 Acta Mater. 48 279–306 ISSN 1359-6454 URL https://www.sciencedirect.com/science/article/pii/S1359645499003006
  * [6] Wang Q, Liu C T, Yang Y, Dong Y D and Lu J 2011 Phys. Rev. Lett. 106(21) 215505 URL https://link.aps.org/doi/10.1103/PhysRevLett.106.215505
  * [7] Liu X, Chen G, Hou H, Hui X, Yao K, Lu Z and Liu C 2008 Acta Mater. 56 2760–2769 ISSN 1359-6454 URL https://www.sciencedirect.com/science/article/pii/S1359645408001444
  * [8] Wang X L, Almer J, Liu C T, Wang Y D, Zhao J K, Stoica A D, Haeffner D R and Wang W H 2003 Phys. Rev. Lett. 91(26) 265501 URL https://link.aps.org/doi/10.1103/PhysRevLett.91.265501
  * [9] Luo W K, Sheng H W, Alamgir F M, Bai J M, He J H and Ma E 2004 Phys. Rev. Lett. 92(14) 145502 URL https://link.aps.org/doi/10.1103/PhysRevLett.92.145502
  * [10] Chen G L, Liu X J, Hui X D, Hou H Y, Yao K F, Liu C T and Wadsworth J 2006 Appl. Phys. Lett. 88 203115 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.2198015/14660781/203115_1_online.pdf) URL https://doi.org/10.1063/1.2198015
  * [11] Martin I, Ohkubo T, Ohnuma M, Deconihout B and Hono K 2004 Acta Mater. 52 4427–4435 ISSN 1359-6454 URL https://www.sciencedirect.com/science/article/pii/S1359645404003118
  * [12] Della Valle R G, Gazzillo D, Frattini R and Pastore G 1994 Phys. Rev. B 49(18) 12625–12632 URL https://link.aps.org/doi/10.1103/PhysRevB.49.12625
  * [13] Kelton K F, Lee G W, Gangopadhyay A K, Hyers R W, Rathz T J, Rogers J R, Robinson M B and Robinson D S 2003 Phys. Rev. Lett. 90(19) 195504 URL https://link.aps.org/doi/10.1103/PhysRevLett.90.195504
  * [14] Cheng Y and Ma E 2011 Prog. Mater. Sci. 56 379–473 ISSN 0079-6425 URL https://www.sciencedirect.com/science/article/pii/S0079642510000691
  * [15] Li M, Wang C Z, Hao S G, Kramer M J and Ho K M 2009 Phys. Rev. B 80(18) 184201 URL https://link.aps.org/doi/10.1103/PhysRevB.80.184201
  * [16] Wu Z W, Li M Z, Wang W H, Song W J and Liu K X 2013 J. Chem. Phys. 138 074502 ISSN 0021-9606 (Preprint https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4792067/14047051/074502_1_online.pdf) URL https://doi.org/10.1063/1.4792067
  * [17] Jiang S Q, Wu Z W and Li M Z 2016 J. Chem. Phys. 144 154502 ISSN 0021-9606 (Preprint https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4946866/15511205/154502_1_online.pdf) URL https://doi.org/10.1063/1.4946866
  * [18] Mokshin A V and Barrat J L 2008 Phys. Rev. E 77(2) 021505 URL https://link.aps.org/doi/10.1103/PhysRevE.77.021505
  * [19] Mokshin A V and Barrat J L 2010 Phys. Rev. E 82(2) 021505 URL https://link.aps.org/doi/10.1103/PhysRevE.82.021505
  * [20] Wu Y, Xu B, Sun Y and Guan P 2021 Chin. Phys. B 30 57103 URL https://cpb.iphy.ac.cn/EN/abstract/article_123471.shtml
  * [21] Cui S, Liu H and Peng H 2022 Chin. Phys. B 31 86108 URL https://cpb.iphy.ac.cn/EN/abstract/article_124977.shtml
  * [22] Liu K X, Liu W D, Wang J T, Yan H H, Li X J, Huang Y J, Wei X S and Shen J 2008 Appl. Phys. Lett. 93 081918 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.2976667/13708787/081918_1_online.pdf) URL https://doi.org/10.1063/1.2976667
  * [23] Gao J A, Huang H S and Lv Y J 2021 Chin. Phys. B 30 66301 URL https://cpb.iphy.ac.cn/EN/abstract/article_123575.shtml
  * [24] Ye L M, Wu Z W, Liu K X, Tang X Z and Xiong X M 2016 Chin. Phys. B 25 068104 URL https://cpb.iphy.ac.cn/EN/abstract/article_118688.shtml
  * [25] Chen S Y, Wu Z W and Liu K X 2014 Chin. Phys. B 23 066802 URL https://dx.doi.org/10.1088/1674-1056/23/6/066802
  * [26] Cheng Y Q and Ma E 2008 Appl. Phys. Lett. 93 051910 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.2966154/14400695/051910_1_online.pdf) URL https://doi.org/10.1063/1.2966154
  * [27] Wu Z W, Li M Z, Wang W H and Liu K X 2013 Phys. Rev. B 88(5) 054202 URL https://link.aps.org/doi/10.1103/PhysRevB.88.054202
  * [28] Li M Z, Peng H L, Hu Y C, Li F X, Zhang H P and Wang W H 2017 Chin. Phys. B 26 016104 URL https://cpb.iphy.ac.cn/EN/abstract/article_119282.shtml
  * [29] Qiao Q, Wang J, Cai Z, Feng S, Song Z, Huo B, Li Z and Wang L M 2023 Chin. Phys. B 32 116401 URL https://cpb.iphy.ac.cn/EN/abstract/article_126217.shtml
  * [30] Peng H L, Li M Z and Wang W H 2011 Phys. Rev. Lett. 106(13) 135503 URL https://link.aps.org/doi/10.1103/PhysRevLett.106.135503
  * [31] Plimpton S 1995 J. Comput. Phys. 117 1–19 ISSN 0021-9991 URL https://www.sciencedirect.com/science/article/pii/S002199918571039X
  * [32] Zhou X W, Johnson R A and Wadley H N G 2004 Phys. Rev. B 69(14) 144113 URL https://link.aps.org/doi/10.1103/PhysRevB.69.144113
  * [33] Yang L, Xia J H, Wang Q, Dong C, Chen L Y, Ou X, Liu J F, Jiang J Z, Klementiev K, Saksl K, Franz H, Schneider J R and Gerward L 2006 Appl. Phys. Lett. 88 241913 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.2213020/13179151/241913_1_online.pdf) URL https://doi.org/10.1063/1.2213020
  * [34] Mendelev M I, Sordelet D J and Kramer M J 2007 J. Appl. Phys. 102 043501 ISSN 0021-8979 (Preprint https://pubs.aip.org/aip/jap/article-pdf/doi/10.1063/1.2769157/14865432/043501_1_online.pdf) URL https://doi.org/10.1063/1.2769157
  * [35] Finney J L and Bernal J D 1970 Proc. R. Soc. A: Math. Phys. Eng. Sci. 319 479–493 (Preprint https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1970.0189) URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1970.0189
  * [36] Borodin V A 1999 Phil. Mag. A 79 1887–1907 (Preprint https://doi.org/10.1080/01418619908210398) URL https://doi.org/10.1080/01418619908210398
  * [37] Schroers J, Busch R, Masuhr A and Johnson W L 1999 Appl. Phys. Lett. 74 2806–2808 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/74/19/2806/7809734/2806_1_online.pdf) URL https://doi.org/10.1063/1.124020
  * [38] Schroers J, Johnson W L and Busch R 2000 Appl. Phys. Lett. 77 1158–1160 ISSN 0003-6951 (Preprint https://pubs.aip.org/aip/apl/article-pdf/77/8/1158/7816107/1158_1_online.pdf) URL https://doi.org/10.1063/1.1289033
  * [39] Pauly S, Liu G, Gorantla S, Wang G, Kühn U, Kim D and Eckert J 2010 Acta Mater. 58 4883–4890 ISSN 1359-6454 URL https://www.sciencedirect.com/science/article/pii/S1359645410003010
  * [40] Wu Z, Li F, Huo C, Li M, Wang W and Liu K 2016 Sci. Rep. 6 35967 URL https://doi.org/10.1038/srep35967
  * [41] Sun G, Giovambattista N, Wang E and Xu L 2013 Soft Matter 9(47) 11374–11382 URL http://dx.doi.org/10.1039/C3SM52206G
  * [42] Desgranges C and Delhommelle J 2018 Phys. Rev. Lett. 120(11) 115701 URL https://link.aps.org/doi/10.1103/PhysRevLett.120.115701
  * [43] Desgranges C and Delhommelle J 2019 Phys. Rev. Lett. 123(19) 195701 URL https://link.aps.org/doi/10.1103/PhysRevLett.123.195701
  * [44] Li F X and Li M Z 2017 J. Appl. Phys. 122 225103 ISSN 0021-8979 (Preprint https://pubs.aip.org/aip/jap/article-pdf/doi/10.1063/1.5004262/15205851/225103_1_online.pdf) URL https://doi.org/10.1063/1.5004262
  * [45] Wu Z W, Kob W, Wang W H and Xu L 2018 Nat. Commun. 9(1) 5334 URL https://doi.org/10.1038/s41467-018-07759-w