# Physics-Informed Convolutional Autoencoder for Cyber Anomaly Detection in
Power Distribution Grids

Mehdi Jabbari Zideh, , Sarika Khushalani Solanki, Authors are with the Lane
Department of Computer Science and Electrical Engineering, West Virginia
University, Morgantown, WV 26505 USA (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The growing trend toward the modernization of power distribution systems has
facilitated the installation of advanced measurement units and promotion of
the cyber communication systems. However, these infrastructures are still
prone to stealth cyber attacks. The existing data-driven anomaly detection
methods suffer from a lack of knowledge about the system’s physics, lack of
interpretability, and scalability issues hindering their practical
applications in real-world scenarios. To address these concerns, physics-
informed neural networks (PINNs) were introduced. This paper proposes a
multivariate physics-informed convolutional autoencoder (PIConvAE) to detect
stealthy cyber-attacks in power distribution grids. The proposed model
integrates the physical principles into the loss function of the neural
network by applying Kirchhoff’s law. Simulations are performed on the modified
IEEE 13-bus and 123-bus systems using OpenDSS software to validate the
efficacy of the proposed model for stealth attacks. The numerical results
prove the superior performance of the proposed PIConvAE in three aspects: a)
it provides more accurate results compared to the data-driven ConvAE model, b)
it requires less training time to converge c) the model excels in effectively
detecting a wide range of attack magnitudes making it powerful in detecting
stealth attacks.

###### Index Terms:

Physics-informed neural networks, cyber anomaly detection, power distribution
grids, unsupervised data-driven methods, convolutional autoencoder, false data
injection attacks.

##  I Introduction

The recent advancement in smart grid monitoring technology has introduced the
extensive installation of phasor measurement units (PMUs) and $\mu$PMUs to
increase the situational awareness of power system operators [1, 2]. This
improvement in the system observability has exposed smart grids to various
types of cyber-attacks threatening the grids’ cyber-physical infrastructures
[3]. As the behavior and patterns of the cyber anomalies are unknown and
unpredictable, their detection is a challenging task to ensure secure and
reliable operation of the system [4].

Machine learning methods as powerful techniques are applied to various aspects
of power systems including forecasting [5], power market operation [6, 7],
power flow calculations [8], and cyber-physical security analysis [9]. More
specifically, in the task of anomaly detection, several research
investigations have applied deep neural networks (DNNs) to distinguish
abnormal data from the normal data patterns within the sensor data. Authors in
[10] applied discrete wavelet transform (DWT) and a recurrent neural network
(RNN) to detect cyber attacks. A similar methodology was proposed in [11] to
detect and classify high-impedance faults in power grids. The work in [12]
proposed a generative adversarial network (GAN) model for event detection in
$\mu$PMU measurements. A multivariate long short-term memory (LSTM)
autoencoder-based model was employed in [13] for cyber anomaly detection in
power distribution systems.

Despite their impressive performance, ML models mainly suffer from poor
extrapolation capabilities, scalability constraints, and lack of understanding
of the system’s physical laws. Consequently, this reduces their reliability as
the system operators often struggle to interpret their implications. To
address the aforementioned challenges and enhance the performance of the
learning algorithms with lower computational effort, physics-informed neural
networks (PINNs) were introduced integrating the fundamental physical
principles of the systems into NN models [14]. The works in [15] and [16]
extensively reviewed the applications of PINN methods in the cyber-power
domain and power systems, respectively. Recently a growing number of research
studies have applied PINN methods for anomaly detection in power systems.
Authors in [17] applied the physics-based inequality constraints to the
isolation forest algorithm to detect falsification attacks on wind generation
signals. The work in [18] proposed a combined autoencoder with a graph neural
network (GNN) model which learns the spatial structure of $\mu$PMU data
through a weighted adjacency matrix to detect events in power distribution
grids. A new framework for combining model-based and NN methods is proposed in
[19] where the prediction residuals of variational autoencoder (VAE), LSTM,
and state estimation are fed to a DNN model for attack detection. A physics-
informed method is employed in [20] for attack detection utilizing the
information of transmission system topology structure and node features in
graph edge-conditioned convolutional networks (GECCN). In [21], the physical
relationships are integrated into the NN loss function through dynamical swing
equations for the detection of load-altering attacks.

In this paper, we propose a novel PINN model for cyber-anomaly detection in
power distribution systems. More specifically, the contributions of this work
are threefold: 1) We present a physics-informed convolutional autoencoder
(PIConvAE) model by encoding the physics-based equations of $\mu$PMU
measurements into the loss function of the model to detect cyber attacks; 2)
The proposed model works in an unsupervised manner and does not require any
prior labeling of the cyber-attacks; and 3) The efficacy of the proposed
PIConvAE model is evaluated for cyber-attack detection on IEEE 13-bus and
123-bus systems.

The remainder of this paper is organized as follows. Section II introduces
data falsification attack mechanisms and data falsification functions. Section
III covers the principles and structure of the PINN model. Section IV
discusses the simulation results of the proposed PIConvAE model under
falsification attacks. Finally, Section V concludes the paper.

##  II Development of Cyber-attacks on Power Distribution Measurements

The real-time data measured by $\mu$PMU and supervisory control and data
acquisition (SCADA) systems including system frequency, voltage and current
magnitudes and angles, are essential for power system control and stability
analysis. The data transferred by the communication channels is prone to
falsified attacks by the attackers to tackle data confidentiality, system
reliability, and integrity. Assume that the measured data by the measurement
units at time slot $t$ is $Z^{i}_{act}(t)$. The reported measurement
$Z^{i}_{rep}(t)$ at time slot $t$ is equal to the measured data when the data
is not falsified. However, if the data from measurement unit $i$ is
compromised by an adversary, $Z^{i}_{act}(t)\neq Z^{i}_{rep}(t)$. We assume
that the adversary can compromise the voltage and current magnitudes and
angles. We apply the attack strategy employed in [22] where the falsified data
is defined as

$Z^{i}_{comp}(t)=(1+\alpha(t))Z^{i}_{act}(t)$ (1)

where $\alpha(t)$ ($\alpha_{min}\leq\alpha(t)\leq\alpha_{max}$) describes the
attack magnitude at time slot $t$, $\alpha(t)>0$ for additive attacks and
$\alpha(t)<0$ for deductive attacks. For combined attacks, the compromised
data consists of additive and deductive attacks in equal proportions. The
objective of this type of attack is to report the falsified data as normal
even when the data is above or below a certain threshold. Therefore, the
control systems are deceived to not function properly. In this work, various
values for $\alpha$ ($-5\%\leq\alpha(t)\leq 5\%$) are chosen to not only
demonstrate the variability of attack intensities but also to make them less
readily detectable.

##  III Physics-informed Neural Network for FDIA Detection

In this section, we propose the physics-informed data-driven method by
integrating the physics-based relationships of $\mu$PMU measurements into the
ConvAE model to achieve more accurate and robust detection of falsified
attacks in power distribution grids. Suppose that we have a system with $N$
buses comprising distributed energy resources (DERs) such as PV and wind
generation units, and power loads. Assuming the local measurements of the
buses by measurement units, the following quantities are measured for each bus
$i$

$x_{i}=[V_{i},I_{i},\theta_{i},\delta_{i},P_{i},Q_{i}]$ (2)

where $V_{i}$ and $I_{i}$ are voltage and current magnitudes, $\theta_{i}$ and
$\delta_{i}$ are voltage and current phase angles, and $P_{i}$ and $Q_{i}$ are
active and reactive power measurements.

###  III-A Autoencoder-based Attack Detector

Attack detection is a challenging task due to the unpredictable nature and
patterns of attacks and their evolution over time. This motivates the
applications of unsupervised learning methods, which can distinguish abnormal
data from normal data samples without relying on data labels. Autoencoders
(AE) as a category of unsupervised neural networks are trained on only normal
training data and identify data that deviate from the distribution of the
normal samples as anomalies.

Autoencoders learn the hidden patterns of the training data by encoding them
to a latent space (bottleneck) and reconstructing the embedded data using a
decoder. If the input data $x$ is passed through an AE model, the
reconstructed output $\hat{x}$ is generated by learning the following
representation:

${\hat{x}(W,b,x)}=\sigma(W^{(n)}\sigma(W^{(n-1)}\sigma{(W^{(n-2)}...})+b^{(n-1)})+b^{(n)})$
(3)

where $W$ and $b$ represent the AE weight and bias vectors, $n$ is the number
of hidden layers, and $\sigma$ represents the non-linear activation function.
By learning this representation, the primary objective of the training
procedure is to generate the reconstructed data as similar as possible to the
normal training data by minimizing the following loss function

$\mathcal{L_{AE}}\left(W,b,x\right)=\frac{1}{m}\sum_{k=1}^{m}\left|\hat{x}(W,b,x^{k})-x^{k}\right|^{2}+\lambda\|W\|^{2}$
(4)

where $m$ is the training data size and $\lambda$ is a hyper-parameter
preventing the weights from overgrowing. The second term is the regularization
to decrease the magnitudes of the weights.

To extract the hidden patterns of the time series data, we use a convolutional
neural network (CNN) in the structure of the encoder and decoder. The CNN
model consists of a non-linear activation function and a convolution operator
that implements a dot product between the input data and a filter described
as:

$y^{l+1}=\sigma(W^{l}*x^{l}+b^{l})$ (5)

where $y^{l+1}$ is the output of the ($l+1$)-th layer, $W^{l}$ and $b^{l}$
represent weights and biases of the $l$-th layer and $*$ is the convolution
operator.

Figure 1: Overview of the training and detection of the developed PIConvAE
model.

###  III-B Physics-informed Convolutional Autoencoder (PIConvAE)

The physics-based principles can be integrated into NN models by various
strategies [15, 16]. In this work, an additional loss term is added to the
ConvAE loss function to constrain the reconstructed outputs to satisfy the
physics-driven laws. For a series of measurements from a measurement unit
(vector $x$ in (2)), the reconstructed outputs should satisfy Kirchhoff’s law
at each node at each time slot $t$. This ensures the physics-based
relationships between measurements for power flow calculations. Following this
principle, the physics-driven loss functions for the reconstructed
measurements at bus $i$ are defined as

$\mathcal{L}_{Phy_{P}}\left(\hat{x}\right)=\frac{1}{m}\sum_{k=1}^{m}\left|\hat{P}_{i}^{k}-\hat{V}_{i}^{k}\hat{I}_{i}^{k}\cos{(\hat{\theta}_{i}^{k}-\hat{\delta}_{i}^{k})}\right|^{2}$
(6)
$\mathcal{L}_{Phy_{Q}}\left(\hat{x}\right)=\frac{1}{m}\sum_{k=1}^{m}\left|\hat{Q}_{i}^{k}-\hat{V}_{i}^{k}\hat{I}_{i}^{k}\sin{(\hat{\theta}_{i}^{k}-\hat{\delta}_{i}^{k})}\right|^{2}$
(7)

where $\mathcal{L}_{Phy_{P}}$ and $\mathcal{L}_{Phy_{Q}}$ represent the
deviations of the reconstructed outputs from active and reactive powers in the
Kirchhoff’s law, respectively. The total function for training the PIConvAE
model is created by combining (4), (6), and (7) as follows

$\mathcal{L}_{T}\left(W,b,x,\hat{x}\right)=\alpha_{d}\mathcal{L_{AE}}+\alpha_{phy}(\mathcal{L}_{Phy_{P}}+\mathcal{L}_{Phy_{Q}})$
(8)

where $\alpha_{d}$ and $\alpha_{phy}$ are data-driven and physics-based loss
coefficients to make a balance between the loss terms. Fig. 1 shows the data
generation, training, and detection framework of the proposed PIConvAE model.
The pseudo-code for the PIConvAE model is shown in Algorithm 1. In this study,
we set $\alpha_{d}=\alpha_{phy}=1$ as it produces better results.

Input Data : Training, Validation, and Test Data

Required : $Ep$, Number of epochs

$b$, batch size

$N_{w}$, window size

$\alpha$, learning rate

Initialization : Weights $W$, $b$

Training and Validation:

for  _each $Ep$_ do

Sample $(x_{i}^{1...N_{w}})_{i=1}^{b}$ from training data

Sample $(v_{i}^{1...N_{w}})_{i=1}^{b}$ from validation data

Encoder:

Fix weights $W_{\mathcal{D}}$

$g_{w_{\mathcal{E}}}=\nabla_{w_{\mathcal{E}}}[\alpha_{d}\mathcal{L_{AE}}\left(W,b,x\right)+\alpha_{phy}(\mathcal{L}_{Phy_{P}}(\hat{x})+\mathcal{L}_{Phy_{Q}}(\hat{x}))]$

Decoder:

Fix weights $W_{\mathcal{E}}$

$g_{w_{\mathcal{D}}}=\nabla_{w_{\mathcal{D}}}[\alpha_{d}\mathcal{L_{AE}}\left(W,b,x\right)+\alpha_{phy}(\mathcal{L}_{Phy_{P}}(\hat{x})+\mathcal{L}_{Phy_{Q}}(\hat{x}))]$

Compute training and validation losses

end for

Testing:

$X_{test}$=$(x_{n}^{1...N_{w}})_{n=1}^{K}$

for  _$n=1,...,K$_ do

$\hat{x}^{n}=\mathcal{D}{(\mathcal{E}(x^{n}))}$

for  _each $l$_ do

${a}^{n}\left(x^{n},\hat{x}^{n}\right)=|x^{n}-\hat{x}^{n}|+{a}_{phy}^{n}$

end for

end for

Algorithm 1 PIConvAE Training, Validation, and Testing

##  IV Experimental Results

The performance of the proposed PIConvAE model is evaluated on the IEEE 13-bus
and 123-bus systems. Training, validation, and test data are generated using
OpenDSS software. To simulate these systems with DER resources, PV and wind
generation units are installed at different nodes to simulate the randomness
and uncertainty of the real-world cases. Solar and wind generation data of San
Diego city including ambient temperature, solar irradiance, and wind speed for
the first week of 2021 are employed [23]. Historical load data of San Diego
Gas and Electric (SDGE) is used to create load data profiles for both systems.
In addition, measurement units are installed at various locations to monitor
the measurements. The monitor and DER placements are summarized in Table I.
The total generated data has 10080 data points where 75% is used for training,
10% for validation, and 15% for testing. Data samples are normalized using the
MinMax function and are divided into subsequences with window size 16 and step
size 1 before feeding to the PIConvAE model. Attacks are implemented on the
current magnitude at bus 671 in IEEE 13-bus and on the voltage magnitude at
bus 35 in IEEE 123-bus systems. All the simulations are implemented with
Tensorflow in Python on a computer with Nvidia RTX3070 8GB and a 2.1 GHz Intel
Core i7 CPU with 32 GB of RAM.

Table I: DGs and Monitors Locations in the Test Systems.

###  IV-A PIConvAE Settings

We employ a two-layer convolutional NN in the encoder design with 64 and 32
filters, with kernel sizes of 5 and 3, respectively. These layers are followed
by a fully connected (FC) layer. The decoder has a reverse structure and an FC
layer after the second convolutional layer. The activation functions are
LeakyReLU and dropout with a rate of 0.2 is applied to the convolutional
layers. Batch normalization is applied to all layers and a decaying learning
rate with an initial rate of 0.001 is used.

###  IV-B Performance Metrics

In the testing stage, both reconstruction and physics-based errors are
employed to find the anomaly score $a$ at each time step $n$ for measurement
$l\in\\{V,I,\theta,\delta,P,Q\\}$ as follows.

$\begin{split}{a}^{l,n}\left(x^{l,n},\hat{x}^{l,n}\right)&=|x^{l,n}-\hat{x}^{l,n}|+{a}_{phy}^{l,n}\end{split}$
(9)

The first term is the point-wise error and ${a}_{phy}^{l,n}$ represents the
physics-based mismatch in the test data. This physics-driven error is defined
for each measurement based on Kirchhoff’s law. For instance, ${a}_{phy}^{l,n}$
for voltage measurements are defined as

$\displaystyle{a}_{phy}^{V,n}$
$\displaystyle=|V^{n}-\frac{P^{n}}{I^{n}\cos(\theta^{n}-\delta^{n})}|^{2}+|V^{n}-\frac{Q^{n}}{I^{n}\sin(\theta^{n}-\delta^{n})}|^{2}\vspace{+1mm}$
(10)

To evaluate the performance of the proposed model under attack scenarios, true
positive (TP), true negative (TN), false positive (FP), and false negative
(FN) are used to find accuracy (ACC=(TP+TN)/(TP+TN+FP+FN), precision
(Prec=TP/(TP+FP)), recall (Rec=TP/(TP+FN)), and F-1 score
(F1=2$\times$(Prec$\times$ Rec)/(Prec+Rec)).

###  IV-C Training and Detection Results

The proposed PIConvAE model offers advantages such as faster training time and
higher detection results. Fig. 2(a) compares the validation loss of the
proposed PIConvAE and ConvAE models for the IEEE 13-bus system. As shown in
this figure, the loss of the PIConvAE model increases after 200 epochs
indicating that the model will be overfit after it while the ConvAE model
still requires additional training. This phenomenon of early stop in the
training of PIConvAE is clearly illustrated in Fig. 2(b) where more training
leads to a reduction in the model performance. In other words, the PIConvAE
model needs 200 epochs to reconstruct the data whereas the ConvAE model is
trained for 1000 epochs. Figures 2(c) and 2(d) show the validation loss and
reconstruction outputs of the PIConvAE model for the IEEE 123-bus system. As
shown in Fig. 2(c), the validation loss of the PIConvAE model reaches a
plateau in around 500 epochs and increases afterward. The performance of the
model will degrade when the number of training epochs increases to 1000 as
shown in Fig. 2(d). Table II also represents the training time of the PIConvAE
and ConvAE models for both systems noting that the integration of physics-
based equations into the loss function significantly reduces the computational
time of PIConvAE compared to the data-driven model.

(a) Validation loss

(b) Current magnitude

(c) Validation loss

(d) Voltage magnitude

Figure 2: Simulation results for IEEE 13-bus (a and b), and IEEE 123-bus (c
and d). Table II: Training time of PIConvAE and ConvAE. Table III: Performance
Comparison for FDIA detection in IEEE 13-bus and 123-bus systems.

In addition to the training time, the PIConvAE model has shown more accurate
results for attack detection. Table III shows the comparison of the detection
performance for the two models. As can be seen from this table, in addition to
higher accuracy rates, the PIConvAE model has perfect precision for both
systems showing the model’s capability to correctly detect all the attacked
data. It also outperforms the ConvAE model with recall and F-1 score at 97.5%
and 98.73%, respectively. Although the ConvAE model’s performance improves for
123-bus systems, it still lags behind the physics-informed model. The PIConvAE
model achieves significantly higher recall and F-1 score compared to the
ConvAE model with a recall of 95% and an F-1 score of 97.44%. The comparison
of the physics-informed model with the data-driven model unequivocally
demonstrates that by incorporating the physical principles into the NN model,
the proposed PIConvAE model achieves higher performances in terms of detection
metrics for both systems. This substantiates the outstanding performance of
the PIConvAE model in comparison to the data-driven ConvAE model.

###  IV-D Evaluation Under Different Attack Magnitudes

In the previous detection phase, we randomly generated attack magnitudes to be
within 5% of the actual measurements. In this section, the sensitivity of the
proposed PIConvAE model against different attack magnitudes is evaluated. Fig.
3(a) shows performance metrics for the IEEE 13-bus system under different
attack magnitudes. As shown in this figure, under very low magnitudes, the
model performance is low, however, as the higher magnitudes of attack can be
more concerning for the system operation, the proposed model presents
remarkable results by the increase in the attack magnitudes. The model
achieves superior performance, particularly under attack magnitudes of 3% and
above. The proposed model also exhibits high detection performance for the
IEEE 123-bus system. As can be seen in Fig. 3(b), the PIConvAE model
demonstrates outstanding results for attack detection when the attack
magnitude increases.

##  V Conclusion

This work proposed a multivariate PIConvAE model to detect cyber attacks in
power distribution grids. The proposed model uses the power flow measurements
to embed the physical principles into the loss function of the ConvAE model
using Kirchhoff’s law. The performance of the proposed PIConvAE model was
evaluated on the modified IEEE 13-bus and 123-bus systems. The encoding of
physic-based relationships enabled the PIConvAE model to achieve higher
detection performance compared to the corresponding data-driven model with
faster training time. The detection results on both systems also showed the
effectiveness of the proposed model in detecting stealth cyber attacks under
different attack magnitudes.

(a)

(b)

Figure 3: Sensitivity of the proposed PIConvAE model to different attack
magnitudes for (a) IEEE 13-bus system, (b) IEEE 123-bus system.

## References

  * [1] A. Shahsavari, M. Farajollahi, E. Stewart, E. Cortez, and H. Mohsenian-Rad, “Situational awareness in distribution grid using micro-pmu data: A machine learning approach,” IEEE Transactions on Smart Grid, vol. 10, no. 6, pp. 6167–6177, 2019.
  * [2] M. Farajollahi, A. Shahsavari, E. Stewart, and H. Mohsenian-Rad, “Locating the source of events in power distribution systems using micro-pmu data,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6343–6354, 2018.
  * [3] R. Liu, H. Mustafa, Z. Nie, and A. K. Srivastava, “Reachability-based false data injection attacks and defence mechanisms for cyberpower system,” Energies, vol. 15, no. 5, p. 1754, 2022.
  * [4] A. A. Habib, M. K. Hasan, A. Alkhayyat, S. Islam, R. Sharma, and L. M. Alkwai, “False data injection attack in smart grid cyber physical system: Issues, challenges, and future direction,” Computers and Electrical Engineering, vol. 107, p. 108638, 2023.
  * [5] R. Nematirad and A. Pahwa, “Solar radiation forecasting using artificial neural networks considering feature selection,” in 2022 IEEE Kansas Power and Energy Conference (KPEC), pp. 1–4, IEEE, 2022.
  * [6] M. Jabbari Zideh and S. S. Mohtavipour, “Two-sided tacit collusion: Another step towards the role of demand-side,” Energies, vol. 10, no. 12, p. 2045, 2017.
  * [7] S. S. Mohtavipour and M. J. Zideh, “An iterative method for detection of the collusive strategy in prisoner’s dilemma game of electricity market,” EEJ Transactions on Electrical and Electronic Engineering, vol. 14, no. 2, pp. 252–260, 2019.
  * [8] X. Hu, H. Hu, S. Verma, and Z.-L. Zhang, “Physics-guided deep neural networks for power flow analysis,” IEEE Transactions on Power Systems, vol. 36, no. 3, pp. 2082–2092, 2020.
  * [9] J. Zhang, L. Pan, Q.-L. Han, C. Chen, S. Wen, and Y. Xiang, “Deep learning based attack detection for cyber-physical system cybersecurity: A survey,” IEEE/CAA Journal of Automatica Sinica, vol. 9, no. 3, pp. 377–391, 2021.
  * [10] J. James, Y. Hou, and V. Li, “Online false data injection attack detection with wavelet transform and deep neural networks,” IEEE Transactions on Industrial Informatics, vol. 14, no. 7, pp. 3271–3280, 2018.
  * [11] V. Veerasamy, N. I. A. Wahab, M. L. Othman, S. Padmanaban, K. Sekar, R. Ramachandran, H. Hizam, A. Vinayagam, and M. Z. Islam, “Lstm recurrent neural network classifier for high impedance fault detection in solar pv integrated power system,” IEEE Access, vol. 9, pp. 32672–32687, 2021.
  * [12] A. Aligholian, A. Shahsavari, E. M. Stewart, E. Cortez, and H. Mohsenian-Rad, “Unsupervised event detection, clustering, and use case exposition in micro-pmu measurements,” IEEE Transactions on Smart Grid, vol. 12, no. 4, pp. 3624–3636, 2021.
  * [13] P. S. Sarker, M. F. Rafy, A. K. Srivastava, and R. K. Singh, “Cyber anomaly-aware distributed voltage control with active power curtailment and ders,” IEEE Transactions on Industry Applications, 2023, doi: 10.1109/TIA.2023.3328850.
  * [14] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational physics, vol. 378, pp. 686–707, 2019.
  * [15] M. J. Zideh, P. Chatterjee, and A. K. Srivastava, “Physics-informed machine learning for data anomaly detection, classification, localization, and mitigation: A review, challenges, and path forward,” arXiv preprint arXiv:2309.10788, 2023. https://doi.org/10.48550/arXiv.2309.10788.
  * [16] B. Huang and J. Wang, “Applications of physics-informed neural networks in power systems-a review,” IEEE Transactions on Power Systems, vol. 38, no. 1, pp. 572–588, 2022.
  * [17] F. Alotibi and D. Tipper, “Physics-informed cyber-attack detection in wind farms,” in GLOBECOM 2022-2022 IEEE Global Communications Conference, pp. 1–6, 2022.
  * [18] A. Ahmed, S. Basumallik, A. Gholami, S. K. Sadanandan, M. H. Namaki, A. K. Srivastava, and Y. Wu, “Spatio-temporal deep graph network for event detection, localization and classification in cyber-physical electric distribution system,” IEEE Transactions on Industrial Informatics, 2023.
  * [19] N. Aftabi, D. Li, and P. Ramanan, “A variational autoencoder framework for robust, physics-informed cyberattack recognition in industrial cyber-physical systems,” arXiv preprint arXiv:2310.06948, 2023. https://doi.org/10.48550/arXiv.2310.06948.
  * [20] B. Chen, Q. Wu, M. Li, and K. Xiahou, “Detection of false data injection attacks on power systems using graph edge-conditioned convolutional networks,” Protection and Control of Modern Power Systems, vol. 8, no. 1, pp. 1–12, 2023.
  * [21] B. Chen, Q. Wu, M. Li, and K. Xiahou, “Data‐driven detection and identification of iot‐enabled load‐altering attacks in power grids,” IET Smart Grid, vol. 5, no. 3, pp. 203–218, 2022.
  * [22] S. Bhattacharjee, V. P. K. Madhavarapu, S. Silvestri, and S. K. Das, “Attack context embedded data driven trust diagnostics in smart metering infrastructure,” ACM Transactions on Privacy and Security (TOPS), vol. 24, no. 2, pp. 1–36, 2021.
  * [23] “Weather Data,” National Renewable Energy Laboratory, [Online]. Available: https://sam.nrel.gov/weather-data.html. [Accessed 5 Dec 2022].


# Kicking Co and Rh atoms on a row-wise antiferromagnet

Felix Zahner  Institute of Nanostructure and Solid State Physics (INF),
University of Hamburg, Jungiusstraße 11, 20355 Hamburg, Germany  Soumyajyoti
Haldar<EMAIL_ADDRESS>Institute of Theoretical Physics and
Astrophysics, University of Kiel, Leibnizstrasse 15, 24098 Kiel, Germany
Roland Wiesendanger  Institute of Nanostructure and Solid State Physics (INF),
University of Hamburg, Jungiusstraße 11, 20355 Hamburg, Germany  Stefan Heinze
Institute of Theoretical Physics and Astrophysics, University of Kiel,
Leibnizstrasse 15, 24098 Kiel, Germany  Kiel Nano, Surface, and Interface
Science (KiNSIS), University of Kiel, 24118 Kiel, Germany  Kirsten von
Bergmann  Institute of Nanostructure and Solid State Physics (INF), University
of Hamburg, Jungiusstraße 11, 20355 Hamburg, Germany  André Kubetzka
<EMAIL_ADDRESS>Institute of Nanostructure and Solid State
Physics (INF), University of Hamburg, Jungiusstraße 11, 20355 Hamburg, Germany

###### Abstract

Diffusion on surfaces is a fundamental process in surface science, governing
nanostructure and film growth, molecular self-assembly, and chemical
reactions. Atom motion on non-magnetic surfaces has been studied extensively
both theoretically and by real-space imaging techniques. For magnetic surfaces
density functional theory (DFT) calculations have predicted strong effects of
the magnetic state onto adatom diffusion, but to date no corresponding
experimental data exists. Here, we investigate Co and Rh atoms on a hexagonal
magnetic layer, using scanning tunneling microscopy (STM) and DFT
calculations. Experimentally, we “kick” atoms by local voltage pulses and
thereby initiate strictly one-dimensional motion which is dictated by the row-
wise antiferromagnetic (AFM) state. Our calculations show that the one-
dimensional motion of Co and Rh atoms results from conserving the Co spin
direction during movement and avoiding high induced Rh spin moments,
respectively. These findings demonstrate that magnetism can be a means to
control adatom mobility.

The movement of atoms and molecules on surfaces is a fundamental and well-
studied process in surface science. Adatoms typically move thermally driven by
random jumps between neighboring binding sites, while longer jumps and
exchange processes with the surface are activated with increasing temperature
[1]. The surface symmetry plays a key role for these processes, see Fig. 1a,b:
whereas diffusion is effectively isotropic on hexagonal surfaces, lower
symmetry surfaces like fcc(110) can show strong directional anisotropies [2,
3], or even strictly one-dimensional movement [4]. Magnetic effects are
usually neglected, either because diffusion is studied above the magnetic
ordering temperature, $T_{\rm C}$, or because their strengths are considered
negligible compared to binding energies and diffusion barriers.

For bulk systems, however, experiments indicate that magnetic states can
affect diffusion. In bcc iron, for example, the self-diffusion coefficient
systematically deviates from a linear Arrhenius relation, i.e. below $T_{\rm
C}$ Fe atoms diffuse more slowly through the Fe crystal than extrapolated from
higher temperatures [5]. In some alloy systems, the structural composition can
be influenced by applying external magnetic fields during annealing; for the
antiferromagnet PdMn it has been argued that excess Pd atoms preferentially
diffuse to the magnetic sublattice which is favored by the external field [6].

On magnetic surfaces, so far no diffusion experiments have been performed, but
STM manipulation experiments have reported alternating forces, when a magnetic
atom is moved across an antiferromagnetic surface with a magnetic tip [7, 8].
On the theory side, almost 20 years ago DFT calculations for self-diffusion on
Fe and Mn monolayers on W(110) and W(001) have predicted a strong influence of
magnetism onto adsorption sites and diffusion properties with a general trend
of enhanced diffusion on ferromagnetic surfaces [9, 10], but respective
experiments are still missing.

Here, we investigate how Co and Rh atoms move at $T=4.2$ K on a hexagonal Mn
layer hosting a row-wise AFM state, when they are “kicked”, i.e. when motion
is initiated by local voltage pulses from an STM tip. For both adatom species
we observe strictly one-dimensional movement along the rows of Mn atoms with
parallel spins ($\uparrow\uparrow$-rows), see Fig. 1c. DFT calculations show
that the pseudomorphic Mn layer is shifted laterally by 15 pm off the fcc
hollow site positions on Re(0001), i.e. the Mn layer and adatom $C_{3}$
(threefold) symmetry is broken by the magnetic state electronically as well as
structurally. For both adatom species the energy barriers differ by about 50
meV for diffusion across Mn bridge sites with the two Mn atoms having the same
($\uparrow\uparrow$-bridge) compared to opposite spin orientation
($\uparrow\downarrow$-bridge), see Fig. 2c, making the potential landscape
highly anisotropic.

##  I Experimental Results

The row-wise AFM state within the hexagonal Mn layer on Re(0001) [11] allows
for three equivalent rotational magnetic domains, two of which are seen in the
spin-polarized (SP) STM data of Fig. 2a. Due to spin-polarized tunneling,
neighboring atomic rows display different signal strengths [12] and the left
domain exhibits a stronger magnetic contrast for this particular tip
magnetization direction [8]. At lower tunnel voltage and smaller tip–sample
distance (higher current) in Fig. 2b the hexagonal atomic lattice becomes
visible and an additional stripe pattern appears with half the period of the
magnetic state. Since all Mn atoms are electronically equivalent, this atomic
stripe pattern must stem from the inequivalent $\uparrow\uparrow$\- and
$\uparrow\downarrow$-bridge sites, see Fig. 2c, indicating that the $C_{3}$
symmetry of the Mn hollow sites is broken. For the same reason many native
defects display a boomerang-like shape which is correlated with the magnetic
state, see Fig. 2c,d and Methods. This defect shape, which has only one mirror
plane, reflects a loss of electronic $C_{3}$ surface symmetry and allows to
determine the orientation of the magnetic domains with a non-magnetic STM tip.
Fig. 2d shows different kinds of these defects with high resolution as well as
the $\uparrow\uparrow$-rows of the magnetic state. In this non-magnetic data,
neighboring $\uparrow\uparrow$-rows appear identical.

With a different non-magnetic tip in Fig. 2e, the Mn layer itself appears
featureless, but the direction of the $\uparrow\uparrow$-rows can still be
inferred from the defect shapes. Co adatoms on fcc-Mn/Re(0001) are easily
moved during standard STM imaging: in Fig. 2e a single Co atom jumps along the
$\uparrow\uparrow$-rows, following the slow scan direction of the tip upward,
and it therefore appears like an atomic chain. Only once does the atom jump
one atomic row to the right (see arrow in Fig. 2e), flipping its magnetic
moment in the process, see Fig. 1c. This data is reminiscent of the
experiments performed by Li et al. with Ag atoms on Ag(110) [3, 13] and is a
first hint that the magnetic state dictates the movement direction of Co
atoms. The line section in Fig. 2f shows the period of the atomic lattice,
indicating that only one type of hollow site is a stable position for Co, in
agreement with no zigzag movement being observed. On all other positions, the
residence time is too short to be detectable in this measurement. We will see
in section II that the Co atoms prefer the so-called hcp hollow site above the
surface Re atom, see red crosses in Fig. 2c.

To minimize the influence of the tip, instead of dragging the adatom along
during imaging, in the following we “kick” the adatoms with a voltage pulse
from a stationary tip and determine the new position by subsequent imaging.
The kicking is performed by a sudden voltage increase at constant tip height
above the atom, which also leads to an increased tunnel current, see inset of
Fig. 3a and Methods. By performing these experiments with one and the same
micro-tip on different rotational domains, we can exclude a possible tip
asymmetry being responsible for the movement direction. Figure 3a shows a
surface area with two rotational domains like in Fig. 2a, this time imaged
with vanishing spin contrast. The domain wall (DW), now visible via a reduced
electron density [8], is marked by a white line. We find that the Co movement
initiated by kicking is always toward the upper left on the left domain and
toward the top of the image on the right domain, in both cases following the
$\uparrow\uparrow$-rows of the magnetic state (Supplementary Video 1).
Examples are given in Fig. 3b,c (left domain) and Fig. 3d,e (right domain).
Independent of the micro-tip used for kicking, we always observe a strictly
one-dimensional movement of Co atoms along the $\uparrow\uparrow$-rows, within
an accuracy of one atomic lattice site. Occasionally, Co atoms which are
further away from the directly kicked atom also move 1–2 lattice sites, see
Fig. 3b–d, indicating a weak long range effect. These short jumps also occur
exclusively along the $\uparrow\uparrow$-rows.

The Co atoms are surprisingly mobile on fcc-Mn/Re(0001). Depending on the
micro-tip, we find a threshold voltage for single-site jumps at rather low
voltage, in the range of 8–50 mV at $I=1$ nA. In Fig. 3b–e Co atoms move 10–15
lattice sites, when kicked by voltage pulses of $U_{\text{P}}=+200$ mV. A
further measurement series of 19 kicking events, using both voltage polarities
(see also section IV.4) is shown in Fig. 3f: the average travel distance is
similar and the distribution is roughly uniform with an accumulation at the
far end at 7–8 nm. We observe a maximum kicking distance of about 10 nm at
$U_{\text{P}}=500$ mV. A random walk type motion of uncorrelated short jumps
would lead to a normal distribution with a maximum at zero and an exponential
tail [1], and can therefore be ruled out. Instead, these results indicate low
diffusion barriers along the $\uparrow\uparrow$-rows and an effectively low
damping, allowing the Co atoms to perform long jumps [14, 15].

To make a more general case we repeated the above experiments with a second
adatom species: Rh is in the same chemical group and thus iso-electronic to
Co, but almost 75% heavier and typically non-magnetic on metal surfaces. We
again use a surface area with two rotational domains to exclude tip artifacts,
see Fig. 4a. We similarly find that movement is always along the
$\uparrow\uparrow$-rows (Supplementary Video 2), as can be seen in Fig. 4b,c
(upper domain) and in Fig. 4d,e (lower domain). The Rh atoms show no spherical
symmetry but instead reflect the broken $C_{3}$ symmetry, like many of the
native defects. In comparison, the Co shows an asymmetry in the extended
electron density surrounding the atom, see Fig. 3a–e. Apparently, Rh is much
less mobile than Co: at $U_{\text{P}}=+200$ mV Rh atoms do not move more than
one atomic site and at $U_{\text{P}}=+1$ V they travel only 1.0–2.5 nm, see
Fig. 4b–e. We observe a maximum travel distance of 3 nm at $U_{\text{P}}=+2$
V.

For both atom species, magnetic Co and non-magnetic Rh, we thus find strictly
one-dimensional movement along the $\uparrow\uparrow$-rows of the AFM state,
within an accuracy of one atomic site. The critical parameter to move an atom
seems to be voltage magnitude $|U_{\text{P}}|$. The threshold voltage is tip-
dependent and element specific with $U_{\text{P}}=8$–50 mV for Co and
$U_{\text{P}}=200$–300 mV for Rh.

##  II First-principles calculations

In order to explain the experimental observations, we studied Co and Rh
adatoms on fcc-Mn/Re(0001) using density functional theory (DFT). First we
discuss the structural and magnetic properties of fcc-Mn/Re(0001). From our
DFT total energy calculations we find that the row-wise AFM state is by 300
meV/Mn atom lower than the ferromagnetic (FM) state. In the structural
relaxations we have taken into account the symmetry breaking due to the row-
wise AFM state highlighted by the experiments. Surprisingly, we find in
addition to a relaxed Mn-Re interlayer distance, a lateral shift by $\sim$15
pm perpendicular to the $\uparrow\uparrow$-rows. Thereby, Re surface atoms are
closer to the two adjacent nearest-neighbor Mn atoms with parallel magnetic
moments ($m_{\rm Mn}\approx 3.4$ $\mu_{\rm B}$) with respect to the nearest-
neighbor Mn atom with an opposite spin moment (Fig. 5a). As a result the
induced magnetic moments of the Re surface atoms are enhanced, strengthening
the antiferromagnetic Mn-Re interaction 111The lateral shift of the Mn layer
results in a total energy gain of 18 meV/Mn atom for the row-wise AFM state.
Thereby, it becomes the magnetic ground state of the system since it is lower
than all spin spiral states as well as the triple-Q state [11].. Note, that
the magnetic moment direction of the Mn atoms is in-plane along the rows of
parallel moments (Fig. 5a) due to the easy-plane magnetocrystalline anisotropy
and the anisotropic symmetric exchange interaction [11].

To understand the experimentally found one-dimensional motion of Rh and Co
atoms we calculate the minimum energy paths and the related energy barriers
for different directions of atom movement using the climbing image nudged
elastic band method (NEB) [17, 18] (see Methods). We begin with an adatom in
its preferred adsorption site (I1), and move it to an equivalent site both
along an $\uparrow\uparrow$-row and along an $\uparrow\downarrow$-row (Fig.
5a) where the final states are indicated as F1 and F3, respectively, and
exhibit opposite spin directions.

We start our discussion with Rh atom movement along the
$\uparrow\uparrow$-rows, i.e. from I1 to F3 (blue lines and blue transparent
spheres in Figs. 5a,b). We find a symmetric minimum energy path with two
identical saddle points at the $\uparrow\downarrow$-bridge sites. The
intermediate energy minimum along this path is the fcc hollow site which is by
only $\approx 4$ meV unfavorable with respect to the adsorption site I1 for a
Rh atom. Note, that the induced magnetic moment of Rh deviates on the path
only little from its value of $0.04$ $\mu_{\rm B}$ at I1. The energy at the
saddle points of $\sim$180 meV corresponds to the energy barrier for Rh atom
movement along an $\uparrow\uparrow$-row.

In contrast, the minimum energy path for Rh exhibits two saddle points of
different barrier height along an $\uparrow\downarrow$-row, i.e. from I1 to F1
(red lines and red transparent spheres in Fig. 5a,b). The first saddle point
results in a barrier of $\sim$240 meV and is due to the movement of the atom
across the $\uparrow\uparrow$-bridge site. At this saddle point the induced Rh
magnetic moment is much enhanced (0.3 $\mu_{\rm B}$). A second smaller energy
barrier ($\sim$180 meV) occurs when the Rh atom starts from the intermediate
state at the fcc hollow site and crosses the second saddle point, the
$\uparrow\downarrow$-bridge site, to reach the final state F1. Therefore, the
overall energy barrier to move a Rh atom along an $\uparrow\downarrow$-row
amounts to $\sim$240 meV. These calculations are consistent with the favored
movement of Rh atoms along $\uparrow\uparrow$-rows observed in our
experiments. The higher barrier for Rh movement along
$\uparrow\downarrow$-rows is related to a larger induced magnetic moment at
the $\uparrow\uparrow$-bridge site and is thus of magnetic origin.

For the Co atom, we also find a lower energy barrier for the path along the
$\uparrow\uparrow$-rows, i.e. from I1 to F3 (blue transparent spheres and blue
lines in Fig. 5a,c). Due to exchange coupling, the Co magnetic moment is
aligned parallel to the nearest-neighbor Mn atoms both at position I1 and F3.
Along the minimum energy path the Co moment remains along this direction,
however, its magnitude varies slightly. The Co atom crosses two identical
saddle points given by $\uparrow\downarrow$-bridge sites resulting in an
energy barrier height of $\sim$150 meV.

For the minimum energy path of a Co atom along an $\uparrow\downarrow$-row
(red transparent spheres and red lines in Fig. 5a,c), its magnetic moment must
flip from $+1.6$ $\mu_{\rm B}$ at I1 to $-1.6$ $\mu_{\rm B}$ at F1 in order to
be in the energetically favorable spin configuration in both initial and final
state. As for Rh, we observe two different energy barriers along this path
with a lower barrier of $\sim$110 meV at the $\uparrow\uparrow$-bridge site
and a more than twice as high barrier of $\sim$260 meV at the
$\uparrow\downarrow$-bridge site. Note, that at the latter saddle point the Co
magnetic moment is reduced to $0.4$ $\mu_{\rm B}$ which results from the
required spin-flip along the path and the constraint to collinear magnetic
moments in our calculation. With an overall barrier height of $\sim$260 meV,
this path is clearly unfavorable with respect to Co movement along
$\uparrow\uparrow$-rows consistent with our experimental observations.

To reveal the role of the Co spin-flip on the energy barrier for the path
along $\uparrow\downarrow$-rows, we have performed a minimum energy path
calculation from I2 to F2 (turquoise transparent spheres and line in Fig.
5a,c). Along this path the magnetic moment of the Co atom stays in the same
direction and a large reduction of the magnetic Co moment on the
$\uparrow\downarrow$-bridge site is avoided. Since the initial state I2 of
this path coincides with the intermediate energy minimum along the path from
I1 to F1 (red line in Fig. 5c), we can combine the first half of the path
I1-F1 with the path I2-F2 to obtain a minimum energy path from I1 to F2. On
this path along an $\uparrow\downarrow$-row the magnetic moment of the Co atom
does not flip and it amounts to 1.5 $\mu_{\rm B}$ at the
$\uparrow\downarrow$-bridge site. The spin of the Co atom may flip after the
metastable state F2 has been reached such that a relaxation into the favorable
state F1 is achieved. Nevertheless, the energy along the path from I1 to F2
without a spin-flip is asymmetric with two saddle points of different energy
as for the I1-F1 path with a spin flip (red line in Fig. 5c). Since the energy
barrier due to the saddle point between I2 and F2 is higher than that for the
path along the $\uparrow\uparrow$-rows (turquoise vs. blue line in Fig. 5c)
the latter remains the favorable path for Co atom motion.

##  III Discussion

The first-principles calculations nicely explain the experimentally observed
moving direction of Co and Rh atoms by anisotropic potential landscapes
arising from the row-wise AFM state. The spontaneous magnetic ordering causes
an electronic anisotropy, which drives the $\sim$15 pm lateral shift of the Mn
layer relative to the Re(0001) surface. The resulting similarity to diffusion
on surfaces with low structural symmetry such as Ag(110) [13] or Cu(211) [19]
is striking. However, there are two important differences: antiferromagnetic
domains can be switched by lateral currents [20], which in principle allows to
externally control diffusion directions on antiferromagnetic surfaces.
Furthermore, since the Mn layer itself remains structurally hexagonal, the
$\uparrow\uparrow$-rows are not necessarily the easy diffusion direction for
all atoms or molecules. We have performed additional DFT calculations which
show for instance that Cu atoms behave inversely to Rh atoms, i.e. the energy
barrier at the $\uparrow\uparrow$-bridge is lower than at the
$\uparrow\downarrow$-bridge, which hinders Cu diffusion along the
$\uparrow\uparrow$-rows.

The experimentally found higher mobility, i.e. longer travel distances of Co
atoms, agrees qualitatively with a lower barrier of $\sim$150 meV compared to
$\sim$180 meV for Rh. Additional effects beyond this adiabatic analysis might
contribute to the high mobility of Co. Firstly, our DFT calculations show a
total binding energy of Co atoms on Mn/Re(0001) which is 1.5 eV lower than for
Rh atoms. This means that a ballistic movement at an increased distance from
the surface [21] requires less energy for Co atoms than for Rh atoms. Since in
the experiment some of the energy of the tunneling electrons is used for the
excitation, the excess energy could be larger for the Co atoms. This is in
accordance with the larger threshold voltage observed for Rh atoms. Secondly,
the $3d$ orbitals of Co have a smaller extent than the $4d$ orbitals of Rh,
which may result in a reduced interaction with the surface and a reduced
damping for a moving Co atom. Thirdly, Rh is about 75% heavier. The lighter
and faster moving Co atom can therefore cover more distance before it comes to
rest. This argument only holds for damping processes which scale with travel
time [22]. Finally, the spin degree of freedom might contribute to the travel
distance. Long jumps have been found experimentally for molecules [23, 24] and
one reason is their ability to store energy during movement in form of
vibrations. Spin excitations of the Co atom and within the Mn layer might play
a similar role here.

In conclusion, the observed one-dimensional movement demonstrates that
magnetism can play a decisive role in controlling atom movement, with possible
consequences for related phenomena such as nanostructure growth, molecular
self-assembly and catalysis.

##  IV Methods

###  IV.1 Experimental details

The experiments were performed in a multi-chamber ultra-high vacuum system
with different chambers for substrate cleaning, metal deposition, and STM
measurements. The Re(0001) single crystal surface was cleaned by cycles of
annealing in an oxygen atmosphere of $10^{-7}$–$10^{-8}$ mbar at temperatures
of up to $T=1400$ K; before metal deposition a final flash to $T=1800$ K was
performed. The Mn was evaporated from a pyrolytic boron nitride (PBN) Knudsen
cell of volume 1 cm3 and held at $T=670$ ∘C, resulting in a flux of
approximately 0.1 atomic layers per minute; during Mn deposition the Re single
crystal was still at elevated temperature ($\approx 100$ ∘C) from the final
flash.

The Mn/Re(0001) samples were then cooled to $T=4.2$ K within the STM. After a
quick transfer to the e-beam evaporator ($\Delta t<5$ s) Co or Rh atoms were
evaporated from a rod of 2 and 1 mm diameter, respectively, at a rate of about
0.1 atomic layers per minute. We estimate a surface temperature of $T=25$ K
$\pm 5$ K during adatom deposition.

We use a Cr bulk STM tip, etched in 1M HCl solution. Within the STM the tip
was cleaned by field emission on a W(110) surface and sharpened by voltage
pulses of $U=4$–6 V and gentle surface collisions on the Mn/Re(0001) sample.
Except for the measurement shown in Fig. 2a,b we did not optimize the tip for
magnetic contrast. Consequently, in all other data the tip exhibits no
significant spin sensitivity at the used tunnel voltages, thereby simplifying
the data analysis. Magnetic domain walls were imaged by making use of their
(spin-averaged) electronic signal [8] and the orientation of rotational
domains was determined via the symmetry of native defects: Extended Data Fig.
1 shows a surface area with all three rotational domains, demonstrating the
strict correlation of domain orientation and defect shape orientation.

For atom “kicking”, the Co or Rh atom was first moved to a defect-free area
and the tip was then stabilized precisely above the atom, using an atom-
tracking feedback by steepest ascent [25] and tunnel parameters on the order
of $U=+10$ mV and $I=200$ pA. The atom-tracking and height-control feedback
were then switched off and the voltage raised abruptly at constant tip height
for $\Delta t=0.4$–1.0 s to $U_{\text{P}}=\pm 200$ mV and $U_{\text{P}}=+1$ V,
respectively. A voltage pulse is thus automatically accompanied by an
increased current through the adatom, but only as long as the atom does not
move. For the measurement series displayed in Fig. 3f, the tip was stabilized
above the atoms at $U=+10$ mV and $I=200$ pA, and after the kick the atoms
were moved back to their original positions using constant current pulling
mode at approximately $U=4$ mV and $I=60$ nA. Since micro-tip changes are more
likely during manipulation than during imaging, the necessity to move an atom
back into position, due to surrounding defects, limits the number of kicks
which can be performed with the same micro-tip. Which of the two equivalent
directions along the $\uparrow\uparrow$-rows is taken, depends on the specific
micro-tip: for some tips, like the one used in Fig. 3a–e, one movement
direction is favored; other tips show a more symmetric distribution. The data
of Fig. 3 has been smoothed by a 3-pixel Gauss filter; all other figures show
raw data.

###  IV.2 Computational details

Our first-principles calculations are based on density-functional theory (DFT)
within the projector augmented wave (PAW) [26] method as implemented in the
vasp code [27, 28]. We have performed spin-polarized calculations using
Wigner-Seitz radii for the elements, i.e. $R_{\rm Mn}^{\rm WS}=1.32$ Å for Mn
and $R_{\rm Re}^{\rm WS}=1.43$ Å for Re. As lattice parameters we used
$a_{NN}=2.78$ Å , $c=4.49$ Å, which were taken from Ref. 29 obtained via DFT
within the generalized gradient approximation (GGA). Films with fcc-stacking
of the Mn monolayer were structurally relaxed using the GGA exchange-
correlation (xc) potential [30] in the row-wise AFM state using an asymmetric
film consisting of six Re(0001) layers with a Mn layer on top. The Mn layer
and the uppermost Re layer were relaxed in $x$, $y$, and $z$-direction, where
the bulk reference is $c/2=2.24$ Å. We used $5\times 5\times 1$ $k$ points in
the irreducible wedge of the two-dimensional (2D) Brillouin zone (BZ) and a
plane wave basis set cutoff of $E_{max}=350$ eV. The structures were optimized
using the conjugate gradient method with forces calculated from the Hellman-
Feynman theorem. Structures were considered to have been optimized when all
the forces were smaller than 0.01 eV/Å. We have used a $(5\times 6)$ super
cell in order to design the row-wise AFM state.

The Co atom preferentially adsorbs in the hcp hollow site (indicated by I1),
which is about 50 meV lower in energy compared to the fcc hollow site (I2). In
both cases a small lateral shift relative to the center of the Mn atom
triangle is found. The configuration of the three nearest Mn magnetic moments
results in a net ferromagnetic coupling to the Co magnetic moment ($1.6$
$\mu_{\rm B}$), where two of them are parallel and one is antiparallel to the
Co moment (Fig. 5a). For the Rh atom the two hollow sites are almost
energetically degenerate, with a difference of only $\sim$4 meV, and in the
hcp site, the adjacent Mn atoms induce a Rh magnetic moment of $\sim$0.04
$\mu_{\rm B}$.

The magnetic moment of the Co adatom at hcp hollow site F1 (Fig. 5a) is $-1.6$
$\mu_{\rm B}$ due the magnetic moment direction of the two nearest-neighbor Mn
moments (net ferromagnetic coupling). In order to obtain a net
antiferromagnetic coupling to the Co magnetic moment with the three nearest
neighbor Mn atoms, we perturbed the position of the Co adatom by a small
amount ($\sim$0.05 Å) and relaxed the structure with an initial high value of
positive magnetic moment for the Co adatom (net antiferromagnetic coupling
with the three nearest neighbor Mn atoms). Upon structural relaxation, we find
that the Co atom moved to position F2 which is shifted by 15 pm with respect
to F1 in the direction perpendicular to the $\uparrow\uparrow$-rows. At this
site the magnetic moment of the Co atom is $+1.2$ $\mu_{\rm B}$ with a net
antiferromagnetic coupling to the three nearest neighbor Mn atoms. Our
calculation shows that F2 is a higher energy local minima state by $\sim$40
meV compared to the F1 state with an opposite spin moment.

We calculate the energy barriers and the minimum energy paths between two
known adatom adsorption sites by using the climbing image nudged elastic band
method (NEB) [17, 18]. In the NEB method we optimize a number of intermediate
images to their lowest energy possible while maintaining equal spacing to
neighboring images. This is obtained by adding a spring force along the band
between the images. The climbing image method is a small modification to the
NEB method where the highest energy image is pushed up to the saddle point to
get the exact position of the saddle point. This is done by inverting the true
force for this image along the tangent. Thus this image does not feel the
spring force. In our NEB calculations we have not relaxed the surface atoms
for simplicity. We have used a mixture of quick-min and LBFGS optimizers [31]
and $-$5.0 eV/Å2 spring force for our calculations.

For the simulated STM and SP-STM images of Mn/Re(0001), we have used the full-
potential linearized augmented plane wave method (FLAPW) as implemented in the
FLEUR code [32]. We have performed spin-polarized calculations using radii for
the muffin tin spheres (in all FLEUR calculations), i.e. $R^{\rm Mn}_{\rm MT}$
= 2.35 a.u. for Mn and $R^{\rm Re}_{\rm MT}$ = 2.40 a.u. for Re. For the FLEUR
calculations of Mn/Re(0001) we have used asymmetric films with 6 layers of Re
substrate. Here we used 1444 $k$ points in the full Brillouin zone and an
energy cutoff of $4.3\,\rm{a.u.}^{-1}$.

###  IV.3 Simulated STM/SP-STM images of Mn/Re(0001)

A comparison of simulated STM and SP-STM images with and without a shift of
the pseudomorphic Mn layer in the row-wise AFM state is shown in Extended Data
Fig. 2 and 3 for a small energy window below and above the Fermi level,
respectively. The STM and SP-STM images are calculated based on the electronic
structure from DFT using the Tersoff-Hamann model [33] and its extension to
spin-polarized STM [12].

The simulated SP-STM images (Extended Data Fig. 2b,d and 3b,d) show the
expected stripe period of the row-wise AFM state observed experimentally in
Fig. 2a. This magnetic contrast is the same for the unshifted and for the
laterally shifted Mn layer. In the simulated STM images for a positive bias
voltage (Extended Data Fig. 3c), we find for the laterally shifted Mn layer
also the experimentally observed stripes of atomic period shown in Fig. 2d.
This stripe contrast is of electronic origin and it is absent in the simulated
STM images of the unshifted Mn layer for both positive and negative bias
voltages (Extended Data Fig. 2a and 3a). This supports that the lateral shift
of the Mn layer predicted by our DFT total energy calculations occurs in the
experiment.

Note, that in the experimental SP-STM image shown in Fig. 2b, the electronic
and magnetic stripes are not in phase as in Extended Data Fig. 3c,d. However,
in the experiment of Fig. 2b the tunneling current was extremely large, i.e.
the tip-sample distance was very small, which might lead to additional effects
due to tip-sample interaction. Another possible reason for the difference with
this particular image may be specific micro-tip properties.

###  IV.4 Co electric dipole and voltage polarity

Concerning voltage polarity, the data of Fig. 3f shows a trend to larger
kicking distances for negative polarity, a much weaker effect than observed
for other systems [34, 35, 24]. In further experiments with $U_{\text{P}}=\pm
500$ mV we found no statistically significant effect of polarity on travel
distance, indicating that the local electric field does not dominate the Co
atom movement. In addition, other effects can cause a voltage polarity
dependence, for instance an asymmetric density of states with respect to the
Fermi level. Indeed, the vacuum density of states of Mn/Re(0001) is about four
times higher 200 meV below compared to 200 meV above the Fermi level (see Fig.
S9 in Ref. 11). This leads to a higher current flowing at pulses of
$U_{\text{P}}=-200$ mV compared to $U_{\text{P}}=+200$ mV, which might explain
the polarity trend observed in Fig. 3f.

To further investigate the impact of the tip’s local electric field onto the
atom movement we have calculated the Co electric dipole using the Bader charge
analysis [36]. It turns out that the value is sizable with $\mu=0.61$ eÅ or
2.9 D; for comparison, H2O has an electric dipole of $\mu=1.85$ D. When a
charge-neutral Co atom is placed on the Mn layer, electrons are transferred
from the surface toward the Co atom, producing an electric dipole pointing
down. Neglecting work function differences between tip and sample, at positive
(negative) voltage the tip’s electric field points up (down), decaying with
lateral distance from the tip position. Thus, a repulse (attractive) force
should act on the Co atom at positive (negative) voltage. This means that the
weak polarity trend observed in Fig. 3f cannot be explained by the tip’s
electric field. Thus, despite the sizable Co dipole, electric fields do not
play a role in the experiments.

##  V Data availability

The STM and DFT data are available upon reasonable request.

##  VI Code availability

The STM data was analyzed with the open access software Gwyddion
(http://gwyddion.net).

##  VII Acknowledgements

We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under project numbers 408119516, 418425860, and
445697818. This work was performed using HPC resources from the North-German
Supercomputing Alliance (HLRN).

##  VIII Author contributions

A.K. devised the experiments. A.K. and F.Z. prepared the samples, performed
the STM experiments, and together with K.v.B. analyzed the experimental data.
S.H. performed and analyzed the calculations with St.H. The figures were
prepared by A.K., F.Z., and S.H. The manuscript was written by A.K., S.H., and
St.H., with all authors contributing.

## References

  * Antczak and Ehrlich [2007] G. Antczak and G. Ehrlich, Jump processes in surface diffusion, Surf. Sci. Rep. 62, 39 (2007).
  * Frenken  _et al._ [1990] J. W. M. Frenken, B. J. Hinch, J. P. Toennies, and C. Wöll, Anisotropic diffusion at a melting surface studied with He-atom scattering, Phys. Rev. B 41, 938 (1990).
  * Li  _et al._ [1996] J. Li, R. Berndt, and W.-D. Schneider, Tip-Assisted Diffusion on Ag(110) in Scanning Tunneling Microscopy, Phys. Rev. Lett. 76, 1888 (1996).
  * Swartzentruber [1996] B. S. Swartzentruber, Direct Measurement of Surface Diffusion Using Atom-Tracking Scanning Tunneling Microscopy, Phys. Rev. Lett. 76, 459 (1996).
  * Iijima  _et al._ [1988] Y. Iijima, K. Kimura, and K. Hirano, Self-diffusion and isotope effect in $\alpha$-iron, Acta Metall. 36, 2811 (1988).
  * Josten  _et al._ [2023] N. Josten, O. Miroshkina, M. Acet, M. E. Gruner, and M. Farle, Annealing time, temperature, and field dependence of pinned magnetic moments in the collinear antiferromagnet PdMn, Phys. Rev. B 108, 064417 (2023).
  * Wolter  _et al._ [2012] B. Wolter, Y. Yoshida, A. Kubetzka, S.-W. Hla, K. von Bergmann, and R. Wiesendanger, Spin Friction Observed on the Atomic Scale, Phys. Rev. Lett. 109, 116102 (2012).
  * Spethmann  _et al._ [2021] J. Spethmann, M. Grünebohm, R. Wiesendanger, K. von Bergmann, and A. Kubetzka, Discovery and characterization of a new type of domain wall in a row-wise antiferromagnet, Nat. Commun. 12, 3488 (2021).
  * Spišák and Hafner [2004] D. Spišák and J. Hafner, Diffusion of Fe atoms on W surfaces and Fe/W films and along surface steps, Phys. Rev. B 70, 195426 (2004).
  * Dennler and Hafner [2005] S. Dennler and J. Hafner, First-principles study of ultrathin magnetic Mn films on W surfaces. II. Surface diffusion, Phys. Rev. B 72, 214414 (2005).
  * Spethmann  _et al._ [2020] J. Spethmann, S. Meyer, K. von Bergmann, R. Wiesendanger, S. Heinze, and A. Kubetzka, Discovery of Magnetic Single- and Triple-$\mathbf{q}$ States in Mn/Re(0001), Phys. Rev. Lett. 124, 227203 (2020).
  * Wortmann  _et al._ [2001] D. Wortmann, S. Heinze, P. Kurz, G. Bihlmayer, and S. Blügel, Resolving Complex Atomic-Scale Spin Structures by Spin-Polarized Scanning Tunneling Microscopy, Phys. Rev. Lett. 86, 4132 (2001).
  * Li  _et al._ [1998] J. Li, W.-D. Schneider, and R. Berndt, Low-temperature manipulation of Ag atoms and clusters on a Ag(110) surface, Appl. Phys. A 66, S675 (1998).
  * Pollak  _et al._ [1993] E. Pollak, J. Bader, B. J. Berne, and P. Talkner, Theory of correlated hops in surface diffusion, Phys. Rev. Lett. 70, 3299 (1993).
  * Ferrando  _et al._ [1994] R. Ferrando, R. Spadacini, G. Tommei, and G. Caratti, Correlation functions in surface diffusion: the multiple-jump regime, Surf. Sci. 311, 411 (1994).
  * Note [1] The lateral shift of the Mn layer results in a total energy gain of 18 meV/Mn atom for the row-wise AFM state. Thereby, it becomes the magnetic ground state of the system since it is lower than all spin spiral states as well as the triple-Q state [11].
  * Henkelman  _et al._ [2000] G. Henkelman, B. P. Uberuaga, and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113, 9901 (2000).
  * Henkelman and Jónsson [2000] G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113, 9978 (2000).
  * Meyer  _et al._ [1195] G. Meyer, B. Neu, and K. Rieder, Controlled lateral manipulation of single molecules with the scanning tunneling microscope, Appl. Phys. A 60, 343 (1195).
  * Baltz  _et al._ [2018] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Antiferromagnetic spintronics, Rev. Mod. Phys. 90, 015005 (2018).
  * Guo  _et al._ [2015] S. Y. Guo, S. J. Jenkins, W. Ji, Z. Ning, J. C. Polanyi, M. Sacchi, and C.-G. Wang, Repulsion-Induced Surface-Migration by Ballistics and Bounce, J. Phys. Chem. Lett. 6, 4093 (2015).
  * Brune  _et al._ [1992] H. Brune, J. Wintterlin, R. J. Behm, and G. Ertl, Surface migration of “hot” adatoms in the course of dissociative chemisorption of oxygen on Al(111), Phys. Rev. Lett. 68, 624 (1992).
  * Hla  _et al._ [2004] S.-W. Hla, K.-F. Braun, B. Wassermann, and K.-H. Rieder, Controlled Low-Temperature Molecular Manipulation of Sexiphenyl Molecules on Ag(111) Using Scanning Tunneling Microscopy, Phys. Rev. Lett. 93, 208302 (2004).
  * Civita  _et al._ [2020] D. Civita, M. Kolmer, G. J. Simpson, A.-P. Li, S. Hecht, and L. Grill, Control of long-distance motion of single molecules on a surface, Science 370, 957 (2020).
  * [25] See https://www.specs-group.com/nc/nanonis/products/detail/atom-tracking/ for details.
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Ji  _et al._ [2016] D.-P. Ji, Q. Zhu, and S.-Q. Wang, Detailed first-principles studies on surface energy and work function of hexagonal metals, Surf. Sci. 651, 137 (2016).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Sheppard  _et al._ [2008] D. Sheppard, R. Terrell, and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys. 128, 134106 (2008).
  * [32] See https://www.flapw.de.
  * Tersoff and Hamann [1983] J. Tersoff and D. R. Hamann, Theory and Application for the Scanning Tunneling Microscope, Phys. Rev. Lett. 50, 1998 (1983).
  * Whitman  _et al._ [1991] L. J. Whitman, J. A. Stroscio, R. A. Dragoset, and R. J. Celotta, Manipulation of Adsorbed Atoms and Creation of New Structures on Room-Temperature Surfaces with a Scanning Tunneling Microscope, Science 251, 1206 (1991).
  * Stroscio and Eigler [1991] J. A. Stroscio and D. M. Eigler, Atomic and Molecular Manipulation with the Scanning Tunneling Microscope, Science 254, 1319 (1991).
  * Tang  _et al._ [2009] W. Tang, E. Sanville, and G. Henkelman, A grid-based Bader analysis algorithm without lattice bias, J. Phys.: Condens. Matter 21, 084204 (2009).

Figure 1: — Structural and magnetic surface symmetry. a, Quasi-isotropic
diffusion of an adatom on a hexagonal surface with six structurally equivalent
directions. b, Fcc(110) surface with typically preferred diffusion along
close-packed atomic rows [3, 13]. c, Adatom on a row-wise antiferromagnet,
with the $\uparrow\uparrow$-rows breaking the hexagonal surface symmetry. Its
magnetic moment is assumed to be parallel to the majority of neighboring Mn
surface atoms. Figure 2: — Broken surface symmetry and one-dimensional Co atom
movement. a, Two rotational domains of the row-wise AFM state in fcc-
Mn/Re(0001), imaged by constant current SP-STM ($U=+14$ mV, $I=1$ nA). b, With
the tip closer to the surface ($U=+5$ mV, $I=70$ nA), additional stripes of
atomic period appear. c, Sketch of the magnetic state indicating different
hollow sites (red and black crosses) and inequivalent $\uparrow\uparrow$\- and
$\uparrow\downarrow$-bridge sites. d, High resolution non-magnetic STM image
($U=+30$ mV, $I=1$ nA) of the $\uparrow\uparrow$-rows with native defects
reflecting the broken symmetry of the hollow site, see inset. e, A single Co
atom is imaged while it is moving to the top of the image, guided by the
$\uparrow\uparrow$-rows (slow scan direction: bottom$\rightarrow$top, fast
scan direction: left$\rightarrow$right). f, The line section shows that the Co
atom is imaged with atomic periodicity of 2.7 Å, see red crosses in c. Figure
3: — Kicking Co atoms on a row-wise antiferromagnet. a, Overview STM image of
Co atoms on fcc-Mn/Re(0001) with two rotational domains; their orientation is
apparent by the defect shapes and the domain wall (DW) is imaged via a reduced
electron density. Kicking is performed by a sudden voltage raise at constant
tip height, see inset. b, The upper left Co atom has moved along the
$\uparrow\uparrow$-rows after kicking. c, The bottom Co atom has moved in the
same direction after kicking. d, The bottom Co atom has moved along the
$\uparrow\uparrow$-rows after kicking. e, The left Co atom has moved in the
same direction after kicking. (All images $U=+8$ mV, $I=1$ nA, all voltage
pulses: $U_{\text{P}}=+200$ mV, $\Delta t=1$ s, tip kicking locations are
marked by crosses in the previous image, previous atom positions are marked by
circles.) f, Kicking distances for both voltage polarities (color-coded),
$|U_{\text{P}}|=200$ mV, $\Delta t=0.5$ s. Figure 4: — Kicking Rh atoms on a
row-wise antiferromagnet. a, Overview STM image of Rh atoms on fcc-Mn/Re(0001)
with two rotational domains. A domain wall (DW) is marked by a black line.
b–e, When kicked, the Rh atoms also move according to the rotational domain
state, but travel shorter distances compared to Co, despite the use of pulses
with higher voltage magnitude. The electron density of native defects and also
the Rh atoms is asymmetric due to the AFM state. (All images $U=+25$ mV, $I=2$
nA, $U_{\text{P}}=+1$ V, $\Delta t=400$ ms). Figure 5: — DFT minimum energy
paths and energy barriers for Rh and Co atoms on Mn/Re(0001). a, Minimum
energy path for the quasi one-dimensional movement of Rh and Co adatoms
between close-packed atomic rows of the fcc-Mn/Re(0001) surface in the row-
wise AFM state obtained via DFT nudged-elastic band calculations. Red and blue
spheres with arrows represent Mn atoms with opposite spin directions. Yellow
spheres denote Re surface atoms, while small green and cyan spheres represent
adatoms. Note that the lateral shift of the Mn layer with respect to the Re
surface obtained from DFT is included in the sketch. Red and blue transparent
spheres indicate the calculated paths from the hcp hollow adsorption site
(initial state 1, I1) to two similar adsorption sites with opposite spin
directions (final states F1 and F3). For the Co atom the minimum energy path
is also given from the initial state I2 to the final state F2 (cyan spheres).
b, c Energy along the minimum energy paths for Co and Rh atoms from the
initial state I1 to the final states F1 and F3 along the red and blue paths
indicated in panel a, respectively. For the Co atom the energy is also given
along the path from I2 to F2 (cyan symbols and line). The magnetic moment of
the Co and Rh atom is given in $\mu_{\rm B}$ at every image along the paths.
Extended Data Fig. 1: — Correlation of row-wise AFM state and defect shape.
Intersection of three rotational domains, two of them showing an almost
vanishing spin contrast ($U=-30$ mV, $I=7$ nA). Some defect types display a
broken $C_{3}$ symmetry, correlated with the direction of the
$\uparrow\uparrow$-rows, compare zoom-ins and yellow boomerangs. This
correlation allows to determine the orientation of row-wise AFM domains in
non-magnetic STM measurements. Extended Data Fig. 2: — STM simulations for
Mn/Re(0001) at a bias voltage of $U=-$50 mV. Comparison of simulated STM and
SP-STM images of fcc-Mn/Re(0001) in the row-wise AFM state based on DFT
calculations at a bias voltage of $-$50 meV and a distance of 3 Å above the
surface without and with a lateral shift of the Mn monolayer. a,b Mn atoms are
at perfect fcc positions. c,d Mn atom rows shifted from the perfect fcc sites
into the energetically favorable relaxed positions. Red and blue spheres
represent Mn atoms of opposite spin directions, while yellow spheres denote Re
atoms. Left panels show STM images and right panels show SP-STM images
assuming a tip spin-polarization of 0.5. Extended Data Fig. 3: — STM
simulations for Mn/Re(0001) at a bias voltage of $U=+$50 mV. Comparison of
simulated STM and SP-STM images of fcc-Mn/Re(0001) in the row-wise AFM state
based on DFT calculations at a bias voltage of $+$50 meV and a distance of 3 Å
above the surface without and with a lateral shift of the Mn monolayer. a,b Mn
atoms are at perfect fcc positions. c,d Mn atom rows shifted from the perfect
fcc sites into the energetically favorable relaxed positions. Red and blue
spheres represent Mn atoms of opposite spin directions, while yellow spheres
denote Re atoms. Left panels show STM images and right panels show SP-STM
images assuming a tip spin-polarization of 0.5.