# Virtual excitations and quantum correlations in ultra-strongly coupled
harmonic oscillators under intrinsic decoherence

Radouan Hab-arrih<EMAIL_ADDRESS>Laboratory of Theoretical Physics,
Faculty of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El Jadida,
Morocco  Ahmed Jellal<EMAIL_ADDRESS>Laboratory of Theoretical Physics,
Faculty of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El Jadida,
Morocco  Canadian Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7,
Canada  El Hassan El Kinani  Mathematical Modeling and Scientific Computing
Team, Department of Mathematics, Faculty of Sciences, Moulay Ismail
University, Meknes, Morocco

###### Abstract

We study the intrinsic decoherence of coupled harmonic oscillators. The
Milburn master equation is solved exactly, and the dynamics of virtual ground
state excitations are investigated. The interaction of quantum correlations
and virtual excitation was then studied. The following is a summary of our
major findings. (i) The damped oscillatory profile of all three quantities is
the same. (ii) Ultra-strong coupling combined with huge anisotropy values
results in the reemergence of entanglement and steering. (iii) To sustain
entanglement and steering, virtual excitations are required. (iv) The quantum
correlations are amplified in the quantum synchronous regime. (v) Ultra-strong
couplings cause inherent decoherence to be avoided.

###### pacs:

###### pacs:

03.65.Fd, 03.65.Ge, 03.65.Ud, 03.67.Hk
Keywords: Harmonic oscillator, intrinsic decoherence, entanglement, steering,
virtual excitations, ground state, ultra-strong coupling, anisotropy.

##  I Introduction

Entanglement is one of the amazing resource of quantum mechanics that have not
any analog in classical arena. In quantum information processing, such
resource must be preserved during all scenarios. Unfortunately, because of the
coherence losses (decoherence), entanglement dies out rapidly in finite time,
such phenomenon is called entanglement sudden death (ESD) [1]. Recently, it
was shown by engineering the physical parameters entanglement revives after
its ESD, the issue that brought some hope on quantum computing and the
promised quantum computer[2, 3]. In 1935, Einstein, Podolsky and Rosen (EPR),
proposed an intellectual experience to refute quantum entanglement for the
completeness of quantum theory (QT) [4]. Later on, Schrödinger introduced
another quantum correlation, the so-called quantum steering to support QT [5].
In the modern quantum information science (QIS), the quantum steering denotes
a quantum correlation situated between entanglement and Bell non locality [6].
It denotes the impossibility to describe the conditional states at one party
by a local hidden state model [7]. Moreover, recently, the steering is widely
discussed in several systems, including continous variables [8, 9, 10] and
discrete systems [11, 12].

Ultrastrong coupling (USC) physics plays a paramount role in several fields,
ranging from quantum optics to condensed matter [13, 14]. USC requires that
the coupling $J$ between the system parts be of the same order of magnitude of
the transition frequencies $\omega$ of the system
$(0.1<\frac{J}{\omega^{2}}<1)$. Furthermore, the USC regime has achieved
important records in several experimental setting, let’s quote, for example,
intersubband polaritons [15], superconducting circuits [16], Landau polaritons
[17] and optomechanics [18]. The importance of USC regime lies in its drastic
effects on some standard fundamental physical effects, including, for
instance, the Purcell effect [19], Zeno effect [20], photon blockade [21] and
excitations in the ground state [22].

The ground state excitations are one of the most amazing predictions of modern
quantum theory. In USC, the ground state is not empty, but populated with a
sea of virtual particles. These short-lived fluctuations are originating from
the couter-rotating terms of the Hamiltonian, because in USC the wave rotating
approximation (WRA) in not valid and anti-resonant terms should be kept in the
Hamiltonian [23, 24]. Those excitations are the background of some of the most
important physical processes in the universe, namely, the Casimir effect [25]
and Hawking radiation [26]. The main open theoretical challenge in USC is to
distinguish between virtual (unobservable) and physical (observable)
excitations [27]. Recently, ground state excitations have been extensively
studied, and it has been demonstrated that they cannot be detected [28, 29].
However, some authors have shown the possibility of indirect (without being
absorbed) detection, for instance, by measuring the change that they produce
in the Lamb shift of an ancillary probe qubit coupled to the cavity [30], or
by probing the radiation pressure that they generate in optomechanics systems
[31]. Although these excitations are not directly detectable, they have the
potential to spontaneously change into real excitations [29].

Intrinsic decoherence (ID) is a mechanism by which the system loses its
coherence because of its intrinsic degrees of freedom. Several models of ID
have been proposed [32, 33, 34], and one of them, is Milburn decoherence (MD)
model. Milburn proposed a simple modification to the Schrödinger equation
based on the assumption that the change in the state of the system, i.e.,
$\rho(t)\rightarrow\rho(t+\tau)$, is uncertain and occurs with a probability
of $0<p(\tau)<1$ on a sufficiently small time scale [35]. The change of state
is always certain in standard quantum mechanic, and $p(\tau)=1$. In the last
decade, the ID of quantum coupled discrete systems has been widely studied
[36, 37, 38, 39], but the majority of those works dealt with diffusion
approximation, i.e., first order correction of the von Neumann equation.

One of the important continuous systems is the quantum-coupled harmonic
oscillator (QCHO). It models, for example, coupled ions in ion traps [40],
coupled nano-sized electromechanical devices arrayed [41], and light
propagation in inhomogeneous media [42]. In the last decades, the
environmental decoherence of QCHO has been well studied in both regimes of
coupling weak and ultrastrong [43, 22], as well as Markovian and non-Markovian
evolution [44, 45]. Fortunately, it was shown that avoiding environmental
decoherence is possible by engineering only the frequencies and coupling
between oscillators [46]. The effect of MD on coupled oscillators beyond WRA
is not yet studied. In this paper we shed light on the effect of MD on the
virtual excitations of the ground state and their interconnection with quantum
entanglement and steering. Subsequently, we discuss how one can avoid
intrinsic decoherence by appropriately choosing the experimentally accessible
parameters.

The following is a breakdown of the current paper’s structure. The Hamiltonian
is diagonalized in the creation and annihilation operators in Sec. II. We
construct the covariance matrix and quantum correlations with virtual
excitations in Sec. III. In Sec. IV, the numerical results and discussions
will be presented. Finally, we present a quick summary of our findings

##  II Diagonalization and Milburn dynamics

Consider a system consisting of two coupled harmonic oscillators that can be
described by the Hamiltonian (in the unit $\hbar=m=1$)

$\displaystyle\hat{H}(\hat{x}_{1},\hat{x}_{2})=\frac{\hat{p}_{1}^{2}}{2}+\frac{\hat{p}_{2}^{2}}{2}+\frac{1}{2}\omega_{1}^{2}\hat{x}_{1}^{2}+\frac{1}{2}\omega_{2}^{2}\hat{x}_{2}^{2}-J\hat{x}_{1}\hat{x}_{2}$
(1)

Figure 1: (color online) Two harmonic oscillators coupled via position-
position type interaction $-J\hat{x}_{1}\hat{x}_{2}$, the system losses its
coherence via Milburn intrinsic decoherence characterized by the factor
$\Gamma$. The excitations in the ground state are virtual (not being absorbed)
and the intrinsic decoherence affects their generation.

To diagonalize the Hamiltonian, we perform a unitary transformation around the
$z$-axis performing a rotation $\mathcal{R}(\theta)$ a round $z-$axis, such
that

$\displaystyle\mathcal{R}(\theta)=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}$ (2)

of angle given by

$\displaystyle\theta=\frac{1}{2}\arctan\left(\frac{2g}{R^{4}-1}\right)$ (3)

where the dimensionless parameters $R=\sqrt{\omega_{1}/\omega_{2}}$ and
$g=J/\omega_{2}^{2}$ have been defined. It is worth noting that when
$\omega_{1}=\omega_{2}$, the rotation angle becomes $\theta=\frac{\pi}{4}$.
Now the Hamiltonian takes the following form

$\displaystyle\hat{H}_{D}(\hat{x}_{1},\hat{x}_{2})=\frac{\hat{p}_{1}^{2}}{2}+\frac{1}{2}\Omega_{1}^{2}\hat{x}_{1}^{2}+\frac{\hat{p}_{2}^{2}}{2}+\frac{1}{2}\Omega_{2}^{2}\hat{x}_{2}^{2}$
(4)

and the normal frequencies are (we assume that $\omega_{1}>\omega_{2}$ without
loss of generality)

$\displaystyle\Omega_{1}^{2}/\omega_{1}^{2}$
$\displaystyle=(1+R^{-4})/2+\frac{1}{2}\sqrt{(1-R^{-4})^{2}+4g^{2}R^{-8}}$ (5)
$\displaystyle\Omega_{2}^{2}/\omega_{2}^{2}$
$\displaystyle=(1+R^{4})/2-\frac{1}{2}\sqrt{(1-R^{4})^{2}+4g^{2}}.$ (6)

We introduce a pair of commuting operators

$\displaystyle
A=\tfrac{1}{2\sqrt{\Omega_{1}}}(\Omega_{1}\hat{x}_{1}+i\hat{p}_{1}),\qquad
B=\tfrac{1}{2\sqrt{\Omega_{2}}}(\Omega_{2}\hat{x}_{2}+i\hat{p}_{2})$ (7)

and then write Eq. (4) as

$\displaystyle\hat{H}_{D}(\hat{x}_{1},\hat{x}_{2})=\Omega_{1}A^{+}A+\Omega_{2}B^{+}B+\frac{\Omega_{1}+\Omega_{2}}{2}.$
(8)

The following quantum squeezers, expressed in terms the original commuting
operators $a$ and $b$

$\displaystyle\mathcal{S}_{w}(s_{j})=\exp[s_{j}(w^{2}-w^{+2})],\qquad w=a,b$
(9)

are used to diagonalize the Hamiltonian in the original creation and
annihilation operators, with $j=1,2$ and the squeezing parameters
$s_{j}=1/2\ln(\Omega_{j}/\omega_{j})$. As a result, we obtain

$\displaystyle\hat{H}_{d}=\Omega_{1}a^{+}a+\Omega_{2}b^{+}b+\frac{\Omega_{1}+\Omega_{2}}{2}.$
(10)

The Milburn master equation will be used to characterize the system’s
intrinsic decoherence [35]. This is

$\dot{\rho}(t)=\Gamma\left\\{\exp\left[-\frac{i}{\Gamma}\hat{H}\right]\rho(0)\exp\left[\frac{i}{\Gamma}\hat{H}\right]-\rho\right\\}$
(11)

where $\Gamma$ is the intrinsic damping factor. Its formal solution is given
by [47, 48]

$\displaystyle\rho(t,\Gamma)=\exp\left[-\Gamma
t\right]\sum_{k=0}^{\infty}\frac{(\Gamma
t)^{k}}{k!}\exp\left[-i\frac{k}{\Gamma}\hat{H}\right]\rho(0)\exp\left[i\frac{k}{\Gamma}\hat{H}\right].$
(12)

We note here that the decoherence is rather high for $\Gamma\rightarrow 0$ and
extremely weak for $\Gamma\rightarrow+\infty$ [35]. The goal of our work is to
investigate a fundamental feature of the ground state, namely virtual
excitation. Then, consider that oscillators are initially in their vacuum
state given by

$\displaystyle\rho(0)=|0\rangle_{1}\langle 0|\times|0\rangle_{2}\langle
0|=|00\rangle\langle 00|.$ (13)

As a matter of convenience, we set

$\displaystyle|\Phi_{k}(\Gamma)\rangle=\exp\left[-i\frac{k}{\Gamma}\hat{H}\right]|00\rangle$
(14) $\displaystyle 0\leq\varOmega_{k}(t,\Gamma)=\exp[-\Gamma t]\frac{(\Gamma
t)^{k}}{k!}\leq 1.$ (15)

Thus, Eq. (12) becomes

$\displaystyle\rho(t,\Gamma)=\sum_{k=0}^{\infty}\varOmega_{k}(t,\Gamma)|\Phi_{k}(\Gamma)\rangle\langle\Phi_{k}(\Gamma)|.$
(16)

Finally, the time evolution of the observable $\mathcal{O}$ is

$\displaystyle\langle{\mathcal{O}}\rangle(t)=\sum\limits_{k=0}^{+\infty}\varOmega_{k}(t,\Gamma)\langle\Phi_{k}(\Gamma)|\mathcal{O}|\Phi_{k}(\Gamma)\rangle.$
(17)

##  III Quantum correlations and virtual excitations

###  III.1 Covariance formalism

In our case, we consider two bosonic modes, $a$ and $b$, which satisfy the
bosonic algebra commutator $[\mathfrak{o},\mathfrak{o}^{+}]=\mathbb{I}$
$(\mathfrak{o}=a,b)$. We collect these operators in the vector
$\mathbb{A}=(\mathbb{A}_{n})=(a,a^{+},b,b^{+})^{T}$, while $T$ stands for the
transpose. The above commutation relation between the creation and
annihilation operators can be rewritten in the matrix form
$[\mathbb{A}_{n},\mathbb{A}_{m}]=i\mathbb{J}_{n,m}$. It is simple to write
$\mathbb{J}$ in the symplectic form shown below

$\displaystyle i\mathbb{J}_{4\times 4}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\tilde{\mathbb{{I}}}_{2\times 2}&0\\\
0&\tilde{\mathbb{{I}}}_{2\times
2}\end{array}\right),\qquad\tilde{\mathbb{{I}}}_{2\times
2}=\left(\begin{array}[]{cc}0&1\\\ -1&0\end{array}\right).$ (22)

A symplectic matrix $\mathbb{S}$ is one that verifies the following
properties: $\mathbb{S}^{+}\mathbb{J}\mathbb{S}=\mathbb{J}$ and
$\det(\mathbb{S})=1$. It is worth noting that the corresponding symplectic
matrix for a given unitary operator $\mathbb{U}$ can be derived using the fact
that $\mathbb{U}\mathbb{A}\mathbb{U}^{+}=\mathbb{S}\mathbb{A}$. To quantify
quantum correlations together with virtual excitations, we use the covariance
matrix $\Sigma$ defined by the compact form

$\displaystyle\Sigma_{n,m}$ $\displaystyle=$
$\displaystyle\langle\\{\mathbb{A}_{n},\mathbb{A}^{+}_{m}\\}\rangle_{\rho}-2\langle\mathbb{A}_{n}\rangle_{\rho}\langle\mathbb{A}^{+}_{m}\rangle_{\rho}.$
(23)

We mention that when
$\langle\mathbb{A}_{n}\rangle_{\rho}=\langle\mathbb{A}^{+}_{m}\rangle_{\rho}=0$,
the covariance matrix reduces to
$\Sigma_{n,m}=2\langle\mathbb{A}_{n}^{+}\mathbb{A}_{m}\rangle+1$. Here, all
the averages will be computed by using Eq. (17), and $\\{x,y\\}:=xy+yx$. It is
known that every quantifier of entanglement for two mode symmetric Gaussian
states is a function of the smallest symplectic eigenvalue of the partial
transpose [49].

###  III.2 Symplectic form of unitary transformations

In the annihilation and creation representation, the rotation unitary operator
Eq. (27) can be reduced to

$\displaystyle\mathcal{R}(\theta)=\exp\left[\tfrac{\theta}{2}\tfrac{1-R^{2}}{R}(a^{+}b^{+}-ab)-\tfrac{\theta}{2}\tfrac{1+R^{2}}{R}(a^{+}b-ab^{+})\right].$
(24)

We show the results

$\displaystyle\mathcal{R}(\theta)a\mathcal{R}^{-1}(\theta)$
$\displaystyle=a\cos\theta+\tfrac{1}{2}b\sin(\theta)(R+R^{-1})+\tfrac{1}{2}b^{+}\sin(\theta)(R-R^{-1})$
(25) $\displaystyle\mathcal{R}(\theta)b\mathcal{R}^{-1}(\theta)$
$\displaystyle=b\cos\theta-\tfrac{1}{2}a\sin(\theta)(R+R^{-1})+\tfrac{1}{2}a^{+}\sin(\theta)(R-R^{-1}).$
(26)

As a result, the symplectic form $\mathbb{S}_{\mathcal{R}}$ of
$\mathcal{R}(\theta)$ is obtained

$\displaystyle\mathbb{S}_{\mathcal{R}}(\theta)=\begin{pmatrix}\cos\theta&0&\frac{R+R^{-1}}{2}\sin\theta&\frac{R-R^{-1}}{2}\sin\theta\\\
0&\cos\theta&\frac{R-R^{-1}}{2}\sin\theta&\frac{R+R^{-1}}{2}\sin\theta\\\
-\frac{R+R^{-1}}{2}\sin\theta&\frac{R-R^{-1}}{2}\sin\theta&\cos\theta&0\\\
\frac{R-R^{-1}}{2}\sin\theta&-\frac{R+R^{-1}}{2}\sin\theta&0&\cos\theta\\\
\end{pmatrix}.$ (27)

Subsequently, one can rewrite the squeezors defined in Eq. (9) in their
symplectic form

$\displaystyle\mathbb{S}_{\mathcal{S}}(s_{1},s_{2})=\begin{pmatrix}\cosh
s_{1}&\sinh s_{1}&0&0\\\ \sinh s_{1}&\cosh s_{1}&0&0\\\ 0&0&\cosh s_{2}&\sinh
s_{2}\\\ 0&0&\sinh s_{2}&\cosh s_{2}\end{pmatrix}.$ (28)

Thus, we show that the symplectic transformation corresponding to
$\exp\left[-i\frac{k}{\Gamma}H_{d}\right]$ takes the following form

$\displaystyle\mathbb{S}_{H}(\Gamma)=\text{diag}\left(\exp\left[-i\frac{k}{\Gamma}\Omega_{1}\right],\exp\left[i\frac{k}{\Gamma}\Omega_{1}\right],\exp\left[-i\frac{k}{\Gamma}\Omega_{2}\right],\exp\left[i\frac{k}{\Gamma}\Omega_{2}\right]\right)$
(29)

###  III.3 Virtual excitations and quantum correlations

To compute virtual excitations and quantum correlations encoded in our state,
we use the covariance matrix formalism. It is not difficult to show that the
covariance matrix $\Sigma$ is

$\displaystyle\Sigma(t,\Gamma)=\exp[-\Gamma
t]\sum\limits_{k=0}^{\infty}\frac{(\Gamma
t)^{k}}{k!}\mathcal{H}(\theta,\Gamma,R)\Sigma(0)\mathcal{H}^{\dagger}(\theta,\Gamma,R)$
(30)

where the symplectic matrix $\mathcal{H}$ is defined as

$\displaystyle\mathcal{H}=\mathbb{S}_{\mathcal{R}}(-\theta,R)\mathbb{S}_{\mathcal{S}}(s_{1},s_{2})\mathbb{S}_{H}(\Gamma)\mathbb{S}_{\mathcal{S}}(-s_{1},-s_{2})\mathbb{S}_{\mathcal{R}}(\theta,R)$
(31)

and $\Sigma(0)$ is the covariance matrix corresponding to vacuum state
$\rho(t=0)=|00\rangle\langle 00|$

$\displaystyle\Sigma(0)=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix}$ (32)

with $x^{\dagger}$ denotes the hermitian conjugate of $x$. After computation,
we end up with

$\displaystyle\Sigma(t,\Gamma)=\begin{pmatrix}\sigma_{11}&\sigma_{12}&\sigma_{13}&\sigma_{23}^{\ast}\\\
\sigma_{12}^{\ast}&\sigma_{11}&\sigma_{23}&\sigma_{13}^{\ast}\\\
\sigma_{13}^{\ast}&\sigma_{23}^{\ast}&\sigma_{33}&\sigma_{34}\\\
\sigma_{23}&\sigma_{13}&\sigma_{34}^{\ast}&\sigma_{33}\end{pmatrix}.$ (33)

It is worthwhile to mention that the inputs of the covariance matrix are too
long, and we avoid writing them out here. Note that they can be obtained by
performing the matrix product of the matrices Eqs. (27-30). Furthermore, for
isotropic oscillators $R=1$, $\Sigma(t,\Gamma)$ reduces to

$\displaystyle\Sigma(t,\Gamma)=\begin{pmatrix}\sigma_{11}&\sigma_{12}&\sigma_{13}&\sigma_{23}^{\ast}\\\
\sigma_{12}^{\ast}&\sigma_{11}&\sigma_{23}&\sigma_{13}\\\
\sigma_{13}&\sigma_{23}^{\ast}&\sigma_{11}&\sigma_{12}\\\
\sigma_{23}&\sigma_{13}&\sigma_{12}^{\ast}&\sigma_{11}\end{pmatrix}$ (34)

where the explicit expressions of the covariance inputs are

$\displaystyle\sigma_{11}=$
$\displaystyle\frac{1}{2}{\sf{ch_{1}}}^{2}{\sf{sh_{1}}}^{2}\left(2-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)$
(35)
$\displaystyle+\frac{1}{2}{\sf{ch_{2}}}^{2}{\sf{sh_{2}}}^{2}\left(2-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)+1$
$\displaystyle\sigma_{12}=$
$\displaystyle{\sf{ch_{1}}}^{3}{\sf{sh_{1}}}\left(-2+\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma})\right)-1\right)\right]+\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)$
(36) $\displaystyle+{\sf{ch_{2}}}^{3}{\sf{sh_{2}}}\left(-2+\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]+\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)$
$\displaystyle+{\sf{ch_{1}}}{\sf{sh_{1}}}\left(1-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)+{\sf{ch_{2}}}{\sf{sh_{2}}}\left(1-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)$
$\displaystyle\sigma_{13}=$
$\displaystyle{\sf{ch_{2}}^{2}}{\sf{sh_{2}}^{2}}\left(2-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)$
(37)
$\displaystyle+\sf{{ch_{1}}^{2}}{\sf{sh_{1}}^{2}}\left(-2+\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]+\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)$
$\displaystyle\sigma_{23}=$
$\displaystyle{\sf{ch_{1}}}^{3}{\sf{sh_{1}}}\left(2-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)$
(38) $\displaystyle+{\sf{ch_{2}}}^{3}{\sf{sh_{2}}}\left(-2+\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]+\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)$
$\displaystyle-{\sf{ch_{1}}}{\sf{sh_{1}}}\left(1-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)+{\sf{ch_{2}}}{\sf{sh_{2}}}\left(1-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right)$

where $\sf{ch_{j}}=\cosh(s_{j})$ and $\sf{sh_{j}}=\sinh(s_{j})$ have been
defined. The virtual excitations in both modes are

$\displaystyle\langle N_{1}\rangle=\langle
a^{+}a\rangle=\tfrac{1}{2}\left(\sigma_{11}-1\right),\qquad\langle
N_{2}\rangle=\langle b^{+}b\rangle=\tfrac{1}{2}\left(\sigma_{22}-1\right).$
(39)

It is straightforward to show that the virtual excitations in both oscillators
are equal and reduce to

$\displaystyle\langle a^{+}a\rangle(t)=\langle b^{+}b\rangle(t)=$
$\displaystyle\frac{1}{2}{\sf{ch_{1}}}^{2}{\sf{sh_{1}}}^{2}\left(2-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{1}}{\Gamma}\right)-1\right)\right]\right)$
(40)
$\displaystyle+\frac{1}{2}{\sf{ch_{2}}}^{2}{\sf{sh_{2}}}^{2}\left(2-\exp\left[\Gamma
t\left(\exp\left(\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]-\exp\left[\Gamma
t\left(\exp\left(-\frac{2i\Omega_{2}}{\Gamma}\right)-1\right)\right]\right).$

Both the single reduced modes and the correlation sub-matrix will be derived
by using the trace-out prescription and, as a consequence, we obtain

$\displaystyle\Sigma^{(a)}=\Sigma^{(b)}=\left(\begin{array}[]{cc}\sigma_{11}&\sigma_{12}\\\
\sigma_{12}^{\ast}&\sigma_{11}\end{array}\right),\qquad\Sigma^{(ab)}=\left(\begin{array}[]{cc}\sigma_{13}&\sigma_{23}^{\ast}\\\
\sigma_{23}&\sigma_{13}\end{array}\right).$ (45)

To evaluate the entanglement, we compute the partial transposition
$\tilde{\Sigma}$ and obtain the four eigenvalues

$\displaystyle\tilde{\nu}_{m,M}^{2}=\frac{1}{2}\left(\delta\mp\sqrt{\delta^{2}-4\det\Sigma}\right)$
(46)

where the seralian $\delta$ is defined by

$\displaystyle\delta=\det\Sigma^{(a)}+\det\Sigma^{(b)}-2\det\Sigma^{(ab)}.$
(47)

Consequently, the logarithmic negativity becomes

$\displaystyle E_{N}=\max\left(0,-\ln(\tilde{\nu}_{m})\right).$ (48)

Distinct from entanglement, quantum steering is a magic quantum correlation
that is generally asymmetric. Thus, one needs to define the steering in both
ways [50]

$\displaystyle S^{A\rightarrow
B}=\max\left(0,\log_{2}\sqrt{\frac{\det\Sigma^{(a)}}{4\det\Sigma}}\right),\qquad
S^{B\rightarrow
A}=\max\left(0,\log_{2}\sqrt{\frac{\det\Sigma^{(b)}}{4\det\Sigma}}\right)$
(49)

and the steering asymmetry is

$\displaystyle\Delta S=|S^{A\rightarrow B}-S^{B\rightarrow A}|.$ (50)

We mention that the steering is symmetric in the isotropic case since
$S^{A\rightarrow B}=S^{B\rightarrow A}$.

##  IV Numerical results

###  IV.1 Avoiding intrinsic decoherence

In Fig. 2, we show in a comparative study the relationship between Milburn and
von Neumann dynamics. For this reason, we plot the evolution of the virtual
excitations for different values of ultra-strongly coupling $J$. As expected,
the dynamics of the virtual excitations under von Neumann dynamics exhibit an
unquenched oscillatory behavior. We set the decoherence rate $\Gamma=100$, the
excitations undergo a quenched oscillatory behavior that becomes very slow
with the increase of ultra-strongly coupling. Subsequently, near the
hermitianity point where the ultra-strong coupling reaches its maximum, the
dynamics overlap each other. Therefore, we conclude that the intrinsic
decoherence of a quantum system can be removed by ultra-strongly coupling its
internal degrees of freedom.

Figure 2: (color online) Milburn and von Neumann dynamics under an ultra-
strongly coupling effect with the blue line represents von Neumann dynamics
and the red line, Milburn dynamics. We choose $\Gamma=100$ and
$\omega_{1}=\omega_{2}=1$.

###  IV.2 Role of anisotropy

In Fig. 3, we investigate the behaviour of quantum entanglement $E_{N}$,
steering ($S^{A\rightarrow B}$ and $S^{B\rightarrow A}$) and virtual
excitations ($\langle N\rangle_{1}$ and $\langle N\rangle_{2}$) versus scaled
time for $\Gamma=100$, $J=0.2$. Before analyzing our results, we recall that
for weak coupling, the ground state is empty and since the coupling takes
important values, namely, the ultra-strong coupling regime, the ground state
becomes populated. According to our previous results [51], these excitations
are linked to the classical instability of classical oscillators and maintain
entanglement between oscillators. Now, in the first line of figures we address
the dynamics in the resonant case $\omega_{1}=\omega_{2}=1$. Excitations and
one-way steering will be the same, as expected. Subsequently, we observe that
dynamics rapidly revives excitations that will maintain entanglement and
steering between oscillators. In addition, we mention that due to intrinsic
decoherence, the three quantities undergo an exponential decay to a steady
value except for steering which suddenly dies out. This shows that
entanglement and excitations are robust against decoherence. This profile is
mathematically due to the intrinsic decoherence factor $\exp[-\Gamma t]$.
Subsequently, we slightly move away from resonance by taking $\omega_{2}=0.95$
and maintaining all other parameters. At first glance, we amazingly notice the
redistribution of excitations as well as their distinguishability $D=\langle
N_{1}\rangle-\langle N_{2}\rangle\neq 0$. As a result, the asymmetry in
quantum steering $\Delta S\neq 0$ and a soft decrease in entanglement compared
to the former case. Thirdly, we vary $\omega_{2}$ to $0.7$ and $0.3$, the
distinguishability $D$ increases, and as a consequence, the profile of
entanglement and steering shows enhancement immediately after their
generation, but because of decoherence, they fastly die out. For completeness,
reduce $\omega_{2}$ to $0.21$ and observe that excitations result in a quantum
synchronisation regime [52]. As a result, the quantum correlations are
significantly enhanced and revived for a long time. Finally, it is worthwhile
to mention the hierarchy of quantum correlations [8] and excitations. From the
above results, we have $\langle N_{1}\rangle=\langle N_{2}\rangle=0$ showing
$E_{N}=0$ and therefore $S^{A\rightarrow B}=S^{B\rightarrow A}=0$. This
hierarchy relation means, for instance, that the establishment of correlations
necessarily entails the appearance of excitations. This will allow a new route
to construct new quantifiers based on excitations [22, 51]. Another issue that
deserves attention is the dependence of the one-way steerability on the
angular velocity $\omega_{1,2}$. More clearly, the steerabilty of a rapid
oscillator from a slow one always exceeds the inverted steerabilty. This may
be due to the obvious distinguishability of excitations. To conclude, moving
away from resonance (increasing the anisotropy
$R=\sqrt{{\omega_{1}}/{\omega_{2}}}$) has the following implications: (i) A
generation of excitations followed by a redistribution phenomenon. A similar
phenomenon was observed in a system made of two oscillators interacting with a
qubit [53]. (ii) The redistribution of excitations implies the redistribution
of quantum steering. For instance, $S^{A\rightarrow B}<S^{B\rightarrow A}$ if
only and if $\langle N_{1}\rangle<\langle N_{2}\rangle$, this is consistent
with recent results of magnon-photon steering in ferromagnetic and anti-
ferromagnetic sub-lattices [9]. (iii) The activation of virtual excitation
generation, which results in the preservation and enhancement of quantum
correlations between oscillators. (iv) The system enters its quantum
synchronous regime, at which point excitations and correlations become
significant. This is similar to what has been discovered in other works, such
as the qubit-qubit system [54, 55] and oscillators [56, 57].

Figure 3: (color online) In the first column, we plot the dynamics of
excitations in both oscillators. In the second column, we plot the dynamics of
quantum steering for both cases. In the third one, we present the dynamics of
entanglement. We choose $\Gamma=100$, $\omega_{1}=1$ and $J=0.2$.

###  IV.3 Effect of ultra-strong coupling

In Fig. 4, we investigate the dynamics of the three quantities and their
interconnection by varying the ultra-strong coupling $J$ for $\Gamma=100$,
$\omega_{1}=1$, $\omega_{2}=0.5$. We start with $J=0.1$, the
distinguishability is evident and $\langle N_{1}\rangle<\langle N_{2}\rangle$
for all time. As a result, the one-way steering $S^{A\rightarrow B}$ is more
fragile than the second steering $S^{B\rightarrow A}$. By gradually increasing
the ultra-strong coupling, we observe that the excitations are significantly
generated, as well as a revival of one-way steering $S^{A\rightarrow B}$. The
two oscillators will become synchronous and the quantum correlations will be
greatly enhanced with a large coupling $J\rightarrow
J_{max}=\sqrt{\omega_{1}\omega_{2}}$. To summarize, a large coupling constant
$J$ combined with a large anisotropy constant $R$ can result in a significant
number of quantum correlations even in the presence of decoherence and over
long simulation times.

Figure 4: (color online) In the first column, we plot the dynamics of
excitations in both oscillators. In the second column, we plot the dynamics of
quantum steering for both cases. In the third one, we present the dynamics of
entanglement. We choose $\Gamma=100$, $\omega_{1}=1$ and $\omega_{2}=0.5$.

##  V Conclusion

We have investigated the dynamics of quantum entanglement, quantum steering,
and their interconnection with virtual excitations generated by the counter-
rotating terms of the Hamiltonian. The physical system of interest is a system
made of two coupled harmonic oscillators coupled via
$J\hat{x}_{1}\hat{x}_{2}$. The system was decoupled and the diagonalisation of
the Hamiltonian was obtained after appropriate transformations. The initial
density matrix is considered to be a separable ground state, namely
$|00\rangle\langle 00|$. Subsequently, we have exactly resolved the Milburn
master equation to end up with the covariance matrix.

Our findings show that intrinsic decoherence dramatically alters the profile
of the oscillatory activity experienced by the three quantities. We’ve
established that the system has a steady state of entanglement and virtual
excitations due to dissipation. It is shown that steering is less stable than
entanglement. Then, we have demonstrated that coupling ultra-strongly with
large anisotropy synchronizes the excitations between the harmonic
oscillators. As a result, correlations were found to survive for long periods
of time in simulations. Virtual excitations are shown to have a large
influence on quantum correlations, which is an issue that can be employed in
quantum information science to quantify, regulate, and process quantum
resources utilizing virtual excitations. We have demonstrated the feasibility
of preventing inherent decoherence in quantum systems by modifying the
experimentally available ultra-strong coupling value.

## References

  * [1] T. Yu and J. H. Eberly, Phys. Rev. B 66, 193306 (2002); ibid, Phys. Rev. Lett. 93, 140404 (2004).
  * [2] J. C. Gonzalez-Henao, E. Pugliese, S. Euzzor, S. F. Abdalah, R. Meucci, and J. A. Roversi, Sci. Rep. 5, 13152 (2015).
  * [3] J. C. Gonzalez-Henao, E. Pugliese, S. Euzzor, R. Meucci, J. A. Roversi, and F. T. Arecchi, Sci. Rep. 7, 9957 (2017).
  * [4] A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. 47, 777 (1935).
  * [5] E. Schrödinger, Math. Proc. Camb. Philos. Soc. 31, 555 (1935).
  * [6] R. Uola, A. C. S. Costa, H. C. Nguyen, and O. Gühne, Rev. Mod. Phys. 92, 15001 (2020).
  * [7] H. M. Wiseman, S. J. Jones, and A. C. Doherty, Phys. Rev. Lett. 98, 140402 (2007).
  * [8] H. S. Qureshi, S. Ullah, and F. Ghafoor, Sci. Rep. 8, 16288 (2018).
  * [9] S. S. Zheng, F. X. Sun, H. Y. Yuan, Z. Ficek, Q. H. Gong, and Q. Y. He, Sci. China Phys. Mech. Astron. 64, 210311 (2021).
  * [10] V. Le Duc, J. K. Kalaga, W. Leoński, M. Nowotarski, K. Gruszka, and M. Kostrzewa, Symmetry 13, 2201 (2021).
  * [11] D. Yao. and C. Ren, Phys. Rev. A 103, 052207 (2021).
  * [12] W. Y. Sun, D. Wang, J. D. Shi, and L. Ye, Sci. Rep. 7, 39651 (2017).
  * [13] S. De Liberato, Nat. Commun. 8, 1465 (2017).
  * [14] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * [15] B. Askenazi et al., ACS Photonics 4, 2550 (2017).
  * [16] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Phys. Rev. A 95, 053824 (2017).
  * [17] C. Maissen et al., Phys. Rev. B 90, 205309 (2014).
  * [18] F. Benz et al., Science 354, 726 (2016).
  * [19] S. De Liberato, Phys. Rev. Lett. 112, 016401 (2014).
  * [20] I. Lizuain, J. Casanova, J. J. García-Ripoll, J. G. Muga, and E. Solano, Phys. Rev. A 81, 062131 (2010).
  * [21] A. Le Boité, M. J. Hwang, H. Nha, and M. B. Plenio, Phys. Rev. A 94, 033827 (2016).
  * [22] J.-Y. Zhou, Y.-H. Zhou, X.-L. Yin, J.-F. Huang, and J.-Q. Liao, Sci. Rep. 10, 12557 (2020).
  * [23] S. De Liberato, C. Ciuti, and I. Carusotto, Phys. Rev. Lett. 98, 103602 (2007).
  * [24] T. Niemczyk et al., Nat. Phys. 6, 772 (2010).
  * [25] J. R. Johansson, G. Johansson, C. M. Wilson, and F. Nori, Phys. Rev. Lett. 103, 147003 (2009).
  * [26] J. Steinhauser, Nat. Phys. 12, 959 (2016).
  * [27] O. Di Stefano, R. Stassi, L. Garziano, A. F. Kockum, S. Savasta, and Franco Nori, New J. Phys. 19, 053010 (2017).
  * [28] A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett. 109, 193602 (2012).
  * [29] R. Stassi, A. Ridolfo, O. Di Stefano, M. J. Hartmann, and S. Savasta, Phys. Rev. Lett. 110, 243601 (2013).
  * [30] J. Lolli, A. Baksic, D. Nagy, V. E. Manucharyan, and C. Ciuti, Phys. Rev. Lett. 114, 183601 (2015).
  * [31] M. Cirio, K. Debnath, N. Lambert, and F. Nori, Phys. Rev. Lett. 119, 053601 (2017).
  * [32] G. C. Ghirardi, A. Rimini, and T. Weber, Phys. Rev. D 34, 470 (1986).
  * [33] L. Diosi, Phys. Rev. A 40, 1165 (1989).
  * [34] J. Ellis, S. Mohanty, and D. V. Nanopoulos, Phys. Lett. B 235, 305 (1990).
  * [35] G. J. Milburn, Phys. Rev. A 44, 9 (1991).
  * [36] S. J. Anwar, M. Ramzan, and M. K. Khan, Quantum Inf. Process. 16, 142 (2017).
  * [37] L. Zheng and G. F. Zhang, Eur. Phys. J. D 71, 288 (2017).
  * [38] Y. Y. Liao, S. R. Jian, and J. R. Lee, Eur. Phys. J. D 73, 47 (2019).
  * [39] A. S. F. Obada and A. B. A. Mohamed, Opt. Commun. 309, 236 (2013).
  * [40] K. R. Brown, C. Ospelkaus, Y. Colombe, A. C. Wilson, D. Leibfried, and D. J. Wineland, Nature 471, 196 (2011).
  * [41] J. Eisert, M. B. Plenio, S. Bose, and J. Hartley, Phys. Rev. Lett. 93, 190402 (2004).
  * [42] A. R. Urzúa, I. Ramos-Prieto, F. Soto-Eguibar, V. Arrizón, and Héctor M. Moya-Cessa, Sci. Rep 9, 16800 (2019).
  * [43] C. Joshi, P. Öhberg, J. D. Cresser, and E. Andersson, Phys. Rev. A 90, 063815 (2014).
  * [44] S. Xue and I. R. Petersen, Quantum Inf. Process. 15, 1001 (2016).
  * [45] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, New J. Phys. 18, 073007 (2016).
  * [46] G. Manzano, F. Galve, and R. Zambrini, Phys. Rev. A 87, 032114 (2013).
  * [47] X. Jing-Bo, Z. Xu-Bo, G. Xiao-Chun, and F. Jian, Commun. Theor. Phys. 37, 733 (2002).
  * [48] A. R. Urzúa and H. M. Moya-Cessa, Pramana 96, 72 (2022).
  * [49] G. Adesso, S. Ragy, and A. R. Lee, Open Syst. Inf. Dyn. 21, 1440001 (2014).
  * [50] I. Kogias, A. R. Lee, S. Ragy, and G. Adesso, Phys. Rev. Lett. 114, 060403 (2015).
  * [51] R. Hab-arrih, A. Jellal, D. Stefanatos, and A. Merdaci, Quantum Rep. 3, 684 (2021).
  * [52] Gj. Qiao, Hx. Gao, Hd. Liu, and X. X. Yi, Sci. Rep. 8, 15614 (2018).
  * [53] B. Militello, H. Nakazato, and A. Napoli, Phys. Rev. A 96, 023862 (2017).
  * [54] A. Roulet and C. Bruder, Phys. Rev. Lett. 121, 063601 (2018).
  * [55] O. V. Zhirov and D. L. Shepelyansky, Phys. Rev. B 80, 014519 (2009).
  * [56] C. G. Liao, R. X. Chen, H. Xie, M. Y. He, and X. M. Lin, Phys. Rev. A 99, 033818 (2019).
  * [57] S. Lorenzo, B. Militello, A. Napoli1, R. Zambrini, and G. M.Palma, New J. Phys. 24, 023030 (2022).


# Quantifying Disorder One Atom at a Time Using an Interpretable Graph Neural
Network Paradigm

James Chapman Corresponding Author<EMAIL_ADDRESS>Materials Science
Division, Lawrence Livermore National Laboratory, Livermore, CA, USA  Tim Hsu
Corresponding Author<EMAIL_ADDRESS>Center for Applied Scientific Computing,
Lawrence Livermore National Laboratory, Livermore, CA, USA  Xiao Chen  Center
for Applied Scientific Computing, Lawrence Livermore National Laboratory,
Livermore, CA, USA  Tae Wook Heo  Materials Science Division, Lawrence
Livermore National Laboratory, Livermore, CA, USA  Brandon C. Wood
Corresponding Author<EMAIL_ADDRESS>Materials Science Division, Lawrence
Livermore National Laboratory, Livermore, CA, USA

###### Abstract

Quantifying the level of atomic disorder within materials is critical to
understanding how evolving local structural environments dictate performance
and durability. Here, we leverage graph neural networks to define a physically
interpretable metric for local disorder. This metric encodes the diversity of
the local atomic configurations as a continuous spectrum between the solid and
liquid phases, quantified against a distribution of thermal perturbations. We
apply this novel methodology to three prototypical examples with varying
levels of disorder: (1) solid-liquid interfaces, (2) polycrystalline
microstructures, and (3) grain boundaries. Using elemental aluminum as a case
study, we show how our paradigm can track the spatio-temporal evolution of
interfaces, incorporating a mathematically defined description of the spatial
boundary between order and disorder. We further show how to extract physics-
preserved gradients from our continuous disorder fields, which may be used to
understand and predict materials performance and failure. Overall, our
framework provides an intuitive and generalizable pathway to quantify the
relationship between complex local atomic structure and coarse-grained
materials phenomena.

##  1 Introduction

Figure 1: General workflow for calculating SODAS values. Atomic structures are
converted into graph representations, which explicitly encode all necessary
geometric information. These atomic graphs are then fed into a graph neural
network, which has been trained to distinguish between the unique local
geometries in different material phases. The graph neural network then gives
each atomic environment a SODAS value, which classifies where in the phase
space, between phases, that local structure is most likely to occur.

Understanding how a material’s structure affects its properties is one of the
most fundamental principles in materials science. At the center of this
paradigm is the fact that macroscopic materials behavior begins at the atomic
scale, with local atomic arrangements ultimately coming together to form
structural features observed at larger length scales [1, 2]. Characterizing
the nature and propagation of these local environments is therefore vital to
understanding macroscale structure-property relationships and their evolution
[3]. Complicating this endeavor is the fact that the long-range features often
depend on structurally disordered atomic environments, which tend to dictate
materials functionality [4]. For instance, transport, chemical reactivity, and
phase nucleation are all profoundly affected by the presence of interfaces,
interphases, and grain boundaries [5, 6, 7, 8, 9, 10]. These processes in turn
are intricately connected to performance-durability trade-offs in both
functional [11] and structural [12] materials. Examples include temperature-
dependent microstructure evolution [13, 14], hotspot formation [15, 16], and
the nucleation and growth of new material phases [17, 18].

However, quantifying local atomic disorder in a physically motivated way in
practice is extraordinarily difficult [19]. Although a number of methods have
been proposed to characterize local atomic environments, these methods are
often not optimized to magnify the subtle differences present in disordered
environments. Existing methods can typically be grouped into three general
classes, each of which carries distinct trade-offs: (1) semi-empirical
structure factors such as Adaptive Common Neighbor Analysis (CNA) [20],
Steinhardt order parameters (SOP) [21], Ackland-Jones order parameters (AJ)
[22], atomic excess volume [23], atomic ”Smoothness” metrics (ASM) [24], and
the local atomic environment metric (LAE) [25]; (2) parameterized symmetry
functions such as the Smooth Overlap of Atomic Positions (SOAP) [26], Behler-
Parinnello functions (BP) [27], Moment Tensor Representations (MTR) [28]; and
the Adaptive Generalizable Neighborhood Informed functions (AGNI) [29]; and
(3) unsupervised machine learning methods which include graph-based [30] and
image-based [31] representations.

In general, it is highly desirable to develop a methodology that is by
construction specifically designed to distinguish, quantify, and physically
interpret regions with varying degrees of atomic disorder. Such a capability
would enable more accurate predictions of how disordered atomic environments
translate to higher-level features and functionality. For instance, mapping
between discretized atomistic models and continuous field representations,
such as phase-field [32, 33, 34] and finite-element [35] models, forces the
use of ill-defined and arbitrary approximations [36], particularly when
disorder is present. Moreover, continuous field representations propagate via
local gradients [37], the evaluation of which amplifies inaccuracies
associated with disordered regions. Addressing these shortcomings is therefore
a critical priority.

To this end we introduce a physics-aware workflow composed of two stages,
which can be seen in Fig. 1. First, we use graph neural networks (GNN) to
explicitly encode local atomic structural information. Next, we apply this
encoding to map the local atomic structure to a novel order parameter that
characterizes local disorder. This order parameter, henceforth referred to as
the Structural Orderness Degree for Atomic Structures (SODAS), $\lambda_{i}$,
quantifies an atom’s local structure in terms of the “closeness” to likely
environments encountered between two limiting cases: a perfect crystal
($\lambda_{i}=$ 1) and a melt ($\lambda_{i}=$ 0). Our approach offers three
distinct advantages: (1) the graph representation accurately encodes the
topology of the connected network of atoms; (2) SODAS is well-defined and
bounded, facilitating universal and intuitive applicability across all atomic
structures of a material; and (3) the compact GNN encoding allows for physical
interpretability. The power of this workflow is demonstrated by application to
several examples of disordered aluminum systems including solid-liquid
interfaces, polycrystalline microstructures, and grain boundaries.

##  2 Results

###  2.1 Definition of SODAS

In principle, the level of configurational disorder can be mapped onto an
equivalent level of thermal disorder in a finite-temperature ensemble. To this
end, we can introduce a fictitious temperature ($T^{\prime}$) that
mathematically represents this configurational disorder. In practice,
$T^{\prime}$ can be parameterized for a given system using explicit MD
simulations as discussed in the Methods section. To physically bound
$T^{\prime}$, we introduce $T_{d}$ as the limit of full disorder (nominally
the melting temperature). The value of $T^{\prime}$ is then confined to the
range between 0 and $T_{d}$. We next define $\gamma$ as a thermodynamic order
parameter:

$\gamma(T^{\prime};T_{d},s)=\mathcal{N}\frac{1}{1+\exp(-(T_{d}/T^{\prime})^{s})}$
(1)

where $\mathcal{N}$ normalizes $\gamma$ between 0 (absolute disorder) and 1
(absolute order), and $s$ is an empirical scaling metric that determines where
to begin the decay of $\gamma$ from ordered to disordered. The introduction of
$s$ makes the definition of $\gamma$ universal, as one can simply tailor its
value for any unique material system. A plot of the relationship between
$\gamma$ and $T^{\prime}$ can be found in the supplemental information.

While $\gamma$ describes the level of disorder of a macroscopic, homogeneously
disordered system, we are primarily interested in local atomic disorder within
a heterogeneous system. To establish this connection, we map the likelihood of
finding a given local atomic environment within an ensemble of configurations,
to a local order parameter (SODAS), $\lambda(n)$, where $n$ indexes an atom.
In practice, due to ergodic constraints, we assume that this ensemble can be
sampled discretely from MD simulations. Note that this is analogous to
determining an atomically resolved configurational entropy density. We
represent this mapping as:

$f(\\{\gamma\\}_{i})\mapsto\lambda(n)$ (2)

where $\\{\gamma\\}_{i}=\\{\gamma_{1},\gamma_{2},...,\gamma_{k}\\}$ represents
the set of $\gamma$ values associated with a given local atomic structural
motif $i$ across a discrete set of $k$ ensembles, and $f$ is a function that
maps $\\{\gamma\\}_{i}$ to $\lambda(n)$. While the function $f$ is unknown, it
can be approximated. In this work we use a graph neural network scheme to
facilitate this approximation, while retaining physical interpretability. Fig.
1 outlines the key steps in this process and is discussed in further detail in
the methods section.

###  2.2 Validation of SODAS

We first validate SODAS by applying our procedure to molecular dynamics
simulations of pristine bulk FCC Al at temperatures up to 1200K (melt). Fig. 2
provides a visual depiction of the $\lambda(n)$ mapping with respect to the
theoretical prediction of $\gamma$. Here we see that all atoms in the
structure at 0K is uniformly predicted to have $\lambda=$ 1, which is
indicative of the perfect crystal. In contrast, at 1200K, all atoms indicate
$\lambda$ to be close to 0 due to the structure existing as a melt. On
average, structures between these limits yield intermediate values of
$\lambda$, as expected.

Figure 2: SODAS calculations on bulk structures taken during a superheating MD
simulation. Values along the y-axis represent the average SODAS value for each
shown structure, whose atoms are colored according to each atom’s SODAS value.
The dashed line indicates the theoretical values of $\gamma$ while the plotted
SODAS values represent the accuracy of the GNN mapping.

At the same time, detailed visualization of intermediate-temperature
configurations reveals a spectrum of atomic environments covering a range of
$\lambda(n)$ in lieu of homogeneously distributed disorder. For example, the
second structure in Fig. 2, which represents a structure at roughly 200K, has
$\lambda$ values ranging from 0.5 to 0.95. Accordingly, as described
previously, similar atomic environments exist at a range of temperatures but
with different degrees of expression according to the average overall level of
disorder. Intuitively this makes sense, as the goal of the SODAS metric is to
judge the likelihood of an atomic environment existing at an arbitrary point
along the thermodynamic spectrum between fully ordered and fully disordered
variants. If a unique atomic environment occurs at multiple temperatures, one
would expect its $\lambda$ to be a weighted combination of the individual
occurrences of the environment along the temperature spectrum.

Importantly, Fig. 2 also compares the GNN-learned SODAS values with respect to
the theoretical values for $\gamma$. Here, we can see excellent agreement
between the GNN mapping and $\gamma$ up to around $T=$1000K. This indicates
that our MD simulations are sufficient to capture the configurational entropy
present within the material at these temperatures. However, deviations exist
above $T=$1000K, with a nearly identical average SODAS value predicted within
this temperature range. This can be explained by our choice of $T_{d}=$ 1200K
(slightly above the EAM melting temperature $\sim$1050K) to improve sampling
in the high-temperature limit; further details can be found in the
supplemental information.

###  2.3 Boundary Identification in Solid-Liquid Interfaces

While crystalline interfaces are easy to identify, disordered interface
boundaries are intrinsically more complex, making classification of the
transition into the boundary region ill-defined. Bond-angle methods such as
CNA and AJ, and more complete methods such as SOAP, often fail to distinguish
between perturbed crystalline and disordered atomic environments. These
difficulties ultimately make defining interface boundaries challenging. In
contrast, the SODAS formalism accomplishes this by providing a continuous
metric that allows for a physically justifiable and mathematically rigorous
definition of the interface boundary transition.

To this end, we have performed two-phase crystal/liquid CMD simulations at
several temperatures (200K, 1000K, and 1500K), to observe how SODAS classifies
the unique structural environments present in each scenario. Further details
regarding the simulation setup can be found in the Methods section. Fig. 3
shows the SODAS characterization of the atomic environments present for the
three temperatures. For the CMD simulation at 200K SODAS correctly identifies
the collapse of the interface boundary. It also locates a small cluster of
local disorder, indicating the presence of a point defect-like region, as a
vestige of solidification of the liquid region. Analogously, at 1500K SODAS
correctly identifies the existence of a single liquid phase and the loss of
the boundary region. SODAS does however, predict varying levels of disorder as
well as the presence of pockets within the liquid that show moderate-to-low
disorder.

Figure 3: Snapshots taken from MD 2-phase simulations. The reference structure
is taken at the moment the 2 phases are allowed to coexist ($t=0$). The 200K,
1000K, and 1500K snapshots are taken at the end of each MD simulation. The
atoms in each system are color-coded based on their SODAS values, as reference
din their colorbar.

On the other hand, at 1000K SODAS correctly captures the interface between the
two phases, as well as moderate perturbations relative to the reference
structure. More importantly, we also see that the interface between the two
phases is not a straight line, instead reflecting the more jagged, realistic
topology of the boundary. SODAS further captures the volumetric nature of the
entire interface boundary region’s complex morphology, which features spatio-
temporally varying 3D shape and thickness.

###  2.4 Autonomous Microstructural Feature Extraction

The ability to define boundary transitions also enables the identification of
larger-scale microstructural features. To demonstrate this capability, we
showcase how a combination of SODAS and a graph characterization algorithm can
autonomously quantify complex polycrystalline microstructures as they evolve
through both time and temperature domains. As described in the Methods
section, we performed three MD simulations of a 1.6 million atom FCC aluminum
system containing 5, 50, and 250 initial grains, respectively.

Figure 4 (a) provides a visual depiction of how one can identify the changes
in microstructure as a function of time, with examples of polycrystalline Al
systems containing 5 and 250 initial grains. In Fig. 4 (b), we provide a
visual workflow that shows how autonomous microstructural feature extraction
is performed using three steps. First, SODAS calculations are performed in
which each atom is assigned a SODAS value as described earlier. Second, an
atomic graph is constructed, where all atoms below a SODAS threshold are
removed and the remaining atoms are converted to an atomic graph [38] in order
to remove grain boundary atoms and leave grain atoms. Third, grain
characterization is performed through subgraph clustering, where an
exhaustive, recursive subgraph search is performed on the atomic graph to
identify the unique grains. Further details can be found in the Methods.
Importantly, this procedure leverages a unique feature of SODAS—namely, the
ability to provide a physically-motivated and mathematically well-defined
threshold for identifying which atoms belong to the grain and which belong to
the grain boundary. This process is extremely challenging using existing
methods due to their lack of a continuously and bounded metric.

Figure 4: (a) SODAS predictions on the initial and final configurations of the
polycrystalline MD simulations. Each column block represents the initial grain
structure present in the system, with the right-side representing the 5
initial grains, and the left side indicative of the 250 initial grains. Each
row, color-coded in green and yellow, represents the temperature of the MD
simulation, with green representing 200K and yellow being 600K. (b) Workflow
of the unsupervised graph-based grain detection algorithm, visualized using
simulated polycrystalline Al. (Left) SODAS values for each atomic environment
present in the system. (Middle) Atoms-to-graph mapping (after SODAS
thresholding), where node colors represent the connectivity of a given atom.
(Right) Autonomous grain detection using recursive subgraph clustering.

Fig. 5 provides a quantitative understanding of how the grain structure
changes as both a function of the initial grain morphology, as well as the
temperature. To quantify this relationship, we invoke our recently developed
graph order parameter (SGOP) [38]. Here, we emphasize that SGOP provides an
approximate measure to a grain’s atomic-level connectivity, which in turn
encodes the ratio between near-boundary and interior atoms, and can therefore
serve as a proxy to that grain’s shape and size. As a result, the distribution
of SGOP values provides a physically-intuitive and robust measure of the
microstructure. This featurization is a more accurate and unique
representation compared to the number of atoms in a grain, the grain radius,
or the grain density.

Figure 5 (a) provides histograms which are generated using the SGOP values
calculated on the final structure from the MD simulations performed at 200K,
400K, and 600K. Here, one can observe different interplay between temperature
and grain size for two examples containing 5 and 250 initial grains. For the
former case (top), at 200K there is an abundance of larger grains, as little-
to-no grain coarsening occurred. At 400K, one can see a broadening of the
distribution, indicating the coexistence of grains both smaller and larger
than those existing at 200K. The same trend continues at 600K, yielding more
medium sized grains. For the latter case of 250 initial grains (bottom), the
distributions are more consistent, indicating a more homogeneous spread of
grain shapes and sizes. Nevertheless, temperature effects continue to broaden
the distributions.

To interpret the results of (a) we track the time-evolution of the
microstructure. To do so, we use a normal distribution to fit the SGOP
distributions, whose mean and standard deviation can be used as a
microstructural feature vector, as shown for the 250-grain case in Figure 5
(b). At 200K there is initially an increase in the mean, shown in Fig. 5
(b,top-left), implying that larger, more connected grains exist throughout the
structure. However, a plateau is reached within 2 ns, and no further increase
in the mean can be seen. At 400K a similar trend is observed, though some
further growth in the grains exists between 10 and 14 ns. At 600K one can
observe a constant trend, indicating that grain growth occurs over the course
of the entire simulation.

Changes in the standard deviation, seen in Fig. 5 (b,top-right), can be
attributed to grain coarsening, because as one grain grows in size, another
neighboring grain must reduce in size as it is swallowed by the growing grain.
This process ultimately causes a fluctuation in the spread of graph order
parameter values, which is captured by the standard deviation of the fitted
distribution. At 200K there is initially an increase in the standard
deviation, indicating the presence of grain coarsening. However, in agreement
with the mean values, a plateau is reached within 2 ns, and no further
increase in the standard deviation can be seen. At 400K, one can observe a
constant increase in the standard deviation, indicating that grain coarsening
is occurring throughout the course of the trajectory. The slope of this line
can be qualitatively linked to the rate of coarsening. At 600K the trend is
similar to 400K but with a larger slope, indicating that grain coarsening at
600K is occurring more rapidly.

Figure 5: Quantitative metrics characterizing microstructural changes over
time. In each case the colors correspond to the temperature of the MD
simulation. (a) Changes in the graph order parameter distribution’s (top,
left) mean, (top, right) standard deviation, (bottom, left) multiplication of
the mean and standard deviation, and (bottom, right) number of grains, all as
a function of time, starting from the 250 grain structure. (b) SGOP
microstructure characterization performed the SODAS predicted values for the
case of (top) 5 initial grains and (bottom) 250 initial grains. Histograms are
determined using the resulting SGOP values, calculated on the final structure
from each MD simulation (the last point in each subplot in (a)). Inserted
images, and corresponding arrow, show the final structure obtained from the
600K MD simulations. All histograms are normalized and show each x-value’s
probability.

Figure 5 (b,bottom-left) combines the mean and standard deviation into a
single metric. Here, the product of the two terms is used to gauge both the
rate of grain growth and the magnitude of the grains themselves. At 200K we
see a trend that confirms our previous observation that grain growth initially
occurs rapidly, but tapers off quickly. At 400K we observe a more complete
picture, where slow grain growth occurs throughout the trajectory. For the
case of 600K, faster growth can be seen, as well as the existence of larger
grains.

This picture can also be complemented by observing how the number of grains
changes as a function of time, as seen in Fig. 5 (b,bottom-right). We can see
that at 200K, there is an initial reduction in the number of grains from 250
to approximately 75, where that number holds for the remainder of the
simulation. At 400K there is a constant reduction in the number of grains from
250 roughly 30 by the end of the MD simulation. At 600K the number of grains
is reduced even further by the end of the simulation, going from 250 to only
12 grains. The picture painted in Fig. 5 (b,bottom-right) aligns well with the
values in Fig. 5 (b,bottom-left), indicating that our product metric can be
used to quantify both the size, shape, and the growth rate of the grains.

###  2.5 Mapping Atoms to a Continuous Field

Lastly, since $\lambda$ is continuously valued over the discrete atoms, it can
be interpolated to a continuous field. This mapping allows for its integration
into continuum models. For instance, we note the similarity between such a
continuous field representation and the phase order parameter used in phase
field models [8, 39]. We showcase this concept for the example of two grain
boundary regions with varying levels of interfacial complexity. Nevertheless,
we note that this method can be used for other classes of crystalline
interfaces, such as symmetric tilt and twin boundaries, and edge/screw
dislocations.

Figure 6: Continuous fields (and its gradient norm) of the originally
discrete, per-particle SODAS value $\lambda$. The discrete-to-continuum
conversion is done by interpolating the discrete $\lambda$ onto a uniform,
fine grid. The gradient information can then be computed over the uniform
grid.

From Fig. 6, one can see the intuitive nature of SODAS, cleanly characterizing
the grain regions with a $\lambda$ close to 1, smoothly transitioning to
higher degrees of disorder present near the boundary. For boundaries that show
higher degrees of crystallinity, such as those in $\sum 5(110)[120]$, the
disorder present at the interface is minimal, as is expected, though is still
clearly present. Likewise, for more disordered boundaries, such as those in
$\sum 9(110)[110]$, a greater degree of disorder is detected within the
interface region. As in the previous section, these characterizations
exemplify the ability of SODAS to intuitively determine where the grain begins
and ends.

Figure 6 shows the continuous fields derived from the originally discrete,
per-particle SODAS value $\lambda$. Additionally, the gradient norm
$||\nabla\lambda||$ was calculated and visualized. This discrete-to-continuum
conversion was done by interpolating the discrete $\lambda$ values onto a
uniform grid using PyVista [40]. When calculating the gradient of this field
we observe areas of the structure where there are sharp changes in the SODAS
values. Notably, the gradient is maximized not at the center of the grain
boundary, but rather at the transition to the boundary region, because these
are locations within the structure where there is an abrupt change in the
level of disorder present.

The sensitivity of this detection can be seen in Fig. 6, where the gradient of
the scalar field predicts two regions where there is an abrupt change in the
SODAS values. As we move from the crystalline regions towards the interface
normal to the boundary region, we first encounter a crystal-to-boundary
region, followed by the boundary itself, and finally a boundary-to-crystal
region as we move away from the interface. Therefore, the gradient predictions
in Fig. 6 highlight the fact that a degree of homogeneity can exist in both
the ordered interior of the grain as well as the disordered interior of the
grain boundary.

##  3 Discussion and Conclusion

In summary, characterizing the nature of local atomic disorder is critical and
necessary to understand how structure-property relationships evolve. SODAS is
a new mathematical framework in which local atomic environments are
transformed into graph representations, encoded via a graph neural network
paradigm, and finally mapped onto a local order parameter. This order
parameter, $\lambda$, is a physically intuitive, continuous, and
mathematically bounded scalar which represents the level of disorder present
within an atomic environment, and is analogous to an atomically resolved
configurational entropy density. In addition to the examples shown throughout
this work, these advantages allow for the universal quantification of a
multitude of complex and heterogeneous materials properties and phenomena.

We also envision our proposed methodology as a novel tool for multiscale model
integration. In particular, SODAS provides an atomistically derived,
physically motivated continuous scalar field representation for phase field
and continuum models. This mapping can be likewise leveraged to output field
quantities such as phase order, grain distribution, concentration,
stress/strain, and so on. Such an approach offers a new perspective and
valuable technique for bridging scales in multiscale models, both between
atomistic and microscale descriptions, as well as between discrete and
continuous representations. We further emphasize that although this work
focuses on single-element systems, our method is generally applicable to
multi-component systems and their corresponding microstructural features.

The advantages of SODAS also become clear for extraction of physical
properties that relate to materials performance or degradation. For instance,
we showed that by interpolating the discrete representation to a continuum
representation, we can analyze or differentiate $\lambda$ to deduce spatially
resolved changes in structural homogeneity. In practice, these structural
changes often map to changes in key response properties, including
diffusivity, dielectric response, electrical conductivity, and elastic
compliance. [41] In cases where such properties can be computed locally or
measured using local probes, SODAS offers a way to extract analytical
relationships between structure and function. Moreover, sharp gradients from
abrupt changes in response functions can concentrate electrical, chemical, or
mechanical potential, creating “hotspots” that can initiate key
electrochemomechanical failure modes. We therefore propose that gradients in
the continuous representation of $\lambda$ may provide a unique way to
identify such “hotspots”, with direct connection to early prediction of
propensity for deleterious outcomes such as fracture, corrosion, and thermal
runaway.

##  4 Methods

###  4.1 Training Data Preparation

Classical molecular dynamics (CMD), using the LAMMPS software package [42],
was used to generate training data for the GNN model. Starting from bulk FCC
aluminum (containing 1024 atoms), CMD was performed in the NVT ensemble for 1
nanosecond using the Zhou et al EAM potential [43], increasing the temperature
linearly with respect to time. The initial and final tempeartures for the CMD
simulation was 100 and 1200, respectively. 1000 equidistant snapshots were
then taken along the trajectory and labelled according to their respective
temperatures. Further information regarding the training data preparation can
be found in the supplemental information.

During training, all atoms at a given temperature are assigned the same value
of $\gamma$. In this way, the value of $\gamma$ is not directly connected to
an atom’s local structural geometry. However, at each temperature local atomic
geometries bounce around an equilibrium point, which in this case can be
thought of as the 0K structure. As the temperature increases the magnitude of
the displacement from this structure increases, though some atomic
neighborhoods may resemble low temperature structural motifs, even at higher
overall temperatures.

###  4.2 Graph Neural Network Implementation

####  4.2.1 Conversion to Graph

Prior to GNN operation, we converted the atomic systems into graphs by a
simple cutoff radius-based neighbor list search (implemented using Atomic
Simulation Environment [44]), with the cutoff $R_{c}=3.5$ Å. Each node of the
converted graph corresponds to the atom type $z$, and each edge the bond
distance $d$.

####  4.2.2 GNN Operation

The GNN model used in this work consists of three components: the initial
embedding, the atom-bond interactions, and the final output layers (Fig. 1).

In the initial embedding, each atom type $z$ is transformed into a feature
vector by an `Embedding` layer (PyTorch [45]). Each bond distance $d$ is
expanded into a $D$-dimensional feature vector by the Radial Bessel basis
functions (RBF) [46]

$\text{RBF}_{n}(d)=\sqrt{\frac{2}{R_{c}}}\frac{\sin(\frac{n\pi}{R_{c}}d)}{d},$
(3)

where $n\in[1..D]$ and $R_{c}$ is the cutoff value. Both atom and bond feature
vectors have the same length $D=100$.

The atom-bond interactions are also known as GNN convolution, aggregation, or
message-passing. There are many variants of GNN convolution operations that
can be adopted from the literature. In this work, we choose the edge-gated
graph convolution [47, 48]. The term _atom-bond interaction_ is based on the
fact that the nodes and the edges exchange information during the convolution
operation. Specifically, the node features $\mathbf{h}^{l+1}_{i}$ of node $i$
at the $(l+1)$th layer is updated as

$\mathbf{h}^{l+1}_{i}=\mathbf{h}^{l}_{i}+\mathrm{SiLU}\left(\mathrm{LayerNorm}\left(\mathbf{W}^{l}_{s}\mathbf{h}^{l}_{i}+\sum_{j\in\mathcal{N}(i)}\hat{\mathbf{e}}^{l}_{ij}\odot\mathbf{W}^{l}_{d}\mathbf{h}^{l}_{j}\right)\right),$
(4)

where SiLU is the Sigmoid Linear Unit activation function [49]; LayerNorm is
the Layer Normalization operation [50]; $\mathbf{W}_{s}$ and $\mathbf{W}_{d}$
are weight matrices; the index $j$ denotes the neighbor node of node $i$;
$\hat{\mathbf{e}}_{ij}$ is the edge gate vector for the edge from node $i$ to
node $j$; and $\odot$ denotes element-wise multiplication. The edge gate
$\hat{\mathbf{e}}^{l}_{ij}$ at the $l$th layer is defined as

$\hat{\mathbf{e}}^{l}_{ij}=\frac{\sigma(\mathbf{e}^{l}_{ij})}{\sum_{j^{\prime}\in\mathcal{N}(i)}\sigma(\mathbf{e}^{l}_{ij^{\prime}})+\epsilon},$
(5)

where $\sigma$ is the sigmoid function, $\mathbf{e}^{l}_{ij}$ is the original
edge feature, and $\epsilon$ is a small constant for numerical stability. The
edge features $\mathbf{e}^{l}_{ij}$ is updated by

$\mathbf{e}^{l+1}_{ij}=\mathbf{e}^{l}_{ij}+\mathrm{SiLU}\left(\mathrm{LayerNorm}\left(\mathbf{W}^{l}_{g}\mathbf{z}^{l}_{ij}\right)\right),$
(6)

where $\mathbf{W}_{g}$ is a weight matrix, and $\mathbf{z}_{ij}$ is the
concatenated vector from the node features $\mathbf{h}_{i}$, $\mathbf{h}_{j}$,
and the edge features $\mathbf{e}_{ij}$:

$\mathbf{z}_{ij}=\mathbf{h}_{i}\oplus\mathbf{h}_{j}\oplus\mathbf{e}_{ij}.$ (7)

Lastly, in the final output layers, each node feature is eventually
transformed into a scalar output $y$ ranging from 0 to 1. Effectively, the GNN
predicts the SODAS metric for every atom.

Further details regarding model training are described in Supporting
Information.

###  4.3 Atomistic Simulation Details

####  4.3.1 2 Phase Simulations

2 phase CMD simulations were performed on a orthogonal block of aluminum
containing 23,040 atoms. All CMD simulations were performed in the NVT
ensemble. Initially, half the box was superheated to 4000K to assure complete
melting, while the remaining half of the box was set at 200K. This initial CMD
simulation was performed for 1 nanosecond with separate NVT thermostats
driving each region. The 4000K and 200K regions were then allowed to coexist
using a single NVT thermostat to drive the CMD simulation. For the purpose of
this analysis upon the introduction of the shared thermostat the temperature
was fixed to 3 cases: one at 200K, another at 1000K, and a final case of
1500K. The purpose of these temperatures is to observe how our proposed
methodology predicts the unique levels of atomic disorder present in each
case. Further information regarding these simulations can be found in the
Results section and in Fig. 3.

####  4.3.2 Grain Coarsening Simulations

CMD simulations in the NVT ensemble were performed for 6 polycrystalline
cases, each with a varying number of initial grains. An initial bulk aluminum
system containing roughly 1.6 million atoms was used to construct 6
polycrystalline structures, using the Atomsk software package [51], containing
5 and 250 initial grains. CMD simulations were performed on each case at 200K,
400K, and 600K. NVT simulations were run for approximately 1.5 nanoseconds for
each combination of initial structure and temperature. All CMD simulations
were done using LAMMPS with the Zhou et. al EAM potential [43]. Further
details regarding the polycrystalline structures can be found in the Results
section as well as the Supplemental Information.

###  4.4 Microstructure Characterization

Microstructure characterization occurs in 4 stages: (1) calculation of SODAS
for all atoms in the system, (2) thresholding of the atomic configuration,
based on an atom’s SODAS value and subsequent removal of all atoms below the
threshold value, (3) conversion of the remaining atoms to a graph
representation for the discovery of subgraphs within the graph, and (4)
characterization of the grains through a graph order parameter. Fig. 4 (b)
depicts this workflow visually. While step (1) requires little-to-no input
from the user, step (2) requires one to define the level of disorder that
needs to be captured when defining the interface regions. This choice
highlights the intuitive nature of our proposed methodology, as the threshold
value defines the structural properties of the interface region itself, with a
near-zero threshold indicating grain boundaries which are extremely disordered
and a value close to 1 representing highly crystalline boundaries. In
principle, both classes of interfaces can exist within the same structure,
which would require a more complex thresholding system, though for this work
we assume a uniform local atomic environment amongst all grain boundaries.

For all microstructure characterization tests in this work we employ a
$\lambda$ threshold of 0.7, implying that we are defining grain boundaries as
local atomic environments showing disorder equivalent to the atomic
perturbation observed at around 900K. $\lambda=$ 0.7 was chosen due to the
maximum temperatures described in the previous section. As we should expect
the grain regions to experience perturbations no greater than those
encountered at 600K, those observed at 900K should provide a good
approximation for atoms that do not belong to the grains. Once the system has
been thresholded, all atoms below the threshold value are removed, leaving
only the atoms belonging to the grains.

This system is then mapped onto a graph, $G$, shown in Fig. 4 (b), where edges
are represented by ij pairwise interactions within a 4Åcutoff radius. A
recursive subgraph search algorithm is employed to discover all connected
subgraphs, $S_{G}$, within the complete graph $G$. This algorithm is extremely
efficient, discovering all subgraphs within a 1.6 million atom system in 1.2
seconds. As all interface atoms were removed prior to the graph construction,
all subgraphs in $G$ represent the resulting grains contained within the
structure. A 2D slice of the discovered grains in the 3D system is seen in
Fig. 4 (b). As each grain is represented by a subgraph in $G$, they can be
characterized as graphs. This provides more detail into the grain’s shape and
connectivity, rather than relying on metrics such as the number of atoms,
diameter, etc, which are inherently not unique and do not provide a true
quantitative measure of the underlying properties of the structure. Here we
use a recently developed graph topology metric [38] to characterize the
subgraphs in $G$. Each edge in the subgraph is labelled according to:

$\omega_{i,j}=\frac{1}{d_{i,j}}\ni d_{i,j}\leq R_{c}$ (8)

$i$ and $j$ are the atomic indices of the atoms (note that self-interaction
terms are not allowed). $d_{i,j}$ is defined as the $l^{2}$-norm between two
atoms. $R_{c}$ is the cutoff radius, which was chosen as 4Å. $\omega_{i,j}$
represents the weight of a given edge for a specific pair of adjacent nodes in
the graph. The degree of each node is finally defined as the sum of the
elements in a node’s edge set, $d_{i}=\sum_{j}^{J}\omega_{i,j}$, where $J$ is
the set of all neighbors of atom $i$. The degree sets are then fed into the
scalar graph order parameter (SGOP) [38] scheme for the final characterization
of the subgraph. The SGOP functional form is defined as:

$\theta_{S_{G},R_{c}}=\sum_{m}^{D_{s}}P(d_{m})\log_{b}P(d_{m})+d_{m}P(d_{m})$
(9)

where $D_{s}$ is the set of unique node degrees in a subgraph, with
$P_{d_{m}}$ being the probability of a given degree, $d_{m}$, occurring in the
subgraph. As discussed in previous works, the SGOP value provides a measure of
the shape and connectivity of the graph at a structural level. In this way we
can observe, quantitatively, the difference between grains with a similar
number of atoms, diameter, circumference, density, etc.

## Supporting Information

This work contains supplemental information which can be found online.

## Data Availability

All data required to reproduce this work can be requested by contacting the
corresponding author.

## Acknowledgements

J. Chapman, T. Hsu, X. Chen, T. W. Heo, and B. C. Wood are partially supported
by the Laboratory Directed Research and Development (LDRD) program (20-SI-004)
at Lawrence Livermore National Laboratory. This work was performed under the
auspices of the US Department of Energy by Lawrence Livermore National
Laboratory under contract No. DE-AC52-07NA27344.

## Author Contributions

X. Chen and B. C. Wood supervised the research. J. Chapman performed all MD
simulations, and devised/implemented the autonomous microstructure feature
extraction methodology. T. Hsu trained the GNN and performed all GNN-related
predictions. J. Chapman and T. Hsu devised the theoretical SODAS framework. B.
Wood, T. Hsu, and J. Chapman devised the atoms-to-field mapping, while T. Hsu
implemented it. T. W. Heo provided insight into the connection between
atomistic and phase field modelling, and helped guide discussions surrounding
the atoms to continuous field methodology. J. Chapman and T. Hsu wrote the
manuscript with inputs from all authors.

## Competing Interests

The authors declare no competing financial or non-financial interests.

## References

  * [1] Jacob Fish, Gregory J. Wagner and Sinan Keten  “Mesoscopic and multiscale modelling in materials”  In _Nature Materials_ 20.6, 2021, pp. 774–786
  * [2] Ronald E Miller and E B Tadmor  “A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods”  In _Modelling and Simulation in Materials Science and Engineering_ 17.5 IOP Publishing, 2009, pp. 053001  DOI: 10.1088/0965-0393/17/5/053001
  * [3] Alain Karma and Damien Tourret  “Atomistic to continuum modeling of solidification microstructures” Recent Advances in Solidification Microstructure- Experiments and Computational Analysis  In _Current Opinion in Solid State and Materials Science_ 20.1, 2016, pp. 25–36  DOI: https://doi.org/10.1016/j.cossms.2015.09.001
  * [4] S.F. Guo et al.  “Fe-based amorphous coating for corrosion protection of magnesium alloy”  In _Materials and Design_ 108, 2016, pp. 624–631
  * [5] Y. Lin et al.  “Nanoparticle Assembly and Transport at Liquid-Liquid Interfaces”  In _Science_ 299.5604, 2003, pp. 226–229
  * [6] Brent A. Wacaser et al.  “Preferential Interface Nucleation: An Expansion of the VLS Growth Mechanism for Nanowires”  In _Advanced Materials_ 21.2, 2009, pp. 153–165
  * [7] Itamar Willner and Eugenii Katz  “Controlling Chemical Reactivity at Solid-Solution Interfaces by Means of Hydrophobic Magnetic Nanoparticles”  In _Langmuir_ 22.4, 2006, pp. 1409–1419
  * [8] Tae Wook Heo et al.  “Microstructural impacts on ionic conductivity of oxide solid electrolytes from a combined atomistic-mesoscale approach”  In _npj Computational Materials_ 7.1, 2021, pp. 214  DOI: 10.1038/s41524-021-00681-8
  * [9] Tae Wook Heo and Long-Qing Chen  “Phase-field modeling of displacive phase transformations in elastically anisotropic and inhomogeneous polycrystals”  In _Acta Materialia_ 76, 2014, pp. 68–81  DOI: https://doi.org/10.1016/j.actamat.2014.05.014
  * [10] Tae Wook Heo, Kimberly B. Colas, Arthur T. Motta and Long-Qing Chen  “A phase-field model for hydride formation in polycrystalline metals: Application to $\gamma$-hydride in zirconium alloys”  In _Acta Materialia_ 181, 2019, pp. 262–277  DOI: https://doi.org/10.1016/j.actamat.2019.09.047
  * [11] Craig S. Gittleman, Anusorn Kongkanand, David Masten and Wenbin Gu  “Materials research and development focus areas for low cost automotive proton-exchange membrane fuel cells” Energy Storage Energy Transformation  In _Current Opinion in Electrochemistry_ 18, 2019, pp. 81–89
  * [12] Reza Abbasi et al.  “A Roadmap to Low-Cost Hydrogen with Hydroxide Exchange Membrane Electrolyzers”  In _Advanced Materials_ 31.31, 2019, pp. 1805876
  * [13] S.. Srinivasan, M.. Baskes and G.. Wagner  “Atomistic simulations of shock induced microstructural evolution and spallation in single crystal nickel”  In _Journal of Applied Physics_ 101.4, 2007, pp. 043504  DOI: 10.1063/1.2423084
  * [14] Shan Liu et al.  “Atomistic simulation of microstructure evolution of NiTi single crystals in bending deformation”  In _Computational Materials Science_ 199, 2021, pp. 110733  DOI: https://doi.org/10.1016/j.commatsci.2021.110733
  * [15] Michael Simon and Edson L. Meyer  “Detection and analysis of hot-spot formation in solar cells”  In _Solar Energy Materials and Solar Cells_ 94.2, 2010, pp. 106–113  DOI: https://doi.org/10.1016/j.solmat.2009.09.016
  * [16] A. Tokmakoff, M.. Fayer and Dana D. Dlott  “Chemical reaction initiation and hot-spot formation in shocked energetic molecular materials”  In _The Journal of Physical Chemistry_ 97.9, 1993, pp. 1901–1913  DOI: 10.1021/j100111a031
  * [17] E. Budevski, G. Staikov and W.J. Lorenz  “Electrocrystallization: Nucleation and growth phenomena”  In _Electrochimica Acta_ 45.15, 2000, pp. 2559–2574  DOI: https://doi.org/10.1016/S0013-4686(00)00353-4
  * [18] Renyi Zhang et al.  “Nucleation and Growth of Nanoparticles in the Atmosphere” PMID: 22044487  In _Chemical Reviews_ 112.3, 2012, pp. 1957–2011  DOI: 10.1021/cr2001756
  * [19] Noam Bernstein et al.  “Quantifying Chemical Structure and Machine-Learned Atomic Energies in Amorphous and Liquid Silicon”  In _Angewandte Chemie International Edition_ 58.21, 2019, pp. 7057–7061
  * [20] Alexander Stukowski  “Structure identification methods for atomistic simulations of crystalline materials”  In _Model Simul. Mat. Sci. Eng._ 20.4 IOP Publishing, 2012, pp. 045021  DOI: 10.1088/0965-0393/20/4/045021
  * [21] Paul J. Steinhardt and P. Chaudhari  “Point and line defects in glasses”  In _Philos. Mag. A_ 44.6 Taylor & Francis, 1981, pp. 1375–1381  DOI: 10.1080/01418618108235816
  * [22] G.. Ackland and A.. Jones  “Applications of local crystal structure measures in experiment and simulation”  In _Phys. Rev. B_ 73 American Physical Society, 2006, pp. 054104  DOI: 10.1103/PhysRevB.73.054104
  * [23] Yasir Mahmood et al.  “Atomistic and machine learning studies of solute segregation in metastable grain boundaries”  In _Scientific Reports_ 12.1, 2022, pp. 6673  DOI: 10.1038/s41598-022-10566-5
  * [24] Tristan A. Sharp et al.  “Machine learning determination of atomic dynamics at grain boundaries”  In _Proceedings of the National Academy of Sciences_ 115.43, 2018, pp. 10943–10947  DOI: 10.1073/pnas.1807176115
  * [25] Conrad W. Rosenbrock, Eric R. Homer, Gábor Csányi and Gus L.. Hart  “Discovering the building blocks of atomic systems using machine learning: application to grain boundaries”  In _npj Computational Materials_ 3.1, 2017, pp. 29
  * [26] Sandip De, Albert P. Bartók, Gábor Csányi and Michele Ceriotti  “Comparing molecules and solids across structural and alchemical space”  In _Phys. Chem. Chem. Phys._ 18 The Royal Society of Chemistry, 2016, pp. 13754–13769  DOI: 10.1039/C6CP00415F
  * [27] Jörg Behler  “Atom-centered symmetry functions for constructing high-dimensional neural network potentials”  In _J. Chem. Phys._ 134.7, 2011, pp. 074106  DOI: 10.1063/1.3553717
  * [28] Alexander V. Shapeev  “Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials”  In _Multiscale Modeling & Simulation_ 14.3, 2016, pp. 1153–1173  DOI: 10.1137/15M1054183
  * [29] J. Chapman, R. Batra and R. Ramprasad  “Machine learning models for the prediction of energy, forces, and stresses for Platinum”  In _Computational Materials Science_ 174, 2020, pp. 109483
  * [30] Cheol Woo Park et al.  “Accurate and scalable graph neural network force field and molecular dynamics with direct force architecture”  In _npj Computational Materials_ 7.1, 2021, pp. 73
  * [31] Henry Chan et al.  “Machine learning enabled autonomous microstructural characterization in 3D samples”  In _npj Computational Materials_ 6.1, 2020, pp. 1
  * [32] Ingo Steinbach  “Phase-field models in materials science”  In _Modelling and Simulation in Materials Science and Engineering_ 17.7 IOP Publishing, 2009, pp. 073001  DOI: 10.1088/0965-0393/17/7/073001
  * [33] Long-Qing Chen  “Phase-Field Models for Microstructure Evolution”  In _Annual Review of Materials Research_ 32.1, 2002, pp. 113–140  DOI: 10.1146/annurev.matsci.32.112001.132041
  * [34] W.. Boettinger, J.. Warren, C. Beckermann and A. Karma  “Phase-Field Simulation of Solidification”  In _Annual Review of Materials Research_ 32.1, 2002, pp. 163–194  DOI: 10.1146/annurev.matsci.32.101901.155803
  * [35] Brian G. Richmond et al.  “Finite element analysis in functional morphology”  In _The Anatomical Record Part A: Discoveries in Molecular, Cellular, and Evolutionary Biology_ 283A.2, 2005, pp. 259–274  DOI: https://doi.org/10.1002/ar.a.20169
  * [36] “A temperature-dependent atomistic-informed phase-field model to study dendritic growth”  In _Journal of Crystal Growth_ 579, 2022, pp. 126461  DOI: https://doi.org/10.1016/j.jcrysgro.2021.126461
  * [37] Heike Emmerich  “Advances of and by phase-field modelling in condensed-matter physics”  In _Advances in Physics_ 57.1 Taylor & Francis, 2008, pp. 1–87  DOI: 10.1080/00018730701822522
  * [38] James Chapman, Nir Goldman and Brandon C. Wood  “Efficient and universal characterization of atomic structures through a topological graph order parameter”  In _npj Computational Materials_ 8.1, 2022, pp. 37  DOI: 10.1038/s41524-022-00717-7
  * [39] Tae Wook Heo et al.  “A mesoscopic digital twin that bridges length and time scales for control of additively manufactured metal microstructures”  In _Journal of Physics: Materials_ 4.3 IOP Publishing, 2021, pp. 034012  DOI: 10.1088/2515-7639/abeef8
  * [40] C. Sullivan and Alexander Kaszynski  “PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK)”  In _Journal of Open Source Software_ 4.37 The Open Journal, 2019, pp. 1450  DOI: 10.21105/joss.01450
  * [41] Tae Wook Heo et al.  “Microstructural impacts on ionic conductivity of oxide solid electrolytes from a combined atomistic-mesoscale approach”  In _npj Computational Materials_ 7.1 Nature Publishing Group, 2021, pp. 1–15
  * [42] Steve Plimpton  “Fast Parallel Algorithms for Short-Range Molecular Dynamics”  In _Journal of Computational Physics_ 117.1, 1995, pp. 1–19  DOI: https://doi.org/10.1006/jcph.1995.1039
  * [43] X.W. Zhou et al.  “Atomic scale structure of sputtered metal multilayers”  In _Acta Materialia_ 49.19, 2001, pp. 4005–4015  DOI: https://doi.org/10.1016/S1359-6454(01)00287-7
  * [44] Ask Hjorth Larsen et al.  “The atomic simulation environment—a Python library for working with atoms”  In _Journal of Physics: Condensed Matter_ 29.27, 2017, pp. 273002  URL: http://stacks.iop.org/0953-8984/29/i=27/a=273002
  * [45] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035
  * [46] Johannes Klicpera, Janek Gross and Stephan Gunnemann  “Directional message passing for molecular graphs”  In _arXiv preprint arXiv:2003.03123_ , 2020
  * [47] Xavier Bresson and Thomas Laurent  “Residual gated graph convnets”  In _arXiv preprint arXiv:1711.07553_ , 2017
  * [48] Vijay Prakash Dwivedi et al.  “Benchmarking graph neural networks”  In _arXiv preprint arXiv:2003.00982_ , 2020
  * [49] Stefan Elfwing, Eiji Uchibe and Kenji Doya  “Sigmoid-weighted linear units for neural network function approximation in reinforcement learning” Special issue on deep reinforcement learning  In _Neural Networks_ 107, 2018, pp. 3–11  DOI: https://doi.org/10.1016/j.neunet.2017.12.012
  * [50] Jingjing Xu et al.  “Understanding and Improving Layer Normalization”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper/2019/file/2f4fe03d77724a7217006e5d16728874-Paper.pdf
  * [51] Pierre Hirel  “Atomsk: A tool for manipulating and converting atomic data files”  In _Computer Physics Communications_ 197, 2015, pp. 212–219  DOI: https://doi.org/10.1016/j.cpc.2015.07.012