# Exact time-dependent density functional theory for non-perturbative dynamics
of helium atom

Davood Dar  Department of Physics, Rutgers University, Newark 07102, New
Jersey USA  Lionel Lacombe  Department of Physics, Rutgers University, Newark
07102, New Jersey USA  Johannes Feist  Departamento de Física Teórica de la
Materia Condensada and Condensed Matter Physics Center (IFIMAC), Universidad
Autónoma de Madrid, Madrid, Spain  Neepa T. Maitra  Department of Physics,
Rutgers University, Newark 07102, New Jersey USA<EMAIL_ADDRESS>
###### Abstract

By inverting the time-dependent Kohn-Sham equation for a numerically exact
dynamics of the helium atom, we show that the dynamical step and peak features
of the exact correlation potential found previously in one-dimensional models
persist for real three-dimensional systems. We demonstrate that the Kohn-Sham
and true current-densities differ by a rotational component. The results have
direct implications for approximate TDDFT calculations of atoms and molecules
in strong fields, emphasizing the need to go beyond the adiabatic
approximation, and highlighting caution in quantitative use of the Kohn-Sham
current.

##  I Introduction

For simulating dynamics of electrons in non-perturbative fields, time-
dependent density functional theory Runge and Gross (1984); Maitra (2016);
Marques  _et al._ (2012); Ullrich (2011) (TDDFT) has emerged as a key
approach, due to its favorable system-size scaling. In theory, TDDFT is an
exact formulation of quantum mechanics which provides a computationally
tractable approach for tackling calculations involving many-body problems in
external time-dependent fields. Mapping to a non-interacting system, the Kohn-
Sham (KS) system, that exactly reproduces the one-body density allows the
computation of much larger systems than with traditional wavefunction methods,
with no restriction on the strength of the applied fields nor on how far the
system is driven from equilibrium, see Refs. Draeger  _et al._ (2017); Sato
(2021); Guandalini  _et al._ (2021); Singh  _et al._ (2021); Folorunso  _et
al._ (2021); Yost  _et al._ (2017) for examples in a range of recent
applications.

TDDFT does not, however, come without its own difficulties; in particular, the
exchange-correlation (xc) potential in which the many-body effects are
“hidden” needs to be approximated as a functional of density, and the exact xc
potential depends on the density in a spatially- and time- non-local way. This
dependence is neglected in adiabatic approximations used in calculations
today, where the instantaneous density is inserted into a ground-state xc
approximation. A crucial question is how well these approximations accurately
capture the true dynamics. The lack of memory-dependence is believed to be
responsible for errors in its predictions e.g. Refs. Raghunathan and Nest
(2011); Ramakrishnan and Nest (2012); Habenicht  _et al._ (2014); Wijewardane
and Ullrich (2008); Gao  _et al._ (2017); Boström  _et al._ (2018); Krumland
_et al._ (2020); Quashie  _et al._ (2017); Gao  _et al._ (2014); Henkel  _et
al._ (2009); Da  _et al._ (2017); Yao  _et al._ (2021), including sometimes
qualitative failures. Still, the approximations are often accurate enough to
be useful, and some characterization of when to expect the adiabatic
approximation to work well has also been done Lacombe and Maitra (2020).
Studies of the exact xc potential have been made, to compare against
approximate potentials, and also to study the impact of its features on the
resulting dynamics. Such studies require two ingredients: first, an exact
calculation of dynamics of interacting electrons, and second, an inversion of
the TDKS equations to find the exact potential. Because of challenges in
obtaining these ingredients, the studies have so far been limited to model
systems Lacombe and Maitra (2020); Luo  _et al._ (2014); Suzuki  _et al._
(2017); Elliott  _et al._ (2012); Fuks  _et al._ (2013); Lacombe and Maitra
(2019); Lacombe  _et al._ (2018); Ramsden and Godby (2012); Fuks  _et al._
(2016, 2018); Hodgson  _et al._ (2016, 2013); D’Agosta and Vignale (2005);
Ruggenthaler  _et al._ (2013); Elliott and Maitra (2012); Thiele  _et al._
(2008); Ullrich (2006) involving either one-dimension (1D) and/or two
electrons, or involving only small perturbations away from the ground-state.

In this work, the exact time-dependent KS (TDKS) potential is found for the
first time, for a real three-dimensional (3D) multi-electron atom in the non-
perturbative regime. We find dynamical step and peak features that have a non-
local-in-space and -in-time dependence on the density. The results have direct
implications for TDDFT calculations of atoms and molecules driven far from
their ground-state, as these features are missing in adiabatic approximations.
They justify the relevance of the previous 1D studies, where similar dynamical
step and peak features are found in the correlation potential. Moreover, the
example explicitly demonstrates that the KS current-density differs from the
true current-density by a rotational component. Although this has been
recognized before to be theoretically possible D’Agosta and Vignale (2005);
Maitra  _et al._ (2002a); Schaffhauser and Kümmel (2016); Thiele and Kümmel
(2009); Gross and Maitra (2012), not only is this difference neglected in
applications today where typically the current calculated from the KS orbitals
is assumed to represent the true current Andrade  _et al._ (2018); Takeuchi
_et al._ (2019), but the difference has not been demonstrated for systems
beyond the linear response regime.

##  II Dynamics in the Helium Atom

The system we study is the field-free evolution of a superposition state of
the Helium atom, as might be reached, for example, by applying a field that is
turned off after some time. The lowest few eigenstates of this atom were found
using the time-dependent close-coupling method, making a partial wave
expansion in coupled spherical harmonics, and using the finite element
discrete variable representation to discretize the radial degrees of freedom
Feist (2009); Feist  _et al._ (2008). We consider here linear superpositions
of the singlet ground state $1^{1}S_{0}$, denoted $\Psi_{0}$, and singlet
first excited state $2^{1}P_{1}$ that has angular quantum numbers $L=1$ and
$M=0$, denoted here $\Psi_{1}$, so the exact time-evolution of the two-
electron state is

$|\Psi(t)\rangle=\frac{1}{\sqrt{1+|a|^{2}}}\left(|\Psi_{0}\rangle+ae^{-i\omega
t}|\Psi_{1}\rangle\right)$ (1)

where $\omega=E_{2^{1}P}-E_{1^{1}S}=0.77980$ in atomic units (a.u.), is the
frequency with which the system oscillates.The parameter $a$ gives the
relative fraction of the excited state,for example $a=1$ in the case of a
50:50 superposition. We aim then to find the time-dependent KS potential which
reproduces the exact density of the interacting state Eq. (1):

$n({\bf r},t)=\frac{1}{1+|a|^{2}}\left(n_{0}({\bf r})+|a|^{2}n_{1}({\bf
r})+2an_{01}({\bf r})\cos(\omega t)\right)$ (2)

where $n_{q}({\bf r})=2\int|\Psi_{q}({\bf r},{\bf r}_{2})|^{2}d^{3}{\bf
r}_{2},q=0,1$ and $n_{01}({\bf r})=2\int\Psi_{0}({\bf r},{\bf
r}_{2})\Psi_{1}({\bf r},{\bf r}_{2})d^{3}{\bf r}_{2}$ .

We note here that the results we find for the xc potential apply to far more
general dynamical situations than the field-free superposition state dynamics:
due to an exact condition Maitra  _et al._ (2002b), the xc potential applies
to any situation where the instantaneous interacting state is given by Eq. (1)
at some time $t$, and the KS state is a Slater determinant (see Appendix A).

Now, the TDDFT xc potential depends on the choice of the initial KS state van
Leeuwen (1999); Runge and Gross (1984); Elliott and Maitra (2012); the 1-1
density-potential mapping holds only for a given initial state, which endows
$v_{\scriptscriptstyle\rm XC}({\bf r},t)$ with a functional dependence on both
the true and KS states, $v_{\scriptscriptstyle\rm XC}[n;\Psi(0),\Phi(0)]({\bf
r},t)$. In principle, one can begin in any initial KS state that reproduces
the density of the initial interacting state and its first time-derivative;
the structure of the exact xc potential has a strong dependence on this choice
Elliott and Maitra (2012); Fuks  _et al._ (2016, 2018); Lacombe and Maitra
(2019); Suzuki  _et al._ (2017). The choice we make here is a Slater
determinant: this is the natural choice if the state Eq. (1) is reached from
applying an external field to a ground-state and then turning the field off.
The Slater determinant is the natural choice in most physical situations,
because they begin in the ground-state (see also discussion in Supplementary
Material). One would use ground-state DFT to find the initial KS orbitals, and
by the ground-state theorems, this is a Slater determinant. Since the KS
evolution involves a one-body Hamiltonian, the state remains a single Slater
determinant. For our two-electron spin-singlet system, this means that we
always have a single spatial KS orbital that is doubly-occupied, and must have
the form:

$\varphi({\bf r},t)=\sqrt{n({\bf r},t)/2}e^{i\alpha({\bf r},t)}$ (3)

to reproduce the exact interacting density of Eq. 2 with the phase $\alpha$
related to the current ${\bf j}$ through the equation of continuity,

$\nabla\cdot{\bf j}=\nabla\cdot\left(n({\bf r},t)\nabla\alpha({\bf
r},t)\right)=-\frac{\partial}{\partial{t}}n({\bf r},t)$ (4)

Inverting the TDKS equation yields the exact KS potential:

$v_{s}({\bf r},t)=\frac{\nabla^{2}\sqrt{n({\bf r},t)}}{2\sqrt{n({\bf
r},t)}}-\frac{|\nabla\alpha({\bf r},t)|^{2}}{2}-\frac{\partial\alpha({\bf
r},t)}{\partial t}$ (5)

The exact xc potential is then obtained from

$v_{\scriptscriptstyle\rm XC}({\bf r},t)=v_{\scriptscriptstyle\rm S}({\bf
r},t)-v_{\scriptscriptstyle\rm H}({\bf r},t)-v_{\rm ext}({\bf r},t)$ (6)

with the Hartree potential $v_{\scriptscriptstyle\rm H}({\bf
r},t)=\int\frac{n({\bf r}^{\prime},t)}{|{\bf r}-{\bf r}^{\prime}|}d^{3}{\bf
r}^{\prime}$ and external potential $v_{\rm ext}({\bf r},t)=-2/|{\bf r}|$.
Further, one can isolate the correlation component by noting that for our
choice of KS state, $v_{\scriptscriptstyle\rm X}({\bf
r},t)=-v_{\scriptscriptstyle\rm H}({\bf r},t)/2$.

Thus, finding the exact xc potential reduces to solving Eq. (4) for
$\alpha({\bf r},t)$. We note here that for a different choice of initial KS
state, e.g. using a two-configuration state that is more similar to that of
the actual interacting state, the inversion to find $v_{\scriptscriptstyle\rm
XC}$ involves an iterative numerical procedure Ruggenthaler and van Leeuwen R.
(2011); Nielsen  _et al._ (2013); Ruggenthaler  _et al._ (2015); some examples
for the 1D analog of the dynamics here can be found in Refs. Fuks  _et al._
(2016, 2018); Lacombe and Maitra (2019). This could be a more natural state to
begin the KS calculation in some situations, e.g. if the state was prepared in
such a superposition at the initial time, however it is inaccessible in a KS
evolution that begins in the ground state, as discussed earlier. The
importance of judiciously choosing the KS initial state when using an
adiabatic approximation has been realized and exploited in strong-field
charge-migration simulations Bruner  _et al._ (2017); Chen and Lopata (2020);
Folorunso  _et al._ (2021).

Eq. (4) has the form of a Sturm-Liouville equation, which has a unique
solution for $\alpha({\bf r},t)$ for a given boundary condition. Thanks to the
azimuthal symmetry of our density ( $M=0$ at all times), we need solve this in
effectively two dimensions. We construct an explicit matrix representation of
the operator $\nabla\cdot n({\bf r},t)\nabla$ in spherical coordinates using
the fourth order finite difference scheme subject to the following boundary
conditions:

$\alpha({\bf
r}\to\infty,t)=0\;\;\;\;\text{and}\;\;\;\;\frac{\partial}{\partial{\theta}}\alpha({\bf
r},t)|_{\theta=\pi,0}=0\,.$ (7)

Choosing this boundary condition at $t=0$ yields $\alpha({\bf r},0)=0$ since
initially the current is zero, and fixes our initial state as $\phi({\bf
r},0)=\sqrt{n({\bf r},0)/2}$. The Runge-Gross theorem then ensures that there
is a unique $v_{\scriptscriptstyle\rm XC}({\bf r},t)$ that reproduces the
exact $n({\bf r},t)$ and yields a unique $\alpha({\bf r},t)$ at later times
111It is not a priori obvious that the unique solution to Eq. (4) with the
boundary-condition Eq. (7) applied with time $t$ as a parameter is compatible
with the TDKS evolution, but as the results evolve smoothly in time, it is..

Subject to the boundary conditions Eq. 7, the numerical inversion of the
matrix operator $\nabla\cdot n({\bf r},t)\nabla$ results in the solution of
Eq. (4) for $\alpha({\bf r},t)$ which in turn when inserted into Eq. (5)
yields the KS potential, $v_{\scriptscriptstyle\rm S}({\bf r},t)$ (some
details in Appendix B).

##  III Results

Several symmetry features of the dynamics of our system simplify the analysis.
The azimuthal symmetry mentioned earlier together with the fact that the
chosen superposition is one of an $L=0$ and $L=1$ state, mean that the
density, current, and potentials in the lower half-plane ($\pi/2<\theta<\pi)$
exactly follow those in the upper half-plane ($0<\theta<\pi/2)$ a half-cycle
out of phase, $O(r,\pi-\theta,t)=O(r,\theta,t+T/2)$. (See also movies of the
density, current-density, and correlation potentials in Supplementary
Material). Further the simple form of the superposition means that $O({\bf
r},T-t)=O({\bf r},t)$. Thus we show time-snapshots only over a half cycle in
the lower octant.

Figure 1: Correlation potential, $v_{\scriptscriptstyle\rm C}({\bf r},t)$ at
$t=0$ for the 50:50 superposition case ($\alpha=1$) in the range
$\pi/2<\theta<\pi$ at times $t=0,T/8,T/4,3T/4$.

Figure 1 shows snapshots of the correlation potential indicated by fractions
of the period of oscillation, $T=2\pi/\omega=8.057$ a.u. One immediately
notices the unmistakable presence of the step and peak features in the exact
correlation potential that have been shown to arise in many 1D model systems
Lacombe and Maitra (2019, 2020); Luo  _et al._ (2014); Suzuki  _et al._
(2017); Elliott  _et al._ (2012); Fuks  _et al._ (2013); Lacombe  _et al._
(2018); Ramsden and Godby (2012); Fuks  _et al._ (2016, 2018); Hodgson  _et
al._ (2016, 2013). The step and peak feature is initially most prominent in
the region swept by $\pi/2<\theta<\pi$, and then decreases in magnitude,
gliding out of this region and appearing on the other side of $\theta=\pi/2$
at $t=T/2$. As in the 1D case, this time-dependent step has a spatially
nonlocal and non-adiabatic dependence on the density and is completely
unaccounted for in the adiabatic approximations: They are missing even in the
exact adiabatic approximation, i.e. evaluating the exact ground-state xc
potential on the instantaneous density Elliott  _et al._ (2012); Fuks  _et
al._ (2013); Luo  _et al._ (2014); Suzuki  _et al._ (2017); Lacombe  _et al._
(2018). These features often dominate the KS potential (see Fig. 2) and have
been shown to be responsible for various errors in simulations using adiabatic
approximations in 1D e.g. Suzuki  _et al._ (2017); Lacombe  _et al._ (2018).
Here we find they persist just as vigorously, with the same order of
magnitude, in real 3D atoms driven far from their ground-state. This justifies
the relevance for real systems of the conclusions drawn from the 1D studies,
and shows that such strong correlation effects are not a consequence of
reduced dimensionality, as might have been assumed from ground-state systems
Brus (2014). These dynamical steps are distinct to those arising from
fractional charges Thiele  _et al._ (2008); Kraisler  _et al._ (2021); Maitra
(2017), or in response situations van Gisbergen  _et al._ (1999), as in the 1D
case, and we expect they generically appear when a system is driven far from
its ground-state.

For the dynamics of this particular superposition state, at any instant of
time, the correlation potential asymptotes to the same (time-dependent) value
in every direction in the lower half plane, while asymptoting to a different
value in the upper half plane (recall $O({\bf r},T-t)=O({\bf r},t)$). At
$\theta=\pi/2$ there is a step and peak in $v_{\scriptscriptstyle\rm C}$ in
the $\theta-$direction which gives a force that ensures KS currents, like the
true currents, do not cross the $xy$-plane. This complex structure makes the
inversion numerically unreliable right at $\theta=\pi/2$. Along the
$\theta=\pi/2$ plane, the density of the $P$ state vanishes, and the large
change in the potential may be somewhat reminiscent of the abnormal divergent
behavior along the HOMO nodal plane found in the ground-state potential Gori-
Giorgi  _et al._ (2016). Unlike the ground-state case, however, our density
does decay differently along that plane than in other directions, and moreover
cannot be captured by any adiabatic approximation.

Decomposing the terms in Eq. (5), we find that the peak tends to arise from
the second term while the third term results in the step. Because the KS
current ${\bf j}_{\scriptscriptstyle\rm S}=n\nabla\alpha$, the second term and
the peak are related to the local velocity ${\bf j}_{\scriptscriptstyle\rm
S}/n$, while the step when a cut is taken across a fixed $\theta$ is related
to the radial integral of the local acceleration,
$\dot{\alpha}(r,t)=\int^{r}\nabla\dot{\alpha}(r^{\prime},t).\hat{r}^{\prime}dr^{\prime}$.

We note that the appearance of such dominating steps in the correlation
potential is fundamentally linked to the difference in configurations of the
interacting and KS states: the interacting system is a superposition of a
ground and excited state, quite distinct from the KS Slater determinant
structure. Tuning down the electron-interaction dampens the peak structure but
the step remains.

Fig. 2 shows the components of the exact KS potential. We observe that in the
central region where most of the density is localized, the force from the
correlation potential is much smaller than that from the exchange and Hartree
terms. It is in fact of similar magnitude to that in the ground-state Umrigar
and Gonze (1994). Near the density-minimum, where the excited state begins to
dominate over the ground-state, the correlation potential rises, and then
falls before leveling out. In this region, the slopes are such that the step
appears to be keeping different parts of the density separate, while the peak
corrects for dynamical Coulombic electron-interaction effects. The lack of
these features in the adiabatic approximations suggests that the resulting
densities will not be as structured, and will tend to underestimate
oscillation amplitudes in the dipole moment (as seen in the 1D case Fuks  _et
al._ (2018)).

Figure 2: Snapshots of potentials $v_{\scriptscriptstyle\rm
S},v_{\scriptscriptstyle\rm HXC},v_{\scriptscriptstyle\rm
XC},v_{\scriptscriptstyle\rm C}$ and density $4\pi r^{2}n$ at $\theta=3\pi/4$,
at times $t=0,T/8,T/4,3T/8$.

Taking different superpositions of the ground and excited states shows that
the step and peak features are universally present in real 3D systems. Figure
3 shows the KS and correlation potentials at the initial time, when $a$ in Eq.
1 is changed through $0,0.5,1,2,\infty$ . We see that, for finite values of
$a$, as the fraction of excited state is increased the step and peak decrease
in magnitude but extend over a larger region and move inward where more of the
density is. The very sharp peak and large step seen when $a=0.5$ (note it is
scaled to fit on the plot) occur at a sharp minimum of the density and has a
smaller impact on the ensuing dynamics than the softer but still prominent
structures at large $a$ occurring in regions of greater density. When the
excited state is fully occupied ($a=\infty$) the KS potential is such to
maintain the constant excited ${}^{1}P$ density at all times with a non-
interacting doubly-occupied orbital, and the structure is not unlike that seen
in the corresponding 1D excited helium atom of Ref. Elliott and Maitra (2012),
in both magnitude and shape; again, even the adiabatically-exact potential
would have a completely different structure.

Figure 3: Top panels: Left, initial densities for superposition states, 20:80
($a=2$), 50:50 ($a=1)$ and 80:20 ($a=1/2)$. Right, ground state ($a=0$) and
excited state ($a=\infty$) densities. Lower panels: The corresponding KS and
correlation potentials.

##  IV True and Kohn-Sham Current

Finally, we ask, how closely does the exact KS system reproduce the exact
interacting current in this case? It was recognized in the early days of
TDDFT, that the exact KS current could differ from the true current by a
rotational component Thiele and Kümmel (2009); Maitra  _et al._ (2002a); Gross
and Maitra (2012); D’Agosta and Vignale (2005); Schaffhauser and Kümmel
(2016), but how large this difference could be for realistic systems in the
non-perturbative regime was unknown. The KS and true currents are equal in
their longitudinal component, thanks to the equation of continuity,
$\nabla\cdot{\bf j}=-\frac{\partial n}{\partial t}$, but they can differ in
their rotational component. Indeed, for the two-electron singlet case with the
KS system represented by a Slater determinant, it follows from Eq. 3 that
${\bf j}_{s}({\bf r},t)=n({\bf r},t)\nabla\alpha({\bf r},t)$. This implies
that the true current-density would need to satisfy $\nabla\times({\bf
j}/n)=0$ in order for the KS current to possibly be equal to it. For our
chosen $M=0$ superposition, the rotational component of the current (which has
only an azimuthal component) is comparable to the longitudinal component (a
movie is given in the Supplementary Material), and differs from the KS
current. The fractional difference in the rotational component increases from
about $10\%$ to $20\%$ as we increase the proportion of the excited state in
Eq. (1) from $a=1$ to $2$. This is shown in Fig. 4, where the fractional
difference in the curl, $\delta(\nabla\times{\bf j})=\frac{\nabla\times({\bf
j}_{\scriptscriptstyle\rm S}-{\bf j}))}{\nabla\times{\bf j}}$, is contrasted
with that in the divergence $\delta(\nabla\cdot{\bf j})=\frac{\nabla\cdot({\bf
j}_{\scriptscriptstyle\rm S}-{\bf j})}{\partial n/\partial
t}=\frac{\nabla\cdot({\bf j}_{\scriptscriptstyle\rm S}+\partial n/\partial
t)}{\partial n/\partial t}$; the latter comes only from numerical error, and
is negligible except near the origin and at large $r$ where the denominator is
very small. We note that the curl of the current is only non-zero in the
azimuthal direction, and this is the component of the curl that is plotted in
the figure.

##  V Conclusion

Figure 4: Top panel: Divergence and the azimuthal component of the curl of the
true and KS currents for $a=1$ and $a=2$. Lower panel: Fractional difference
$\delta$ in the azimuthal component of curls and divergences of the current-
densities along $\theta=3\pi/4,t=T/8$

In summary, we have shown that the non-adiabatic dynamical features of the
exact correlation potential, previously seen to arise in 1D model systems
persist with comparable magnitudes in real 3D systems, and are not a
consequence of reduced dimensionality. The results inform the on-going
development of more accurate functionals in TDDFT that capture these features
Lacombe and Maitra (2019, 2020), pressing the case to go beyond adiabatic
functionals. Hybrid functionals, including range-separated ones, where non-
local density-dependence arises from the orbital dependence in exact exchange,
do not capture these features; this is particularly evident in the present
two-electron case, where exact-exchange simply cancels the self-interaction in
the Hartree potential. Furthermore, we have demonstrated that the true
interacting current differs from its KS counterpart, with the difference
depending on the relative proportions of ground and excited state composing
the state. The results thus advise caution when computing the current-density
from the KS orbitals; this will be inherently approximate even if the exact xc
functional was somehow known and used.

###### Acknowledgements.

Financial support from the National Science Foundation Award CHE-1940333 (DD)
and from the Department of Energy, Office of Basic Energy Sciences, Division
of Chemical Sciences, Geosciences and Biosciences under Award No. DESC0020044
(NTM, LL) are gratefully acknowledged. J.F. acknowledges financial support
from the European Research Council through grant ERC-2016-StG-714870, and by
the Spanish Ministry for Science, Innovation, and Universities – Agencia
Estatal de Investigación through grant RTI2018-099737-B-I00

## References

  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Maitra (2016) N. T. Maitra, J. Chem. Phys. 144, 220901 (2016).
  * Marques  _et al._ (2012) M. A. Marques, N. T. Maitra, F. M. Nogueira, E. K. Gross, and A. Rubio, eds.,  _Fundamentals of time-dependent density functional theory_ , Vol. 837 (Springer, 2012).
  * Ullrich (2011) C. A. Ullrich,  _Time-dependent density-functional theory: concepts and applications_ (Oxford University Press, 2011).
  * Draeger  _et al._ (2017) E. W. Draeger, X. Andrade, J. A. Gunnels, A. Bhatele, A. Schleife, and A. A. Correa, Journal of Parallel and Distributed Computing 106, 205 (2017).
  * Sato (2021) S. A. Sato, Computational Materials Science , 110274 (2021).
  * Guandalini  _et al._ (2021) A. Guandalini, C. Cocchi, S. Pittalis, A. Ruini, and C. A. Rozzi, Phys. Chem. Chem. Phys. , (2021).
  * Singh  _et al._ (2021) N. Singh, P. Elliott, J. K. Dewhurst, and S. Sharma, Phys. Rev. B 103, 134402 (2021).
  * Folorunso  _et al._ (2021) A. S. Folorunso, A. Bruner, F. m. c. Mauger, K. A. Hamer, S. Hernandez, R. R. Jones, L. F. DiMauro, M. B. Gaarde, K. J. Schafer, and K. Lopata, Phys. Rev. Lett. 126, 133002 (2021).
  * Yost  _et al._ (2017) D. C. Yost, Y. Yao, and Y. Kanai, Phys. Rev. B 96, 115134 (2017).
  * Raghunathan and Nest (2011) S. Raghunathan and M. Nest, J. Chem. Theory and Comput. 7, 2492 (2011).
  * Ramakrishnan and Nest (2012) R. Ramakrishnan and M. Nest, Phys. Rev. A 85, 054501 (2012).
  * Habenicht  _et al._ (2014) B. F. Habenicht, N. P. Tani, M. R. Provorse, and C. M. Isborn, J. Chem. Phys. 141, 184112 (2014).
  * Wijewardane and Ullrich (2008) H. O. Wijewardane and C. A. Ullrich, Phys. Rev. Lett. 100, 056404 (2008).
  * Gao  _et al._ (2017) C.-Z. Gao, P. M. Dinh, P.-G. Reinhard, and E. Suraud, Phys. Chem. Chem. Phys. 19, 19784 (2017).
  * Boström  _et al._ (2018) E. V. Boström, A. Mikkelsen, C. Verdozzi, E. Perfetto, and G. Stefanucci, Nano Letters 18, 785 (2018).
  * Krumland  _et al._ (2020) J. Krumland, A. M. Valencia, S. Pittalis, C. A. Rozzi, and C. Cocchi, The Journal of Chemical Physics 153, 054106 (2020), https://doi.org/10.1063/5.0008194 .
  * Quashie  _et al._ (2017) E. E. Quashie, B. C. Saha, X. Andrade, and A. A. Correa, Phys. Rev. A 95, 042517 (2017).
  * Gao  _et al._ (2014) C.-Z. Gao, J. Wang, F. Wang, and F.-S. Zhang, The Journal of Chemical Physics 140, 054308 (2014).
  * Henkel  _et al._ (2009) N. Henkel, M. Keim, H. J. Lüdde, and T. Kirchner, Phys. Rev. A 80, 032704 (2009).
  * Da  _et al._ (2017) B. Da, J. Liu, M. Yamamoto, Y. Ueda, K. Watanabe, N. T. Cuong, S. Li, K. Tsukagoshi, H. Yoshikawa, H. Iwai, S. Tanuma, H. Guo, Z. Gao, X. Sun, and Z. Ding, Nat. Commun. 8, 15629 (2017).
  * Yao  _et al._ (2021) X. Yao, Y. Lee, D. Ceresoli, and K. Cho, The Journal of Physical Chemistry A 125, 4524 (2021), pMID: 34019398, https://doi.org/10.1021/acs.jpca.0c11309 .
  * Lacombe and Maitra (2020) L. Lacombe and N. T. Maitra, Faraday Discuss. 224, 382 (2020).
  * Luo  _et al._ (2014) K. Luo, J. I. Fuks, E. D. Sandoval, P. Elliott, and N. T. Maitra, J. Chem. Phys 140, 18A515 (2014).
  * Suzuki  _et al._ (2017) Y. Suzuki, L. Lacombe, K. Watanabe, and N. T. Maitra, Phys. Rev. Lett. 119, 263401 (2017).
  * Elliott  _et al._ (2012) P. Elliott, J. I. Fuks, A. Rubio, and N. T. Maitra, Phys. Rev. Lett. 109, 266404 (2012).
  * Fuks  _et al._ (2013) J. I. Fuks, P. Elliott, A. Rubio, and N. T. Maitra, J. Phys. Chem. Lett. 4, 735 (2013).
  * Lacombe and Maitra (2019) L. Lacombe and N. T. Maitra, Journal of Chemical Theory and Computation 15, 1672 (2019), https://doi.org/10.1021/acs.jctc.8b01159 .
  * Lacombe  _et al._ (2018) L. Lacombe, Y. Suzuki, K. Watanabe, and N. T. Maitra, The European Physical Journal B 91, 96 (2018).
  * Ramsden and Godby (2012) J. Ramsden and R. Godby, Phys. Rev. Lett. 109, 036402 (2012).
  * Fuks  _et al._ (2016) J. I. Fuks, S. Nielsen, M. Ruggenthaler, and N. T. Maitra, Phys. Chem. Chem. Phys. 18, 20976 (2016).
  * Fuks  _et al._ (2018) J. I. Fuks, L. Lacombe, S. E. B. Nielsen, and N. T. Maitra, Phys. Chem. Chem. Phys. 20, 26145 (2018).
  * Hodgson  _et al._ (2016) M. J. P. Hodgson, J. D. Ramsden, and R. W. Godby, Phys. Rev. B 93, 155146 (2016).
  * Hodgson  _et al._ (2013) M. J. P. Hodgson, J. D. Ramsden, J. B. J. Chapman, P. Lillystone, and R. W. Godby, Phys. Rev. B 88, 241102 (2013).
  * D’Agosta and Vignale (2005) R. D’Agosta and G. Vignale, Phys. Rev. B 71, 245103 (2005).
  * Ruggenthaler  _et al._ (2013) M. Ruggenthaler, S. E. B. Nielsen, and R. van Leeuwen, Phys. Rev. A 88, 022512 (2013).
  * Elliott and Maitra (2012) P. Elliott and N. T. Maitra, Phys. Rev. A 85, 052510 (2012).
  * Thiele  _et al._ (2008) M. Thiele, E. K. U. Gross, and S. Kümmel, Phys. Rev. Lett. 100, 153004 (2008).
  * Ullrich (2006) C. A. Ullrich, The Journal of Chemical Physics 125, 234108 (2006), https://doi.org/10.1063/1.2406069 .
  * Maitra  _et al._ (2002a) N. T. Maitra, K. Burke, H. Appel, and E. Gross, Reviews of Modern Quantum Chemistry 2, 1186 (2002a).
  * Schaffhauser and Kümmel (2016) P. Schaffhauser and S. Kümmel, Phys. Rev. B 93, 035115 (2016).
  * Thiele and Kümmel (2009) M. Thiele and S. Kümmel, Phys. Rev. A 79, 052503 (2009).
  * Gross and Maitra (2012) E. Gross and N. Maitra, in  _Fundamentals of Time-Dependent Density Functional Theory_ , edited by M. A. Marques, N. T. Maitra, F. M. Nogueira, E. Gross, and A. Rubio (Springer Berlin Heidelberg, 2012) pp. 53–99.
  * Andrade  _et al._ (2018) X. Andrade, S. Hamel, and A. Correa, Eur. Phys. J. B. 91, 229 (2018).
  * Takeuchi  _et al._ (2019) T. Takeuchi, M. Noda, and K. Yabana, ACS Photonics 6, 2517 (2019), https://doi.org/10.1021/acsphotonics.9b00889 .
  * Feist (2009) J. Feist,  _Two-photon double ionization of helium, PhD Thesis_ (2009).
  * Feist  _et al._ (2008) J. Feist, S. Nagele, R. Pazourek, E. Persson, B. I. Schneider, L. A. Collins, and J. Burgdörfer, Phys. Rev. A 77, 043420 (2008).
  * Maitra  _et al._ (2002b) N. T. Maitra, K. Burke, and C. Woodward, Phys. Rev. Lett. 89, 023002 (2002b).
  * van Leeuwen (1999) R. van Leeuwen, Phys. Rev. Lett. 82, 3863 (1999).
  * Ruggenthaler and van Leeuwen R. (2011) M. Ruggenthaler and van Leeuwen R., Europhys. Lett. 95, 13001 (2011).
  * Nielsen  _et al._ (2013) S. E. B. Nielsen, M. Ruggenthaler, and R. van Leeuwen, EPL (Europhysics Letters) 101, 33001 (2013).
  * Ruggenthaler  _et al._ (2015) M. Ruggenthaler, M. Penz, and R. van Leeuwen, J. Phys. Condens. Matter 27, 203202 (2015).
  * Bruner  _et al._ (2017) A. Bruner, S. Hernandez, F. Mauger, P. M. Abanador, D. J. LaMaster, M. B. Gaarde, K. J. Schafer, and K. Lopata, The Journal of Physical Chemistry Letters 8, 3991 (2017).
  * Chen and Lopata (2020) M. Chen and K. Lopata, Journal of Chemical Theory and Computation 16, 4470 (2020), pMID: 32470295, https://doi.org/10.1021/acs.jctc.0c00122 .
  * Note (1) It is not a priori obvious that the unique solution to Eq. (4) with the boundary-condition Eq. (7) applied with time $t$ as a parameter is compatible with the TDKS evolution, but as the results evolve smoothly in time, it is.
  * Brus (2014) L. Brus, Accounts of Chemical Research 47, 2951 (2014), pMID: 25120173, https://doi.org/10.1021/ar500175h .
  * Kraisler  _et al._ (2021) E. Kraisler, M. J. P. Hodgson, and E. K. U. Gross, Journal of Chemical Theory and Computation 17, 1390 (2021), pMID: 33595312, https://doi.org/10.1021/acs.jctc.0c01093 .
  * Maitra (2017) N. T. Maitra, Journal of Physics: Condensed Matter 29, 423001 (2017).
  * van Gisbergen  _et al._ (1999) S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Gritsenko, E. J. Baerends, J. G. Snijders, B. Champagne, and B. Kirtman, Phys. Rev. Lett. 83, 694 (1999).
  * Gori-Giorgi  _et al._ (2016) P. Gori-Giorgi, T. Gál, and E. J. Baerends, Molecular Physics 114, 1086 (2016), https://doi.org/10.1080/00268976.2015.1137643 .
  * Umrigar and Gonze (1994) C. J. Umrigar and X. Gonze, Phys. Rev. A 50, 3827 (1994).


11institutetext: Imperial College London

# Automated Enriched Medical Concept Generation for Chest X-ray Images

Aydan Gasimova

###### Abstract

Decision support tools that rely on supervised learning require large amounts
of expert annotations. Using past radiological reports obtained from hospital
archiving systems has many advantages as training data above manual single-
class labels: they are expert annotations available in large quantities,
covering a population-representative variety of pathologies, and they provide
additional context to pathology diagnoses, such as anatomical location and
severity. Learning to auto-generate such reports from images present many
challenges such as the difficulty in representing and generating long,
unstructured textual information, accounting for spelling errors and
repetition/redundancy, and the inconsistency across different annotators. We
therefore propose to first learn visually-informative medical concepts from
raw reports, and, using the concept predictions as image annotations, learn to
auto-generate structured reports directly from images. We validate our
approach on the OpenI [2] chest x-ray dataset, which consists of frontal and
lateral views of chest x-ray images, their corresponding raw textual reports
and manual medical subject heading (MeSH®) annotations made by radiologists.

###### Keywords:

nlp medical imaging deep learning

##  1 Introduction

Radiologists are faced daily with the very time-consuming and repetitive task
of looking at hundreds of radiography images and writing up radiological
reports. The fast turn-arounds they are expected to produce leads to fatigue
that can negatively affect diagnostic accuracy [19]. Supervised learning for
automated pathology detection from images has the potential for clinical-
decision support, however, such image segmentation and classification learning
tasks require detailed annotations covering a large distribution of input data
for the algorithms to be able to make robust predictions. Such annotations
must be made by qualified radiologists, which, for the detail and breadth of
annotation required, will be an equally if not more time-consuming task than
manually creating the reports. In addition, classification and semantic
segmentation tasks only solve for the prediction of presence of pathologies,
and not the generation of reports which contain additional information such as
severity, location, and absence of pathologies.

Recently, we have seen supervised learning approaches that aim to take
advantage of past radiological exams containing reports in order to either
auto-generate the reports [17, 9, 23], or to assist in classification tasks
[16, 21, 20, 24, 22]. The noise present in medical reports in addition to the
presence of non-visually significant information, such as the negation of
pathologies, make it difficult to learn from them directly as done in natural
image captioning frameworks. Additionally, high recall/precision of
pathologies is more crucial in the medical domain where the risk of mis-
labelling is much higher.

We therefore propose to use a limited number of manual medical concept
annotations in order to first learn to extract them from raw reports, and then
take advantage of the model predictions as image annotations, thus providing a
method for augmenting an image-annotation dataset. We then demonstrate how
these image-concept annotations can be learned through sequence models
conditioned on image features, and generate a more readable context for the
diagnosis that can be used as part of a clinical decision support system, thus
greatly alleviating the burden on radiologists. Our approach can be summarised
in the following steps:

  1. 1.

We propose a network that learns to extract visually-significant medical
concepts from raw reports. To our knowledge, this is a first attempt to that
goes beyond simple pathology detection to include concepts such as anatomical
position and severity.

  2. 2.

We explore several sequence-learning networks that aim to condition the
sequence generation process on image features in order to learn to auto-
generate structured reports from radiological images.

  3. 3.

We use the predictions made by the structured report generation process in
step 1 to demonstrate how they can be used to create an image-report training
set for step 2.

##  2 Related Work

###  2.1 Data-mining Image Labels

There are two common approaches to extracting image labels from raw reports:
statistical and tool based. Radiological text mapping tools such as DNorm [13]
and MetaMap [3] have been used to extract labels for multi-label
classification [22] and in weakly supervised localisation learning frameworks
[21]. However, other biological concepts in the reports, such as location,
severity, and other visually descriptive features of the pathology are not
taken advantage of. Unsupervised, statistical methods such as latent Dirichlet
allocation [16] and clustering [20] have been used to implicitly define topics
and cluster groups containing key words and propose classification into these
topics and groups. These approaches are heavily dependent on the number of
topics/groups providing the lowest perplexity score, which can be a range of
values. In addition, these are not generative models, therefore reports can
only be selected based on nearest-neighbour methods. To this end, we propose
instead to learn to generate reports comprised of medical concept from images,
in a similar style to natural image captioning.

###  2.2 Radiology Report Generation

Closest to our work, Shin et al. proposed a cascaded learning framework to
auto-generate MeSH annotations from chest X-rays [17] whereby image embeddings
are first extracted from a pre-trained classification network, and then used
to initialise a sequence prediction model to auto-generate MeSH sequences.
Zhang et al. [24, 23] leverage manually created structured reports in a dual-
attention framework to improve features used for classifying histopathology
images and to provide interpretability to the classification. The reports used
in both cases are far more structured than their raw counterparts and so this
approach cannot be directly translated to hospital data. Training on raw
hospital reports, Jing et al. [9] demonstrated how they can be generated by
first training a multi-label CNN on the images and the Medical Text Indexer
(MTI) tags identified in the original raw reports of the Openi chest x-ray
dataset. However, reports can be very long and heterogeneous, and the authors
do not evaluate the model’s ability to determine whether visually and
clinically-relevant medical concepts have been identified. To address the
challenges of learning from raw reports directly, we first learn to generate
structured reports made up of only visually-significant medical concepts that
correspond directly to features seen in the images. Being shorter and
vocabulary-controlled, the generation process is easier to evaluate for
correct identification of pathologies.

##  3 Method

###  3.1 Enriched Concept Extraction from Raw Reports

We approach learning structured reports from raw textual reports as a multi-
label classification task since the vocabulary of MeSH terms is consistent
across annotators, and limited. We modify the shallow CNN first introduced by
Kim [10] for multi-class text classification and later adapted for multi-label
text classification by Liu at al. [14] by introducing a learn-able embedding
layer as we do not have the advantage of pre-trained word embeddings for
medical text, and by introducing dropout followed by a fully-connected layer
to each convolutional output prior to the concatenation to aid regularisation.

Let $\bm{x}_{i}\in\mathbb{R}^{d}$ be the d-dimensional word vector for the
$i$-th word of report $\bm{p}$. The textual report is thus represented as a
concatenation of word embeddings:
$\bm{p}=[\bm{x}_{0},...\bm{x}_{i},...\bm{x}_{M}]\in\mathbb{R}^{M\times d}$
where $M$ is the maximum length of the reports. A filter
$\bm{m}\in\mathbb{R}^{hd}$ is convolved with a window of $h$ words to produce
a new feature $c_{i}$:

$c_{j}=f(\bm{m}*\bm{x}_{i:i+h-1}+b)$ (1)

where $f$ is a non-linear activation function and $b$ is a bias term. The
filter is applied consecutively to every $h$-word window in the sentence,
resulting in a feature map
$\bm{c}=[c_{0},....c_{i},...c_{M-h+1}]\in\mathbb{R}^{M-h+1}$. Max-over-time
pooling [4] is applied over each feature map to capture the most important
feature $\hat{c}=max(\bm{c})$. In this way we apply many filter operations
with varying window widths in order to obtain multiple features that are able
to capture semantic information of reports with varying word lengths. We use
the sigmoid activation function as we require an independent prediction for
each class and train by minimising the multi-class sigmoid cross-entropy (SCE)
loss. In addition, we add terms to balance maximising the true positive class
prediction with true negative class predictions as the positive label space is
very sparse:

$\begin{split}\widehat{SCE}_{i}=&-\lambda_{1}\sum_{j=1}^{K}\left(y_{j}log(f(s_{ij}))+(1-y_{j}log(1-f(s_{j}))\right)\\\
&-\lambda_{2}\sum_{j=1}^{K}y_{j}f(s_{j})/\sum_{j=1}^{K}(y_{j}f(s_{j})+y_{j}(1-f(s_{j})))\\\
&-\lambda_{3}\sum_{j=1}^{K}(1-y_{j})(1-f(s_{j}))/\sum_{j=1}^{K}((1-y_{j})(1-f(s_{j}))+(1-y_{j})(f(s_{j})))\end{split}$
(2)

where $K$ is the number of classes, $y_{j}$ is presence/absence of class label
$j$ for instance $i$, $f(s_{j})$ is the prediction for instance $i$ on label
$j$ made through a sigmoid activation:

$f(s_{i})=\frac{1}{1+e^{-s_{i}}}$ (3)

The weights of each loss term, $\lambda_{1},\lambda_{2},\lambda_{3}$ are non-
negative, sum to 1 and chosen through cross-validation. Finally, the modified
SCE loss is averaged over batches.

###  3.2 Report Generation from Images

Given that it is possible to learn structured report outputs from raw reports,
we propose a method of learning to auto-generate structured reports directly
from images. We explore multiple ways of conditioning the MeSH sequence
learning on the image embeddings that aims to maintain the dependency between
the word generation process and the image embedding at every time-step. The
MeSH sequence is modelled using an RNN, specifically the Long Short-Term
Memory (LSTM) implementation proposed in [8]. Each LSTM unit has three sigmoid
gates to control the internal state: ‘input’, ‘output’ and ‘forget’. At each
time step, the gates control how much of the previous time steps is propagated
through to determine the output. For an input word sequence
$\\{x_{1},\dots,x_{n}\\}$ where $x_{i}\in\mathbb{R}^{d}$, the internal hidden
state $h_{t}\in\mathbb{R}^{h}$ and memory state $m_{t}\in\mathbb{R}^{m}$ are
updated as follows:

$\begin{split}h_{t}&=f_{t}\odot
h_{t-1}+i_{t}\odot\tanh(W^{(hx)}x_{t}+W^{(hm)}m_{t-1})\\\
m_{t}&=o_{t}\odot\tanh(h_{t})\end{split}$ (4)

where $x_{t}\in\mathbb{R}^{D}$ is the word embedding, $W^{(hx)}$ and
$W^{(hm)}$ are the trainable weight parameters, and $i_{t}$, $o_{t}$ and
$f_{t}$ are the input, output and forget gates respectively. Bias terms are
left out for readability.

The image embedding, $\bm{im}_{i}=\text{CNN}(I)$ where
$\bm{im}_{i}\in\mathbb{R}^{g}$ is extracted from the final spatial-average
pooling layer of the pre-trained CNN. We explore three ways of conditioning
the sequence learning process:

  1. 1.

RNN0: The image embedding is projected into the same embedding space as the
word embeddings via a dense transition layer:
$\bm{im}=\text{relu}(W^{(dg)}\text{CNN}(I))$. The image embedding is
concatenated with the word sequence and thus treated as the first ‘word’ in
the MeSH sequence.

  2. 2.

RNN1: The image embedding is projected via a dense transition layer into a
fixed embedding width and combined with the output of the recurrent layer
through either concatenation or summation operation, and passed to the decoder
$dec$:

$\bm{dec}_{t}=\text{relu}(W^{(z)}(o_{t}*\text{relu}(W^{(dg)}\text{CNN}(I)))$
(5)

where $*$ represents concatenation or summation and $W^{z}$ are the weights of
the decoder.

  3. 3.

RNN2: The image embedding is projected via a dense transition layer into a
fixed embedding width and combined with the input of the recurrent layer
through either concatenation or summation operation, and passed to the encoder
$enc$:

$\bm{enc}_{t}=\text{relu}(W^{(a)}(x_{t}*\text{relu}(W^{(dg)}\text{CNN}(I)))$
(6)

where $W^{a}$ are the weights of the encoder.

The model architectures are illustrated in Figure 1. For all models, the
decoder outputs are passed to the prediction layer $s(t)=f(W^{T}x_{t})$ where
$f$ is the softmax function. The models are all trained by minimising the
cross-entropy loss between the output and true sequence:

$L(S,I)=-\sum_{t=0}^{T}\log p(P_{t}=T_{t}|\text{CNN(I)},P_{0}\ldots P_{t-1})$
(7)

where $p$ is the probability that the predicted word $P_{t}$ equals the true
word $T_{t}$ at time step $t$ given image features $\text{CNN}(I)$ and
previous words $P_{0}\ldots P_{t-1}$, and $T$ is the LSTM sequence length.

Figure 1: Image-MeSH sequence learning model architectures.

##  4 Experiments

###  4.1 Dataset

We evaluated our models on the OpenI [2] Indiana U. Chest X-ray Collection.
This dataset consists of 7,470 frontal and lateral chest x-ray images and
3,955 associated radiological reports from the hospital’s picture archiving
systems. They have all been fully anonymised to remove patient names. In
addition to the raw text reports, each exam has MeSH annotations made by
qualified radiologists. MeSH annotations are (with some exceptions) formatted
as [pathology/description,… pathology/description] where _description_ is a
combination of _anatomy/position/severity_. The number of captions per image
is on average 2.33, with an average of 2.68 MeSH terms per caption.

####  4.1.1 Preprocessing

This involved lower-casing, punctuation and non-alpha-numeric character
removal from reports and MeSH. We limit the MeSH annotations to just one
_pathology/description_ pair by selecting the caption with the most common
pathology. Additionally, as the negation of pathologies was generally standard
across reports, we performed negation removal using regex. Finally, the text
reports were cropped/padded to 32 words based on the average of 20.23 +1 std
of 11.9. MeSH captions were cropped/padded to length 5 based on average+1 std.
Empty reports were removed. This resulted in 3,023 unique report-MeSH term
pairs, of which 300 were randomly selected for validation and 300 for test.

###  4.2 Experimental Settings

We first investigate whether the structured reports can be learned from raw
reports by creating a sub-set of size=1000 of the MeSH annotated reports, and
training the text CNN on the report-MeSH pairs in the sub-set. The trained
text CNN is then used to make MeSH prediction on the remaining set of raw
reports, and these (together with the gold-standard annotated sub-set from the
previous step) are used to train the image-MeSH sequence model. We compare
this to training the image-MeSH sequence model on the entire gold-standard
annotated set of 3,023.

####  4.2.1 Text CNN

For the text CNN model, we use rectified linear units as activation function
on the convolutional layers, one-dimensional convolutional filters of width
3,4,5 with 512 feature maps for each filter, dropout rate [5] of p=0.5, with
254 hidden units for the dense layer, and with
$\lambda_{1}=0.5,\lambda_{2}=0.2,\lambda_{3}=0.3$ for the loss terms. The
model was trained through batch backpropagation, batch size=128 and using Adam
optimisation [12] with learning rate=0.001 for 100 epochs with early-stopping.
To compensate for the class imbalance of ‘normal’ vs. diseased cases, we
select batches with uniform distribution over the classes, augmenting the
instances by sentence-shuffling.

####  4.2.2 Sequence Models

Image embeddings are extracted from the last average pooling layers of Vgg16
[18] and Resnet50 [7] models, pre-trained on ImageNet [6] to extract
$\bm{im}\in\mathbb{R}^{4096}$ and $\bm{im}\in\mathbb{R}^{2048}$ respectively.
For RNN0, the joint image-word embedding dim is set to 2048 for the Vgg input,
and 1024 for Resnet. For RNN1 and RNN2, the dense transition layer dimension
is set to 1024. For all the sequence models, the LSTM hidden state is set to
dim 512, and the LSTM units are unrolled up to 6 time steps (1 for the start
token, and 5 for MeSH sequence). All models are trained with batch size 128,
using Adam optimisation [11], learning rate=0.001 and early-stopping.

Table 1: Text CNN classification metrics for sub-sampled and full gold-
standard annotated data. Metrics reported on test data. All Classes
---
Training sample size | Acc. | R | R-OC | R-OS | P | P-OC | P-OS
1000 | 98.26 | 67.82 | 40.49 | 65.79 | 70.15 | 45.15 | 67.18
3023 (all) | 99.48 | 92.07 | 84.50 | 91.77 | 89.90 | 84.82 | 89.47
Pathology Classes
Training sample | Acc. | R | R-OC | R-OS | P | P-OC | P-OS
1000 | 98.64 | 67.41 | 44.72 | 60.00 | 69.73 | 44.57 | 56.22
3023 (all) | 99.54 | 90.74 | 86.94 | 80.67 | 88.45 | 86.35 | 78.72
Table 2: BLEU1-4 score comparisons on model in [17] and our RNN0, RNN1 and RNN2 trained on all gold-standard MeSH annotations, and trained on 1000 gold-standard MeSH annotations+textCNN predictions. | Train | Val | Test
---|---|---|---
Model | B-1/B-2/B-3/B-4 | B-1/B-2/B-3/B-4 | B-1/B-2/B-3/B-4
Learning to read [17] |  97.2 / 67.1 / 14.9 / 2.8 |  68.1 / 30.1 / 5.2 / 1.1 |  79.3 / 9.1 / 0.0 / 0.0
RNN0+vgg16+all | 8.8 / 1.8 / 0.7 / 0.2 | 7.8 / 2.6 / 1.1 / 1.4 | 6.9 / 2.3 / 0.7 / 0.1
RNN0+resnet50+all | 16.5 / 8.6 / 4.6 / 2.4 | 16.7 / 8.7 / 3.9 / 1.0 | 18.8 / 10.4 / 3.8 / 1.9
RNN1+resnet50+all | 77.9 / 45.6 / 29.4 / 18.9 | 65.7 / 51.6 / 30.2 / 17.1 | 66.7 / 47.1 / 26.8 / 15.9
RNN2+resnet50+all | 74.1 / 42.0 / 26.8 / 17.3 | 63.2 / 47.5 / 27.3 / 15.8 | 63.2 / 43.9 / 25.5 / 13.6
RNN0+resnet50+pred | 22.9 / 15.5 / 7.8 / 4.0 | 13.6 / 8.3 / 4.0 / 0.9 | 14.7 / 9.3 / 2.7 / 1.5
RNN1+resnet50+pred | 73.6 / 50.0 / 30.9 / 17.8 | 41.5 / 29.7 / 15.9 / 7.2 | 41.6 / 28.2 / 13.2 / 8.1
RNN2+resnet50+pred | 69.4 / 47.6 / 29.6 / 16.7 | 39.4 / 28.0 / 14.6 / 6.7 | 39.8 / 26.4 / 12.7 / 8.0

##  5 Results

####  5.0.1 Enriched Concept Extraction from Reports

We evaluate the MeSH term prediction from the text reports by calculating the
total binary accuracy (Acc), precision (P) and recall (R), and the mean-over-
class (P-OC, R-OC) and mean-over-samples (P-OS, R-OS) precision and recall of
the 102 classes. In addition, we report metrics of the ‘pathology’ classes
separately by manual allocation based on the definitions on the MeSH term
online library [1]. Complete metrics are compared in Table 1.

####  5.0.2 Report Generation from Images

During inference, the first word is sampled from the LSTM, concatenated to the
input, and used to predict consequent words. The quality of the generated
reports was evaluated by measuring BLUE [15] scores averaged over all the
reports, which are a form of $n$-gram precision commonly used for evaluating
image captioning as they maintain high correlation with human judgement. BLEU
scores of RNN0, RNN1 and RNN2 trained on all gold-standard annotations and on
the predictions made by the text CNN are presented in Table 2. RNN0 is the
same framework used in [17], however, they additionally train their model in a
cascaded fashion which significantly improves the model’s ability to predict
the first word, but struggles to maintain visual correspondence in generating
subsequent words, hence the steep reduction in higher $n$-gram precision.
Additionally, cascaded models suffer from error propagation during test time,
hence the poor performance on test data. RNN1 and RNN2 solve both problems by
conditioning the word generation process on the images at every time-step and
by being trained end-to-end, hence achieving higher $n$-gram scores on the
test data. In addition, we have shown that we can achieve comparably high BLEU
metrics when training on the predicted MeSH terms made by the text CNN.

##  6 Conclusion

We demonstrate how, given a small amount of manual annotations, clinically and
visually-important concepts can be learned from raw textual radiology reports.
We then demonstrate how these concepts can be used as radiological image
annotations and used in an image-sequence learning model to auto-generate
reports as part of a clinical decision support system.

## References

  * [1] Medical subject headings - national library of medicine, https://meshb.nlm.nih.gov/
  * [2] Open-i: An open access biomedical search engine., https://openi.nlm.nih.gov/
  * [3] Aronson, A.R.: Effective mapping of biomedical text to the umls metathesaurus: the metamap program. In: Proceedings of the AMIA Symposium. p. 17. American Medical Informatics Association (2001)
  * [4] Collobert, R., Weston, J., Bottou, L., Karlen, M., Kavukcuoglu, K., Kuksa, P.: Natural language processing (almost) from scratch. Journal of machine learning research 12(Aug), 2493–2537 (2011)
  * [5] Dahl, G.E., Sainath, T.N., Hinton, G.E.: Improving deep neural networks for lvcsr using rectified linear units and dropout. In: 2013 IEEE international conference on acoustics, speech and signal processing. pp. 8609–8613. IEEE (2013)
  * [6] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [7] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [8] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural computation 9(8), 1735–1780 (1997)
  * [9] Jing, B., Xie, P., Xing, E.: On the automatic generation of medical imaging reports. arXiv preprint arXiv:1711.08195 (2017)
  * [10] Kim, Y.: Convolutional neural networks for sentence classification. arXiv preprint arXiv:1408.5882 (2014)
  * [11] Kingma, D., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [12] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [13] Leaman, R., Khare, R., Lu, Z.: Challenges in clinical natural language processing for automated disorder normalization. Journal of biomedical informatics 57, 28–37 (2015)
  * [14] Liu, J., Chang, W.C., Wu, Y., Yang, Y.: Deep learning for extreme multi-label text classification. In: Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval. pp. 115–124. ACM (2017)
  * [15] Papineni, K., Roukos, S., Ward, T., Zhu, W.J.: Bleu: a method for automatic evaluation of machine translation. In: Proceedings of the 40th annual meeting on association for computational linguistics. pp. 311–318. Association for Computational Linguistics (2002)
  * [16] Shin, H.C., Lu, L., Kim, L., Seff, A., Yao, J., Summers, R.M.: Interleaved text/image deep mining on a large-scale radiology database for automated image interpretation. Journal of Machine Learning Research 17(1-31), 2 (2016)
  * [17] Shin, H.C., Roberts, K., Lu, L., Demner-Fushman, D., Yao, J., Summers, R.M.: Learning to read chest x-rays: Recurrent neural cascade model for automated image annotation. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2497–2506 (2016)
  * [18] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556 (2014)
  * [19] Stec, N., Arje, D., Moody, A.R., Krupinski, E.A., Tyrrell, P.N.: A systematic review of fatigue in radiology: is it a problem? American Journal of Roentgenology 210(4), 799–806 (2018)
  * [20] Wang, X., Lu, L., Shin, H.C., Kim, L., Bagheri, M., Nogues, I., Yao, J., Summers, R.M.: Unsupervised joint mining of deep features and image labels for large-scale radiology image categorization and scene recognition. In: 2017 IEEE Winter Conference on Applications of Computer Vision (WACV). pp. 998–1007. IEEE (2017)
  * [21] Wang, X., Peng, Y., Lu, L., Lu, Z., Bagheri, M., Summers, R.M.: Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2097–2106 (2017)
  * [22] Yan, K., Peng, Y., Sandfort, V., Bagheri, M., Lu, Z., Summers, R.M.: Holistic and comprehensive annotation of clinically significant findings on diverse ct images: Learning from radiology reports and label ontology. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 8523–8532 (2019)
  * [23] Zhang, Z., Chen, P., Sapkota, M., Yang, L.: Tandemnet: Distilling knowledge from medical images using diagnostic reports as optional semantic references. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 320–328. Springer (2017)
  * [24] Zhang, Z., Xie, Y., Xing, F., McGough, M., Yang, L.: Mdnet: A semantically and visually interpretable medical image diagnosis network. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 6428–6436 (2017)


# Unidirectional propagation of zero-momentum magnons

Ondřej Wojewoda<EMAIL_ADDRESS>CEITEC BUT, Brno University of
Technology, Purkyňova 123, Brno, 612 00, Czech Republic
Jakub Holobrádek  CEITEC BUT, Brno University of Technology, Purkyňova 123,
Brno, 612 00, Czech Republic
Dominik Pavelka  Faculty of Mechanical Engineering, Institute of Physical
Engineering, Brno University of Technology, Technická 2, Brno, 616 69, Czech
Republic  Ekaterina Pribytova  Faculty of Mechanical Engineering, Institute of
Physical Engineering, Brno University of Technology, Technická 2, Brno, 616
69, Czech Republic  Jakub Krčma  Faculty of Mechanical Engineering, Institute
of Physical Engineering, Brno University of Technology, Technická 2, Brno, 616
69, Czech Republic  Jan Klíma  Faculty of Mechanical Engineering, Institute of
Physical Engineering, Brno University of Technology, Technická 2, Brno, 616
69, Czech Republic  Jan Michalička  CEITEC BUT, Brno University of Technology,
Purkyňova 123, Brno, 612 00, Czech Republic
Tomáš Lednický  CEITEC BUT, Brno University of Technology, Purkyňova 123,
Brno, 612 00, Czech Republic
Andrii V. Chumak  Faculty of Physics, University of Vienna, A-1090 Wien,
Austria  Michal Urbánek<EMAIL_ADDRESS>CEITEC BUT, Brno
University of Technology, Purkyňova 123, Brno, 612 00, Czech Republic
Faculty of Mechanical Engineering, Institute of Physical Engineering, Brno
University of Technology, Technická 2, Brno, 616 69, Czech Republic

###### Abstract

We report on experimental observation of unidirectional propagation of zero-
momentum magnons in synthetic antiferromagnet consisting of strained
CoFeB/Ru/CoFeB trilayer. Inherent non-reciprocity of spin waves in synthetic
antiferromagnets with uniaxial anisotropy results in smooth and monotonous
dispersion relation around Gamma point, where the direction of the phase
velocity is reversed, while the group velocity direction is conserved. The
experimental observation of this phenomenon by intensity-, phase-, and time-
resolved Brillouin light scattering microscopy is corroborated by analytical
models and micromagnetic simulations.

††preprint: APS/123-QED

Synthetic antiferromagnets (SAFs) recently attracted a lot of attention in the
magnonics community, as they can provide additional degrees of freedom for
spin wave dispersion engineering compared to conventional ferromagnetic
systems [1, 2, 3, 4, 5, 6, 7]. One of the most prominent features of
dispersion relation in SAFs is non-reciprocity caused by the dipolar field
interaction between the layers [8, 9, 10, 11, 12, 13, 14]. Spin wave non-
reciprocity and symmetry breaking in material systems are of fundamental
interest and can also be exploited in practical magnonic applications, such as
diodes or circulators. Non-reciprocal dispersion can be achieved also in
ferromagnets, e.g., by exploiting interfacial effects [15, 16, 17, 18], or
current induced Doppler shift [19, 20]. However, the non-reciprocity in SAFs
is much more pronounced. Another way to create non-reciprocal spin-wave system
is to use a grating coupler [13, 21], however, these systems can be tuned only
to a specific $k$-vector and frequency.

In this Letter, we report on the observation of smooth and monotonous spin
wave dispersion in the vicinity of $\Gamma$ point ($k$=0) in SAFs. This
property leads to unidirectional spin wave propagation over a wide range of
frequencies. The group velocity points in the same direction at both sides of
the $\Gamma$ point, whereas the phase velocity is reversed upon the $\Gamma$
point crossing (see Fig. 1a), resulting in a situation where we can observe
unidirectional energy flow with two opposing phase velocities (wavevectors).
The smooth and monotonous dispersion causes that even zero-momentum spin waves
(with $k$=0) propagate with nonzero group velocity. This phenomenon exhibits
itself as a change of the sign of phase evolution for the same direction of
the spin wave propagation, and is directly observed by phase-resolved
Brillouin light scattering (BLS) microscopy (Fig. 1a).

To study this behavior we prepared a sample consisting of a multilayer stack
of two 15 nm thick CoFeB layers separated by a 0.6 nm thick Ru spacer layer
which were grown by magnetron-sputtering on Si(100) substrate covered by 285
nm thick thermal SiO2 layer (Fig. 1b). On top of the SAF layer, we patterned a
500 nm wide microwave antenna for spin wave excitation using electron beam
lithography and lift-off process. The antenna consists of a stack of Ti 5
nm/SiO2 10 nm/Ti 5 nm/Au 50 nm thin films deposited by e-beam evaporation. The
SAF layer was imaged by scanning transmission electron microscopy with a high-
angle annular dark-field detector (STEM-HAADF) showing individual layers with
Z-contrast and its chemical composition was investigated by energy-dispersive
X-ray spectroscopy (STEM-EDX), see Fig. 1c. The STEM image shows good quality
of the interfaces and measured distances of 15 nm between Ru interlayers
revealed by the EDX chemical composition profile well correspond with the
nominal thicknesses of the layers [$t_{\mathrm{CoFeB}}=(15\pm 1)\,$nm].

Figure 1:  Experimental geometry. (a) Schematics of the experiment and
unidirectional excitation of spin waves with positive and negative
wavevectors. Due to the linear dispersion, the group velocity direction is not
changed upon the reversal of wavevector direction. (b) Schematics of the
studied system. The two CoFeB layers are coupled antiferromagnetically through
RKKY interaction. (c) STEM image of the lamella prepared under the antenna
layer (right panel), and atomic fraction of individual elements obtained by
EDX analysis.

In-plane (IP) and out-of-plane (OOP) hysteresis loops of the sample were
acquired by vibrating sample magnetometry (VSM), see Fig. 2a and supplementary
Fig. SUnidirectional propagation of zero-momentum magnons. To obtain static
magnetic parameters of the SAF layer, we model the hysteresis loop using
minimization of the energy in a similar manner as it is done in Stoner-
Wohlfarth model [22, 23]. The total energy reads as [24, 25]

$E=E_{\mathrm{dip}}+E_{\mathrm{Zeeman}}+E_{\mathrm{ex}}+E_{\mathrm{anis}}+E_{\mathrm{J}}+E_{\mathrm{inter}},$
(1)

where $E_{\mathrm{dip}}$ is dipolar energy, $E_{\mathrm{Zeeman}}$ is energy
caused by external field, $E_{\mathrm{ex}}$ is intralayer exchange,
$E_{\mathrm{anis}}$ is energy of uniaxial anisotropy, $E_{\mathrm{inter}}$ is
dipolar energy between the two layers, and interlayer exchange
$E_{\mathrm{J}}=J_{\mathrm{bl}}(1-m_{1}\cdot
m_{2})+J_{\mathrm{bq}}[1-(m_{1}\cdot m_{2})]^{2}$ , where $J_{\mathrm{bl}}$
and $J_{\mathrm{bq}}$ are bilinear and biquadratic exchange constants,
respectively [26, 27]. Using this model, we obtained saturation magnetization
($M_{\mathrm{s}}$) and interlayer coupling constants from the simultaneous
fitting of the experimental IP and OOP hysteresis loops. The fitting results
are summarized in Table 1.

We also implemented bilinear and biquadratic coupling into MuMax3 code [28]
using custom fields. For this, we have to derive the effective field caused by
the coupling and switch off the standard exchange interaction between the
layers. The discretized effective coupling field $B_{\mathrm{J}}$ then reads
as

$B_{\mathrm{J}}=\frac{\partial
E_{\mathrm{J}}}{\partial\bm{M}}=\frac{\sum_{i}\bm{m}_{i}-\bm{m}}{c_{z}M_{\mathrm{s}}}\left(J_{\mathrm{bl}}+2J_{\mathrm{bq}}\sum_{i}\bm{m_{i}}\cdot\bm{m}\right),$
(2)

where $c_{z}$ is cell size in the OOP direction. The results (see Fig. 2a,
circles) are in perfect agreement with the experiment and results from Eq. 1,
which proves that the implementation of the coupling constants in Mumax3 is
valid, and the code can also be used for simulations of dispersion relations.

Figure 2:  Static characterization of SAF sample. (a) In-plane VSM hysteresis
curve, analytic modeling based on energy minimalization, and MuMax3 simulation
results. (b) Hysteresis loops measured by wide-field Kerr microscope. The
external field was applied perpendicular to uniaxial anisotropy axis. (c) VNA-
SAFMR measurement with external field parallel and perpendicular to the
uniaxial anisotropy axis and relevant analytical modeling.

To get further insight into the sample static behavior, we performed wide-
field Kerr microscopy with the external magnetic field applied at different
angles to the easy axis of the SAF. The magnetization components perpendicular
and parallel to the applied field direction are shown in supplementary Fig.
S6, together with corresponding loops obtained from the energy minimalization.
When the external magnetic field is applied perpendicularly to the easy axis,
it is possible to stabilize the single domain state, with the Neél vector
aligned with the easy axis, see Fig. 2b. This property of the studied system
is essential for later BLS experiments which are performed in zero external
magnetic field.

Angle-resolved Kerr microscopy experiments revealed that the switching is
deterministic and depends on the magnetic field sweep direction (see
supplementary Fig. S7). This deterministic switching is unexpected in
perfectly symmetric SAFs. Here, it is probably caused by slightly different
magnetic properties of the individual layers or interfaces [29], and thus
slightly different total energies for the two opposite directions of
magnetization. It could be also caused by different layer thicknesses; however
our STEM data (Fig. 1c) do not show any difference above the measurement
uncertainty of $\pm 1\,$nm.

Table 1: Fitted and fixed material parameters. $\mathrm{t}_{\mathrm{CoFeB}}$
is thickness of the individual CoFeB layers, $\mathrm{t}_{\mathrm{Ru}}$ is
thickness of the Ru spacer, $M_{\mathrm{s}}$ is saturation magnetization,
$J_{\mathrm{bl}}$ ($J_{\mathrm{bq}}$)is bilinear (biquadratic) interlayer
coupling, $H_{\mathrm{u}}$ is uniaxial anisotropy field, $\gamma$ is
gyromagnetic ratio, $A_{\mathrm{ex}}$ is exchange constant, $\alpha$ is
dimensionless damping constant.

Parameter |  $t_{\mathrm{CoFeB}}$ (nm) |  $t_{\mathrm{Ru}}$ (nm) |  $M_{\mathrm{s}}$ (kA/m) |  $J_{\mathrm{bl}}$ $(\mathrm{mJ}/\mathrm{m}^{2})$ |  $J_{\mathrm{bq}}$ $(\mathrm{mJ}/\mathrm{m}^{2})$ |  $H_{\mathrm{u}}$ (kA/m) |  $\gamma$ (GHz/T) |  $A_{\mathrm{ex}}$ (pJ/m) | $\alpha$
---|---|---|---|---|---|---|---|---|---
Value  | $15\pm 1^{*}$ | $0.6\pm 0.2^{*}$ | $1299\pm 14^{\dagger}$ | $-0.67\pm 0.02^{\dagger}$ | $-0.28\pm 0.02^{\dagger}$ | $1.5\pm 0.2^{\ddagger}$ | $30^{\mathsection}$ | $15^{\mathsection}$ | $0.004^{\mathsection}$

∗ Measured from TEM images.
---
† Obtained from the simultaneous fit of OOP and IP hysteresis loops.
‡ Estimated from the offset of the acoustic mode synthetic antiferromagnetic
resonance.
§ These parameters were fixed [17].

Flip-chip vector network analyzer (VNA) synthetic antiferromagnetic resonance
(SAFMR) [30] was performed in two geometries: with the external magnetic field
perpendicular and parallel to the easy axis of the SAF, see supplementary Fig.
S8. The positions of SAFMR, determined by fitting the Lorentz function for
each applied external magnetic field, are shown in Fig. 2c. The experimental
data agree well with the analytic model of SAFMR (shown in Fig. 2c as black
and red lines) developed by Gallardo et al. [24, 25]. This model is based on
solving the following eigenproblem (assuming $k=0$):

$\hat{A}\bm{m}=i\left(\frac{\omega}{\mu_{0}\gamma}\bm{m}\right)$ (3)

where $\hat{A}$ is a system matrix, $\omega$ is a spin-wave frequency,
$\gamma$ is a gyromagnetic ratio, and $\bm{m}$ is dynamic normalized
magnetization vector in both layers. The system matrix depends on the sample
geometry (CoFeB and Ru thicknesses), magnetization orientation of the
individual CoFeB layers [obtained by energy minimization (see Eq. 1)], and on
the magnetic properties (saturation magnetization $M_{\mathrm{s}}$, exchange
constant $A_{\mathrm{ex}}$, uniaxial anisotropy field $H_{\mathrm{u}}$ and on
interlayer exchange coupling $J_{\mathrm{bl}}$ and $J_{\mathrm{bq}}$). To
calculate results shown in Figs. 2c, 3a, and 4c, we used parameters from Table
1.

Figure 3:  Unidirectional excitation of the spin waves in SAF. (a) Dispersion
relation obtained by analytic model and MuMax3 simulation. The red curve shows
the excitation efficiency of the used antenna. (b) BLS intensity measured
approximately $1.5\,\mu$m away from the antenna in the positive (red symbols)
and the negative (black symbols) $z$-direction. (c) 2-dimensional mapping of
the spin-wave intensity at frequency of the $2.5\,$GHz. (d) Extracted spin-
wave phase from five BLS measurements at $2.5\,$GHz.

To obtain SAF dispersion, we can use the same analytic model (Eq. 3), which
can be extended for
$k\in\left[-20\,\mathrm{rad}/\mu\mathrm{m},20\,\mathrm{rad}/\mu\mathrm{m}\right]$
(gray line in Fig.‘3a). To further confirm the results of the analytic model
we performed micromagnetic simulations in MuMax3 (see colormap in Fig.‘3a).
The simulation cell size was set to ($10\times 10$) nm2 in the IP direction
and 15 nm in the OOP direction, i.e., we neglected the Ru spacer and assumed
only one cell per one CoFeB layer. The simulation size was $2048\times
2048\times 2$ cells. The dispersions were calculated using the approach
described in [31].

Figure 4:  Phase velocity reversal upon $\Gamma$ point crossing. (a,b) Spin-
wave intensity (blue squares) and spin-wave phase (black squares) at 1.8 GHz
(a) and 2.6 GHz (b). The red lines are fits of the linear evolution of phase
and exponential decay of spin waves, respectively. (c) Dispersion relation and
phase velocity obtained experimentally (black and red squares) and by analytic
modeling (black and red solid lines) Group velocity measured by time-resolved
BLS (red circles) and calculated from the analytical model (red dashed line).
(d) Simulated space-time map of the zero-momentum spin-wave evolution. Black
dashed line shows the group velocity. The appearance of the spin wave in the
left part of the space-time map is caused by periodic boundary conditions used
in the simulation.

Before BLS experiments, we stabilized the single-domain magnetization state in
wide-field Kerr microscope. The Neél vector pointed perpendicular to the
stripline antenna (connected to RF generator), which was used to excite spin
waves in the SAF layer. As the first BLS experiment, we swept the excitation
frequency and compared the BLS spectra acquired at a distance of $1.5\,\mu m$
from the antenna on both sides, see Fig. 3b. Higher excitation efficiency on
the positive side of the antenna is clearly visible. On the negative side, the
spin wave signal above 2 GHz is completely suppressed. At this region, the
antenna cannot efficiently excite spin waves due to the small excitation
efficiency of the antenna in the range of available wavenumbers, see Fig. 3a,
red line. The excitation efficiency was calculated using fast Fourier
transform of the field calculated by FEMM solver [32]. The signal below 2 GHz
visible on the negative side of the antenna originates from the dispersion
region with negative group velocity ($k<-5\,\mathrm{rad}/\mu\mathrm{m}$).
Note, that this signal can be entirely suppressed by modifying the excitation
efficiency, i.e. by using a slightly wider antenna, thus the spin wave
emission can be completely unidirectional.

To further investigate this non-reciprocal behavior, we performed a
2-dimensional mapping of BLS signal at the frequency of 2.5 GHz, shown in Fig.
3c. At this frequency, we can observe unidirectional excitation of uniform
spin waves. We have performed full phase-resolved mapping of spin waves with
use of five different measurements [33]. The results (Fig. 3d) show clear
spatially uniform evolution of the phase in the direction perpendicular to the
antenna in the positive $z$-direction. On the other side of the antenna, no
phase evolution was observed.

As the next experiment, we performed two full phase-resolved linescans in the
direction perpendicular to the antenna at the frequencies of 1.8 GHz and 2.6
GHz, and extracted the intensity and phase of spin waves (Fig. 4a, b). For 1.8
GHz, we can observe a negative phase evolution in the positive values of $z$
(upper side of the antenna in Fig. 3c,d). On the other hand, at the frequency
of 2.6 GHz, we can observe a positive phase evolution on the same side of the
antenna. The intensity data show prominent non-reciprocity, where exponential
decay is visible only towards the positive direction. By repeating this
experiment at different frequencies we obtained full dispersion relation (gray
squares) and frequency-dependent phase velocity (red squares), see Fig. 4c. We
also performed time-resolved BLS experiments and extracted group velocities
(see red circles in Fig. 4c, and supplementary Fig. S9). The experimental data
agree well with our analytical model (see solid and dashed lines in Fig. 4c).

Spin waves at the $\Gamma$ point ($k=0$) are shifted up by approx. 1.9 GHz,
which agrees well with the data in Fig. 2c, and is caused by the uniaxial
anisotropy. The anisotropy and the dipolar interaction induce the monotonous
and smooth dispersion through the $\Gamma$ point and leads to unidirectional
group velocity (see dashed line in Fig. 4c). The group velocity is conserved
despite the change of the wavenumber (phase velocity) from positive to
negative.

This property leads to interesting phenomena, where the wave with zero
momentum ($k=0$), can still propagate in space. To demonstrate this, we
performed micromagnetic simulation where the dynamics was excited by OOP sinc
pulse centered at $z=0$ $\mu$m. The evolution of the whole excited spin wave
packet in space and time is shown in supplementary Fig. S10. If we filter only
wavevectors with $k\to 0$ we can clearly see the propagation of this ”wave”
towards the positive value of $z$ (see Fig. 4d). The phase velocity
$v_{\mathrm{p}}\to\infty$, while the group velocity
$v_{\mathrm{g}}=2.3\,\mu\mathrm{m}/\mathrm{ns}$ is given by the dispersion
relation (see Fig. 4c). The orientation of the unidirectional propagation can
be reversed by switching the SAF magnetization (see supplementary Fig. S10).

In conclusion, we observed unidirectional propagation of zero-momentum magnons
in synthetic antiferromagnet. Inherent non-reciprocity of spin waves in
synthetic antiferromagnets with uniaxial anisotropy shapes the dispersion
relation into a smooth and monotonous function in the vicinity of the $\Gamma$
point. This situation leads to unidirectional propagation of spin waves with
negative, zero and positive wavenumbers. This phenomenon was observed by
intensity-, phase-, and time-resolved BLS microscopy and further corroborated
by micromagnetic simulations and analytical modelling. The unusual dispersion
relation offers many possibilities for further investigations of e.g.,
propagation of wavepackets with bidirectional phase velocities, or spin wave
accumulation on domain boundaries. In practical applications such
unidirectional system can pave the way towards designing new non-reciprocal
elements for signal processing and unconventional spin-wave computing.

###### Acknowledgements.

The authors thank to Roman Verba for discussion on analytical models for SAFs.
This research was supported by the Horizon Europe EIC Pathfinder program,
project ID 101098651 (METASPIN). CzechNanoLab project LM2023051 is
acknowledged for the financial support of the measurements and sample
fabrication at CEITEC Nano Research Infrastructure. Specific research projects
IDs CEITEC VUT-J-23-8435 and CEITEC VUT-J-23-8411 are also acknowledged. O.W.
was supported by Brno PhD talent scholarship. A.V.C. acknowledges support by
the Austrian Science Fund (FWF) under Grant No. I 4917-N (MagFunc).

## References

  * Ishibashi  _et al._ [2020] M. Ishibashi, Y. Shiota, T. Li, S. Funada, T. Moriyama, and T. Ono, Switchable giant nonreciprocal frequency shift of propagating spin waves in synthetic antiferromagnets, Sci. Adv. 6, eaaz6931 (2020).
  * Shiota  _et al._ [2020] Y. Shiota, T. Taniguchi, M. Ishibashi, T. Moriyama, and T. Ono, Tunable magnon-magnon coupling mediated by dynamic dipolar interaction in synthetic antiferromagnets, Phys. Rev. Lett. 125, 017203 (2020).
  * Mouhoub  _et al._ [2023] A. Mouhoub, F. Millo, C. Chappert, J.-V. Kim, J. Létang, A. Solignac, and T. Devolder, Exchange energies in CoFeB/Ru/CoFeB synthetic antiferromagnets, Phys. Rev. Mater. 7, 044404 (2023).
  * Chumak  _et al._ [2022] A. V. Chumak _et al._ , Advances in magnetics roadmap on spin-wave computing, IEEE Trans. Magn. 58, 1 (2022).
  * Pirro  _et al._ [2021] P. Pirro, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Advances in coherent magnonics, Nat Rev. Mater. 6, 1114 (2021).
  * Gerevenkov  _et al._ [2023] P. I. Gerevenkov, V. D. Bessonov, V. S. Teplov, A. V. Telegin, A. M. Kalashnikova, and N. E. Khokhlov, Nonreciprocal collective magnetostatic wave modes in geometrically asymmetric bilayer structure with nonmagnetic spacer, Nanoscale 15, 6785 (2023).
  * Wintz  _et al._ [2016] S. Wintz, V. Tiberkevich, M. Weigand, J. Raabe, J. Lindner, A. Erbe, A. Slavin, and J. Fassbender, Magnetic vortex cores as tunable spin-wave emitters, Nat Nanotechnol. 11, 948 (2016).
  * Grünberg [1981] P. Grünberg, Magnetostatic spin-wave modes of a heterogeneous ferromagnetic double layer, J. Appl. Phys. 52, 6824 (1981).
  * Grünberg  _et al._ [1986] P. Grünberg, R. Schreiber, Y. Pang, M. B. Brodsky, and H. Sowers, Layered magnetic structures: Evidence for antiferromagnetic coupling of Fe layers across Cr interlayers, Phys. Rev. Lett. 57, 2442 (1986).
  * Cochran  _et al._ [1990] J. F. Cochran, J. Rudd, W. B. Muir, B. Heinrich, and Z. Celinski, Brillouin light-scattering experiments on exchange-coupled ultrathin bilayers of iron separated by epitaxial copper (001), Phys. Rev. B 42, 508 (1990).
  * Vavassori  _et al._ [2000] P. Vavassori, M. Grimsditch, E. Fullerton, L. Giovannini, R. Zivieri, and F. Nizzoli, Brillouin light scattering study of an exchange coupled asymmetric trilayer of Fe/Cr, Surf. Sci. 454, 880 (2000).
  * Millo  _et al._ [2023] F. Millo, J.-P. Adam, C. Chappert, J.-V. Kim, A. Mouhoub, A. Solignac, and T. Devolder, Unidirectionality of spin waves in synthetic antiferromagnets, arXiv preprint 2306.05259 (2023).
  * Chen  _et al._ [2019] J. Chen, T. Yu, C. Liu, T. Liu, M. Madami, K. Shen, J. Zhang, S. Tu, M. S. Alam, K. Xia, M. Wu, G. Gubbiotti, Y. M. Blanter, G. E. W. Bauer, and H. Yu, Excitation of unidirectional exchange spin waves by a nanoscale magnetic grating, Phys. Rev. B 100, 104427 (2019).
  * Gallardo  _et al._ [2021] R. A. Gallardo, P. Alvarado-Seguel, A. Kákay, J. Lindner, and P. Landeros, Spin-wave focusing induced by dipole-dipole interaction in synthetic antiferromagnets, Phys. Rev. B 104, 174417 (2021).
  * Zakeri  _et al._ [2010] K. Zakeri, Y. Zhang, J. Prokop, T.-H. Chuang, N. Sakr, W. X. Tang, and J. Kirschner, Asymmetric spin-wave dispersion on Fe(110): Direct evidence of the Dzyaloshinskii-Moriya interaction, Phys. Rev. Lett. 104, 137203 (2010).
  * Di  _et al._ [2015] K. Di, V. L. Zhang, H. S. Lim, S. C. Ng, M. H. Kuok, J. Yu, J. Yoon, X. Qiu, and H. Yang, Direct observation of the Dzyaloshinskii-Moriya interaction in a Pt/Co/Ni film, Phys. Rev. Lett. 114, 047201 (2015).
  * Vaňatka  _et al._ [2021] M. Vaňatka, K. Szulc, O. c. v. Wojewoda, C. Dubs, A. V. Chumak, M. Krawczyk, O. V. Dobrovolskiy, J. W. Kłos, and M. Urbánek, Spin-wave dispersion measurement by variable-gap propagating spin-wave spectroscopy, Phys. Rev. Appl. 16, 054033 (2021).
  * Gladii  _et al._ [2016] O. Gladii, M. Haidar, Y. Henry, M. Kostylev, and M. Bailleul, Frequency nonreciprocity of surface spin wave in permalloy thin films, Phys. Rev. B 93, 054430 (2016).
  * Vlaminck and Bailleul [2008] V. Vlaminck and M. Bailleul, Current-induced spin-wave Doppler shift, Science 322, 410 (2008).
  * Kim  _et al._ [2021] D.-H. Kim, S.-H. Oh, D.-K. Lee, S. K. Kim, and K.-J. Lee, Current-induced spin-wave doppler shift and attenuation in compensated ferrimagnets, Phys. Rev. B 103, 014433 (2021).
  * Temdie  _et al._ [2023] L. Temdie, V. Castel, C. Dubs, G. Pradhan, J. Solano, H. Majjad, R. Bernard, Y. Henry, M. Bailleul, and V. Vlaminck, High wave vector non-reciprocal spin wave beams, AIP Adv. 13, 025207 (2023).
  * Stoner and Wohlfarth [1948] E. C. Stoner and E. P. Wohlfarth, A mechanism of magnetic hysteresis in heterogeneous alloys, Philos. Trans. Phys. Sci. Eng. 240, 599 (1948).
  * Skomski [2008] R. Skomski,  _Simple models of magnetism_ (Oxford university press, 2008).
  * Gallardo  _et al._ [2019] R. Gallardo, T. Schneider, A. Chaurasiya, A. Oelschlägel, S. Arekapudi, A. Roldán-Molina, R. Hübner, K. Lenz, A. Barman, J. Fassbender, J. Lindner, O. Hellwig, and P. Landeros, Reconfigurable spin-wave nonreciprocity induced by dipolar interaction in a coupled ferromagnetic bilayer, Phys. Rev. Appl. 12, 034012 (2019).
  * Wojewoda [2023] O. Wojewoda, Spin-wave toolkit, GitHub (2023).
  * Slonczewski [1991] J. C. Slonczewski, Fluctuation mechanism for biquadratic exchange coupling in magnetic multilayers, Phys. Rev. Lett. 67, 3172 (1991).
  * Guslienko [1994] K. Y. Guslienko, Bilinear and biquadratic interlayer exchange coupling in ferromagnet-nonmagnetic metal multilayer film, IEEE Trans. Magn. 30, 757 (1994).
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, The design and verification of Mumax3, AIP Adv. 4, 107133 (2014).
  * Jang  _et al._ [2010] S. Y. Jang, S. H. Lim, and S. R. Lee, Magnetic dead layer in amorphous CoFeB layers with various top and bottom structures, J. Appl. Phys. 107, 09C707 (2010).
  * Adams  _et al._ [2018] D. J. Adams, D. Cimpoesu, M. A. Khan, P. Poudyal, and L. Spinu, Critical switching curves of FeCoB synthetic antiferromagnets, J. Phys. D Appl. Phys. 51, 055005 (2018).
  * Wojewoda  _et al._ [2023a] O. Wojewoda, F. Ligmajer, M. Hrtoň, J. Klíma, M. Dhankhar, K. Davídková, M. Staňo, J. Holobrádek, J. Krčma, J. Zlámal, _et al._ , Observing high-k magnons with Mie-resonance-enhanced Brillouin light scattering, Commun. Phys. 6, 1 (2023a).
  * Meeker [2023] D. Meeker, Finite element method magnetics, www.femm.info (2023), accessed: 2023-10-16.
  * Wojewoda  _et al._ [2023b] O. Wojewoda, M. Hrtoň, M. Dhankhar, J. Krčma, K. Davídková, J. Klíma, J. Holobrádek, F. Ligmajer, T. Šikola, and M. Urbánek, Phase-resolved optical characterization of nanoscale spin waves, Appl. Phys. Lett. 122, 202405 (2023b).

Supplemental Material: Unidirectional propagation of zero-momentum magnons

FIG. S 5:  The out-of-plane hysteresis loop measured with VSM. FIG. S 6:
Hysteresis loops acquired by wide-field Kerr microscope where the angle
between the external field and the uniaxial anisotropy was 0 deg (a), 30 deg
(b), 60 deg (c), and 90 deg (d). FIG. S 7:  The remanent magnetization
perpendicularly to the direction of applied field for different angles and
sweep directions. FIG. S 8:  The calculated transmitted power from S12
parameter. The uniaxial anisotropy is perpendicular (a), or parallel (b) to
external field. FIG. S 9:  (a,b) The time-resolved BLS signal for different
distances from antenna at 2 GHz (a) and 3 GHz (b). The extracted rise times
from logistic curve fit at 2 GHz (c) and 3 GHz (d). The data are subsequently
fitted with linear curve. FIG. S 10:  (a, b) The spatio-temporal map of the
sinc pulse for two opposite configuration of SAF. We can observe that the
zero-momentum ($k=0$) spin waves propagate to the right (a) or left (b) with
the group velocity of $\left|2.3\right|\,\mu\mathrm{m}/\mathrm{ns}$ (black
dashed lines). The waves with lower (higher) group velocities and negative
(positive) phase velocities, are above (below) the dashed lines. The change in
the phase velocity direction (the sign of wavenumber), exhibits itself in this
representation as the change of the slope of wavefronts. (c) For comparison we
present also simulation of 30 nm thick ferromagnetic slab with the same
material parameters, see Table 1.