# ACCELERATED INTRAVASCULAR ULTRASOUND IMAGING USING DEEP REINFORCEMENT
LEARNING

###### Abstract

Intravascular ultrasound (IVUS) offers a unique perspective in the treatment
of vascular diseases by creating a sequence of ultrasound-slices acquired from
within the vessel. However, unlike conventional hand-held ultrasound, the thin
catheter only provides room for a small number of physical channels for signal
transfer from a transducer-array at the tip. For continued improvement of
image quality and frame rate, we present the use of deep reinforcement
learning to deal with the current physical information bottleneck. Valuable
inspiration has come from the field of magnetic resonance imaging (MRI), where
learned acquisition schemes have brought significant acceleration in image
acquisition at competing image quality. To efficiently accelerate IVUS
imaging, we propose a framework that utilizes deep reinforcement learning for
an optimal adaptive acquisition policy on a per-frame basis enabled by actor-
critic methods and Gumbel top-$K$ sampling.

Index Terms—  Deep reinforcement learning, intravascular ultrasound,
compressed sensing

##  1 Introduction

Minimally invasive vascular interventions are increasingly guided by phased
array intravascular imaging. Intravascular ultrasound (IVUS) is a catheter-
based imaging modality that enables detailed assessment of the lumen,
endothelium and surrounding tissue of blood vessels. Most importantly, it
enables the verification of the treatment, generally the placement of one or
more stents.

Ultrasound (US) is an affordable, safe and overall convenient medical imaging
tool. However, in general, ultrasound images are notoriously difficult to
interpret, due to the characteristics of current methods of ultrasound image
formation. For IVUS images, this means that the artifacts such as clutter and
speckle strongly affect the appearance of small details and of subtle
variations in texture. A good image quality facilitates the important
distinction between tissues but also e.g. between different types of arterial
plaque. It also makes it easier to verify the proper placement of a stent.

Image quality in phased array IVUS is fundamentally bounded by the information
extracted by the adopted pulse-echo imaging sequence. Subsampling of this
acquisition sequence to boost frame rates while maintaining image quality is
an active area of research in ultrasound in general [1, 2], but its benefits
are especially promising for IVUS as its form-factor highly limits the
bandwidth for the raw sensor data. Careful selection of which measurements to
sample will not only reduce the amount of data, but subsequently decrease
acquisition time and potentially simplify the hardware. Moreover, it allows
extra time for dedicated imaging modes, such as spectral Doppler imaging and
high-resolution elastography, potentially providing additional clinical
information. Furthermore, the reconstruction of the blood vessel is captured
in a sequence of frames by the gradual withdrawal of the transducer-bearing
catheter. These _pullbacks_ are often accompanied by motion artifacts [3].
High-frame rate imaging is therefore desirable as it alleviates these adverse
motion effects and allows for faster pullbacks thereby reducing the overall
procedure time. In a clinical setting, this directly translates to quicker
interventions.

Model-based approaches towards subsampling of ultrasound data focus on either
sub-aperture imaging [4] or sparse array designs such as Cantor arrays [5].
These sequences are hand-engineered, static, and not jointly optimized with
downstream image reconstruction methods. In the past decades, digitization of
raw radio-frequency (RF) channel data and the subsequent beamforming required
custom signal processing hardware, to achieve real-time performance. Today,
the availability of affordable, fast general-purpose high-performance
processors opens the way to implement increasingly sophisticated algorithms
for real-time ultrasound image formation. This revolution in high-performance
computing has spurred the adoption of deep learning (DL) in ultrasound imaging
[6, 7, 8]. Naturally, this trend has also led to the use of compressed sensing
(CS) in ultrasound imaging [1, 2].

Therefore, we propose learning optimal IVUS-specific acquisition-sequencing
using state-of-the-art DL methods in a framework we dub as _AiVUS_. The
framework learns to optimally extract information from the acoustic scene by
adaptive sampling of transducer element pairs. AiVUS poses imaging as an
information-gain problem, where the intelligent system measures precisely what
is needed to provide at any time the optimal trade-off between imaging
properties, disturbing artifacts, and constraints on data rates.

To that end, AiVUS relies on deep reinforcement learning (DRL). Researchers
have employed DRL with great success in a variety of applications, from
magnetic resonance imaging (MRI) [9] and computed tomography (CT) [10] to
radar [11]. Adaptive sampling methods for medical imaging that use DRL have
predominantly been developed for MRI [9, 12]. These frameworks seek to
iteratively optimize the measurements taken for a single reconstruction. This
means that the RL loop is executed for every measurement necessary to
construct the frame. However, we argue this approach does not extrapolate well
to IVUS, as the scan time in MRI is orders of magnitude longer compared to
ultrasound acquisition [13]. In contrast to these methods, AiVUS selects all
measurements needed for a single IVUS reconstruction concurrently, such that
the agent is employed on a per-frame basis. This means that the sampling
strategy of the current frame is conditioned on the reconstruction of the
previous frame, which benefits the faster IVUS acquisition process. Granted,
sufficient temporal correlation between successive frames is assumed. In other
words, the current imaging scene is informative enough for the policy to
provide valid sampling for the next frame.

Other deep learning subsample strategies (not using DRL), such as Deep
Probabilistic Subsampling (DPS) [2] and Active-DPS (A-DPS) [14] leverage the
so-called ‘Gumbel-max trick’ [15], a reparameterization ’trick’ to jointly
optimize a sampling- and reconstruction-network. DPS learns a fixed sample
scheme using a generative sampling model, while A-DPS extends this to an
active acquisition framework that conditions the sampling procedure on already
acquired data. Both DPS and A-DPS require fully differentiable signal
processing pipelines for joint optimization with the sample network. This
inherently limits the scope of applications.

##  2 Methods

###  2.1 IVUS Signal Acquisition

Circular-array intravascular ultrasound imaging utilizes uniformly spaced
transducers to generate a tomographic reconstruction of the blood vessel wall.
Multiple elements eliminate the need for rotating mechanical systems by virtue
of an array processing technique known as beamforming, which aims to maximize
spatial resolution. During beamforming, focused virtual beams are constructed
from unfocused ultrasound pulse-echo signals. The fundamental principle of
delay-and-sum beamforming is to arrange the raw data such that the coherent
signals are enhanced while the incoherent signals are suppressed through
temporal adjustment (time-of-flight correction) and spatial adjustment
(apodization). Due to the limited space in the operational domain, severe
hardware constraints are imposed on the catheter-based IVUS probe, leaving
room for only small circuitry and for only a small number of signal wires.
Therefore, each acquisition, corresponding to a transmitting-receiving
transducer-element pair, is time-multiplexed and sent over the catheter one at
a time. This provokes a trade-off between the image resolution and data
acquisition time, which scales linearly with the amount of measured data.

Because of the circular nature of the array, it is advantageous to only
combine pulse-echo responses of neighboring transducers. Only elements that
significantly contribute to the resulting focused beam are used. For
reconstruction of a single IVUS image, the number of measurements is equal to
$N=E\times A$, where $E$ is the number of ultrasound elements and $A$ is the
sub-aperture size. Each of these $N$ measurements is a vector containing the
fast-time RF samples.

###  2.2 Subsampling

Given the RF channel data $\mathbf{x}\in\mathbb{R}^{N}$, we are looking for an
optimal binary subsampling mask $\mathbf{a}\in\left\\{0,1\right\\}^{N}$ that
selects $K\leq N$ measurements from the total of $N$ possible measurements
(i.e. element pairs), resulting in an $N/K$ acceleration of the acquisition
process. The measured data is given by

$\mathbf{y}=\mathbf{a}\odot\mathbf{x},\vspace{-2mm}$ (1)

where $\odot$ denotes an element-wise multiplication. Given the current
partial acquisition $\mathbf{y}\in\mathbb{R}^{N}$, which is zero-filled at the
non-sampled indices, we would like to design an adaptive subsampling mask
$\mathbf{a}$ that maximizes image quality of the next frame. The
reconstruction operation $g(\cdot)$ encompasses all signal processing
operations for B-mode imaging, after which the IVUS image
$\mathbf{s}=g(\mathbf{y})$ follows. A fixed apodization scheme was used for
the sake of simplicity. This may require further investigation in future work.

###  2.3 Reinforcement Learning for Adaptive Subsampling

The model in AiVUS that performs the adaptive subsampling is trained using
reinforcement learning (RL). An RL agent takes actions upon observation of its
interaction with the environment and can learn to update its behavior
accordingly. Transitions in this perception-action-learning loop are modeled
as a Partially Observable Markov Decision Process (POMDP). Each step $t$ in
the POMDP corresponds to the acquisition and reconstruction of a single IVUS
image $\mathbf{s}$. An action $\mathbf{a}_{t}$ represents a partial
acquisition of IVUS data, according to some transmit scheme. The agent
produces a policy $\pi(\mathbf{a}_{t}|\mathbf{s}_{t})$ over actions given the
current reconstruction $\mathbf{s}_{t}$ and is incentivised by a reward signal
$r_{t}(\mathbf{s}_{t},\mathbf{a}_{t})$ that captures the reconstruction
quality of $\mathbf{s}$ with respect to the ground truth image
$g(\mathbf{x})$. The goal of the agent is to maximize the expected discounted
cumulative reward
$J_{\pi}=\sum_{t=0}^{T-1}\mathbb{E}_{(\mathbf{s}_{t},\mathbf{a}_{t})\sim\pi}\left[\gamma^{t}r(\mathbf{s}_{t},\mathbf{a}_{t})\right]$,
where the episodes of length $T$ are indexed by $t\in\left[0,T-1\right]$. The
discount factor $\gamma$ weighs the reward scores such that more emphasis is
placed on immediate rewards. In order to meet the high dimensional state and
action spaces, an actor-critic framework is leveraged for training the RL
policy [16, 17]. Actor-critic methods extend the use of an action-value
function $Q$, modeled by the critic network, with an explicit definition of
the policy, i.e. the actor network $\pi$ parameterized by $\theta$. The critic
network $Q$ parameterized by $\phi$ can be iteratively updated through
stochastic gradient descent using the Bellman error [16, 17]:

$J_{Q}(\phi)=\mathbb{E}_{(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{s}_{t+1})\sim\mathcal{D}}\left[\big{(}Q_{\phi}(\mathbf{s}_{t},\mathbf{a}_{t})-\hat{Q}(\mathbf{s}_{t},\mathbf{a}_{t})\big{)}^{2}\right],\vspace{-1mm}$
(2)

where $\mathcal{D}$ is the Replay Memory, which stores state transitions
sampled from the environment. The target value $\hat{Q}$ given by

$\hat{Q}(\mathbf{s}_{t},\mathbf{a}_{t})=r(\mathbf{s}_{t},\mathbf{a}_{t})+\gamma
Q_{\phi^{\prime}}(\mathbf{s}_{t+1},\pi_{\theta^{\prime}}(\mathbf{s}_{t+1})),\vspace{-1mm}$
(3)

defines the recursive relation of the action-value function through temporal
difference (TD) learning. The target is computed with separate target networks
parameterized by $\phi^{\prime},\theta^{\prime}$ to improve stability of
learning [16]. The actor network is updated through gradient ascent using the
policy gradient

$\nabla_{\theta}J_{\pi}(\theta)\simeq\mathbb{E}_{\mathbf{s}_{t}\sim\mathcal{D},\mathbf{e}_{t}\sim\text{Gumbel}}\left[\nabla_{\theta}Q_{\phi}(\mathbf{s}_{t},\pi_{\theta}(\mathbf{s}_{t}))\right],\vspace{-1mm}$
(4)

such that the expected return, modeled by $Q_{\phi}(\cdot)$, is maximized.
Current actor-critic methods such as Deep Deterministic Policy Gradient (DDPG)
[16] and Soft Actor-Critic (SAC) [17] do not directly support multiple action
selection. In the IVUS setting, we cannot afford running the network
iteratively to generate multiple actions for a single reconstruction, due to
the high pulse-repetition interval. Multiple actions can be modeled by a
larger underlying unstructured distribution such that the policy only has to
sample a single action. However, because the number of parameters of a
distribution over an $N$-choose-$K$ problem scales factorially, modeling the
underlying categorical distribution using the policy is intractable.
Therefore, AiVUS leverages Gumbel top-$K$ sampling, which allows exact
sampling of $K$ actions simultaneously by a single pass of the policy network
[18]. Gumbel top-$K$ sampling is equivalent to sampling $K$ times without
replacement from the same categorical distribution. In this way, the logits of
the policy network represent only a single categorical distribution of size
$N$ from which we take $K$ actions. Gumbel Top-$K$ sampling is a
generalization of the Gumbel-max trick [15] which perturbs the unnormalized
logits with i.i.d. Gumbel noise samples

$\mathbf{e}_{t}=-\log\left(-\log(U)\right)*\sigma\sim\text{Gumbel}(0,\sigma),\vspace{-1mm}$
(5)

where $U\in[0,1]$ is sampled from a continuous uniform distribution. The
Gumbel samples are sampled with scaling term $\sigma$. The sampling procedure
is computed as follows:

$\mathbf{a}_{t}=\mathrm{topK}(f_{\theta}(\mathbf{s}_{t})+\mathbf{e}_{t}),\vspace{-1mm}$
(6)

where $f_{\theta}$ is a deterministic function of the state whose logits
represent the unnormalized log-probabilities over actions, from which
$\pi_{\theta}$ follows after top-$K$ sampling. The non-differentiable
$\text{argmax}(\cdot)$ operation (used to compute top-$K$) can be relaxed with
reparameterizable subset sampling using $\text{softmax}(\cdot)$ [19], to
permit error backpropagation through $f_{\theta}(\cdot)$, similar as done in
[2]. Intuitively, we can scale the Gumbel noise with $\sigma$ to sample more
uniform actions, which in return ensures exploration. During training we use
an annealing scheme for $\sigma$, starting with a relatively high value
($\sigma=2$), and reducing it until reaching a set minimum value
($\sigma=0.2$). In practice, this annealing approach led to more stable
training compared to a fixed setting of $\sigma=1$. During evaluation, we
disable the Gumbel noise and take the top-$K$ largest values of the policy
distribution. A schematic overview of the full signal path is shown in Fig. 1.

Fig. 1: Block diagram of the perception-action-learning loop of the RL agent
in an IVUS imaging setting.

##  3 Experiments

In this work, a prototype capacitive micromachined ultrasonic transducer
(CMUT) IVUS sensor with a 114 element circular phased-array (Philips Research
[20]) is used to provide the raw RF data for training and verification of the
models. We demonstrate our framework on three different datasets, increasing
in complexity. Namely, simulated wire targets, wire phantoms and _in-vivo_
data from a porcine model. The data is recorded in raw channel RF format and
after beamforming each frame contains 456 scanlines with 520 samples along the
penetration depth of $8\text{\,}\mathrm{m}\mathrm{m}$. We consider episodic
environments with fixed length $T=10$.

Fig. 2: Six successive wire phantom frames constructed using AiVUS. The
agent’s action is displayed in the top row, where the circle represents the
elements in the transducer array. Actions are displayed such that black
represents 1 receiving element and white $A$ receiving elements, for each
transmitting element $\in\left\\{1,\ldots,E\right\\}$.

Fig. 3: SSIM performance on the _in-vivo_ test data for both learned (AiVUS)
and random sampling strategies. Average and standard deviations are computed
over five seeds.

###  3.1 Simulated Wire Target

The first dataset is created by simulating moving wire targets that revolve
around the center at a random distance. We test if AiVUS can learn an adaptive
policy and track the targets. In this specific case, it is very clear what
kind of policy is expected, i.e. sample the channels that originate from
transducers that are in close vicinity to the wire targets. The policy network
is trained on 200 frames and evaluated on 100 test examples. A mean squared
error (MSE) is used as reward:

$\vspace{-2mm}r_{t}=\eta_{MSE}=-||\mathbf{s}_{t}-g(\mathbf{x}_{t}))||_{2}^{2}.\vspace{-2mm}$
(7)

###  3.2 Wire Phantom

To bridge the gap from simulated data to more realistic US data, we subject
the DRL agent to images of wire phantoms. It allows to assess the performance
of the agent in a more realistic environment which encompasses interference
and noise as well as more subtle US features such as reverberation and
ringdown. In total 15 sequences with 8767 frames are recorded, of which 20%
are randomly assigned to the test set. For the reward calculation, (7) is
used, however, the images are filtered using straightforward thresholding to
reduce background noise and interference. We found this step was useful for
constructing a consistent reward signal.

###  3.3 In-vivo data

Although the wire target experiment shows proof of concept of AiVUS,
experiments on _in-vivo_ data are necessary to show the framework has
potential for clinical relevance. The prior acquired dataset consists of eight
different pullbacks from a porcine model with 8679 frames in total of which
20% is randomly assigned as test set. Each pullback is done at approximately
$1\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ for a length of
$30\text{\,}\mathrm{m}\mathrm{m}$. The structural similarity index measure
(SSIM) is used as opposed to the MSE for training the policy network.
Furthermore, we observe that the speckle, which is the granular appearance
inherent to the ultrasound imaging, has a significant effect on the image
quality metric. Different sample schemes lead to small shifts in the speckle
texture which translate to the reward score. Slight changes in the speckle
texture do not positively contribute to the visual appearance as the speckle
pattern appears random. To promote sharp reconstruction of local features, we
add an adversarial loss by adopting a discriminator network $D_{\psi}$. The
resulting reward score is a combination of the SSIM score and discriminator
loss $L_{D_{\psi}}$ :

$r_{t}=\eta_{SSIM}-\lambda_{D}\mathcal{L}_{D_{\psi}},\vspace{-2mm}$ (8)

where $\lambda_{D}$ weights the discriminator loss to the image quality metric
and was empirically set to $110-5$. To further reduce the effect of the
speckle, an anisotropic diffusion filter is applied to both state and ground
truth images before the reward calculation.

Table 1: Quantitative results on the test sets of all three experiments with subsampling factor $N/K=4$, comparing a random agent (I) with a trained agent (II) (AiVUS). Dataset | Sim. | Phantom | In-vivo
---|---|---|---
Agent | I | II | I | II | I | II
MSE $\downarrow$ | 3.25 | 1.42 | 0.067 | 0.047 | 0.078 | 0.070
MAE $\downarrow$ | 0.069 | 0.034 | 0.169 | 0.133 | 0.211 | 0.200
PSNR $\uparrow$ | 44.33 | 49.95 | 61.59 | 64.48 | 59.43 | 59.90
SSIM $\uparrow$ | 0.996 | 0.998 | 0.308 | 0.447 | 0.552 | 0.578

##  4 Results

In Table 1 we report test scores of the final reconstructions using four
common image quality metrics. Namely, the mean-squared error (MSE), mean-
absolute error (MAE), peak signal-to-noise ratio (PSNR) and the structural
similarity index measure (SSIM). The trained DRL agent is compared to an agent
which randomly samples from the transducer element pairs with the same
subsampling factor. AiVUS outperforms a random agent on all four metrics. Fig.
3 shows a sequence of successive frames from the wire phantom dataset
alongside the sampling actions that were used to facilitate the
reconstruction. AiVUS is able to show adaptivity by increasing the aperture
size of specific transducers that are located near the target and tracking the
targets from frame to frame. The agent demonstrates it can exploit knowledge
of the current frame to provide an optimal acquisition sequence for the next
reconstruction. For the _in-vivo_ data, a comparison of the SSIM scores for
different subsampling factors, both for random and learned sampling by AiVUS,
is visualized in Fig. 3. For higher subsampling factors, the learned strategy
is preferred, while for a subsampling factor of two there is almost no
performance difference.

##  5 Conclusion

In this paper, we have proposed a deep reinforcement learning framework that
learns subsampling policies for IVUS on a per-frame-basis using actor critic
methods. The framework allows for any non-differentiable reconstruction method
and quality metric. We have demonstrated that AiVUS is able to navigate in
controlled IVUS environments with high dimensional state and action spaces.
AiVUS outperforms a random agent using a learned acquisition strategy. In
future work, we foresee this framework can be put to use for control of other
ultrasound transmit parameters.

## References

  * [1] Noam Wagner, Yonina C Eldar, Arie Feuer, Gilad Danin, and Zvi Friedman,  “Xampling in ultrasound imaging,”  in Medical Imaging 2011: Ultrasonic Imaging, Tomography, and Therapy. International Society for Optics and Photonics, 2011, vol. 7968, p. 796818\.
  * [2] Iris AM Huijben, Bastiaan S Veeling, Kees Janse, Massimo Mischi, and Ruud JG van Sloun,  “Learning sub-sampling and signal recovery with applications in ultrasound imaging,”  IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 3955–3966, 2020.
  * [3] Nima Torbati, Ahmad Ayatollahi, and Parham Sadeghipour,  “Image-based gating of intravascular ultrasound sequences using the phase information of dual-tree complex wavelet transform coefficients,”  IEEE transactions on medical imaging, vol. 38, no. 12, pp. 2785–2795, 2019.
  * [4] Mustafa Karaman and Matthew O’Donnell,  “Subaperture processing for ultrasonic imaging,”  IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 45, no. 1, pp. 126–135, 1998.
  * [5] Chun-Lin Liu and PP Vaidyanathan,  “Maximally economic sparse arrays and cantor arrays,”  in 2017 IEEE 7th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP). IEEE, 2017, pp. 1–5.
  * [6] Ruud JG Van Sloun, Regev Cohen, and Yonina C Eldar,  “Deep learning in ultrasound imaging,”  Proceedings of the IEEE, vol. 108, no. 1, pp. 11–29, 2019.
  * [7] Ben Luijten, Regev Cohen, Frederik J de Bruijn, Harold AW Schmeitz, Massimo Mischi, Yonina C Eldar, and Ruud JG van Sloun,  “Adaptive ultrasound beamforming using deep learning,”  IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 3967–3978, 2020.
  * [8] Hannes Hase, Mohammad Farid Azampour, Maria Tirindelli, Magdalini Paschali, Walter Simson, Emad Fatemizadeh, and Nassir Navab,  “Ultrasound-guided robotic navigation with deep reinforcement learning,”  arXiv preprint arXiv:2003.13321, 2020.
  * [9] Tim Bakker, Herke van Hoof, and Max Welling,  “Experimental design for mri by greedy policy search,”  Advances in Neural Information Processing Systems, vol. 33, 2020\.
  * [10] Florin-Cristian Ghesu, Bogdan Georgescu, Yefeng Zheng, Sasa Grbic, Andreas Maier, Joachim Hornegger, and Dorin Comaniciu,  “Multi-scale deep reinforcement learning for real-time 3d-landmark detection in ct scans,”  IEEE transactions on pattern analysis and machine intelligence, vol. 41, no. 1, pp. 176–189, 2017.
  * [11] Tristan SW Stevens, R Firat Tigrek, Eric S Tammam, and Ruud JG van Sloun,  “Automated gain control through deep reinforcement learning for downstream radar object detection,”  arXiv preprint arXiv:2107.03792, 2021.
  * [12] Kyong Hwan Jin, Michael Unser, and Kwang Moo Yi,  “Self-supervised deep active accelerated mri,”  arXiv preprint arXiv:1901.04547, 2019.
  * [13] Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, et al.,  “fastmri: An open dataset and benchmarks for accelerated mri,”  arXiv preprint arXiv:1811.08839, 2018.
  * [14] Hans Van Gorp, Iris Huijben, Bastiaan S Veeling, Nicola Pezzotti, and Ruud JG Van Sloun,  “Active deep probabilistic subsampling,”  in International Conference on Machine Learning. PMLR, 2021, pp. 10509–10518.
  * [15] Chris J Maddison, Daniel Tarlow, and Tom Minka,  “A* sampling,”  arXiv preprint arXiv:1411.0030, 2014.
  * [16] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra,  “Continuous control with deep reinforcement learning,”  arXiv preprint arXiv:1509.02971, 2015.
  * [17] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine,  “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,”  in International Conference on Machine Learning. PMLR, 2018, pp. 1861–1870.
  * [18] Wouter Kool, Herke Van Hoof, and Max Welling,  “Stochastic beams and where to find them: The gumbel-top-k trick for sampling sequences without replacement,”  in International Conference on Machine Learning. PMLR, 2019, pp. 3499–3508.
  * [19] Sang Michael Xie and Stefano Ermon,  “Reparameterizable subset sampling via continuous relaxations,”  arXiv preprint arXiv:1901.10517, 2019.
  * [20] Reza Pakdaman Zangabad, Johan G Bosch, Frits Mastik, Robert HSH Beurskens, Vincent A Henneken, Johannes W Weekamp, Antonius FW van der Steen, and Gijs van Soest,  “Real-time coded excitation imaging using a cmut-based side looking array for intravascular ultrasound,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 68, no. 6, pp. 2048–2058, 2021.


# Computational Approach of Designing Magnetfree Nonreciprocal Metamaterial

Swadesh Poddar*, 1 Md. Tanvir Hasan  and Md. Ragib Shakil Rafi

###### Abstract

This article aims at discussing computational approach to design magnet-free
nonreciprocal metamaterial. Detailed mathematical derivation on floquet mode
analysis is presented for Faraday and Kerr rotation. Non-reciprocity in the
designed metasurface is achieved in the presence of biased transistor loaded
in the gap of circular ring resonator. Based on the derived mathematical
model, co and cross-polarized components have been extracted, which helps find
Faraday and Kerr rotation and compare/contrast the reciprocal and
nonreciprocal systems.

Vol. x, y–z, 2021

111Received date222* Corresponding author: Swadesh Poddar
(poddarswadesh@gmail.com).3331 Design Engineer, Qorvo, US.

Md. Tanvir Hasan and Md. Ragib Shakil Rafi

##  1 Introduction

A system is called reciprocal when signal transmission between any two ports
is independent with respect to propagation direction. Signal transmission
between transmitter and receiver can be represented in terms of scattering
($S$) as $S^{T}=S$). In our daily life application such as antenna, passive
electrical circuits, and components are reciprocal. Non-reciprocal systems,
the opposite of reciprocal, exhibit diﬀerent received-transmitted field ratios
when the source and observation points interchange. Non-reciprocal components
such as circulators, isolators, and gyrators play a pivotal role in many
microwave wireless applications, including high-power transmitters, FDD
system, emerging quantum computing and readout, defense and satellite
communication, to name a few [8, 1].

Since 1950, there are significant amount of researches on non-reciprocity and
over time different topologies have been developed to achieve nonreciprocity
in the presence of permanent magnet, using nonreciprocal component, time
modulation, parametric amplifier and so on [8, 14, 7, 19, 2]. Advances in
solid state technology and significant research effort on artificially
engineered material, known as metamaterial is to induce customized properties
in a material that originally does not exist, have initiated a prime interest
on achieving “magnet-free” non-reciprocity, where the time reversal symmetry
can be broken in absence of static magnets. Over the past decade, research
into metamaterials has been extended to a search for real-world applications,
leading to the concept of meta-devices, defined as metamaterial-based devices
that can operate in an active manner [20]. This article will mostly aim on the
computational approach of designing nonreciprocal metamaterial leading to
practical realization of co-simulation and machine learning based
optimization. The proposed design in this article is the extension of the
transistor based design first proposed by Toshiro et al. [6, 3] and later in
our previous articles a nonreciprocal metasurface was investigated with a
precise electronic tunability features [14, 12] as well as analytical
modelling of polarization rotation in a nonreciprocal media. [15].

The article is organized as follows. In Section II, an insight to the
mechanism of floquet mode analysis and corresponding analytical approach will
be provided, in Section III, the principle of ferromagnetic resonance, and in
section IV, the design approaches and process flow of our nonreciprocal
devices will be described. At last in section V our discussion will be mostly
focused on simulated performance, benchmarking reciprocal and nonreciprocal
design, highlighting major challenges, limitation, and future works.

##  2 Floquet Mode Analysis and Mathematical Modeling

In this section, based on floquet mode analysis, a mathematical model is
presented to characterize the unit cell of our designed metasurface. The
Floquet port is used with planar-periodic structures. The analysis of the
periodic structure is accomplished by analyzing a unit cell and master and
slave boundaries are used to model a unit cell of a repeating structure. The
fields on the slave surface are constrained to be identical to those on the
master surface, with phase shift. A coordinate system must be identified on
the master and slave boundary to identify point-to-point correspondence and
master and slave surfaces must be of identical shapes and sizes. A set of
modes called floquet modes represent the fields on the port boundary.
Fundamentally, Floquet modes can be assumed as plane waves where propagation
direction is set by the frequency, and geometry of periodic structures. To
simplify the discussion, a 3D structure has been developed and 1st 2 modes
have been taken into consideration to propagate through the metamaterial unit
cell. Floquet boundary at z direction has been calculated accordingly.

Figure 1: Geometry of proposed metamaterial unit cell. Primary and secondary
master and slave boundaries form the side walls of unit cell. In both ports
floquet environment scattered and incident wave is shown with arrow direction

Figure 1 shows the proposed metamaterial unit cell in the floquet port
environment. $a_{mn}$ and $b_{mn}$ are the incident wave from the floquet port
1 and floquet port 2, respectively. On the other hand $c_{mn}$ and $d_{mn}$
are the scattered wave. In a normal incidence case where the 1st 2 modes,
namely TE and TM is evaluated and can be expressed in S matrix. In Ansys HFSS
simulation, floquet mode has been analyzed where transverse electric (TE)
modes are y polarized and transverse magnetic (TM) modes are x polarized.

$\left[\begin{array}[]{c}\mathbf{c}_{mn}^{TE}\\\ \mathbf{c}_{mn}^{TM}\\\
\mathbf{d}_{mn}^{TE}\\\
\mathbf{d}_{mn}^{TM}\end{array}\right]=\left[\begin{array}[]{cccc}\mathbf{S}_{11}^{TE,TE}&\mathbf{S}_{11}^{TE,TM}&\mathbf{S}_{12}^{TE,TE}&\mathbf{S}_{12}^{TE,TM}\\\
\mathbf{S}_{11}^{TM,TE}&\mathbf{S}_{11}^{TM,TM}&\mathbf{S}_{12}^{TM,TE}&\mathbf{S}_{12}^{TM,TM}\\\
\mathbf{S}_{21}^{TE,TE}&\mathbf{S}_{21}^{TE,TM}&\mathbf{S}_{22}^{TE,TE}&\mathbf{S}_{22}^{TE,TM}\\\
\mathbf{S}_{21}^{TM,TE}&\mathbf{S}_{21}^{TM,TM}&\mathbf{S}_{22}^{TM,TE}&\mathbf{S}_{22}^{TM,TM}\end{array}\right]\left[\begin{array}[]{c}\mathbf{a}_{mn}^{TE}\\\
\mathbf{a}_{mn}^{TM}\\\ \mathbf{b}_{mn}^{TE}\\\
\mathbf{b}_{mn}^{TM}\end{array}\right],$ (1)

where, $\mathbf{c}_{mn}^{TE}$ and $\mathbf{a}_{mn}^{TE}$ defines the scattered
and incident field respectively. Besides, $\mathbf{S}_{ij}^{pol1,pol2}$,
Where, i (pol1) and j (pol2) represents incident and received port
(polarization) respectively. For example, $\mathbf{S}_{11}^{TE,TE}$ is the co
polarized reflection case or in other words incident wave is TE(y) polarized
from port 1 and received wave TE(y) polarized to port 1. Similarly,
$\mathbf{S}_{21}^{TM,TE}$ defines a cross polarization transmission case. The
above equation 1 can be written in a short form as below in terms of floquet
modes.

$\left[\begin{array}[]{c}\mathbf{c}_{mn}^{TE}\\\ \mathbf{c}_{mn}^{TM}\\\
\mathbf{d}_{mn}^{TE}\\\
\mathbf{d}_{mn}^{TM}\end{array}\right]=[\mathbf{S}]\cdot\left[\begin{array}[]{c}\mathbf{a}_{mn}^{TE}\\\
\mathbf{a}_{mn}^{TM}\\\ \mathbf{b}_{mn}^{TE}\\\
\mathbf{b}_{mn}^{TM}\end{array}\right],$ (2)

where

$\mathbf{S}=\left[\begin{array}[]{cccc}S(\text{FP}1:1,\text{FP}1:1)&S(\text{FP}1:1,\text{FP}1:2)&S(\text{FP}1:1,\text{FP}2:1)&S(\text{FP}1:1,\text{FP}2:2)\\\
S(\text{FP}1:2,\text{FP}1:1)&S(\text{FP}1:2,\text{FP}1:2)&S(\text{FP}1:2,\text{FP}2:1)&S(\text{FP}1:2,\text{FP}2:2)\\\
S(\text{FP}2:1,\text{FP}1:1)&S(\text{FP}2:1,\text{FP}1:2)&S(\text{FP}2:1,\text{FP}2:1)&S(\text{FP}2:1,\text{FP}2:2)\\\
S(\text{FP}2:2,\text{FP}1:1)&S(\text{FP}2:2,\text{FP}1:2)&S(\text{FP}2:2,\text{FP}2:1)&S(\text{FP}2:2,\text{FP}2:2)\end{array}\right]$
(3)

where FP is the short form of floquet port used here in eq. 3. In summary, FP
$p:u$, FP $q:u$ means $S_{pq}$ for $u$ polarized as incident and received wave
define co polarization. Similarly, FP $p:v$, FP $q:u$ means $S_{pq}$ for
incident mode as u and received wave as v polarized. In this explanation, $p$
and $q$ were used to mention port number and $v$ and $u$ to mention
polarization type.

$\left[\begin{array}[]{c}\mathbf{c}_{mn}^{TE}\\\ \mathbf{c}_{mn}^{TM}\\\
\mathbf{d}_{mn}^{TE}\\\
\mathbf{d}_{mn}^{TM}\end{array}\right]=\left[\begin{array}[]{cccc}\mathbf{R}_{yy}^{(1)}&\mathbf{R}_{yx}^{(1)}&\mathbf{T}_{yy}^{(2)}&\mathbf{T}_{yx}^{(2)}\\\
\mathbf{R}_{xy}^{(1)}&\mathbf{R}_{xx}^{(1)}&\mathbf{T}_{xy}^{(2)}&\mathbf{T}_{xx}^{(2)}\\\
\mathbf{T}_{yy}^{(1)}&\mathbf{T}_{yx}^{(1)}&\mathbf{R}_{yy}^{(2)}&\mathbf{R}_{yx}^{(2)}\\\
\mathbf{T}_{xy}^{(1)}&\mathbf{T}_{xx}^{(1)}&\mathbf{R}_{xy}^{(2)}&\mathbf{R}_{xx}^{(2)}\end{array}\right]\left[\begin{array}[]{c}\mathbf{a}_{mn}^{TE}\\\
\mathbf{a}_{mn}^{TM}\\\ \mathbf{b}_{mn}^{TE}\\\
\mathbf{b}_{mn}^{TM}\end{array}\right]$ (4)

Equation 4 has been represented in terms of reflection and transmission for
the easier understanding of the floquet mode analysis. where for R/T matrix,
superscript describes the received port number and subscript describes the
polarization. e.g. $T_{xx}^{1}$ equals the wave that is being propagated from
port 2 and received at port 1, for generalized understanding, this can be
mapped with $S_{12}$. In general the reflected and transmitted fields consist
of both x- and y- polarized components. Transmission field $T_{xy}$ can be
defined as the ratio of the transmitted field with polarization $x$ to the
incident field with polarization $y$ Based on same concept, $R_{xy}$ can be
defined as the ratio of the reflected field with polarization $x$ to the
incident field with polarization $y$.

The incident and transmitted fields are related through Jones transmission
matrix in cartesian basis as below [11, 15]

$\mathbf{T}=\left[\begin{array}[]{cc}\mathbf{T}_{xx}&\mathbf{T}_{xy}\\\
\mathbf{T}_{yx}&\mathbf{T}_{yy}\end{array}\right]$ (5)

The transmission matrix in linear basis can be mapped to transmission matrix
in circular basis as below

$\mathbf{T}_{CP}=\left[\begin{array}[]{cc}\mathbf{T}_{\circlearrowright\circlearrowright}&\mathbf{T}_{\circlearrowright\circlearrowleft}\\\
\mathbf{T}_{\circlearrowleft\circlearrowright}&\mathbf{T}_{\circlearrowleft\circlearrowleft}\end{array}\right]=\frac{1}{2}\left[\mathbf{L}\right]$
(6)

$\mathbf{T}_{CP}=\left[\begin{array}[]{cc}\mathbf{T}_{\circlearrowright\circlearrowright}&\mathbf{T}_{\circlearrowright\circlearrowleft}\\\
\mathbf{T}_{\circlearrowleft\circlearrowright}&\mathbf{T}_{\circlearrowleft\circlearrowleft}\end{array}\right]=\frac{1}{2}\left[\mathbf{L}\right]=\frac{1}{2}\left[\begin{array}[]{cc}T_{xx}+T_{yy}+i(T_{xy}-T_{yx})&T_{xx}-T_{yy}-i(T_{xy}+T_{yx})\\\
T_{xx}-T_{yy}+i(T_{xy}\text{ }+\text{
}T_{yx})&T_{xx}+T_{yy}-i(T_{xy}-T_{yx})\end{array}\right]$ (7)

Based on our above discussion, $\mathbf{T}$ matrix in circular basis can be
generalized as below

$\displaystyle\mathbf{T}_{\circlearrowright\circlearrowright}$
$\displaystyle=$ Incident RHCP to recieved RHCP (8)
$\displaystyle\mathbf{T}_{\circlearrowright\circlearrowleft}$ $\displaystyle=$
Incident LHCP to recieved RHCP (9)
$\displaystyle\mathbf{T}_{\circlearrowleft\circlearrowright}$ $\displaystyle=$
Incident RHCP to recieved LHCP (10)
$\displaystyle\mathbf{T}_{\circlearrowleft\circlearrowleft}$ $\displaystyle=$
Incident LHCP to recieved LHCP (11)

where RHCP/LHCP stands for right/left handed circular polarization. Once the
transmission matrix in circular or linear basis is achieved, Faraday rotation
and ellipticity can be calculated from there. Ellipticity, $\delta_{F}$ is a
measure of the polarization state of an electromagnetic wave and it varies
from $+1$ to $-1.$In terms of RHCP and LHCP, this can be written as

$\delta_{F}=\frac{\mid T_{\circlearrowright\circlearrowright}\mid-\mid
T_{\circlearrowleft\circlearrowleft}\mid}{\mid
T_{\circlearrowright\circlearrowright}\mid+\mid
T_{\circlearrowleft\circlearrowleft}\mid}$ (12)

Therefore, from above equation, if $\delta_{F}=+1$, wave is RHCP and
$\delta_{F}=-1$, the wave is LHCP, whereas, $\delta_{F}=0$ means Linear
polarized wave wave.

Faraday rotation can be calculated from the equation below [10, 14]

$\theta_{F}=\frac{1}{2}\tan^{-1}(\frac{T_{\circlearrowleft\circlearrowleft}}{T_{\circlearrowright\circlearrowright}})$
(13)

Similarly, for ellipticity for a reflection based metasurface and Kerr
rotation can be calculated similarly following the same steps as described
above.

$\delta_{K}=\frac{\mid R_{\circlearrowright\circlearrowright}\mid-\mid
R_{\circlearrowleft\circlearrowleft}\mid}{\mid
R_{\circlearrowright\circlearrowright}\mid+\mid
R_{\circlearrowleft\circlearrowleft}\mid}$ (14)

$\theta_{K}=\frac{1}{2}\tan^{-1}(\frac{R_{\circlearrowleft\circlearrowleft}}{R_{\circlearrowright\circlearrowright}})$
(15)

##  3 Ferromagnetic Resonance

The concept related to the motion of the magnetic dipole in the presence of
the magnetic field and a self-consistent RF magnetic field represents the
principle idea of ferromagnetic resonance [10]. Microwave magnetism in a
ferrite material is based on the precession of the magnetic dipole moments
which arise from unpaired electron spins along the axis of an externally
applied static magnetic bias. For resonance condition, it is necessary to
apply the RF magnetic field (the microwave field) perpendicular to the dc
field. The quantum mechanical phenomena can be described from Landau-Lifshitz-
Gilbert equation as [10, 9, 12]

$\frac{d\mathbf{M}}{dt}=\gamma(\mathbf{M}\times\mathbf{H})-\frac{\alpha}{\left|M\right|}\mathbf{M}\times\frac{d\mathbf{M}}{dt},$
(16)

where $\alpha$ is the Gilbert damping term, $\gamma$ is the gyromagnetic
ratio, $\mathbf{M}$ = $\textbf{i}_{z}$${M}_{0}$ \+ m$e^{j\omega t}$, and
$\mathbf{H}$ = $\textbf{i}_{z}$${H}_{i}$ \+ h$e^{j\omega t}$, where
$\textbf{i}_{z}$ is the unit vector in the z direction, m, and h are RF
quantities, and ${H}_{i}$ is the total internal dc magnetic field as shown in
Fig. 2. The nearer the frequency of the microwave field to the natural
precession frequency, the greater will be the energy absorbed by the spins
[10]. The longitudinal (z) component does not contribute to precession and
therefore, to magnetism. Since, $\textbf{H}_{0}$ $\parallel$
$\mathbf{\widehat{z}}$, the z-component of M produced by $\textbf{H}_{t}^{RF}$
would lead to
$[M_{z}^{RF}(\mu_{0}{H}+H^{RF})](\widehat{\mathbf{z}}\times\widehat{{\mathbf{z}}})=0.$
Therefore, the only torque being produced by the transverse component at x-y
plane can be stated as
$(M_{t}^{RF}\mu_{0}{H})(\widehat{\mathbf{t}}\times\widehat{\mathbf{z}})\neq
0.$

Figure 2: Magnetic dipole precession, arising from electron spinning in a
ferrite material about the z axis of an externally applied static magnetic
bias field, $H_{0}$, with effective unidirectional current loops around z axis
and magnetic dipole moment $m_{\rho}$

.

The magnetization vector described above has a singularity, when $\omega$ =
$\omega_{0}$ = $\gamma H_{0}$, which can be defined as the resonance
condition. Resonance can be achieved by varying the operating frequency or the
applied field until the precession frequency equals the microwave frequency.
Permeability tensor of ferrite material can be represented as [16]

$\overrightarrow{\mathbf{\mu}}=\mu_{0}\begin{bmatrix}\mu&-jk&0\\\ jk&\mu&0\\\
0&0&1\end{bmatrix},$ (17)

where $\mu$ = $1$ +$\frac{\omega_{0}\omega_{M}}{\omega_{0}^{2}-\omega^{2}}$,
$jk$ = $\frac{j\omega\omega_{M}}{\omega_{0}^{2}-\omega^{2}}$, and $\omega_{M}$
= $\gamma 4\pi M_{0}$.

##  4 Nonreciprocal Device Design

Active devices have been used in the design of magnet-free non-reciprocal
metamaterials. For example, a transistor amplifier can be viewed as an
isolator that amplifies the input signal in one direction and blocks it in the
reverse direction and very popular in RF microwave devices, metamaterial and
antenna design [18, 17]. In a series of publications between 2011 and 2019,
Kodera et al. explored the unidirectional properties of transistors in the
design of non-reciprocal metamaterials [7, 6, 4, 5]. A ring resonator supports
two counter propagating modes with the same resonance frequency, a result of
reciprocity. However, when unidirectional component such as transistor is used
in the ring gap, one of the modes is blocked and the ring loaded with
transistor resonates only for one type of circularly polarized waves [8, 14].

Figure 3: Non-reciprocal gyrotropy. Rotating dipole moment in (a) the proposed
metamaterial structure consisting of slot-ring resonator loaded with active
components (transistor). (b) travelling wave along the ring at four quarter-
period spaced time instants. (c) Rotating electric field in the slot-ring and
the corresponding electric dipole moment due to vectorial magnetic field
pointed outward and inward. (d) Dimensions of the ring resonator (e) Unit
super-cell of a 3 by 3 super-cell board consisting of 4 rings in a
$90^{0}$-symmetric configuration (Top view). (f) Bias connection layout of the
complete board (Bottom view) [14, 12]

Figure. 3 shows travelling wave resonance of a slot ring resonator loaded with
transistor and passive components mimicking the electron spin precession in a
biased ferromagnetic material. The electrical length ($\lambda$) of the ring
resonator should be $2\pi$ due to phase matching. The wave is traveling (as
opposed to standing wave) due to the presence of the unilateral (transistor)
element. Based on the applied electric bias to the active component, the
electric dipole moment (red arrow) rotates along the circular ring at four
quarter period spaced time instants. The direction of the electric field
follows the direction of the dipole moment. Therefore, non-reciprocal
gyrotropy is established in the slot-ring resonator. The circular ring
resonates when the electrical length is equal to an integer multiple of
$2\pi$.

$\beta l+\alpha_{FET}=2m\pi,$ (18)

where $\beta$ is the propagation constant defined by
$\frac{2\pi}{\lambda_{g}}$, $\alpha_{FET}$ is the physical distance of the
ring gap as shown in Fig. 3 (d), and l is the physical length of the slotted
transmission line defined as $l$ = $(2\pi-\alpha_{FET})r_{mean}$. In this
research work, FR4 ($\varepsilon_{r}$ = 4.4) substrate has been used with a
thickness of 0.8mm. The major focus in this research work is the Faraday
rotation, an important key performance indicator to quantify nonreciprocity,
which is obtained based on the transmission properties. Ground plane have been
designed on geometry of the same front side so the wave can propagate through
the metasurface itself. Details design modelling has been discussed in our
previous work [14].

The effective dielectric constant for a microstrip geometry is

$\varepsilon_{e}=\frac{\varepsilon_{r}+1}{2}+\frac{\varepsilon_{r}-1}{2}\left\\{1+12\left(\frac{H}{W}\right)\right\\}^{-\frac{1}{2}},$
(19)

where Eq. 19 follows the condition $\frac{W}{H}\geq 1$. and $Z_{0}$ can be
calculated as

$Z_{0}=\frac{120\pi}{\sqrt{\varepsilon_{e}}\left[\frac{W}{H}+1.393+\frac{2}{3}\ln\left(\frac{W}{H}+1.444\right)\right]}.$
(20)

Based on the above analysis, initial parameters have been calculated
considering a resonance frequency at 5.7 GHz and dielectric properties of the
FR4 substrate.

Figure 4: Simulated performance of ring resonator based metamaterial. (a)
Simple circular ring resonator based unit cell on proposed metasurface design.
(b) The design is reciprocal, hence, no nonreciprocal properties can be
observed. (c) Both co pol and (d) cross pol components exhibit the reciprocal
behaviour, resonance is observed at 5.69 GHz in this design.

Figure 4 (a) shows the simulated performance of a ring resonator in a
reciprocal condition i.e. without any active component connected. Based on our
previous analysis, nonreciprocal properties can’t be achieved as shown in 4
(b), hence no Faraday or Kerr rotation. Propagating wave exhibits a linear
polarization which can be shown from ellipticity. Figure 4 (c) and (d) shows
reflection and transmission behavior in terms of scattering parameter. One
important observation to highlight for reciprocal device is that the co and
cross pol components will be in line with respect to each other. For example,
we observe no difference between $S_{11}^{xy}$ and $S_{11}^{yx}$ which is
cross pol reflection parameter. Co and cross pol transmission parameter
follows the similar pattern.

Figure 5: Electromagnetic co-simulation (a) Circuit configuration that was
used to achieve nonreciprocity, nonreciprocal component such as transistor was
used to bridge the ring resonator gap, other component are for impedance
matching and biasing the transistor in various condition (b) 3D
electromagnetic model of the metamaterial using floquet mode analysis. (c) S
parameter file of the 3D model and circuit model has been connected using
dynamic link and to evaluate co-simulation results.

Figure 5 shows the setup of the co-simulation. The 3D electromagnetic
structure of the metasurface have been solved in Ansys HFSS which is a passive
solver. Besides, the circuit model has been solved independently in a specific
bias condition. Later the S parameters are connected using the dynamic link to
run the co-simulation and optimization.

Figure 6: Process flow of the computational approach of the nonreciprocal
device design

Figure 6 shows the process-flow of design and optimization. At first step, the
frequency of operation along with preferred dielectric material is selected.
After that, the initial dimension has been calculated based on microwave
equation discussed earlier in the section. All dimensions are used as
variables for later optimization and sensitivity analysis. Lumped ports are
used in the gap and floquet based S matrix have been solved for the
metasurface. Besides, in a circuit solver S parameter from our designed
circuit and solved S parameter of the 3D model are cascaded for the co-
simulation. Once the initial design and result are obtained, based on the
targeted goal and machine learning optimization algorithm design is optimized
to achieve nonreciprocity at desired resonance. In this design, Quasi newton
algorithm was used from HFSS. An in house script was developed and used during
the optimization to calculate and optimize the performance at resonance on Key
performance indicators such as Faraday, Kerr rotation, co-cross polarized
components.

##  5 Results and Discussion

In this section the simulation results of the nonreciprocal metasurface will
be discussed. From antenna 1 to antenna 2, the x-polarized wave is transmitted
from region 1 and y-polarized probe (antenna) pick the signal at region 2 and
if the time is reversed, from antenna 2 to antenna 1, the y-polarized wave is
transmitted and x-polarized probe (antenna) pick the signal at region 1. For
x-polarized wave and x-polarized probing, phase of the transmission
coefficients are same when the time is reversed. Therefore, likewise
magnetized ferrites, the co-polarized component does not show any change in
their phase when the time is reversed, however, the cross-polarized components
shows $180^{\circ}$ phase difference at the resonance frequency so that the
designed structure behaves like a magnetized ferrite where the permeability
tensor has cross-components and they are out-of phase such that the non-
reciprocal gyrotropic response is evident. In Fig. 7 (a) the nonreciprocal
behaviour is evident where we see the magnitude of the cross polarized
component is different and at resonance magnitude of co polarized component is
negligible. Almost $90^{0}$ Faraday rotation at resonance has been achieved.
Similarly, from fig 7 (b), at resonance $\delta_{F}$ = 0, i.e. the wave is
linear polarized. Figure 8 shows the Kerr rotation which is a reflection based
measurements and can be explained similarly as we discussed for the Faraday
rotation. The practical implementation and measurement procedures have been
discussed in detail in our previous articles [14, 13, 15]

Figure 7: (a) Faraday rotation and scattering parameter over frequency. co
polarized components are same and at resonance those are almost 0. Cross
polarized components are different and we observe maximum isolation at
resonance. (b) Faraday rotation and ellipticity over frequency. Figure 8: (a)
Kerr rotation and scattering parameter over frequency. co polarized components
are same and at resonance those are almost 0. (b) Faraday rotation and
ellipticity over frequency.

Design steps of a magnet-less nonreciprocal devices such as isolators,
circulators, and leaky wave antennas are same however, design and optimization
to achieve meaningful results are very challenging. For example, Impedance and
phase matching between circular ring resonator and active/passive component at
microwave frequencies are the major challenges to create artificial non-
reciprocity. The excited resonance is accompanied by loss in the presence of
mismatch leading to nonzero standing wave ratio (SWR). For an infinite
standing wave ratio, the device does not exhibit gyrotropy, hence, non-
reciprocity. Therefore, full wave 3D simulation along with optimization cost
functions were used in this work and highly recommended for any non-reciprocal
device design. Lastly, the number of active/passive components required to
implement the artificial gyrotropy can be optimized and reduced. Characterize
nonreciprocal devices are very susceptible where precise bias condition, setup
test bench and careful extraction of data is very important to take into
consideration.

##  6 Conclusion

To summarize, a detailed computational approach has been presented on a fully
electronic non-reciprocal metamaterial. Necessary condition to achieve
nonreciprocity, analytical modelling along with key steps required for the
design has been proposed, discussed and implemented. With the precision level
of engineering, magnet free non-reciprocal devices can be a probable
substitute of conventional ferrite material biased by permanent magnet in
industrial applications such as isolator, circulator, gyrator, leaky wave
antenna, and 5G applications in future. We envision that the computational
approaches in this work will be very helpful in the existing research area of
magnet-free non-reciprocal metamaterials.

## References

  * [1] Christophe Caloz, Andrea Alù, Sergei Tretyakov, Dimitrios Sounas, Karim Achouri, and Zoé-Lise Deck-Léger.  Electromagnetic nonreciprocity.  Physical Review Applied, 10(4), Oct 2018.
  * [2] Alexander M. Holmes, Mohsen Sabbaghi, Swadesh Poddar, Samane Pakniyat, and George W. Hanson.  Experimental realization of topologically protected surface magnon polaritons on ceramic yig ferrites.  In 2021 International Conference on Electromagnetics in Advanced Applications (ICEAA), pages 204–204, 2021.
  * [3] T. Kodera and C. Caloz.  Unidirectional loop metamaterials (ulm) as magnetless artificial ferrimagnetic materials: Principles and applications.  IEEE Antennas and Wireless Propagation Letters, 17(11):1943–1947, 2018.
  * [4] T. Kodera, D. Sounas, and C. Caloz.  Magnetless nonreciprocal metamaterial (mnm) technology: Application to microwave components.  IEEE Transactions on Microwave Theory and Techniques, 61:1030–1042, 2013.
  * [5] T. Kodera, D. L. Sounas, and C. Caloz.  Switchable magnetless nonreciprocal metamaterial (mnm) and its application to a switchable faraday rotation metasurface.  IEEE Antennas and Wireless Propagation Letters, 11:1454–1457, 2012\.
  * [6] Toshiro Kodera and Christophe Caloz.  Unidirectional Loop Metamaterials (ULM) as Magnetless Artificial Ferrimagnetic Materials: Principles and Applications.  IEEE Antennas and Wireless Propagation Letters, 17(11):1943–1947, November 2018.
  * [7] Toshiro Kodera, Dimitrios L. Sounas, and Christophe Caloz.  Artificial faraday rotation using a ring metamaterial structure without static magnetic field.  Applied Physics Letters, 99(3):031114, 2011.
  * [8] A. Kord, D. L. Sounas, and A. Alù.  Microwave nonreciprocity.  Proceedings of the IEEE, 108(10):1728–1758, 2020.
  * [9] L. Landau and E. Lifshitz.  On the theory of the dispersion of magnetic permeability in ferromagnetic bodies.  Phys. Z. Sowjetunion, 8, 01 1992.
  * [10] Benjamin Lax, Kenneth Button, and H. Hagger.  Microwave ferrites and ferrimagnetics.  Physics Today - PHYS TODAY, 16, 01 1963.
  * [11] Christoph Menzel, Carsten Rockstuhl, and Falk Lederer.  Advanced jones calculus for the classification of periodic metamaterials.  Phys. Rev. A, 82:053811, Nov 2010.
  * [12] Swadesh Poddar.  Design and analysis of fully-electronic magnet-free non-reciprocal metamaterial.  Theses and Dissertations. 2578, 2020.  https://dc.uwm.edu/etd/2578.
  * [13] Swadesh Poddar, Alexander M. Holmes, and George W. Hanson.  Automatic measurement technique of electromagnetic rotation in a nonreciprocal medium, 2022.
  * [14] Swadesh Poddar, Alexander M. Holmes, and George W. Hanson.  Design and analysis of an electronically tunable magnet-free non-reciprocal metamaterial.  IEEE Transactions on Antennas and Propagation, pages 1–1, in press, 2022.
  * [15] Swadesh Poddar, Ragib Shakil Rafi, and Md. Tanvir Hasan.  Analytical modeling of electromagnetic rotation in nonreciprocal media, 2022.
  * [16] D. Polder.  Viii. on the theory of ferromagnetic resonance.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 40(300):99–115, 1949.
  * [17] Dominik Gerhard Rabus and Cinzia Sada.  Ring Resonators: Theory and Modeling, pages 3–46.  Springer International Publishing, Cham, 2020.
  * [18] Sunanda Roy, M. Abdus Samad, and Swadesh Podder.  Effect of complementary triangular split ring resonator on microstrip patch antenna.  In 2015 2nd International Conference on Electrical Information and Communication Technologies (EICT), pages 353–358, 2015.
  * [19] Freek Ruesink, Mohammad-Ali Miri, Andrea Alù, and Ewold Verhagen.  Nonreciprocity and magnetic-free isolation based on optomechanical interactions.  Nature Communications, 7(1), Nov 2016.
  * [20] Shuyuan Xiao, Tao Wang, Tingting Liu, Chaobiao Zhou, Xiaoyun Jiang, and Jianfa Zhang.  Active metamaterials and metadevices: a review.  Journal of Physics D: Applied Physics, 53(50):503002, sep 2020.


11institutetext: Institute for Solar Physics, Dept. of Astronomy, Stockholm
University, Albanova University Center, 10691 Stockholm, Sweden
11email<EMAIL_ADDRESS>22institutetext: High Altitude Observatory,
National Center for Atmospheric Research, 80307, Boulder, CO, USA

# Rapid Blue- and Red-shifted Excursions in H$\alpha$ line profiles
synthesized from realistic 3D MHD simulations

S. Danilovic 11 J. P. Bjørgen  J. Leenaarts 11 M. Rempel 22

###### Abstract

Context. Rapid blue- and red-shifted events (RBEs/RREs) may have an important
role in mass-loading and heating the solar corona, but their nature and origin
are still debatable.

Aims. We aim to model these features to learn more about their properties,
formation and origin.

Methods. A realistic three-dimensional (3D) magneto-hydrodynamic (MHD) model
of a solar plage region is created. Synthetic H$\alpha$ spectra are generated
and the spectral signatures of these features are identified. The magnetic
field lines associated with these events are traced and the underlying dynamic
is studied.

Results. The model reproduces well many properties of RBEs and RREs, such as
spatial distribution, lateral movement, length and lifetimes. Synthetic
H$\alpha$ line profiles, similarly to observed ones, show strong blue- or red-
shift and asymmetries. These line profiles are caused by the vertical
component of velocity with magnitudes larger than $30-40$ km/s that appear
mostly in the height range of $2-4$ Mm. By tracing magnetic field lines, we
show that the vertical velocity that causes the appearance of RBE/RREs to
appear is always associated with the component of velocity perpendicular to
the magnetic field line.

Conclusions. The study confirms the hypothesis that RBEs and RREs are signs of
Alfvénic waves with, in some cases, a significant contribution from slow
magneto-acoustic mode.

###### Key Words.:

Sun: atmosphere – Sun: chromosphere – Sun: transition region –
magnetohydrodynamics (MHD) – radiative transfer

##  1 Introduction

The solar chromosphere is shaped from long, dense fibrils (Rutten 2007) that
connect large magnetic features organized at mesogranular to supergranular
scales. These magnetic features are the centers of most chromospheric dynamics
(Tsiropoula et al. 2012). Radially extending from those are dark stable
features called mottles, visible in cores of chromospheric spectral lines.
Moving to the wings of chromospheric lines, features named rapid blue- and
red-shifted events (RBEs/RREs) become dominant. Unlike mottles, these features
are very dynamic, short-lived, long and thin, but also visible near the edges
of the magnetic network. Although it is generally believed that they have an
important role in mass-loading and heating the solar corona, they are still a
puzzle.

The RBEs were first detected in the blue wing of the Ca II $854.2$ nm line as
sudden broadening in $\lambda$t-dataslices (Langangen et al. 2008). Their
lifetimes, spatial extent and location, as well as sudden disappearance,
suggested that they might be the on-disk counterparts of then newly discovered
off-limb spicule dynamics (de Pontieu et al. 2007). Subsequent papers based on
observations with higher spatial and temporal resolution further strengthened
this hypothesis (Rouppe van der Voort et al. 2009; Sekse et al. 2012, 2013).
The spectral imaging data taken in both Ca II $854.2$ nm and H$\alpha$ lines
revealed that the RBEs vary strongly between these lines. The former shows
significantly fewer and shorter features appearing closer to the magnetic
field concentrations and generally before the H$\alpha$ RBEs. Also, the
corresponding Doppler shifts differ, being of order $20-50$ km/s in the
H$\alpha$ and $15-20$ km/s in the Ca II $854.2$ nm line. The differences were
explained by the chromospheric opacity in the lines, with H$\alpha$ sampling
higher layers and naturally faster upflows.

The properties of these features also differ depending on the detection method
used and the temporal cadence of the data. The distribution of the RBEs
lifetimes ranges from $5$ to $60$ s, with a slow cadence data shifting it to
the higher values. This may be due to the high recurrence times around
magnetic elements which gives apparently higher lifetimes at a lower cadence.
Furthermore, the overall behavior of RBEs changed with temporal resolution.
The first studies claimed that RBEs moved away from magnetic concentration
with Doppler shifts and widths increasing from the footpoint to the top. The
higher cadence data, however, revealed much more RBEs that display erratic
behavior. While there were cases that had a distinct rise phase, there were
also features that showed a complicated mixture of longitudinal variations.
There were many cases that seemingly came to a halt or even retreated or just
appeared suddenly along the full length within a few seconds. Some RBEs moved
sideways over a few hundred kilometers, while others hardly got displaced. It
also became clear that RBEs appear in groups with large and dark jet-like
features appearing first, followed by numerous shorter, thinner jets
(Yurchyshyn et al. 2013; Samanta et al. 2019).

Figure 1: Synthetic H$\alpha$ images at different wavelength positions: left -
blue wing at $v_{dopp}=-36$ km s-1, middle - nominal line center and right -
red wing $v_{dopp}=36$ km s-1. Blue vertical line marks the position of the
cut shown in Fig. 4. The movie is available in the supplementary
material.[Movie]

The RREs are similar absorption features but in the red wings of the Ca II
$854.2$ nm and H$\alpha$ lines (Sekse et al. 2013; Kuridze et al. 2015; Shetye
et al. 2016; Bose et al. 2019). RREs are found over the whole solar disk and
are located in the same regions as RBEs. Lengths, widths, lifetimes, and
average Doppler velocity of RREs are similar to RBEs. RREs have similar
occurrence rates as RBEs with the total occurrence rate increasing towards the
limb. Most RREs and RBEs are observed in isolation, but many examples of
parallel and touching RRE/RBE pairs are also found to be part of the same
spicule. Kuridze et al. (2015) detect more RREs and Sekse et al. (2013) more
RBEs. The latter study also shows that the RRE/RBE detection count ratio
increases from the disk center to the limb from $0.52$ to $0.74$. The higher
number of RBEs and the decreased imbalance towards the limb is interpreted as
an indication that field-aligned upflows have a significant contribution to
the net Dopplershift of the structure.

Transverse displacements of RREs and RBEs are also similar. There are examples
of transitions from RRE to RBE and vice versa that sometimes appear to
propagate along the structure. In some cases, even oscillatory behavior can be
found with periods averaging $90$ s and amplitudes of $200$ km. All these
characteristics are in line with a hypothesis that these are an indication of
Alfvénic waves associated with the swaying motion. Namely, the movement of the
’flux tube’ is detected in both image planes as transverse motion and along
the line-of-sight (LOS) as a sudden blue- and red-shifts in the wings of
spectral lines. The formation of such phenomena would also arise naturally in
case of magnetic reconnection, which was indicated by observations that show
RBEs generally appearing when new, mixed or unipolar fields are detected in
close proximity to network fields (Yurchyshyn et al. 2013; Deng et al. 2015;
Samanta et al. 2019). Although this scenario explains all peculiarities of
RBEs and RREs, alternative hypotheses must also be mentioned. One alternative
explanation is that the observed transverse motions are not waves, but just
transverse displacement of the flux tube as a whole from its initial position
as it happens when the interaction timescale between the granule and the flux
tube is comparable to or greater than the chromospheric cutoff period (Hasan &
Kalkofen 1999). Another alternative explanation is that the sudden shifts in
the spectral lines can be an indication of sheet-like structures (Judge et al.
2012; Lipartito et al. 2014).

There is a large number of numerical experiments that try to explain the
formation and origin of spicules (e.g. Uchida 1969; Hollweg 1982; Suematsu et
al. 1982; Murawski & Zaqarashvili 2010; Jess et al. 2012; Iijima & Yokoyama
2017; Martínez-Sykora et al. 2017). However, only two studies tried to
explicitly model RBEs and RREs (Kuridze et al. 2015; Srivastava et al. 2017).
These models fail to reproduce many properties of RBEs and RREs. Still, they
indicate that RBEs/RREs may be signatures of Alfvén waves generated by
torsional driving in the photosphere. The present manuscript describes a
three-dimensional (3D) model that reproduces spatial distribution, lateral
movement, length and lifetimes of RBEs and RREs, as well as strongly blue- or
red-shift asymmetric H$\alpha$ line profiles. We start with identifying the
RBE/RRE signatures in synthetic spectra generated from the numerical models.
We then compare their properties with observations and further study their
formation.

##  2 Model atmospheres and forward synthesis

The MURaM code (Vögler et al. 2005; Rempel 2017) is used to generate the
numerical 3D models of the solar atmosphere. The model is built in phases
starting from the nonmagnetic convection simulation to which a uniform bipolar
field of $200$ G is added. After the field configuration evolved, the
potential field extrapolations are used to extend the computational domain and
include the upper atmosphere. The final extent of the simulation domain is
$40\times 40\times 22$ Mm, spanning vertically from $-8$ Mm to $14$ Mm above
the average $\tau_{500}=1$ height. At the bottom boundary, a horizontal field
at roughly equipartition field strengths is allowed to emerge into the domain
(Rempel 2014). At the upper boundary, the field is set to be vertical. Also,
the upper boundary is open to outflows but closed to inflows. The horizontal
boundaries are periodic. In the last phases, the model is first run at a lower
resolution and then again at the double resolution until a relaxed state is
achieved. The grid spacing of the simulation used in this paper is $39$ km and
$21$ km in the horizontal and vertical directions, respectively. The model
contains only the most basic physics to be called ‘comprehensive’: 3D grey LTE
radiative transfer, a tabulated LTE equation of state, Spitzer heat
conduction, and optically thin radiative losses in the corona based on CHIANTI
(Landi et al. 2012). The numerical diffusive and resistive terms in this run
are set so that an effective magnetic Prandtl number P${}_{\rm m}>1$. This
makes viscous heating to be the largest contributor. The maximum Alfvén
velocity is chosen as max$(2$c${}_{s},3|$v${}_{\rm max}|)$ where vmax is the
maximum velocity and cs the speed of sound (Rempel 2017). This gives the
largest unaffected Alfvén velocity of $900$ km/s (Chen et al. 2021). The
snapshots are outputted every $\leavevmode\nobreak\ 2.5$ s and every fourth
snapshot is used as input for the subsequent radiative transfer computations.

For these computations, we used Multi3D (Leenaarts & Carlsson (2009)) with the
model atmosphere that we cropped in two ways. The vertical span of the input
model is limited to $360$ points which include the range of temperatures
relevant for the formation of the H$\alpha$ line. In the horizontal direction,
we use every second grid point to make computations faster. Tests show that
this does not affect the appearance of synthetic features. H$\alpha$ was
calculated using a 3-level plus continuum H model atom created by Bjørgen et
al. (2019) to increase numerical stability in Multi3D. The synthetic spectra
are generated for $29$ snapshots in total.

Figure 2: Automatically identified synthetic RBEs and RREs. The panels show
summed masks for RBEs and RREs, respectively. The black contours outline the
unsigned photospheric vertical field of $1000$ G.

Figure 3: Lifetime (left) and length (right) of tracked synthetic RBE (black)
and RREs (red).

##  3 Results

###  3.1 Spatial distribution of synthetic RBEs and RREs

Figure 1 shows synthetic H$\alpha$ images generated from one of the snapshots.
The magnetic field configuration is visible from the distribution of long
fibrils in the H$\alpha$ core shown in the middle panel. Panels on the sides
show synthetic RBEs and RREs. The wavelength positions correspond to Doppler
shifts of $\pm 36$ km/s, similar to observations. Thin elongated dark features
are visible in both H$\alpha$, wings, resembling their observed counterparts.
Besides thin features, the red-wing image shows slightly more spread-out
structures especially close to the ’slope’ generated by the inclined magnetic
field around $[7,25]$ Mm. The dark circular structure at $[33,20]$ Mm visible
in both wings is a pore.

Inspection of the movie corresponding to Fig. 1 reveals different kinds of
features. In both RBE and RRE images, there are examples that show apparent
motion or increase in length away and towards magnetic network. There are also
many examples that appear and disappear in situ, seemingly halfway between the
magnetic footpoints. Many features show transverse motion. All this agrees
with findings based on observations.

To study the statistics of the synthetic RBEs and RREs, we apply simple
filtering on the H$\alpha$ wing images, exposing all regions where the
intensity was below $22.3\times 10^{-6}$ erg/cm2/s/Hz/ster and aspect ratio
larger than 15. These numbers are somewhat arbitrarily chosen. The final masks
summed over $29$ snapshots are visible in Fig. 2. The figure shows the
occurrence and feature position with respect to magnetic elements, which are
outlined with black contours. The chosen criteria allow the inclusion of
fairly short features as the one at $[3,30]$ Mm in the top panel. It also
allows many extended structures visible in the H$\alpha$ red wing images to be
counted as RREs, as the bottom panel Fig. 2 shows. Although magnetic elements
are distributed all over the simulation domain, maps nevertheless show areas
with reduced activity at $[25,2.5]$, $[28,15]$ and $[35,30]$ Mm. In other
regions, RREs and RBEs can be found in seemingly equal numbers. The exception
is visible at the base of the ’slope’ where RREs prevail.

###  3.2 Properties of synthetic RBEs and RREs

We cluster masked features that overlap from one snapshot to the other and
count them as one feature. This gives the total number of RBEs and RREs to be
231 and 322, respectively. The distribution of their lifetimes is shown in
Fig. 3. Many features traced this way appear and disappear in one snapshot,
giving them lifetimes of less than $10$ s. Most of the extra RREs is placed
into this bin. The distribution exponentially decreases with the lifetime.
Only around $5$ % of both RRE and RBEs seems to be visible for at least a
minute.

For each case, we track the evolution and include their most extended length
in the distribution shown on the right of Fig. 3. We see that most features
are shorter than $3$ Mm in size, although some cases extend to over $5$ Mm.
There seems to be no preference for the extra short-lived RREs. They seem to
be distributed proportionally to all length bins.

###  3.3 Formation of synthetic RBEs and RREs

Figure 4 shows a vertical cut over several RBE and RRE at the time instant
shown in Fig. 1. Overploted are $\tau=1$ contours for the same wavelength
positions as in Fig. 1. Corrugated black lines thus show that the H$\alpha$
line core is formed in the range $2-4$ Mm, while jumps in green and red
contours indicate the positions of the crossed RBEs and RREs respectively.
Most of them, the four RBEs and two RRE, are formed close to the H$\alpha$
line core.They are formed at similar heights, which results in similar
intensities visible in Fig. 1. The bottom panel of Fig. 4 shows that at each
position where RBE/RREs is formed, there is a strong upflow/downflow in the
vertical velocity map. The field strength in those regions is still higher
than $100$ G, as the third panel in Fig. 4 shows.

Figure 4: Vertical cut along the blue vertical lines shown in Fig. 1. Vertical
cuts in temperature, density, magnetic field strength and vertical velocity
are shown from top to bottom. The black line in all panels shows the $\tau=1$
at the nominal position of the H$\alpha$ line center. Green and red contours
mark $\tau=1$ levels at wavelength positions that correspond to Doppler shifts
of $v_{\rm dopp}=-36$ km s-1 and $v_{\rm dopp}=36$ km s-1 respectively.
Figure 5: The H$\alpha$ four-panel diagram for the RBE at $x=10.7$ Mm in Fig.
4. The plotted quantities, as functions of frequency and height, are indicated
in the upper left corner of each subplot and plotted in a grayscale. The
vertical component of velocity as a function of height is shown in all panels
(solid white line), with upwards velocity being positive. The red line marks
the $\tau=1$ position as a function of frequency and height. Planck and total
source function are plotted in the upper right panel with dotted and dashed
lines, respectively. The line profile is shown with solid white line in lower
right panel.  Figure 6: The H$\alpha$ four-panel diagram for RRE at $x=11.25$
Mm in Fig. 4. The outline is the same as in Fig. 5.

Examples at $x=10.7$ Mm and $x=11.25$ Mm in Fig. 4 are analyzed in more detail
using four-panel schematic H$\alpha$ line formation diagrams (Carlsson & Stein
1997; Leenaarts et al. 2010) shown in Figs. 5 and 6. An example of an RBE is
shown in Fig. 5. The emergent profile in the lower right panel is strongly
blue-shifted with a small reversal at the red wing. Overplotted in the same
panel is the line-of-sight (LOS) velocity profile displaying a shift in the
same direction. The predominant upflows reaching $40$ km/s are spread in the
height range $2-3$ Mm. Because of these upflows, the location with high
$\chi_{\nu}/\tau_{\nu}$ (upper left panel) is strongly shifted at z$=3$ Mm,
exactly where the optical depth reaches unity at those frequencies (lower left
panel). The source function is completely decoupled from the Planck function
and stays fairly flat with two small local maxima higher up, at $3.5$ and $4$
Mm. The opacity has a shape characteristic for the H$\alpha$ line with sharp,
sudden shifts of $\tau=1$ in both wings but strongly shifted to the positive
frequencies. As a result, the contribution function jumps to the height of $3$
Mm and stays there over a wide range of frequencies. A small reversal in the
red wing is formed due to a temperature and density increase that are also
visible in Fig. 4. The formation mechanism of the red wing reversal is the
same as for spectra signatures of the so-called Ellerman bombs (Danilovic
2017).

The four-panel diagram created for the pixel at x$=11.25$ Mm in Fig. 4 is
presented in Fig. 6. This pixel experiences a much more complex opacity
distribution than in the case of Fig. 5. Here, we have numerous bright regions
in the upper left panel at a variety of altitudes and frequencies. Two major
effects on the contribution function come from high $\chi_{\nu}/\tau_{\nu}$
located at z$=3$ and z$=2.3-2.5$ Mm. The latter is the location where a strong
downflow reaches $50$ km/s. At the former location, there is a switch from an
upflow to a downflow, both slightly weaker reaching $20$ km/s. The source
function shows a local maximum in both locations. As a result, line wings and
core intensities are all experiencing contributions from both heights. The
blue-shifted material also coincides with the location where the optical depth
reaches unity at that frequency, so the blue wing gets extended, resulting in
a very asymmetric line profile.

Figure 7: The overview of the RBE/RRE sample that is chosen for magnetic field
line tracing. For each of the $14$ examples, a traced magnetic field line is
plotted over the vertical component of the photospheric magnetic field. The
diamonds mark the position of the starting seed used for line tracing: red
marks RREs and blue RBEs. Yellow stars mark the position of the F1.

Figure 8: The horizontal (y, solid line) and vertical (z, dashed line)
displacement at constant x for $4$ examples of RBE/RRE. The vertical line
marks the time when the seed is ’planted’. The cases F0 and F1 are RREs and F6
and F9 RBEs. The labels are the same as in Fig. 7.

###  3.4 Origin of synthetic RBEs and RREs

We further choose $14$ examples around the center of the simulation domain to
study what causes their appearance. An equal number of RBEs and RREs are
chosen. They appear at different time instances distributed over $9$ different
snapshots. Figure 7 gives an overview of positions of one traced field line
per each chosen feature. We trace the magnetic field lines as curvatures in 3D
space (Leenaarts et al. 2015). The time cadence between snapshots used for the
field tracing is around $2.5$ s. For each of the examples, we pick several
seeds along the feature, visible either in blue or red synthetic H$\alpha$
wing, to advect as corks forward and backward in time. For each seed, a
formation height at that specific wavelength position is taken as the
z-coordinate in the 3D space. The individual cases are labeled with F$0-13$ as
shown in Fig. 7. Cases F12 and F13 are the RBE and RRE visible in Fig. 4 and
further discussed in four-panel diagrams. The yellow stars in Fig. 7 mark the
position of the RRE associated with F1, which appears close to F0 when F0
disappears. The line tracing indicates that these RREs belong to the same
event since their origin can be traced to overlapping field lines111The movie
is available as supplementary material..

Figure 9: The evolution along a field line traced for four cases labeled in
the top left corner. From left to right: total velocity v, vertical velocity
vz, longitudinal vL and perpendicular vP velocity. Diamonds mark the location
of starting seeds. The dotted lines follow a signal propagating with the
Alfvén speed along the loop.

To analyze the displacement of each chosen example, we place the ’slit’ at the
x coordinate of the starting seed and then follow the y and z coordinates of
the traced field line in time at that fixed x position. Figure 8 shows
resulting horizontal and vertical displacements for four cases after a linear
fit is subtracted. The shifts in two directions get in the phase when the
RBE/RRE appears and sometimes show oscillatory behavior over a short time. For
example, F0 shows a sinusoidal wave pattern in y coordinate that lasts for
$\sim 40$ s, while in the case of F9, the same lasts $\sim 100$ s. In most
cases, we see more erratic behavior like for F6. The four cases shown in Fig.
8 illustrate that the range of displacements stays within $200$ km with
amplitudes of a few $10$ km during the oscillatory period.

Figure 10: Proxies of different wave modes: Alfvén (left column), fast (middle
column) and slow mode waves convolved with the contribution function at
different wavelength positions: top row - blue wing at v${}_{d\rm opp}=-36$ km
s-1, middle row - red wing v${}_{\rm dopp}=36$ km s-1 and bottom row - nominal
line center. Blue vertical line marks the position of the cut shown in Fig. 4.

The field tracing reveal two common characteristics for all cases. This is
illustrated in Fig. 9 on four RBE/RRE examples. The figure shows velocities
along a traced field line as a function of time. Diamonds mark the location of
starting seeds and dotted lines follow signals propagating with the Alfvén
speed along the loop. The first common characteristic for all cases is that we
detect a strong flow at one of the footpoints which perturbs the field line.
This is visible in the vertical component of velocity displayed in the second
column. For cases F3, F9 and F13, alternating upflow and downflow are visible
close to the first footpoint. For case F0, the same happens at the second
footpoint. The flow propagates along loops with Alfvén speed of $\sim 300$
km/s for all inspected cases. The second common characteristic is that the
vertical velocity that causes RBE/RREs to appear in H$\alpha$ line wings is
always associated with the component of velocity perpendicular to the magnetic
field line. This is visible in the third and fourth column, where we see high-
velocity pattern only in vP and nothing in vL at the time instant in proximity
to the position of the starting seed. Only in the case F13, strong flow is
visible along the field line vL, although not at the position of this specific
seed.

We can also use instantaneous proxies for three wave modes slow $f_{long}$ and
fast $f_{fast}$ magneto-acoustic and Alfvén waves $f_{alf}$ (Khomenko et al.
2018; Cally 2017):

$\displaystyle
f_{alf}=\hat{e_{\parallel}}\cdot\bigtriangledown\times\mathbf{v}$ (1)
$\displaystyle
f_{fast}=\bigtriangledown\cdot\left(\mathbf{v}-\hat{e_{\parallel}}v_{\parallel}\right)$
(2) $\displaystyle
f_{long}=\hat{e_{\parallel}}\cdot\bigtriangledown\left(\mathbf{v}\cdot\hat{e_{\parallel}}\right)$
(3)

where $\mathbf{v}$ is the velocity vector and $\hat{e_{\parallel}}$ the field-
aligned unit vector. To visualize where the wave power for each mode is
located with respect to the formation of RBEs and RREs, each quantity is
convoluted with a contribution function at the corresponding wavelength. The
resulting maps are shown in Fig. 10.

There is a significant difference between maps corresponding to the line
center and the line wings. The latter are largely contaminated with the
contribution from the photosphere. Only in regions where RBEs and RREs occur
the contribution function shifts to higher layers, as illustrated in the
previous section, and the maps show elongated features. If we compare Fig. 10
with Fig. 1, we can identify many features in the intensity. Most of them are
visible only in $f_{alf}$. In some cases, weak signatures can be also traced
in $f_{fast}$ maps. A few cases can also be identified in all three wave
proxies. Cases F12 and F13 belong in the last group of features. There also
seems to be a trend that $f_{long}$ and $f_{fast}$ signatures are located
closer to the footpoints and $f_{alf}$ proxy extends further along a feature
(see also Danilovic 2022).

##  4 Discussion and conclusions

H$\alpha$ line profiles generated from 3D model of solar plage show RBEs and
RREs signatures. Their properties: spatial distribution, length and lifetimes
are very similar to the properties of their observed counterparts. The
synthetic features appear close to the magnetic network or at the apex of the
loops. Similarly to the observations, their behavior can be both consistent or
erratic. In the first case, they move away or closer to the network. In the
second case, they appear and disappear far from any magnetic concentration.
They show lateral motion of a few $10-100$ km. This agrees with observations
(Sekse et al. 2012, 2013). We identify cases where features exhibit lateral
motion by being offset along the whole length, as reported in observations
(Kuridze et al. 2015). An example of this is the feature labeled F0 and F1.

Synthetic H$\alpha$ line profiles associated with these features are either
completely blue- or red-shifted or they show asymmetric, extended wing. These
line profiles are caused by the vertical component of velocity with magnitudes
larger than $30-40$ km/s. In our model, these velocities appear mostly in the
height range of $2-4$ Mm, at the line core formation height or slightly below
that. By tracing magnetic field lines, we show that the vertical velocity that
causes RBE/RREs to appear in H$\alpha$ line wings is always associated with
the component of velocity perpendicular to the magnetic field line. Further
analysis shows that features mainly outline the proxy location of Alfvén waves
and to a lesser degree of fast magneto-acoustic waves. In some cases, proxies
of all three wave modes can be found coinciding with the RBE/RREs formation.
As most events were associated to Alfven waves, the hypothesis that RBEs and
RREs are signs of Alfvénic waves (Sekse et al. 2013) is confirmed.

Tracing the magnetic field lines also demonstrates that strong flows at one of
the footpoints perturbs field lines, as a result of which RBE/RREs are formed.
This explains why RBEs appear in groups (Yurchyshyn et al. 2013; Samanta et
al. 2019). Strong jets would generate perturbations in all adjacent field
lines which would result in numerous shorter, thinner features like the ones
modeled in this study. Our model, however, does not show many cases with
strong jets. Likewise, the model generates more RREs than RBEs. This is
inconsistent with observations. We find several possible reasons for this. Our
detection method included many features that might not be counted as RRE if
found in the observations. Among these are the ones that resemble dynamic
fibrils. These features, when observed, show $5$ minute oscillations in the
region with a more inclined field (De Pontieu et al. 2007). Our synthetic
H$\alpha$ dataset covers less than that so there is a possibility that we
caught most of them in their descending phase. This will be the subject of our
next study.

The resolution of our model can also be a reason for generating RBE/RRE ratio
different than the one observed. Namely, the higher resolution would produce
more turbulent flows, resulting in footpoint motions that could lead to more
frequent magnetic reconnection. As the reconnection would happen near
footpoints, the upflows would prevail in the height ranges where RBE/RREs form
and thus generate more RBEs than RREs. Also, flux recycling or small-scale
emergence generated at higher resolution would result in fast lateral movement
and would have the same effect. Although most of the flux emerging on granular
scales gets pulled back (Danilovic et al. 2010), we do not exclude the
scenario where these features are generated by the reconnection of emerging
with the pre-existing field, which our simulations do not model. Finally, the
velocity magnitude increases with resolution, hence more events will result in
sufficient line shifts. The second possible reason for the mismatch of
simulations and observations could be in not including the ion-neutral
interaction effects i.e. ambipolar diffusion (Martínez-Sykora et al. 2017).
With ambipolar diffusion included, the plasma is not fully frozen in, so the
jets may be more vertical, not following the field lines. This would lead to a
larger vertical velocity component and the identification of more jets.

Apart from the unresolved cause of the mismatch between simulated and observed
RBE/RREs ratio, a few more questions remain for the following studies. One of
which is what determines the width of these features. In models like
Srivastava et al. (2017) or in a cartoon by Sekse et al. (2013), RBEs and RREs
are generated at the side edges of magnetic flux tubes. Based on this
interpretation, observations (Srivastava et al. 2017; Shetye et al. 2021)
would suggest that these flux tubes are resolved and sometimes as an arcsec
wide. Srivastava et al. (2021) further argue that if the waveguide is not
spatially resolved the opposing Doppler shifts will result in non-thermal
broadening in optically thin lines while the chromospheric observables would
sample only ’the outermost shells of the structure’. We find no distinct
structures as such in our model, although some synthetic RBE/RREs pairs can be
traced to the same source. The model shows that while H$\alpha$ line core
outlines density ridges (Leenaarts et al. 2012; Danilovic 2022), the formation
of the line wings is more complex.

Finally, we did not discuss the photospheric source of these features in the
model. Movies of traced field lines show braiding of field lines close to a
footpoint. This suggest the formation of photospheric vortex flow. The vortex
signatures are present in all areas of the simulation domain where magnetic
field is stronger. They are most prominent in the temperature maps sampled at
upper photosphere, as shown by Moll et al. (2012). Vorticity caused by
torsional motions in the photosphere can excite torsional Alfvén pulses that
propagate along the magnetic field lines to the upper layers in the form of
torsional Alfvén waves (Shelyag et al. 2013; Liu et al. 2019; Battaglia et al.
2021). All this is agreement with the results presented in this study.

###### Acknowledgements.

This project has received funding from Swedish Research Council (2021-05613),
Swedish National Space Agency (2021-00116) and the Knut and Alice Wallenberg
Foundation. This material is based upon work supported by the National Center
for Atmospheric Research, which is a major facility sponsored by the National
Science Foundation under Cooperative Agreement No. 1852977. This research data
leading to the results obtained has been supported by SOLARNET project that
has received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement no 824135. The calculations were
performed on resources provided by the Swedish National Infrastructure for
Computing (SNIC) at the National Supercomputer Centre (NSC) at Linköping
University and the PDC Centre for High Performance Computing (PDC-HPC) at the
Royal Institute of Technology in Stockholm.

##  Appendix A Supplementary material

Besides the movie version of Fig. 1, we share two additional movies that show
magnetic field lines for cases F9, F0 and F1. The field lines of F0 and F1 are
plotted together, as shown in Fig. 11 so that initial overlapping of magnetic
field lines is visible.

Figure 11: Magnetic field lines traced for RREs F0 (blue lines) and F1 (green
lines). Panel on the left shows top view with vertical photospheric magnetic
field in the background. Panels on the right show projection of field lines in
x (bottom) and y (top) planes with density cut at x,y coordinates of one of
the initial seeds.

## References

  * Battaglia et al. (2021) Battaglia, A. F., Canivete Cuissa, J. R., Calvo, F., et al. 2021, A&A, 649, A121
  * Bjørgen et al. (2019) Bjørgen, J. P., Leenaarts, J., Rempel, M., et al. 2019, A&A, 631, A33
  * Bose et al. (2019) Bose, S., Henriques, V. M. J., Joshi, J., et al. 2019, A&A, 631, L5
  * Cally (2017) Cally, P. S. 2017, MNRAS, 466, 413
  * Carlsson & Stein (1997) Carlsson, M. & Stein, R. F. 1997, ApJ, 481, 500
  * Chen et al. (2021) Chen, F., Rempel, M., & Fan, Y. 2021, arXiv:2106.14055
  * Danilovic (2017) Danilovic, S. 2017, A&A, 601, A122
  * Danilovic (2022) Danilovic, S. 2022, arXiv:2208.03744
  * Danilovic et al. (2010) Danilovic, S., Beeck, B., Pietarila, A., et al. 2010, ApJ, 723, L149
  * De Pontieu et al. (2007) De Pontieu, B., Hansteen, V. H., Rouppe van der Voort, L., et al. 2007, ApJ, 655, 624
  * de Pontieu et al. (2007) de Pontieu, B., McIntosh, S., Hansteen, V. H., et al. 2007, PASJ, 59, S655
  * Deng et al. (2015) Deng, N., Chen, X., Liu, C., et al. 2015, ApJ, 799, 219
  * Hasan & Kalkofen (1999) Hasan, S. S. & Kalkofen, W. 1999, ApJ, 519, 899
  * Hollweg (1982) Hollweg, J. V. 1982, ApJ, 257, 345
  * Iijima & Yokoyama (2017) Iijima, H. & Yokoyama, T. 2017, ApJ, 848, 38
  * Jess et al. (2012) Jess, D. B., Pascoe, D. J., Christian, D. J., et al. 2012, ApJ, 744, L5
  * Judge et al. (2012) Judge, P. G., Reardon, K., & Cauzzi, G. 2012, ApJ, 755, L11
  * Khomenko et al. (2018) Khomenko, E., Vitas, N., Collados, M., et al. 2018, A&A, 618, A87
  * Kuridze et al. (2015) Kuridze, D., Henriques, V., Mathioudakis, M., et al. 2015, ApJ, 802, 26
  * Landi et al. (2012) Landi, E., Del Zanna, G., Young, P. R., et al. 2012, ApJ, 744, 99
  * Langangen et al. (2008) Langangen, Ø., De Pontieu, B., Carlsson, M., et al. 2008, ApJ, 679, L167
  * Leenaarts & Carlsson (2009) Leenaarts, J. & Carlsson, M. 2009, The Second Hinode Science Meeting: Beyond Discovery-Toward Understanding, 415, 87
  * Leenaarts et al. (2012) Leenaarts, J., Carlsson, M., & Rouppe van der Voort, L. 2012, ApJ, 749, 136
  * Leenaarts et al. (2015) Leenaarts, J., Carlsson, M., & Rouppe van der Voort, L. 2015, ApJ, 802, 136
  * Leenaarts et al. (2010) Leenaarts, J., Rutten, R. J., Reardon, K., et al. 2010, ApJ, 709, 1362
  * Lipartito et al. (2014) Lipartito, I., Judge, P. G., Reardon, K., et al. 2014, ApJ, 785, 109
  * Liu et al. (2019) Liu, J., Nelson, C. J., Snow, B., et al. 2019, Nature Communications, 10, 3504
  * Martínez-Sykora et al. (2017) Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017, Science, 356, 1269
  * Moll et al. (2012) Moll, R., Cameron, R. H., & Schüssler, M. 2012, A&A, 541, A68
  * Murawski & Zaqarashvili (2010) Murawski, K. & Zaqarashvili, T. V. 2010, A&A, 519, A8
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132
  * Rempel (2017) Rempel, M. 2017, ApJ, 834, 10
  * Rouppe van der Voort et al. (2009) Rouppe van der Voort, L., Leenaarts, J., de Pontieu, B., et al. 2009, ApJ, 705, 272
  * Rutten (2007) Rutten, R. J. 2007, The Physics of Chromospheric Plasmas, 368, 27
  * Samanta et al. (2019) Samanta, T., Tian, H., Yurchyshyn, V., et al. 2019, Science, 366, 890
  * Sekse et al. (2012) Sekse, D. H., Rouppe van der Voort, L., & De Pontieu, B. 2012, ApJ, 752, 108
  * Sekse et al. (2013) Sekse, D. H., Rouppe van der Voort, L., De Pontieu, B., et al. 2013, ApJ, 769, 44
  * Sekse et al. (2013) Sekse, D. H., Rouppe van der Voort, L., & De Pontieu, B. 2013, ApJ, 764, 164
  * Shelyag et al. (2013) Shelyag, S., Cally, P. S., Reid, A., et al. 2013, ApJ, 776, L4
  * Shetye et al. (2016) Shetye, J., Doyle, J. G., Scullion, E., et al. 2016, A&A, 589, A3
  * Shetye et al. (2021) Shetye, J., Verwichte, E., Stangalini, M., et al. 2021, arXiv:2112.14486
  * Srivastava et al. (2017) Srivastava, A. K., Shetye, J., Murawski, K., et al. 2017, Scientific Reports, 7, 43147
  * Srivastava et al. (2021) Srivastava, A. K., Ballester, J. L., Cally, P. S., et al. 2021, Journal of Geophysical Research (Space Physics), 126
  * Suematsu et al. (1982) Suematsu, Y., Shibata, K., Neshikawa, T., et al. 1982, Sol. Phys., 75, 99
  * Tsiropoula et al. (2012) Tsiropoula, G., Tziotziou, K., Kontogiannis, I., et al. 2012, Space Sci. Rev., 169, 181
  * Uchida (1969) Uchida, Y. 1969, PASJ, 21, 128
  * Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
  * Yurchyshyn et al. (2013) Yurchyshyn, V., Abramenko, V., & Goode, P. 2013, ApJ, 767, 17