# A Domain Adaptation Model for Carotid Ultrasound: Image Harmonization, Noise
Reduction, and Impact on Cardiovascular Risk Markers

Mohd Usama∗
<EMAIL_ADDRESS>&Emma Nyman†
<EMAIL_ADDRESS>&Ulf Näslund†
<EMAIL_ADDRESS>&Christer Grönlund∗
<EMAIL_ADDRESS>&

∗Department of Diagnostics and Intervention, and Biomedical Engineering, Umea
University, Umea, Sweden
†Department of Public Health and Clinical Medicine, Umea University, Umea,
Sweden

A preprint -

###### Abstract

Deep learning has been used extensively for medical image analysis
applications, assuming the training and test data adhere to the same
probability distributions. However, a common challenge arises when dealing
with medical images generated by different systems or even the same system
with varying parameter settings. Such images often contain diverse textures
and noise patterns, violating the assumption. Consequently, models trained on
data from one machine or setting usually struggle to perform effectively on
data from another. To address this issue in ultrasound images, we proposed a
Generative Adversarial Network (GAN) based model in this paper. We formulated
image harmonization and denoising tasks as an image-to-image translation task,
wherein we modified the texture pattern and reduced noise in Carotid
ultrasound images while keeping the image content (the anatomy) unchanged. The
performance was evaluated using feature distribution and pixel-space
similarity metrics. In addition, blood-to-tissue contrast and influence on
computed risk markers (Gray scale median, GSM) were evaluated. The results
showed that domain adaptation was achieved in both tasks (histogram
correlation 0.920 and 0.844), as compared to no adaptation (0.890 and 0.707),
and that the anatomy of the images was retained (structure similarity index
measure of the arterial wall 0.71 and 0.80). In addition, the image noise
level (contrast) did not change in the image harmonization task (-34.1 vs 35.2
dB) but was improved in the noise reduction task (-23.5 vs -46.7 dB). The
model outperformed the CycleGAN in both tasks. Finally, the risk marker GSM
increased by 7.6 (p<0.001) in task 1 but not in task 2. We conclude that
domain translation models are powerful tools for ultrasound image improvement
while retaining the underlying anatomy but that downstream calculations of
risk markers may be affected.

_K_ eywords Deep learning $\cdot$ Generative Adversarial Network $\cdot$
Medical image analysis $\cdot$ Cardiovascular disease assessment $\cdot$
Domain Adaptation $\cdot$ Image harmonization $\cdot$ Noise reduction $\cdot$
Carotid ultrasound images

##  1 Introduction

Deep learning has been used extensively for medical image analysis
applications, assuming the training and test data adhere to the same
probability distributions. However, a common challenge arises when dealing
with medical images generated by different systems or even the same system
with varying parameter settings. Such images often contain diverse textures
and noise patterns, violating the assumption. Consequently, models trained on
data from one machine or setting usually struggle to perform effectively on
data from another.

Deep learning for automatic analysis of ultrasound images has attracted much
attention for multiple clinical applications, such as atherosclerosis risk
assessment [1, 2], artery plaque detection and segmentation [3, 4], breast
cancer detection and segmentation [5], etc. However, images from different
ultrasound systems often vary in many aspects due to different imaging
parameter settings, algorithm design, hardware components, etc. Therefore, the
texture pattern of images acquired by different ultrasound systems may not
belong to the same feature population distribution. For example, images taken
from different ultrasound machines or different settings usually differ in
texture pattern and noise. In addition, the images’ texture and noise can vary
depending on the composition and distribution of tissues of the imaged
subject. This may cause generalizability problems when using pre-trained
models for new images from different systems or on images from patients with
varying noise distributions.

To overcome these limitations, several studies have leveraged Generative
Adversarial Networks (GANs) [6] in ultrasound imaging for domain adaptation
and denoising individually. These networks, known for their ability to
generate realistic images, have shown promise in various areas of medical
imaging, such as denoising [7], domain adaptation [8], disease translation
[9], synthetic image generation [10], segmentation[5], and classification
[11]. Several researchers have introduced different methods and techniques in
ultrasound image processing to enhance these networks’ effectiveness and
efficiency [12]. They have explored the potential of GANs, which has resulted
in a surge of publications focusing on quality adaptation and denoising [13].
Zhou et al. [14] proposed a two-stage GAN for translating low-quality
ultrasound images into high-quality images, which are usually degraded due to
the limited size of portable ultrasound devices. Xia et al. [15] propose a
Multilevel Structure-Preserved Generative Adversarial Network (MSP-GAN) for
domain adaptation of ultrasound images while maintaining image intravascular
structures. The CycleGAN network [16] has been widely used for different
applications in domain-to-domain translation tasks that do not require paired
image data. Its application spans various modalities of medical images,
encompassing ultrasound, MRI and CT [17, 18]. Moreover, CycleGAN architecture
has shown promising histopathology results in assisting diagnoses more
precisely by standardizing microscopy staining [19]. Additionally, several
other works used CycleGAN for image quality adaptation and denoising, such as
Athreya et al. [20], who proposed a CycleGAN-based model with perceptual loss
for quality enhancement of ultrasound images. Huang et al. [21] proposed a
stability-enhanced CycleGAN-based model to harmonize ultrasound images from
various medical centres through domain transformation with good preservation
of image details. Tan et al.[22] proposed an unsupervised model using CycleGAN
with selective kernel networks for adapting select features, patchGAN
discriminator and perceptual loss in LDCT imaging to generate CT. Tang et al.
[23] proposed a CycleGAN-based denoising network for low-dose unpaired CT
images with prior image information. Similar work is also done by Zhao et
al.[24], who proposed a CycleGAN-based dual-scale similarity-guided model.

Image harmonization in carotid ultrasound images is crucial for standardizing
measurements and reducing variability across different scanners and imaging
conditions. Various methods have been developed to address this issue. For
instance, a study by Elatrozy et al. [25] demonstrated the standardization of
plaque echodensity in carotid ultrasound images, allowing for comparable
measurements. Furthermore, Kakkos et al. [26] highlighted the importance of
standardizing resolution settings during ultrasound scanning to minimize
variability in texture analysis of carotid plaques, emphasizing the need for
resolution normalization algorithms. Additionally, model-based unsupervised
domain adaptation has been successfully applied in vessel segmentation tasks,
where neural networks are initially trained on simulated ultrasound images and
then fine-tuned on in-vivo data, resulting in significant performance
improvements without the need for extensive manual labelling [27]. These
approaches demonstrate the potential of leveraging advanced techniques to
enhance the generalizability and accuracy of carotid ultrasound image analysis
through domain adaptation and harmonization strategies. Moreover, Cheung et
al. [28] developed an algorithm to correct attenuation artifacts and normalize
intensity in contrast-enhanced ultrasound images of carotid arteries,
enhancing quantification and reducing variability, and Loizou et al. [29]
focus on evaluating the quality of ultrasound images in the carotid artery
after applying normalization and speckle reduction filtering techniques. Image
quality assessment is crucial for tasks like assessing atherosclerotic disease
and telemedicine image transfer. Noise reduction in ultrasound images plays
another crucial role in enhancing the quality of the images. Research has
shown that GANs, such as the wavelet-based GAN proposed in the studies [30,
31], can effectively suppress speckle noise while preserving image features,
leading to high-quality image reconstruction. Additionally, the use of GANs in
ultrasound image processing, as demonstrated in various studies [32, 33, 34],
has shown promising results in improving image quality. By leveraging GANs,
particularly customized models like conditional GANs and wavelet-based GANs,
researchers have been able to address the challenges of noise reduction in
carotid ultrasound images, ensuring accurate diagnostic information is
preserved for clinical applications.

Taken together, previous work shows that, among others, image quality
adaptation and denoising are important challenges in carotid ultrasound
images, and separate models have been developed to solve these tasks
individually [7]. However, this requires building different models for each
task, including data collection from each machine with similar parameter
settings to train the models, which is time-consuming and tedious, requiring
separate resources and computation power, which is expensive [19, 35, 36].

This paper proposes a new GAN-based model that can address image harmonization
and image denoising tasks with the same architecture. The proposed model
contains one generator and two discriminators, and model training is
regularized with content, noise and adversarial losses. The generator aims to
take a source image and generate an image with features of the target domain.
The first discriminator with content loss ensures that the generated image’s
content (anatomy) is preserved, and the second discriminator with noise loss
enforces the noise features in the image.

The main contributions of this paper are:

  * •

We propose a GAN-based model to address the problem of denoising and domain
adaptation of ultrasound images through an unpaired image-to-image translation
task. The model is trained to transfer target domain image features into the
source domain image without changing the underlying anatomical content of the
image.

  * •

We introduce the noise loss, which employs the Wasserstein distance to measure
the dissimilarity between the stylistic features derived from the generator’s
early layer and the style of target domain images.

  * •

We fine-tune the generator architecture from Jun-Yan Zhu et al. [16] by adding
two extra blocks (one convolution block and one deconvolution block) to manage
the feature map extraction from the first three layers.

  * •

We evaluate the performance of the domain adaptation task and preservation of
anatomical content. In addition, we evaluate the domain adaptation’s impact on
clinical risk markers. To our knowledge, this is the first study to do so.

The rest of the paper is organized as follows. Section II explains the method.
Section III describes the results, and section IV discusses the findings. The
last section concludes the paper.

##  2 Methods

This paper presents an unpaired image-to-image translation model such that an
input image with one feature distribution can be converted to an output image
with another feature distribution while preserving the underlying anatomical
content of the images. We target the problems of image feature harmonization
and noise reduction. Generally, paired images from the corresponding domains
(different systems, different noise levels) are often unavailable. Therefore,
we chose to solve this as an unpaired image-to-image translation task. In this
work, we consider the underlying anatomical information of the images as
content and noise pattern as a style and formulate the task as a style
transfer task.

###  2.1 GAN-based Model

To solve the aforementioned problem, we propose a GAN-based model consisting
of one generator and two discriminators. We formulate the task as an unpaired
image-to-image translation task and incorporate a noise loss. This loss is
determined using the Wasserstein distance and the style of feature maps
derived from the early layers of the generator and the target domain images
(see figure 1).

Our approach is inspired by methods [37, 13], demonstrating the independent
processing of feature representations learned by neural networks. Our features
transfer approach combines a domain-to-domain translation concept, which
builds a mapping between two domains for adaptation instead of mapping between
two individual images, unlike [37], with a noise loss employing the
Wasserstein distance and the stylistic features derived from the generator
instead of the discriminator, unlike NAGAN [13].

We define the training source dataset as $S$ and the target dataset as $T$.
Our main objective is to learn a mapping from $S$ to $T$ so that the feature
distribution of image $x\epsilon S$ will be translated to $x^{{}^{\prime}}$
images such that anatomical content is preserved while the style of image
feature distribution is similar to that of $T$. The generator’s job is to
translate the source images that map into the target domain. The $D_{n}$
distinguishes the texture pattern and noise between domains $S$ and $T$. The
$D_{c}$ distinguishes the anatomical content of images. $\mathcal{L}_{an}$ and
$\mathcal{L}_{ac}$ are used for adversarial training. $L_{c}$ and $L_{n}$ are
used to preserve the anatomical contents, transfer texture patterns and reduce
the noise, respectively. Finally, adversarial losses are combined with content
and noise losses to train the model. The computation of all loss functions is
given below.

Figure 1: a) Proposed model. $G$ denotes the generator network and two
discriminators, $D_{c}$ and $D_{n}$, trained by adversarial losses
$\mathcal{L}_{ac}$ and $\mathcal{L}_{an}$, respectively. $x\sim S$ denotes the
image set from the source domain. $y\sim T$ represents the image set from the
target domain. $G(x)$ denotes the image set generated by $G$. $L_{c}$
indicates the content loss calculated from the feature difference of the
source and $G(x)$ image. $L_{n}$ represents the noise loss calculated from the
feature difference between the target image and the early layers of the
generator. b) a residual block.

####  2.1.1 Adversarial losses

The $\mathcal{L}_{an}$ and $\mathcal{L}_{ac}$ are used to train the $D_{n}$
and $D_{c}$, respectively. The $D_{n}$ job distinguishes between the target
image $y\epsilon T$ and the translated image $G(x)$. The $D_{c}$ distinguishes
between the translated image $G(x)$ and the real image $x\epsilon S$. The
Computation of the losses are as follows:

$\mathcal{L}_{an}(G,D_{n})=\mathbb{E}_{y\sim T}[\log
D_{n}(y)]+\mathbb{E}_{x\sim S}[\log(1-D_{n}(G(x)))]$ (1)

$\mathcal{L}_{ac}(G,D_{c})=\mathbb{E}_{x\sim S}[\log
D_{c}(x)]+\mathbb{E}_{x\sim S}[\log(1-D_{c}(G(x)))]$ (2)

####  2.1.2 Content loss

To preserve the anatomical contents of the image, the content loss $L_{c}$ is
computed from the last layer of generator $G$ as follows:

$\mathcal{L}_{c}(G)=\mathbb{E}_{x\sim
S}\left[\left|F(G(x_{i}))^{l}-F(x_{i})^{l}\right|\right]$ (3)

where $F(\cdot)^{l}$ represents the feature maps from the last layer.

####  2.1.3 Noise loss

We formulate noise as style [37]. The style of an image was computed
mathematically by Gram Matrix [38] $G_{r}$ as follows:

$G_{r}(y)^{(G,l)}_{(i,j)}=\text{vec}\left[F(y)^{(G,l)}_{i}\right]\otimes\text{vec}\left[F(y)^{(G,l)}_{j}\right]$
(4)

where $G_{r}(y)^{(G,l)}$ is the Gram matrix of $l^{th}$ layer from generator’s
feature maps. The subscript $i$ and $j$ denotes the $i^{th}$ and $j^{th}$ maps
of the $G_{r}(y)^{(G,l)}$, $vec[\cdot]$ represents the vectorization
operation, and $\otimes$ represents the inner product.

$L_{n}$ is computed as follows:

$\mathcal{L}_{n}(G)=\mathbb{E}_{(x,y)\sim(S,T)}\left[\sum_{l=1}^{3}W\left(G_{r}(G(x)),G_{r}(y)\right)\right]$
(5)

Where $W$ is the Wasserstein distance [39], we trained our model with a
feature map of various early layers such as 1, 2, and 3 and decided to keep 3
for final results based on the outcomes. Thus, we used only the first three
layers of the feature map corresponding to noise. The $l$ denotes the layers
in the generator.

####  2.1.4 Objective function

Full objective function is calculated as follows:

$\arg{\min\max}_{G,D_{n},D_{c}}\mathcal{L}(G,D_{n},D_{c})=\mathcal{L}_{an}(G,D_{n})+\mathcal{L}_{ac}(G,D_{c})+\lambda_{1}\cdot\mathcal{L}_{c}(G)+\lambda_{2}\cdot\mathcal{L}_{n}(G)$
(6)

Where $\lambda_{1}$ and $\lambda_{2}$ are the hyperparameters, our aim is to
optimize the above equation.

###  2.2 Implementation Details

####  2.2.1 Network architecture

We adopted the generator and discriminator networks architecture from Jun-Yan
Zhu et al. [16]. The generator is a Resnet-15 network that we modified by
adding one extra convolution block and one deconvolution block to manage the
feature map extraction from the first three layers. The generator network
contains nine residual blocks (including a reflection padding operation, an
instance normalization layer, a 2-D convolution layer, a ReLu activation
function, an instance normalization layer, a 2-D convolution layer, and a plus
operation sequentially), four convolutional blocks (including a ReLu
activation function, an instance normalization layer, and a 2-D convolution
layer) and four deconvolutional blocks (including a ReLu activation function,
an instance normalization layer, and a 2-D deconvolution layer). The structure
of both discriminators is the same. It contains four convolutional blocks
(including a ReLu activation function, an instance normalization layer, and a
2-D convolution layer) and a fully connected layer.

####  2.2.2 Training details

We used Adam optimizer [40] with batch size 1 to train the model to 200
epochs. The learning rate was kept at 0.0002 for the first 100 epochs, then
linearly decreased to zero over the following 100 epochs. The $\lambda_{1}$
and $\lambda_{2}$ values in equation 6 were set to 10 for all experiments. The
number of filters in the generator and discriminator’s first convolutional
layer was 64. Model training was done using two NVIDIA GeForce RTX 2080TI
GPUs.

####  2.2.3 Datasets

We retrospectively included carotid ultrasound 2D images from the longitudinal
VIPVIZA study cohort of Umea University [41]. The study population are healthy
subjects (50:50 men and women, 40-60 years old at base-line examination), with
subclinical atherosclerosis. Expert operators acquired images in a
longitudinal 2D B-mode of the carotid arteries. In the baseline and 3-year
follow-up of the VIPVIZA study, a Cardio Health Station (CHS) and a linear 7
MHz transducer (Panasonic Healthcare Corporation of North America, Newark, NJ,
USA) ultrasound machine (“US System A”) were used. In the 6-year follow-up of
the study, a GE Vivid IQ with the 9 MHz transducer (GE HealthCare, Chicago,
IL, USA) ultrasound machine (“US System B”) was used.

Table 1: Description of the images used in the experiments. Here, GSM, CHS, GE IQ, and VIPVIZA stand for grey scale median, cardio health station, general electric image quality, and västerbotten intervention program visualization of atherosclerosis, respectively. Characteristics | Experiment 1 | Experiment 2
---|---|---
| Domain A | Domain B | Domain C | Domain D
Anatomical content | Carotid plaques | Carotid plaques | Carotid intima-media | Carotid intima-media
Projection | Longitudinal | Longitudinal | Longitudinal | Longitudinal
US system | CHS (sys A) | GE IQ (sys B) | mixed | mixed
Source | VIPVIZA 3y FU | VIPVIZA 6y FU | mixed | mixed
Noise level | mixed | mixed | high | low
N train | 443 | 546 | 225 | 215
N test | 50 | - | 15 | -
Plaque GSM | 44.6 (22.4) | 48.6 (20.3) | - | -

##### Experiment 1 – Image harmonization:

In this experiment, we included carotid plaque images from the two different
systems (System A and System B). The task for the domain adaptation was to
adapt the system A image features to resemble those of system B images while
retaining the anatomical content. For system A, the plaque acquisition mode
was used, producing a B-mode image with no harmonic imaging and spatial
angular compounding of three projections – here defined as “domain A”. Images
were reconstructed offline from the .cine format, and 493 consecutive images
from different subjects were included from system A, and 546 for system B. For
system B, standard B-mode acquisition mode was used, including harmonic
imaging, and images were stored in dicom screen capture format – here defined
as “domain B”. Images had 8-bit intensity depth, and all images were resampled
to 400x400 pixels, corresponding to a spatial resolution on the order of
0.1mm/px.

##### Experiment 2 – Noise reduction:

In this experiment, we included carotid images of the intima-media complex
(without plaques) from both systems A and B. For system A, the CHS IMT mode,
which uses harmonic imaging for intima-media images and measurement, was used.
For system B, conventional B-mode and harmonic imaging were used, and images
were acquired to measure the intima-media thickness. Out of these images,
noisy and clear images were manually classified by an expert operator
(biomedical scientist and sonographer with +10 years of clinical ultrasound) -
noisy and clear images were defined as “domain C” and “domain D”,
respectively. Noisy images were included based on the criteria of visual noise
in the lumen of the images (a haze of bright pixels in the region of the
blood), and clear images as its counterpart – with no noise in the lumen (dark
pixels in the region of the blood). See figure 2 for examples. The task for
the domain adaptation was to adapt the noisy image features to those of the
clear images while retaining the anatomical content, Images had 8-bit
intensity depth, and all images were resampled to 400x400 pixels,
corresponding to a spatial resolution on the order of 0.1mm/px.

###  2.3 Evaluation

The performance of the translation was evaluated from four perspectives: 1)
their feature distribution similarity, 2) pixel-space similarity, 3) down-
stream estimation of an image-based cardiovascular risk marker, and 4)
(reverberation) noise level. For comparison, we also trained the default
cycleGAN model.

The similarity in feature distributions was computed both without any
translation and with translation, using the Bhattacharya distance (BD) and
Histogram correlation (HC) [42] metrics. Mathematically, they are given by:

$BD(h_{1},h_{2})=\sqrt{1-\frac{1}{\sqrt{\bar{h}_{1}\bar{h}_{2}N^{2}}}\sum_{I}\sqrt{h_{1}(I)h_{2}(I)}}$
(7)

$HC(h_{1},h_{2})=\frac{\sum_{I}(h_{1}(I)-\bar{h}_{1})(h_{2}(I)-\bar{h}_{2})}{\sqrt{\sum_{I}(h_{1}(I)-\bar{h}_{1})^{2}\sum_{I}(h_{2}(I)-\bar{h}_{2})^{2}}}$
(8)

Where $h_{1}$ and $h_{2}$ are the normalized histograms of the images.
$\bar{h}_{1}$ and $\bar{h}_{2}$ represents the mean value of the histograms
$h_{1}$ and $h_{2}$, respectively. $N$ is the total number of bins in the
histogram. $I$ represent the background region pixels. High HC and low BD
between images represent high feature distribution similarity.

The pixel-space similarity was computed using the structure similarity index
measure (SSIM) between the input image and the corresponding translated image.
The purpose was to evaluate the preservation of anatomical content. The SSIM
metric is based on different regional pixel statistics in the image, and a
high value means strong similarity. The SSIM was computed for the whole image
and individual tissues by averaging manually segmented regions of interest
(ROIs) of the plaque/intima-media tissues, the lumen, and adventitia (e.g.
figures 5 and 7).

To assess the impact of the translation on image-based risk markers, we
computed the grey scale median (GSM) and plaque vulnerability re-
classification rate. The GSM is a risk marker that is related to plaque
vulnerability [43], where a low value is related to symptomatic (vulnerable)
plaques and a high value to asymptomatic (stable) plaques. One threshold that
has been used is GSM=25 [11]. Using this threshold, we computed the number of
plaques that were re-classified, i.e. crossed this threshold due to the domain
adaptation.

$GSM=\text{median}\left(190\cdot\left(I_{\text{whole\\_image}}-\text{mean}(I_{\text{ROI\\_Lumen}}\right)\right)/\max(I_{\text{ROI\\_Adventitia}})$
(9)

Finally, the noise level was estimated using the lumen and wall ROIs as

$\text{Contrast}=-20\log_{10}\left(\frac{\text{mean}(I_{\text{ROI\\_Lumen}})}{\text{mean}(I_{\text{ROI\\_Adventitia}})}\right)$
(10)

If there is low noise, the lumen should appear dark, and the
$I_{\text{ROI\\_lumen}}$ should have low values, causing the Contrast to be
strongly negative in the dB scale.

##  3 Results

Figure 2: Examples of carotid ultrasound images used in the training of the
domain adaptation task for the two experiments. a) show examples for
experiment 1 (image harmonization) that we define as domains A and B,
respectively. b) show examples for experiment 2 (noise reduction) that we
define as domains C and D. In the lower panel, c) and d), regions of interests
(ROIs) from the upper images are shown to zoom in on a homogeneous tissue
segment (A1 and B1), carotid plaques (A2 and B2), Lumen/blood (C1 and D1), and
arterial walls (C2 and D2).

###  3.1 Experiment 1 – Image harmonization

Figure 3 shows examples of two domain A images (cropped) and their
corresponding translations using cycleGAN and the proposed GAN. It can be seen
that the anatomical content of the input images is retained in the translated
versions, but the fine-detailed texture pattern is modified. Evaluation of
feature distribution similarity showed that both translation models improved
the similarity, but the proposed model had the highest performance (Table 2
and Figure 4). The proposed model achieved a lower BD and higher HC compared
to CycleGAN and no translation, demonstrating that our model outperformed
CycleGAN.

Figure 5 shows an example of an input image, generated translation and the
corresponding pixel-space similarity comparison using SSIM. Evaluation of
pixel-space similarity showed that SSIM had higher values (higher similarity)
for the proposed model than the GAN model in all tissues except the lumen
(Table 3 and Figure 4).

The difference in plaque GSM between input and output images was 7.6 (6.5) vs
15.0 (13.4) for the proposed and cycleGAN models, respectively (Table 4). In
the training data, the difference in GSM plaque between domains B and A was
4.0 (p<0.006, independent samples t-test). The contrast between the arterial
wall (adventitia) and lumen ROIs was -34.1 (3.8) dB, -35.2 (4.1) dB and -35.8
(4.3) dB for the original image, cycleGAN, and proposed model translations
(Table 4).

Figure 3: Example of results for experiment 1 (Image harmonization). a) and d) shows input images (Domain A). b) and e) shows the proposed model’s generated images. c) and f) shows the corresponding generated images of the cycleGAN model. Figure 4: Summary of results. The upper panel is for experiment 1 (image harmonization), and the lower panel is for experiment 2 (noise reduction). The evaluation was done in four perspectives: 1) feature distributions, a) and d), 2) pixel space b) and e), and 3) image-based risk marker for cardiovascular disease (GSM) and 4) noise level (Contrast), c) and f). Feature distributions were compared using bhattacharyya distance (BD) and histogram correlation (HC), and pixel-space similarity was evaluated using the structure similarity index measure (SSIM). The risk marker and contrast measures are presented as the relative change between the corresponding input and domain-adapted output images. Table 2: Evaluation of feature transfer for Experiment 1 (Image Harmonization). To measure the dissimilarity and correlation between two domain images, bhattacharyya distance (BD) and histogram correlation (HC) are used. The results include mean and standard deviation (SD) values of these metrics. Method | Mean(SD)
---|---
| BD | HC
No Translation A vs A | 0.044 (0.024) | 0.958 (0.022)
No Translation B vs B | 0.033 (0.012) | 0.968 (0.011)
No Translation A vs B | 0.120 (0.095) | 0.890 (0.077)
CycleGAN A to B’ | 0.094 (0.065) | 0.912 (0.057)
Proposed model A to B’ | 0.084 (0.047) | 0.920 (0.043)
Figure 5: Illustration of ROIs used in Experiment 1 to assess translation performance (Lumen, Plaque and Wall/Adventitia). a) shows an example image of domain A, b) the generated image (domain B) and c) the corresponding structure similarity index measure (SSIM) map. Table 3: Performance evaluation for pixel-space similarity for whole image and individual tissues - Experiment 1 (Image harmonization). Model | Tissue segment
---|---
| Plaque | Lumen | Adventitia | Whole Image
Proposed GAN | 0.63 (0.08) | 0.31 (0.15) | 0.71 (0.09) | 0.53 (0.09)
CycleGAN | 0.46 (0.08) | 0.33 (0.18) | 0.59 (0.12) | 0.41 (0.03)
Table 4: Performance evaluation - Impact on risk markers and blood-to-tissue
contrast - Experiment 1 (Image harmonization). GSM stands for grey scale
median.

| A | B’ | Diff | p
---|---|---|---|---
| GSM |  |  |
Proposed GAN | 42.7 (25.0) | 50.3 (26.2) | 7.6 (6.5) | P<0.001 ***
CycleGAN | 42.7 (25.0) | 57.7 (31.8) | 15.0 (13.4) | P<0.001 ***
| Re-classification |  |  |
Proposed GAN | 5 (10%) |  |  |
CycleGAN | 12 (24%) |  |  |
| Contrast (dB) |  |  |
Proposed GAN | -34.1 (3.8) | -35.2 (4.1) |  | P=0.78
CycleGAN | -34.1 (3.8) | -35.8 (4.3) |  | P=0.73

Plaque area = 24.9 (18.0) mm2

###  3.2 Experiment 2 – Noise reduction

Figure 6 shows examples of two noisy domain C images (cropped) and their
corresponding translations using cycleGAN and the proposed GAN, where noise in
the lumen region was reduced.

The feature similarity distribution metrics between input and translated
images showed that both models improved the BD and HC similarities as compared
to no translation and that the proposed model was the best (Table 5, Figure
4). Figure 7 shows example results for noise suppression in carotid ultrasound
images by SSIM map and illustrates that the translation models strongly modify
the lumen ROI while retaining the overall anatomy. The pixel-space similarity
estimations (SSIM) in different tissues also showed that the lumen was
strongly dissimilar with SSIM values close to 0, whereas tissues had strong
similarity 0.7-0.8 for both models.

Evaluation of the influence of translation on risk marker showed that the
intima-media GSM changed from X to Y and Z for the cycleGAN and proposed GAN
models, respectively. The contrast between adventitia and lumen ROIs changed
from -23.5 (3.2) dB to -35.7 (2.8) dB and -46.7 (18.1) dB for the input image
and the cycleGAN and proposed GAN translations, respectively. Note: The
proposed model had the most significant contrast, indicating the highest noise
reduction.

Figure 6: Example of results for experiment 2 (Noise reduction). a) and d) shows input images (Domain C) and b) and e) shows the proposed model’s generated images. c) and f) show the corresponding generated images of the cycleGAN model. Table 5: Evaluation of feature transfer for Experiment 2 (Noise Reduction). To measure the dissimilarity and correlation between two domain images, bhattacharyya distance (BD) and histogram correlation (HC) are used. The results include mean and standard deviation (SD) values of these metrics. Method | Mean(SD)
---|---
| BD | HC
No Translation C vs C | 0.024 (0.008) | 0.976 (0.007)
No Translation D vs D | 0.015 (0.008) | 0.985 (0.007)
No Translation C vs D | 0.357 (0.152) | 0.707 (0.098)
CycleGAN C to D’ | 0.181 (0.112) | 0.839 (0.091)
Proposed model C to D’ | 0.172 (0.075) | 0.844 (0.062)
Figure 7: Illustration of ROIs used in Experiment 2 to assess translation performance (Lumen, inner arterial wall (intima media) and outer wall (adventitia)). a) shows an example image of domain C, b) the generated image (domain D) and c) the corresponding structure similarity index measure (SSIM) map. Table 6: Performance evaluation for pixel-space similarity for whole image and individual tissues - Experiment 2 (Noise reduction). Model | Tissue segment
---|---
| IM | Lumen | Adventitia | Whole Image
Proposed GAN | 0.74 (0.16) | 0.03 (0.03) | 0.80 (0.08) | 0.36 (0.04)
CycleGAN | 0.71 (0.15) | 0.06 (0.02) | 0.72 (0.08) | 0.43 (0.05)
Table 7: Performance evaluation - Impact on riskmarkers and blood-to-tissue contrast - Experiment 2 (Noise reduction). GSM stands for grey scale median. Method | C | D’ | Diff | p
---|---|---|---|---
| GSM |  |  |
Proposed GAN | 51.2 (13.9) | 45.8 (17.7) | -5.3 (11.8) | P=0.10
CycleGAN | 51.2 (13.9) | 42.6 (18.1) | -8.6 (6.2) | P<0.001 ***
| Contrast (dB) |  |  |
Proposed GAN | -23.5 (3.2) | -46.7 (18.1) |  | P<0.001 ***
CycleGAN | -23.5 (3.2) | -35.7 (2.75) |  | P<0.001 ***

##  4 Discussion

In this work, we proposed a GAN-based model for multi-purpose domain
adaptation tasks and tested its performance for image harmonization between
carotid ultrasound images captured with different ultrasound systems and image
noise reduction in carotid ultrasound images with high noise. The main results
show that 1) the model can adapt the image features between two different
ultrasound systems, as well as reduce the noise in noisy images while still
retaining the overall anatomy (content) in the images, and 2) the domain
adaptation influenced critical risk markers for cardiovascular risk computed
on the input and output images.

###  4.1 Image harmonization

Results showed adapted feature distributions compared to those without
adaptation (0.12 vs 0.89) but not as similar to those within domains (0.03 vs
0.97) for BD and HC, respectively (Table 2). The feature similarities were
higher than that of the NAGAN model [13], where they found 0.32 and 0.58 for
BD and HC and 0.36 and 0.58 for the CycleGAN model. However, it should be
noted that they adapted OCT and US images and, therefore, have larger
differences between the image domains. The underlying content (anatomy) of the
images was retained, as shown by the SSIM values (Table 3). For the whole
image, SSIM was 0.53 (0.09). This is lower than those found by Ali et al. [9],
translating US images of healthy and diseased arteries, with an overall SSIM
of 0.78 (0.02). However, those images only included the common carotid artery,
where the anatomy has less variation than the bulb. The bulb is the position
where the carotid artery divides into two arteries, the internal and external
carotids. The latter was often included in the images of this study since this
is where the plaques most commonly form. In addition, it should be noted that
the SSIM were higher in our results for plaque and adventitia arterial wall
segments SSIM were 0.63 and 0.71 as compared to the overall image - and these
tissues are the ones that are used for assessing cardiovascular disease [43].
Thus, a lower preservation of the overall details of the image anatomy is
likely not critical for our application.

###  4.2 Noise reduction

The primary source of noise in the ultrasound images is caused by diffuse
reverberation which is caused by the multiple reflections of the transmitted
ultrasound pulse in the tissue before returning to the transducer. The results
of this work showed adapted feature distributions compared to those without
adaptation (0.36 vs 0.71) but not as similar to those within domains (0.02 vs
0.98) for BD and HC, respectively (Table 5). The noise reduction in the image
can be clearly seen in Figure 6 and is supported by the low SSIM values in the
lumen (0.03), while the arterial wall tissues are unaffected (SSIM 0.74 and
0.80), as shown in Table 6). In addition, the noise reduction was significant
for both models, but the proposed model was better, as seen in the Contrast
measurements changing from -23.5 dB to -46.7 dB and -35.7 dB for the proposed
vs CycleGAN models, respectively (Table 7). Recent work by other groups, e.g.,
by Brickson et al. [44], used a 3D CNN model on the underlying radio-frequency
data of ultrasound images and showed a Contrast reduction from -18.9 (5.8) to
-34.9 (8.1) dB. Thus, noise reduction performance in our work is similar and
potentially better. The proposed GAN preserved the wall and adventitia tissues
the best, as seen in SSIM for wall and adventitia (Table 6). It should be
noted that it slightly distorted other parts of the image (see Figure 7),
e.g., muscle and vein, which is also indicated by a slightly lower overall
SSIM compared to experiment 1. However, this may not be a critical problem for
our application of assessing cardiovascular disease (as commented for
experiment 1).

###  4.3 Influence on image-based computation of cardiovascular risk markers

In experiment 1, the impact of the domain adaptations on image-based risk
markers showed that both models increased the GSM of the plaques, with a
consequence of re-classification of 12 plaques (24%) for cycleGAN and 5 (10%)
plaques for the proposed model. Thus, the proposed model preserves the
computed risk marker integrity better. The SSIM was different in the lumen,
adventitia, and plaque tissues, indicating that the model translated them
differently (Eq 9). Since the GSM risk marker is computed from ROIs in the
lumen, adventitia and arterial plaque, the former finding can likely explain
this change in GSM. In comparison, our group previously found that the blood
pressure variation over the cardiac cycle caused a 16% re-classification in
the same study population [45], and 30% re-classification in another study
population with asymptomatic atherosclerosis (more severe disease) [46]. This
indicates that the re-classification by our model is similar to or lower than
other known modulators of GSM.

In contrast, in experiment 2, the risk marker IM-GSM was used to assess the
impact on clinical assessment, and it showed a slight decrease in the IM-GSM
for both models but with large variability. The IM-GSM risk marker has not
been as extensively studied in the literature as the GSM for plaques, and its
relevance for assessing cardiovascular risk remains to be established.
However, IM-GSM has been shown to be associated with the progression of
atherosclerotic disease in the same study population [47] and to be associated
with plaque GSM in populations with more advanced disease [48]. Since the GSM
computations in experiments 1 and 2 use the same tissue ROIs (Lumen and
adventitial wall), and since the lumen was much more modified in experiment 2
compared to experiment 1, GSM estimations will be influenced differently in
the two experiments. It should be noted that the increase in GSM of the
plaques may be built into the models because the training target domain B had
significantly higher GSM values of the plaques than source domain A (48 vs 44,
see Table 1). However, based on these results, we propose that evaluation of
influence on risk marker estimations should always be done and possibly
adjusted for when applying domain adaptation models to carotid US images since
our results show that different models gave different changes in the risk
marker.

###  4.4 The proposed model

The proposed model addressed both image harmonization and noise reduction, and
results show that it performed better than the state-of-the-art model,
cycleGAN, in both tasks. The proposed model was inspired by the methods [37,
13], demonstrating the independent processing of feature representations
learned by neural networks. In our model, we leverage the potential of
Wasserstein distance to compute the noise loss using the style of feature maps
driven from the early layers of the generator and the target domain images.
The Wasserstein distance reduces dissimilarity by finding the most efficient
way to map the distribution of one image set to another and accounting for the
spatial relationships within the images by considering the geometry of the
data in images[39, 49]. We consider that this makes it more suitable for
comparing images where the spatial relationship between pixels is crucial,
such as in medical images (e.g., carotid ultrasound images). Our features
transfer approach builds a domain-to-domain mapping for both image
harmonization and noise reduction.

We conducted experiments with both Resnet and Unet architectures as generators
and observed that Resnet outperformed Unet in our applications. Therefore, the
final results are reported based on the use of Resnet. Our generator employs a
ResNet-15 architecture, while the discriminator is constructed using a
classical CNN. The ResNet holds an advantage over CNN due to its incorporation
of residual units, enabling deeper layers that directly leverage information
from shallower ones. Consequently, feature maps extracted from the ResNet-15
are likely to exhibit more refined patterns and smoother transitions compared
to those from a standard CNN. In addition, we fine-tuned the generator
architecture from Jun-Yan Zhu et al. [16] by adding two extra blocks (one
convolution block and one deconvolution block) to manage the feature map
extraction from the first three layers. We trained our model for various
layers and decided to use the first three layers of the generator for
stylistic feature extraction based on the training outcomes.

###  4.5 Limitations and future work

The proposed GAN model is designed for unsupervised image-to-image translation
between two specific domains. It may not generalize well to tasks involving
other domains or tasks with highly diverse datasets.

The features of carotid ultrasound images are influenced both by the specific
imaging system (producing e.g. different spatial resolution and noise levels)
as well as the severity of the disease (resulting in changes in the arterial
wall, including plaques of various sizes and composition). While the present
work and previous work by others show that adaptation between two different
imaging systems can be done, the population is subclinical, meaning that the
amount of disease is relatively low, and it is not clear how different amounts
of disease in the arteries affect the translation performance. A potential
solution to this could be to, e.g. train the model with a dataset including a
larger distribution of the amount of disease in the arteries.

The results of the image harmonization task are between two specific
ultrasound systems and will not, by default, generalize to other US systems.
For other pairs of US systems, the model needs to be retained. The noise
reduction task was trained using images from both systems A and B, potentially
allowing the model to learn more robust features (i.e., features not related
to a specific system). Therefore, this model could potentially generalize to
some degree to images of other ultrasound systems. The next step should be to
train it on a larger set of images with low and high image quality and from
multiple systems and use an external cohort for validation.

Moreover, in future work, the models could be trained using the risk marker
data as input to the discriminator to minimize the change in risk markers that
we found in this work.

The implementation of the model is a bit complex compared to traditional GAN
architectures because it involves modifying the generator network to adjust
harmonization features and noise levels between domains adaptively. This
increased complexity can lead to higher computational overhead during training
and inference.

##  5 Conclusion

This paper proposes a GAN-based domain adaptation model for image
harmonization and noise reduction of carotid ultrasound images. The results
show that the model can 1) adapt features of ultrasound images of one system
to images from another and 2) reduce the noise in the lumen (blood), both
without changing the anatomical contents of the images. The proposed model
performed better than the state-of-the-art model, CycleGAN. In addition, the
results showed that the downstream task of computing image-based risk markers
for cardiovascular risk assessment was, in general, influenced by the domain
adaptations. We conclude that domain adaptation is a powerful tool to improve
image quality and harmonization of image features in carotid ultrasound
images, but the influence on computed risk markers should be evaluated for the
specific model used.

##### Acknowledgment:

This work was supported by the Heart Foundation of Northern Sweden, the Kempe
foundations (JCK-3172), and Region Västerbotten.

## References

  * [1] Andrew Ward, Ashish Sarraju, Sukyung Chung, Jiang Li, Robert Harrington, Paul Heidenreich, Latha Palaniappan, David Scheinker, and Fatima Rodriguez.  Machine learning and atherosclerotic cardiovascular disease risk prediction in a multi-ethnic population.  NPJ digital medicine, 3(1):125, 2020.
  * [2] Zihao Fan, Zhi Du, Jinrong Fu, Ying Zhou, Pengyu Zhang, Chuning Shi, and Yingxian Sun.  Comparing the performance of machine learning and conventional models for predicting atherosclerotic cardiovascular disease in a general chinese population.  BMC Medical Informatics and Decision Making, 23(1):134, 2023.
  * [3] Fisnik Jashari, Pranvera Ibrahimi, Gani Bajraktari, Christer Grönlund, Per Wester, and Michael Y Henein.  Carotid plaque echogenicity predicts cerebrovascular symptoms: a systematic review and meta-analysis.  European journal of neurology, 23(7):1241–1247, 2016.
  * [4] Xiao-xiao Chen, Zi-xiang Kong, Shu-fang Wei, Fei Liang, Ting Feng, Shan-shan Wang, and Jian-song Gao.  Ultrasound lmaging-vulnerable plaque diagnostics: Automatic carotid plaque segmentation based on deep learning.  Journal of Radiation Research and Applied Sciences, 16(3):100598, 2023.
  * [5] Mehreen Mubashar, Hazrat Ali, Christer Grönlund, and Shoaib Azmat.  R2u++: a multiscale recurrent residual u-net with dense skip connections for medical image segmentation.  Neural Computing and Applications, 34(20):17723–17739, 2022.
  * [6] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [7] Zhihao Wang, Jian Chen, and Steven CH Hoi.  Deep learning for image super-resolution: A survey.  IEEE transactions on pattern analysis and machine intelligence, 43(10):3365–3387, 2020.
  * [8] Hao Guan and Mingxia Liu.  Domain adaptation for medical image analysis: a survey.  IEEE Transactions on Biomedical Engineering, 69(3):1173–1185, 2021.
  * [9] Hazrat Ali, Emma Nyman, Ulf Näslund, and Christer Grönlund.  Translation of atherosclerotic disease features onto healthy carotid ultrasound images using domain-to-domain translation.  Biomedical Signal Processing and Control, 85:104886, 2023.
  * [10] Tomoyuki Fujioka, Kazunori Kubota, Mio Mori, Leona Katsuta, Yuka Kikuchi, Koichiro Kimura, Mizuki Kimura, Mio Adachi, Goshi Oda, Tsuyoshi Nakagawa, et al.  Virtual interpolation images of tumor development and growth on breast ultrasound image synthesis with deep convolutional generative adversarial networks.  Journal of Ultrasound in Medicine, 40(1):61–69, 2021.
  * [11] Christina I Christodoulou, Constantinos S Pattichis, Marios Pantziaris, and Andrew Nicolaides.  Texture-based classification of atherosclerotic carotid plaques.  IEEE transactions on medical imaging, 22(7):902–912, 2003.
  * [12] Geert Litjens, Thijs Kooi, Babak Ehteshami Bejnordi, Arnaud Arindra Adiyoso Setio, Francesco Ciompi, Mohsen Ghafoorian, Jeroen Awm Van Der Laak, Bram Van Ginneken, and Clara I Sánchez.  A survey on deep learning in medical image analysis.  Medical image analysis, 42:60–88, 2017.
  * [13] Tianyang Zhang, Jun Cheng, Huazhu Fu, Zaiwang Gu, Yuting Xiao, Kang Zhou, Shenghua Gao, Rui Zheng, and Jiang Liu.  Noise adaptation generative adversarial network for medical image analysis.  IEEE transactions on medical imaging, 39(4):1149–1159, 2019.
  * [14] Zixia Zhou, Yuanyuan Wang, Yi Guo, Yanxing Qi, and Jinhua Yu.  Image quality improvement of hand-held ultrasound devices with a two-stage generative adversarial network.  IEEE Transactions on Biomedical Engineering, 67(1):298–311, 2019.
  * [15] Menghua Xia, Hongbo Yang, Yanan Qu, Yi Guo, Guohui Zhou, Feng Zhang, and Yuanyuan Wang.  Multilevel structure-preserved gan for domain adaptation in intravascular ultrasound analysis.  Medical Image Analysis, 82:102614, 2022.
  * [16] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In Proceedings of the IEEE international conference on computer vision, pages 2223–2232, 2017.
  * [17] Long Teng, Zhongliang Fu, and Yu Yao.  Interactive translation in echocardiography training system with enhanced cycle-gan.  IEEE access, 8:106147–106156, 2020.
  * [18] Zizhao Zhang, Lin Yang, and Yefeng Zheng.  Translating and segmenting multimodal medical volumes with cycle-and shape-consistency generative adversarial network.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 9242–9251, 2018.
  * [19] Mingsheng Long, Han Zhu, Jianmin Wang, and Michael I Jordan.  Deep transfer learning with joint adaptation networks.  In International conference on machine learning, pages 2208–2217. PMLR, 2017.
  * [20] Shreeram Athreya, Ashwath Radhachandran, Vedrana Ivezić, Vivek Sant, Corey W Arnold, and William Speier.  Ultrasound image enhancement using cyclegan and perceptual loss.  arXiv preprint arXiv:2312.11748, 2023.
  * [21] Lihong Huang, Zixia Zhou, Yi Guo, and Yuanyuan Wang.  A stability-enhanced cyclegan for effective domain transformation of unpaired ultrasound images.  Biomedical Signal Processing and Control, 77:103831, 2022.
  * [22] Chaoqun Tan, Mingming Yang, Zhisheng You, Hu Chen, and Yi Zhang.  A selective kernel-based cycle-consistent generative adversarial network for unpaired low-dose ct denoising.  Precision Clinical Medicine, 5(2):pbac011, 2022.
  * [23] Chao Tang, Jie Li, Linyuan Wang, Ziheng Li, Lingyun Jiang, Ailong Cai, Wenkun Zhang, Ningning Liang, Lei Li, Bin Yan, et al.  Unpaired low-dose ct denoising network based on cycle-consistent generative adversarial network with prior image information.  Computational and mathematical methods in medicine, 2019, 2019.
  * [24] Feixiang Zhao, Mingzhe Liu, Zhihong Gao, Xin Jiang, Ruili Wang, and Lejun Zhang.  Dual-scale similarity-guided cycle generative adversarial network for unsupervised low-dose ct denoising.  Computers in Biology and Medicine, 161:107029, 2023.
  * [25] T Elatrozy, A Nicolaides, TH Tegos, AZ Zarka, et al.  The effect of b-mode ultrasonic image standardisation on the echodensity of symptomatic and asymptomatic carotid bifurcation plaques.  International Angiology, 17(3):179, 1998.
  * [26] Stavros K Kakkos, Andrew N Nicolaides, Efthyvoulos Kyriacou, Constantinos S Pattichis, and George Geroulakos.  Effect of zooming on texture features of ultrasonic images.  Cardiovascular Ultrasound, 4:1–10, 2006.
  * [27] Mostafa Sharifzadeh, Ali KZ Tehrani, Habib Benali, and Hassan Rivaz.  Ultrasound domain adaptation using frequency domain analysis.  In 2021 IEEE International Ultrasonics Symposium (IUS), pages 1–4. IEEE, 2021.
  * [28] Wing Keung Cheung, Dorothy M Gujral, Benoy N Shah, Navtej S Chahal, Sanjeev Bhattacharyya, David O Cosgrove, Robert J Eckersley, Kevin J Harrington, Roxy Senior, Christopher M Nutting, et al.  Attenuation correction and normalisation for quantification of contrast enhancement in ultrasound images of carotid arteries.  Ultrasound in Medicine & Biology, 41(7):1876–1883, 2015.
  * [29] Christos P Loizou, Constantinos S Pattichis, Marios Pantziaris, T Tyllis, and Andrew Nicolaides.  Quality evaluation of ultrasound imaging in the carotid artery based on normalization and speckle reduction filtering.  Medical and Biological Engineering and Computing, 44:414–426, 2006.
  * [30] Mirko Kemna, Alexander Heinlein, and Cornelis Vuik.  Reduced order fluid modeling with generative adversarial networks.  PAMM, 23(1):e202200241, 2023.
  * [31] Hee Guan Khor, Guochen Ning, Xinran Zhang, and Hongen Liao.  Ultrasound speckle reduction using wavelet-based generative adversarial network.  IEEE Journal of Biomedical and Health Informatics, 26(7):3080–3091, 2022.
  * [32] Mercy N Asiedu, Alex R Benjamin, Vivek K Singh, Shuhang Wang, Kevin Wu, Anthony E Samir, and Viksit S Kumar.  A generative adversarial network for ultrasound signal enhancement by transforming low-voltage beamformed radio frequency data to high-voltage data.  In Medical Imaging 2022: Ultrasonic Imaging and Tomography, volume 12038, pages 246–254. SPIE, 2022.
  * [33] Silvia Seoni, Massimo Salvi, Giulia Matrone, and Kristen M Meiburger.  Ultrasound image beamforming optimization using a generative adversarial network.  In 2022 IEEE International Ultrasonics Symposium (IUS), pages 1–4. IEEE, 2022.
  * [34] Lun Zhang and Junhua Zhang.  Ultrasound image denoising using generative adversarial networks with residual dense connectivity and weighted joint loss.  PeerJ Computer Science, 8:e873, 2022.
  * [35] Min Meng, Mengcheng Lan, Jun Yu, Jigang Wu, and Ligang Liu.  Dual-level adaptive and discriminative knowledge transfer for cross-domain recognition.  IEEE Transactions on Multimedia, 2022.
  * [36] Xiubin Zhu, Dan Wang, Witold Pedrycz, and Zhiwu Li.  Transfer learning realized with nonlinearly transformed input space.  IEEE Transactions on Emerging Topics in Computing, 2022.
  * [37] Leon A Gatys, Alexander S Ecker, and Matthias Bethge.  Image style transfer using convolutional neural networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2414–2423, 2016.
  * [38] Yanghao Li, Naiyan Wang, Jiaying Liu, and Xiaodi Hou.  Demystifying neural style transfer.  arXiv preprint arXiv:1701.01036, 2017.
  * [39] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR, 2017.
  * [40] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [41] Ulf Näslund, Nawi Ng, Anna Lundgren, Eva Fhärm, Christer Grönlund, Helene Johansson, Bernt Lindahl, Bertil Lindahl, Kristina Lindvall, Stefan K Nilsson, et al.  Visualization of asymptomatic atherosclerotic disease for optimum cardiovascular prevention (vipviza): a pragmatic, open-label, randomised controlled trial.  The Lancet, 393(10167):133–142, 2019.
  * [42] Anil Bhattacharyya.  On a measure of divergence between two statistical populations defined by their probability distribution.  Bulletin of the Calcutta Mathematical Society, 35:99–110, 1943.
  * [43] AN Nicolaides, SK Kakkos, E Kyriacou, M Griffin, M Sabetai, DJ Thomas, T Tegos, G Geroulakos, N Labropoulos, CJ Doré, et al.  Asymptomatic carotid stenosis and risk of stroke (acsrs) study group. asymptomatic internal carotid artery stenosis and cerebrovascular risk stratification.  J Vasc Surg, 52(6):1486–1496, 2010.
  * [44] Leandra L. Brickson, Dongwoon Hyun, Marko Jakovljevic, and Jeremy J. Dahl.  Reverberation noise suppression in ultrasound channel signals using a 3d fully convolutional neural network.  IEEE Transactions on Medical Imaging, 40(4):1184–1195, 2021.
  * [45] Emma Nyman, Per Lindqvist, Ulf Näslund, and Christer Grönlund.  Risk marker variability in subclinical carotid plaques based on ultrasound is influenced by cardiac phase, echogenicity and size.  Ultrasound in Medicine & Biology, 44(8):1742–1750, 2018.
  * [46] Isak Stenudd, Elias Sjödin, Emma Nyman, Per Wester, Elias Johansson, and Christer Grönlund.  Ultrasound risk marker variability in symptomatic carotid plaque: impact on risk reclassification and association with temporal variation pattern.  The International Journal of Cardiovascular Imaging, 36:1061–1068, 2020.
  * [47] Emma Nyman, Per Liv, Per Wester, Ulf Näslund, and Christer Grönlund.  Carotid wall echogenicity at baseline associates with accelerated vascular aging in a middle-aged population.  The International Journal of Cardiovascular Imaging, 39(3):575–583, 2023.
  * [48] Pranvera Ibrahimi, Fisnik Jashari, Elias Johansson, Christer Grönlund, Gani Bajraktari, Per Wester, and Michael Y Henein.  Common carotid intima–media features determine distal disease phenotype and vulnerability in asymptomatic patients.  International journal of cardiology, 196:22–28, 2015.
  * [49] Yonatan Dukler, Wuchen Li, Alex Lin, and Guido Montúfar.  Wasserstein of wasserstein loss for learning generative models.  In International conference on machine learning, pages 1716–1725. PMLR, 2019.


# Coherent emission from QED cascades in pulsar polar caps

Fábio Cruz GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior
Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal  Thomas Grismayer
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico,
Universidade de Lisboa, 1049-001 Lisboa, Portugal  Alexander Y. Chen JILA,
University of Colorado Boulder, 440 UCB, Boulder, CO 80309, USA  Anatoly
Spitkovsky Department of Astrophysical Sciences, Princeton University,
Princeton, NJ 08544, USA  Luis O. Silva GoLP/Instituto de Plasmas e Fusão
Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa,
Portugal

###### Abstract

Pulsar magnetospheres are thought to be filled with electron-positron plasma
generated in pair cascades. The driving mechanism of these cascades is the
emission of gamma-ray photons and their conversion into pairs via Quantum
Electrodynamics (QED) processes. In this work, we present 2D particle-in-cell
simulations of pair cascades in pulsar polar caps with realistic magnetic
field geometry that include the relevant QED processes from first principles.
Our results show that, due to variation of magnetic field curvature across the
polar cap, pair production bursts self-consistently develop an inclination
with respect to the local magnetic field that favors the generation of
coherent electromagnetic modes with properties consistent with pulsar radio
emission. We show that this emission is peaked along the magnetic axis and
close to the polar cap edge and may thus offer an explanation for the core and
conal components of pulsar radio emission.

Pulsars (1306) — Radio sources (1358) — Plasma astrophysics (1261)

##  1 Introduction

Cascades of electron-positron pairs are a key source of plasma in pulsar
magnetospheres. They result from a positive feedback loop that develops in
vacuum gaps, regions of unscreened electric field in the magnetosphere. Polar
caps have been proposed to host rotation-induced vacuum gaps (Sturrock, 1971;
Ruderman & Sutherland, 1975), where TeV energy electrons and positrons emit
gamma-ray photons via curvature radiation and these are absorbed in the local
$\sim 10^{12}$G magnetic field to produce new pairs. Cascades cease when the
fresh pair plasma screens the vacuum gap electric field. The plasma is then
advected into the magnetosphere, the gap reopens and a new cascade begins
(e.g., Timokhin, 2010). The time-dependent dynamics of polar cap vacuum gaps
has recently been proposed as a primary ingredient to explain the nature of
pulsar radio emission (Beloborodov, 2008; Timokhin & Arons, 2012; Philippov et
al., 2020; Melrose et al., 2020).

Kinetic plasma simulations are ideally suited to study the highly nonlinear
interplay between time-dependent pair cascades and coherent plasma processes.
Timokhin (2010) presented the first 1D particle-in-cell (PIC) simulations
including curvature radiation and pair production from first principles. These
simulations were performed in a frame corotating with the neutron star, which
embeds the magnetic field twist imposed at macroscopic scales (Bai &
Spitkovsky, 2010) via a background current density
$j_{\text{m}}=\alpha\rho_{\text{GJ}}c$, where $\alpha$ is a constant of order
unity, $\rho_{\text{GJ}}$ is the Goldreich-Julian (Goldreich & Julian, 1969)
charge density and $c$ is the speed of light. They showed that cascades
develop over regular short bursts followed by long quiet phases during which
no pairs are produced. As the gap is locally screened, electric field
oscillations are inductively driven due to collective kinetic-scale plasma
motions (Levinson et al., 2005; Cruz et al., 2021). Similar simulations have
been used by Timokhin & Arons (2012); Timokhin & Harding (2015) to determine
the spectra and multiplicity of the pair plasma created in cascades for a
variety of initial conditions.

Philippov et al. (2020, hereafter PTS20) have recently resorted to an
heuristic description of the emission and pair production processes to perform
the first 2D Cartesian PIC simulations of pair cascades and identified a new
process of coherent radiation emission. The key ingredient for this process is
a finite angle between the pair production front and the background magnetic
field, which cannot be captured in 1D. While PTS20 have been able to
demonstrate this coherent emission mechanism using a simplified configuration,
the shape of the pair production front and the spectrum and Poynting flux
profile of the emitted radiation in realistic field geometries depend on the
microscopic details of the emission and pair production cross sections. First-
principles simulations of polar cap pair cascades are very challenging
computationally, due to i) the rapid and localized creation of a large number
of particles and ii) the large separation between gap and plasma kinetic
scales. While the polar cap vacuum gap extends for $\sim 100$ m (Ruderman &
Sutherland, 1975), the electron skin depth associated with the dense plasma
produced in the cascade is $\sim 1$ cm. This scale separation is a consequence
of the large multiplicity of the cascade process that can be estimated as the
ratio between the energies of primary and secondary particles, respectively
$\varepsilon_{\pm}/m_{e}c^{2}\sim 10^{7}$ and
$\varepsilon_{\pm}^{\prime}/m_{e}c^{2}\sim 10^{2}$, where $m_{e}$ is the
electron mass. In multidimensional configurations, the difficulties of
simulating such scale disparity are heightened, and first-principles numerical
models of the pulsar polar cap have not yet been possible.

In this Letter, we adopt a first-principles rescaling of the Quantum
Electrodynamics (QED) processes responsible for gamma-ray and pair production
processes in 2D axisymmetric PIC simulations of pulsar polar caps with a
realistic field geometry. Using these simulations, we determine the
multidimensional properties of pair cascades and their observational
signatures.

##  2 QED-PIC Simulations

We perform 2D QED-PIC simulations with the code OSIRIS (Fonseca et al., 2002,
2008) in axisymmetric cylindrical coordinates $(R,z)$. The lower $z$ boundary
($z=0$) is a rotating conducting disk of radius $R_{0}$. The reference scale
$R_{0}$ should be interpreted as the polar cap radius, and the $z$ axis as the
magnetic (and rotation) axis of the neutron star. The rotation of the
conducting disk is imposed by forcing the radial electric field at this
boundary to be $E_{R}(R,z=0)=E_{0}\times(R/R_{0})\times g(R/R_{0})$, where
$E_{0}$ is a constant proportional to the angular velocity of the disk
$\Omega$ and $g(x)=0.5\times(1-\tanh((x-1)/0.2))$ is a smooth cutoff function
at $x\simeq 1$. This boundary condition induces a unipolar electric field near
the conductor, forcing plasma in this region into corotation with the disk.
The upper $z$ boundary is open for fields, such that any incident wave escapes
freely from this boundary. Both $z$ boundaries are open for particles.

The simulation domain is also permeated by the externally imposed dipolar
magnetic field $\mathbf{B_{d}}$ with magnitude $B_{0}$ at $R=z=0$. The center
of the magnetic dipole is at $(R,z)=(0,-R_{*})$, with $R_{*}/R_{0}=10$. In
analogy with a realistic polar cap, all magnetic field lines that cross the
lower $z$ boundary at a radius $R<R_{0}$ are open to infinity, and those that
cross this boundary at $R>R_{0}$ are closed on the other hemisphere of the
neutron star. We assume that closed field lines are filled with dense plasma
that we model as an ideal conductor, i.e., we set
$\mathbf{E}_{\parallel}=(\mathbf{E}\cdot\mathbf{B_{d}})\mathbf{B_{d}}/|\mathbf{B_{d}}|^{2}$
in this region to zero. This simulation setup forces a current to be driven
along open field lines while providing a (virtual) return current along the
edge of the conducting disk, mimicking the local conditions in pulsar polar
caps.

The QED processes governing pair cascades in our simulations are photon
emission via nonlinear Compton scattering (Erber, 1966), the QED equivalent of
curvature radiation for classical emission from ultra-relativistic particles
(Kelner et al., 2015; Del Gaudio, 2020), and multiphoton Breit-Wheeler pair
production (Ritus, 1985). In our simulations, these processes are included
using Monte-Carlo methods (Grismayer et al., 2016, 2017). In each time step,
we first compute the quantum parameter $\chi_{\pm,\gamma}$ of each particle
(subscripts $\pm$, $\gamma$ correspond to electrons, positrons and photons,
respectively), defined as
$\chi_{\pm,\gamma}=\sqrt{(p_{\mu}F^{\mu\nu})^{2}}/(B_{Q}m_{e}c^{2})$, where
$p^{\mu}$ is the four-momentum of the particle, $F^{\mu\nu}$ is the
electromagnetic tensor and $B_{Q}\simeq 4.4\times 10^{13}$ G is the Schwinger
field. We then rescale
$\chi_{\pm}\to\chi_{\pm}^{\prime}\equiv\zeta_{\pm}\chi_{\pm}$ and
$\chi_{\gamma}\to\chi_{\gamma}^{\prime}\equiv\zeta_{\gamma}\chi_{\gamma}$,
with $\zeta_{\pm,\gamma}\gg 1$, and use $\chi_{\pm,\gamma}^{\prime}$ to
evaluate the probability of creating new particles. The rescaling is done to
reduce the scale separation described above. A key aspect of our approach is
that it retains the fundamental properties of the QED processes: i) the
spectrum of photons emitted via nonlinear Compton scattering is preserved (in
particular its dependence on the energy and curvature of the trajectory of the
emitting particle), and ii) the probability of pair production critically
depends on the angle between the photon propagation direction and the local
electromagnetic field.

Figure 1:  Electron/positron/photon densities (a-c, respectively) and electric
field component parallel to the background magnetic field (d) at a time
$tc/R_{0}=5.4$, where two bursts are visible. All panels show the closed field
line region displayed in white and the magnetic field lines in solid (white in
a-c, black in d) lines. A visual aid identifying the pair production bursts
(i.e., outlining large density regions) is shown in dashed lines in all
panels.

The surface electric field $E_{0}=\Omega B_{0}R_{*}/c$ and the rescaling
parameters $\zeta_{\pm,\gamma}$ are adjusted such that a lepton is accelerated
to the radiation reaction limited Lorentz factor (Daugherty & Harding, 1982)

$\hat{\gamma}_{\pm}=\left(\frac{3}{2\alpha_{\text{fs}}}\frac{E_{0}}{B_{Q}}\right)^{1/4}\left(\frac{\lambdabar_{\text{C}}}{\rho}\right)^{-1/2}\zeta_{\pm}^{-1/2}$
(1)

over a distance

$\ell_{a}\simeq\hat{\gamma}_{\pm}\left(\frac{B_{Q}}{E_{0}}\right)\lambdabar_{\text{C}}\
,$ (2)

that we take to be $\sim 0.1R_{0}$. In equations (1) and (2),
$\alpha_{\text{fs}}$ is the fine-structure constant, $\lambdabar_{\text{C}}$
is the reduced Compton wavelength and $\rho\simeq R_{*}^{2}/R_{0}$ is the
radius of curvature of the last open field line. We fix
$\zeta_{\pm,\gamma}=10^{3}$, $eB_{0}R_{0}/m_{e}c^{2}=10^{8}$ and
$B_{0}/B_{Q}=0.1$, and vary only $E_{0}$ such that
$\ell_{a}/R_{0}\simeq\\{0.2,0.35,0.5,0.6\\}$. These values were chosen for
simulations to be computationally feasible yet comparable to real pulsars, in
which $\ell_{a}/R_{0}$ varies from $\sim 10^{-4}$ to $\sim 1$ with increasing
rotation period. The mean free path of photons emitted at the critical energy
$\varepsilon_{c}=(3/2)(\lambdabar_{\text{C}}/\rho)\hat{\gamma}_{\pm}^{3}m_{e}c^{2}\zeta_{\pm}$,
defined as (e.g., Timokhin & Harding, 2015)

$\ell_{\gamma}\simeq\frac{8}{3}\frac{B_{Q}}{B}\frac{m_{e}c^{2}}{\varepsilon_{c}}\rho\zeta_{\gamma}^{-1}\
,$ (3)

is always $\sim 10^{-2}~{}R_{0}$ for our choice of parameters. Thus, the ratio
$\ell_{\gamma}/\ell_{a}\sim 0.1$ is also within the range $\sim 10^{-2}-1$
expected in real systems. The electron skin depth associated with a density
$n_{\text{GJ}}=|\rho_{\text{GJ}}|/e\equiv\Omega B_{0}/2\pi ce$ is $\sim
5\times 10^{-3}-10^{-2}~{}R_{0}$. The simulation domain has a size
$L_{R}\times L_{z}=(1.5R_{0})\times(2.5R_{0})$ discretized in $N_{R}\times
N_{z}=6000\times 10000$ cells, and the time step is $\Delta t=10^{-4}R_{0}/c$.
The grid size was chosen to resolve the electron skin depth associated with
the maximum density generated during pair cascades.

The rotating conducting disk is gradually spun up to its maximum angular
velocity in the first 200 time steps of the simulations and kept constant
thereafter. There is initially no plasma in the simulation domain, so this
spin induces a vacuum corotation electric field in the open field line region:
$E_{R}$ increases and $E_{z}$ decreases with $R$, whereas both components
decay with $z$ within a distance $\sim R_{0}$. We let the corotation field
develop for a light-crossing time, $L_{z}/c=2.5~{}R_{0}/c$, and then start
injecting plasma in cells just above the $z=0$ boundary for all $R<R_{0}$. At
every time step, we inject one electron-positron pairs per cell carrying a
density $n_{\text{inj}}=\kappa(E_{\parallel}/eR_{0})$, where $\kappa=0.1$ and
$E_{\parallel}=|\mathbf{E}_{\parallel}|$. Pairs are injected at rest. We
choose a value of $\kappa\ll 1$ to ensure this injection provides only a seed
plasma for the cascades and does not dominate the plasma outflow.

##  3 Results

As pairs are injected and experience the vacuum electric field, only electrons
are able to escape the surface of the conducting disk, whereas positrons are
immediately reabsorbed at the boundary. In their acceleration along the
magnetic field lines, electrons emit curvature photons, which then decay into
pairs, triggering the pair cascade. The accelerating electric field
$E_{\parallel}$ is then locally screened and the cascade stops. When the
plasma flows away from the conducting disk, the vacuum gap reopens and the
process restarts.

In Figure 1, we show the electron, positron and photon densities and
$E_{\parallel}$ for a simulation with $\ell_{a}/R_{0}\simeq 0.2$ when two pair
production bursts are visible in the simulation domain. Due to horizontal
gradients of the magnetic field curvature, photon emission and pair production
do not occur uniformly across the open field lines. Instead, the cascade is
triggered at $R\simeq R_{0}/2$ and $z\simeq\ell_{a}$. For $R\gtrsim R_{0}/2$,
$E_{\parallel}$ decreases, limiting acceleration and consequently photon
emission, whereas for $R/R_{0}\ll 1$ the radius of curvature of magnetic field
lines rapidly diverges and pair production is suppressed (Arons & Scharlemann,
1979). For $R\lesssim R_{0}/2$, the pair cascade develops at $z<\ell_{a}$,
i.e., the pair production front is inclined towards the magnetic axis.

As pair production bursts develop, inductive plasma waves are self-
consistently excited (Levinson et al., 2005; Cruz et al., 2021). Due to the
two-dimensional structure of the pair production front, these waves are not
only emitted in the positive $z$ direction, but also in all other directions —
see Figure 1d. The inclination of the wavevectors relative to $\mathbf{B}$ is
not the result of cross-field particle motion, but rather of field-aligned
currents coherently developed on adjacent field lines.

This inclination is a key ingredient in the self-consistent excitation of
electromagnetic plasma modes. The wave generation process can be understood as
follows (PTS20, Melrose et al., 2020): first, the non-uniformity of
$\mathbf{E}_{\parallel}$ across dipolar field lines gives rise to an
oscillating azimuthal component of the magnetic field $\tilde{B}_{\phi}$ via
$(\partial\mathbf{\tilde{B}}/\partial
t)\sim-c\nabla\times\mathbf{E}_{\parallel}$, where
$\mathbf{\tilde{B}}=\tilde{B}_{\phi}\mathbf{e}_{\boldsymbol{\phi}}$; then, an
oscillatory electric field component $\mathbf{\tilde{E}}$ is excited via
$(\partial\mathbf{\tilde{E}}/\partial t)\sim c\nabla\times\mathbf{\tilde{B}}$.
Gradients in $\mathbf{\tilde{B}}$ occur predominantly along the normal to the
pair production front, which is also the direction of the wavevector
$\mathbf{k}$ of these waves. In general, both $\mathbf{k}$ and
$\mathbf{\tilde{E}}$ have components parallel and perpendicular to
$\mathbf{B}_{\textrm{d}}$. The requirement that the angle between the normal
to the pair production front and $\mathbf{B}_{\textrm{d}}$ is finite is
essential for this process to operate. Hence, modes with
$\mathbf{k}\parallel\mathbf{B}_{\textrm{d}}$ are never produced without
accompanying modes with $\mathbf{k}\perp\mathbf{B}_{\textrm{d}}$. All our
simulations show that these modes are excited at the local plasma frequency
$\omega_{0}$ (i.e., at the frequency of inductive $\mathbf{E}_{\parallel}$
oscillations) and are linearly polarized. As pairs are produced, the local
density increases and the wave frequency extends to $\sim 10-100\omega_{0}$.
The properties of these modes are consistent with the superluminal
O-modes111This mode should not be confused with the O-mode presented in plasma
physics textbooks (e.g., Nicholson, 1983; Stix, 1992). The mode presented here
propagates in a range of directions from purely perpendicular to purely
parallel to $\mathbf{B}_{\textrm{d}}$, while textbook O-modes exist only for
$\mathbf{k}\perp\mathbf{B}_{\textrm{d}}$. identified in previous works (Arons
& Barnard, 1986, PTS20).

Figure 2:  Poynting vector components $S_{z}$ and $S_{R}$ associated with the
fluctuating electromagnetic field components (calculated by subtracting a
local average to the total fields). This snapshot was taken at the same time
as Figure 1. Red (blue) tones identify regions where there is a positive
(negative) electromagnetic energy density flux in the $z$ and $R$ directions
in panels a and b, respectively. The magnetic field lines are shown in solid
black lines.

A snapshot of the components of the Poynting vector
$\mathbf{S}=(c/4\pi)\mathbf{\tilde{E}}\times\mathbf{\tilde{B}}$ is shown in
Figure 2. The wave components $\mathbf{\tilde{E}}$ and $\mathbf{\tilde{B}}$
were calculated by subtracting local time-averaged $\mathbf{E}$ and
$\mathbf{B}$, respectively. Figure 2a shows rings of electromagnetic flux
emitted around the most recently created plasma burst, centered at
$z/R_{0}\simeq 0.4$ and $R/R_{0}\simeq 0.5$. Waves emitted upward carry a
positive $S_{z}$, whereas waves emitted downward have initially a negative
$S_{z}$ that is then reversed after the waves are reflected at the $z=0$
boundary. Electromagnetic waves are emitted in all directions. However, due to
the inclination of the bursts relative to the magnetic field lines, a larger
Poynting flux is generated on the flanks of the bursts. After some altitude,
$S_{z}$ thus exhibits a double-peaked structure, being larger close to the
magnetic axis and the conducting field line boundary and smaller in the center
of the open field line bundle.

Figure 3:  Comparison between $E_{\parallel}$ oscillations induced by cascades
with different $\ell_{a}/R_{0}$, where $\ell_{a}$ is the characteristic
distance required for electrons to be accelerated to energies capable of
emitting pair producing photons. Both the wavelength of the oscillations and
the aspect ratio of the pair production burst (given by the ratio between its
height and width) increase with $\ell_{a}/R_{0}$. The snapshots were taken at
times $t=5.4$, $5.25$, $5.75$ and $6.25R_{0}/c$ for $\ell_{a}/R_{0}=0.2$,
$0.35$, $0.5$ and $0.6$, respectively. A movie showing the temporal evolution
of these quantities is available online.

Figure 2b shows that the part of the electromagnetic flux rings convert into
bands extended in the $z$ direction (see e.g. $R/R_{0}\simeq 0.5$ and
$z/R_{0}\simeq 1$) over time. These bands result from the continuous
reflection of the waves between the magnetic axis and the conducting region of
closed field lines, and they drift across the open field line region over a
time scale $\sim R_{0}/c$.

The electrodynamics of the two-dimensional pair cascades identified above for
$\ell_{a}/R_{0}\simeq 0.2$ holds also for different values of
$\ell_{a}/R_{0}$. However, there are important effects to note: first, the
wavelength of the $E_{\parallel}$ oscillations (and consequently of
$\mathbf{\tilde{E}}$ and $\mathbf{\tilde{B}}$) decreases with decreasing
$\ell_{a}$, due to the larger multiplicity generated in the cascades in this
regime; second, the pair production burst extends to smaller $\theta$ for
smaller $\ell_{a}/R_{0}$; third, the shape of the pair production burst is
flatter in the radial direction for smaller $\ell_{a}/R_{0}$, and more round
for larger $\ell_{a}/R_{0}$ — see Figure 3. The latter effect is responsible
for a more efficient generation of Poynting flux, since it allows for a larger
portion of the pair production bursts to excite waves with $\mathbf{k}$ almost
purely aligned with $z$. This is visible in Figure 4, where we show $\langle
S_{r}\rangle\equiv\langle S_{z}\cos\theta+S_{R}\sin\theta\rangle$, averaged
for all times after the initial vacuum transient at $z/R_{0}=2$ for
simulations performed with different values of $\ell_{a}/R_{0}$. The angle
$\theta$ is normalized to $\theta_{0}\equiv\arcsin(R_{0}/R_{*})$, i.e., the
analogue of the polar cap angle in pulsar magnetospheres. For
$\ell_{a}/R_{0}\simeq 0.2$, we observe the double-peaked structure in this
profile described before (see peaks at $\theta\simeq 0$ and
$\theta\simeq\theta_{0}$). However, for large $\ell_{a}/R_{0}$, the edge
component is absent. We also observe that the peak value of $\langle
S_{r}\rangle$ on the magnetic axis decreases with $\ell_{a}/R_{0}$, but is
always $\sim 10^{-6}-10^{-4}~{}S_{0}$, where
$S_{0}=ce^{2}n_{\text{GJ}}^{2}R_{0}^{2}/16\pi$ is the total Poynting flux
launched by the rotating conductor, a result consistent with the fraction of
pulsar spin down power observed in the radio (e.g, Lorimer & Kramer, 2004).

Figure 4:  Average Poynting flux $\langle S_{r}\rangle$ at $z/R_{0}=2$ for
simulations with $\ell_{a}/R_{0}=0.2$ (blue) and $\ell_{a}/R_{0}=0.6$
(orange). The shaded regions above each curve extend for a standard deviation
above the time average.

##  4 Discussion

In this Letter, we have presented the first 2D simulations of pulsar polar cap
pair cascades including the QED effects from first principles. Our results
show that the gap dynamics has two significant two-dimensional features that
could not be captured in one-dimensional simulations: a) pair production is
inhibited close to the magnetic axis, due to the null curvature of the
magnetic field in this region, and b) gradients of the magnetic field
curvature across the polar cap induce an inclination between the normal to the
pair production front and the background magnetic field.

We have observed that $E_{\parallel}$ oscillations inductively driven in the
inclined pair production bursts can act as a source for coherent
electromagnetic waves. Although the initial oscillations in $E_{\parallel}$ is
longitudinal, the induced waves are oblique and electromagnetic in nature,
which makes them distinct from the L mode waves discussed in previous works
(Rafat et al., 2019; Melrose et al., 2020). Moreover, the electromagnetic
modes observed in our simulations are naturally produced in the oblique pair
production fronts, and do not require the production of intermediary modes
such as the L modes suggested by Melrose et al. (2020). The induced
electromagnetic waves propagate in all directions, but the Poynting flux flows
predominantly outwards, away from the star. We have identified the coherent
electromagnetic modes to be linearly polarized, superluminal O-modes, and
verified that the resulting Poynting flux from this emission has two peaks:
one on the magnetic axis and another on the edge of the open field line
bundle. We interpret this as a consequence of the relative orientation between
the normal to the pair production bursts and the background magnetic field.
Thus, we expect this mechanism to operate also in magnetic field topologies
more complex than a pure dipole (e.g. in multipolar magnetic fields), provided
that the field curvature changes over a scale larger than the gap height
$\sim\ell_{a}$. In such topologies, the configuration of the open field lines,
where cascades develop, may lead to complex-shaped pair production bursts and
thus different emission power profiles.

We have also shown that a fraction of the emitted waves is continuously
reflected at last open field line boundary, giving rise to a drifting
component in the Poynting flux on a time scale $\sim 1$ $\mu$s. We note that
this should not be confused with the observed radio drifting sub-pulses, that
occur on time scales $\sim 0.1-1$ s (e.g., Weltevrede et al., 2007). The
drifts identified in our simulations should be detectable at high time
resolution and could be used to diagnose the height of the emission.

The simulations presented in this work adopt a very dilute, space charge-
limited flow from a rotating conductor, and focus on the role of pair cascades
in providing the current to screen the vacuum gap. In reality, the neutron
star is expected to provide a larger current density $\sim
j_{\text{GJ}}=\rho_{\text{GJ}}c$. In that case, the vacuum gaps are induced by
the inability of the star to match the current density $j>j_{\text{GJ}}$,
required by the global magnetosphere when general relativistic frame-dragging
effects are taken into account (Gralla et al., 2016). It is not expected that
pair cascades operate differently when these effects are included; however,
given the same $B_{0}$ and $\Omega$, the ratio $\ell_{a}/R_{0}$ may differ
slightly from what is presented in this work.

FC, TG and LOS are supported by the European Research Council (ERC-2015-AdG
Grant 695088). FC is also supported by FCT (Portugal) (grant
PD/BD/114307/2016) in the framework of the Advanced Program in Plasma Science
and Engineering (APPLAuSE, FCT grant PD/00505/2012). AC is supported by NSF
grants AST-1806084 and AST-1903335. AS is supported by NASA grant
80NSSC18K1099 and NSF grant PHY-1804048. We acknowledge PRACE for granting
access to MareNostrum (Barcelona Supercomputing Center, Spain) and to TGCC
Joliot Curie (CEA, France), where the simulations presented in this work were
performed.

## References

  * Arons & Barnard (1986) Arons, J., & Barnard, J. J. 1986, Astrophys. J., 302, 120, doi: 10.1086/163978
  * Arons & Scharlemann (1979) Arons, J., & Scharlemann, E. T. 1979, Astrophys. J., 231, 854, doi: 10.1086/157250
  * Bai & Spitkovsky (2010) Bai, X.-N., & Spitkovsky, A. 2010, Astrophys. J., 715, 1282, doi: 10.1088/0004-637X/715/2/1282
  * Beloborodov (2008) Beloborodov, A. M. 2008, Astrophys. J. Lett., 683, L41, doi: 10.1086/590079
  * Cruz et al. (2021) Cruz, F., Grismayer, T., & Silva, L. O. 2021, Astrophys. J., 908, 149, doi: 10.3847/1538-4357/abd2c0
  * Daugherty & Harding (1982) Daugherty, J. K., & Harding, A. K. 1982, Astrophys. J., 252, 337, doi: 10.1086/159561
  * Del Gaudio (2020) Del Gaudio, F. 2020, PhD thesis, Instituto Superior Técnico
  * Erber (1966) Erber, T. 1966, Rev. Mod. Phys., 38, 626, doi: 10.1103/RevModPhys.38.626
  * Fonseca et al. (2008) Fonseca, R. A., Martins, S. F., Silva, L. O., et al. 2008, Plasma Phys. Control. Fusion, 50, 124034, doi: 10.1088/0741-3335/50/12/124034
  * Fonseca et al. (2002) Fonseca, R. A., et al. 2002, in Computational Science — ICCS 2002, ed. P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, & J. J. Dongarra (Berlin, Heidelberg: Springer Berlin Heidelberg), 342–351
  * Goldreich & Julian (1969) Goldreich, P., & Julian, W. H. 1969, Astrophys. J., 157, 869, doi: 10.1086/150119
  * Gralla et al. (2016) Gralla, S. E., Lupsasca, A., & Philippov, A. 2016, ApJ, 833, 258, doi: 10.3847/1538-4357/833/2/258
  * Grismayer et al. (2016) Grismayer, T., Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Phys. Plasmas, 23, 056706, doi: 10.1063/1.4950841
  * Grismayer et al. (2017) —. 2017, Phys. Rev. E, 95, 023210, doi: 10.1103/PhysRevE.95.023210
  * Kelner et al. (2015) Kelner, S. R., Prosekin, A. Y., & Aharonian, F. A. 2015, ApJ, 149, 33, doi: 10.1088/0004-6256/149/1/33
  * Levinson et al. (2005) Levinson, A., Melrose, D., Judge, A., & Luo, Q. 2005, ApJ, 631, 456, doi: 10.1086/432498
  * Lorimer & Kramer (2004) Lorimer, D. R., & Kramer, M. 2004, Handbook of Pulsar Astronomy (Cambridge University Press)
  * Melrose et al. (2020) Melrose, D. B., Rafat, M. Z., & Mastrano, A. 2020, Mon. Not. R. Astron. Soc., 500, 4549, doi: 10.1093/mnras/staa3529
  * Nicholson (1983) Nicholson, D. R. 1983, Introduction to Plasma Theory (New York: John Wiley & Sons)
  * Philippov et al. (2020) Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, Phys. Rev. Lett., 124, 245101, doi: 10.1103/PhysRevLett.124.245101
  * Rafat et al. (2019) Rafat, M. Z., Melrose, D. B., & Mastrano, A. 2019, J. Plasma Phys., 85, 905850305, doi: 10.1017/S0022377819000394
  * Ritus (1985) Ritus, V. I. 1985, J. Sov. Laser Res., 6, 497, doi: 10.1007/BF01120220
  * Ruderman & Sutherland (1975) Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51, doi: 10.1086/153393
  * Stix (1992) Stix, T. H. 1992, Waves in Plasmas (New York: AIP Press)
  * Sturrock (1971) Sturrock, P. A. 1971, ApJ, 164, 529, doi: 10.1086/150865
  * Timokhin (2010) Timokhin, A. N. 2010, Mon. Not. R. Astron. Soc., 408, 2092, doi: 10.1111/j.1365-2966.2010.17286.x
  * Timokhin & Arons (2012) Timokhin, A. N., & Arons, J. 2012, Mon. Not. R. Astron. Soc., 429, 20, doi: 10.1093/mnras/sts298
  * Timokhin & Harding (2015) Timokhin, A. N., & Harding, A. K. 2015, Astrophys. J., 810, 144, doi: 10.1088/0004-637x/810/2/144
  * Weltevrede et al. (2007) Weltevrede, P., Stappers, B. W., & Edwards, R. T. 2007, Astron. Astrophys., 469, 607, doi: 10.1051/0004-6361:20066855