# Explainable Artificial Intelligence driven mask design for self-supervised
seismic denoising

Claire Birnie
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>&Matteo Ravasi
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
###### Abstract

The presence of coherent noise in seismic data leads to errors and
uncertainties, and as such it is paramount to suppress noise as early and
efficiently as possible. Self-supervised denoising circumvents the common
requirement of deep learning procedures of having noisy-clean training pairs.
However, self-supervised coherent noise suppression methods require extensive
knowledge of the noise statistics. We propose the use of explainable
artificial intelligence approaches to “see inside the black box” that is the
denoising network and use the gained knowledge to replace the need for any
prior knowledge of the noise itself. This is achieved in practice by
leveraging bias-free networks and the direct linear link between input and
output provided by the associated Jacobian matrix; we show that a simple
averaging of the Jacobian contributions over a number of randomly selected
input pixels, provides an indication of the most effective mask to suppress
noise present in the data. The proposed method therefore becomes a fully
automated denoising procedure requiring no clean training labels or prior
knowledge. Realistic synthetic examples with noise signals of varying
complexities, ranging from simple time-correlated noise to complex pseudo rig
noise propagating at the velocity of the ocean, are used to validate the
proposed approach. Its automated nature is highlighted further by an
application to two field datasets. Without any substantial pre-processing or
any knowledge of the acquisition environment, the automatically identified
blind-masks are shown to perform well in suppressing both trace-wise noise in
common shot gathers from the Volve marine dataset and colored noise in post
stack seismic images from a land seismic survey.

##  1 Introduction

Machine learning, and in particular Deep Learning (DL), has recently become
part of everyone’s life; geophysics has not been immune to such a trend, and
as geophysicists, we have experienced the great benefits it can bring in terms
of quality and speed of execution. In the last decade, a swath of new
algorithms have been proposed across the entire seismic processing pipeline
that leverage DL-based methodologies. As deep learning is particularly
powerful in finding patterns in complex, high-dimensional data, it has been
successfully adopted for noise suppression/signal enhancement tasks across
many domains, ranging from natural images [e.g., 1] to medical images [e.g.,
2] to synthetic-aperture radar images [e.g., 3]. Within geophysics, different
neural network architectures have been proposed for the suppression of a
variety of noise signals. For example, the Denoising CNN (DnCNN) architecture
has been used for random noise suppression [4], a convolutional auto-encoder
for erratic trace-wise noise suppression [5], and generative adversarial
networks for ground-roll suppression [6], just to name a few. In addition to
this, a number of studies have been focused on developing networks that can
tackle the complex noise field as a whole [e.g., 7, 8, 9].

The majority of the denoising networks proposed so far are trained in a
supervised manner, requiring pairs of noisy-clean training data, often
obtained by either pre-processing of the noisy dataset by means of
conventional (non-DL) denoising procedures or through the creation of
realistic synthetic training data that is highly-similar to the field
application data. Given that creating training data in such a manner is far
from trivial, several researchers have turned their interest towards the
application of self-supervised procedures for denoising tasks. Here, the noisy
data represents both the input to the network and the network’s training
target, often with the input undergoing some minor modifications such that the
input and target are no longer identical. Arguably, AutoEncoders (AEs) were
the first of such approaches, where the input is reduced (encoded) to a latent
space representation of a significantly smaller dimensionality prior to being
up-scaled (decoded) back to the original space. Assuming the seismic signal of
interest is the dominant component in the data, and a suitably small latent
space is used, AEs have been shown to successfully suppress random noise in
seismic data [10]. Similarly, AEs can be used to attenuate coherent noise,
such as high-frequency noise components in marine seismic data [11] or
diffractions from stacked data [12]. Blind-spot networks represent another
family of self-supervised methods recently introduced in computer vision: to
counter the challenge of sourcing noisy-clean training pairs, [13] proposed a
pre-processing methodology to ‘blind’ a central pixel from the network’s
input, which they termed Noise2Void (N2V). Under the assumption that noise is
random and independent among neighbouring pixels, and due to the removal of
information from the central pixel, the network can only learn to recreate the
coherent signal - hence, its ability to act as a natural denoiser. [14]
illustrated the potential of blind-spot networks for suppression of random and
pseudo-random noise from seismic data. Acknowledging that noise is rarely
random, [15] extended the methodology of blind-spot networks to create a
blind-mask network where all pixels within a masked area have their values
replaced; thereby, preventing them from contributing to a central pixel’s
predicted value. [15] showed that when the mask covers noise with a high
correlation to the central pixel, their blind-mask methodology, termed
StructN2V, could efficiently suppress spatially-structured noise in
fluorescence microscopy images.

One of the biggest challenges with StructN2V is selecting a mask that blinds
all pixels that can ‘leak’ noise information into the network’s input, whilst
retaining enough pixels containing signal that the network can use them to
reproduce the signal in the central pixel. In seismic applications, [16] used
prior knowledge of the noise’s spatio-temporal characteristics to design a
trace-wise noise mask. Similarly, [17] and [18] adapted the trace-blinding
procedure for seismic deblending purposes. Whilst, [19] proposed an extended
mask for the StructN2V workflow for the suppression of dipping and hyperbolic
noise instances. In all applications, the purpose of the mask is to ensure
that neighbouring pixels contributing to a central pixel’s predicted value do
not contain coherent noise, as this could be detrimental to the network’s
training process. In order to analyse the influence of neighbouring pixels, we
propose to adapt a methodology from the field of eXplainable Artificial
Intelligence (XAI). More precisely, XAI is used here to remove the black-box
element of machine learning and provide meaning to a prediction by
highlighting how the input has been utilised to compute such an output. To
date, there has been relatively little use of XAI within geophysics. [20] and
[21] both utilise XAI methods for analysing NN predictions in facies
classification, whilst [22] investigates the contribution of different areas
in a seismic shot gather with respect to the prediction of the associated
surface wave dispersion curves.

In this study, we propose a workflow that leverages self-supervised, blind-
spot networks in combination with XAI for the identification and subsequent
application of an optimal noise mask for the training of blind-mask networks;
this ultimately results in an automated denoising scheme for coherent noise
suppression, which requires neither a labelled training dataset nor prior
knowledge of the statistical properties of the noise. The proposed workflow is
initially validated on a synthetic data for the scenarios of colored noise,
trace-wise noise and pseudo-rig noise. Subsequently, the workflow is further
tested on the Volve field dataset to denoise trace-wise noise present in the
common shot gathers, before being utilised for psuedo-random noise suppression
of a post-stack image from a land seismic survey. The successful suppression
of a wide range of noise types, as well as the application to data at varying
stages of processing, highlights the robustness of the proposed procedure.

##  2 Methodology

This study combines self-supervised blind-spot networks as presented in [13]
with [3]’s proposal of using the Jacobian matrix of the network to study the
behaviour of bias-free neural networks. In this section, we detail these two
core concepts before outlining how we combine them to determine an optimal
denoising mask for the suppression of coherent noise in seismic data.

###  2.1 Self-supervised, blind-mask denoisers

Similar to statistical methods like local or non-local means [23], blind-spot
denoising aims to reproduce a central pixel’s value using only the values of
neighbouring pixels. In this paper, we use a variation of the original (N2V)
blind-spot methodology of [13]. First, a number of pixels that will become the
focus of the training procedure are selected; throughout this work we will
refer to these pixels as active pixels. To create the so-called blind-spots in
the input image, the original value of each active pixel is replaced by the
value of a randomly selected pixel within a radius, $r$. This selection and
corruption process is repeated at every epoch, prior to passing the corrupted
image to the network; thereby, the location and values of the active pixels
change on every exposure to the network. As a self-supervised procedure, the
network’s target is the raw noisy data. Due to the fact that the network
focuses on learning only from neighbouring pixels, its performance is only
evaluated at the active pixel locations, i.e., at the blind-spots. This is
implemented by computing the training loss (for example, the $L_{1}$ norm)
across the full image before applying a binary mask such that only the values
at the active pixels are effectively evaluated. As with most deep learning
procedures, the prediction errors are then back-propagated through the network
and used to update its parameters to be used when evaluating the network in
the next training epoch.

As previously mentioned, the blind-spot approach is only applicable for the
scenario of random noise – i.e., noise that is independent between pixels. In
order to tackle noise that presents any degree of coherency, a blind-mask is
required instead of a single blind-spot \- where the mask shape is designed to
remove the influence of any pixels that contain noise correlated to that
within the central (active) pixel. Similar to the original procedure, a number
of active pixels are initially selected. In this instance, not only is the
active pixel’s value replaced, as are the values of all pixels within the
masked area, where the mask is centered on the active pixel. Due to the larger
number of pixels to be replaced, the replaced values are now extracted from a
pre-selected noise distribution, as opposed to using a value from a
neighbouring pixel. The corrupted image is then passed through the network as
the training input, the raw image is used as the training target, and the loss
is calculated only at the active pixel locations.

It is important to note that the implementation of the blind-spots (or masks)
is only performed at the training stage. During inference, the raw, noisy
images are fed directly to the trained network, without any need for pre-
processing or masking of the predictions. Finally, it is worth remembering
that since blind-spot networks are convolutional neural networks, they benefit
from the ability convolutions have that allow them to be applied to inputs of
varying shapes/sizes. As such, whilst training may be performed using patches
of the training data, the inference process is usually carried out on the
entire data to be denoised.

###  2.2 Input Contribution Analysis via the Jacobian Matrix

Feed-forward neural networks with REctified Linear Units (ReLUs) activation
functions, such as the commonly-used UNet architecture, are piecewise affine
transformations: as a consequence, one can express the nonlinear
transformation $y=f(x)$ as the sum of two terms [3]

$y=f(x)=Jx+b$ (1)

where $x$ and $y$ are the input and output of size $N_{h}\times N_{w}$, $J$ is
the Jacobian of size $N_{h}N_{w}\times N_{h}N_{w}$ of $f(\cdot)$, and $b$ is
the net bias of the same size of the input. In the special case of bias-free
networks (i.e., networks whose additive constants are all fixed to zero), the
net bias disappears, and the complex, non-linear relationship between input
and output resorts to the dot product between the Jacobian and the input. Note
that the relationship between the input and output is only locally-linear; in
other words, it is only linear between one given input and its corresponding
output. As such, the benefit of non-linearity introduced by the use of neural
networks remains. [3] used this principle to study the behaviour of a bias-
free, convolutional supervised denoiser.

To be able to apply the same concept to self-supervised denoisers, we utilise
here a bias-free UNet architecture; one can therefore consider the denoised
product as the application of an adaptive linear transform where each pixel in
the denoised image is a weighted sum of the input noisy pixels. Note that this
is nothing more than the non-local means algorithm, with the main difference
that the weights are learned via a neural network instead of being computed
using simple pre-defined functions.

Evaluated during any stage of training or at the inference stage, the Jacobian
matrix identifies the weighting of any given pixel from the input image on the
prediction of a selected pixel in the denoised image. In this sense, the
Jacobian matrix can be utilised in two manners for the analysis of blind-spot
networks:

  * •

firstly, during the training stage the computation of the Jacobian matrix can
highlight whether the network is learning to avoid utilizing the blind-spots;
in other words, utilising the Jacobian matrix we can compute the influence
that a blind-spot pixel’s value has on its own prediction. In theory, this
influence should tend to zero during the training process; and,

  * •

secondly, the Jacobian matrix can also indicate the neighbouring pixels which
heavily influence the the blind-spot’s predicted value.

In this study, we leverage the second opportunity and utilise the Jacobian
matrix to gain knowledge of the neighbouring pixels’ contributions to a
central pixel’s predicted value. In the instance of the application of a
blind-spot network on data contaminated with coherent noise, the neighbouring
pixels’ influence will highlight not only the pixels that contribute to the
recreation of the central pixel’s signal value, but also the pixels that
contribute to the recreation of the central pixel’s noise value - where we
consider the central pixel’s raw value as a summation of a signal component
and noise component. Therefore, by computing the Jacobian matrix over a large
number of pixels at the inference stage of a blind-spot network, a noise mask
can be designed to evade the contribution of pixels that are identified as
high influencers for later use with a self-supervised, blind-mask denoiser.

Figure 1: Proposed denoising workflow that involves: 1) performing a standard,
self-supervised, blind-spot denoising approach; 2) identifying the
contributing pixels through a Jacobian matrix analysis; 3) designing a noise
mask that removes contributing pixels from the training input; and, 4)
training a blind-mask network in a self-supervised fashion. High amplitude
pixels in the training datasets indicate corrupted (i.e., masked) data points.

###  2.3 XAI for the mask creation

The proposed workflow is provided in Figure 1. In order to determine the
contributions from neighbouring pixels, initially a standard blind-spot
network is trained, where only the active pixels have their values replaced
prior to feeding into the network. After training to convergence, the Jacobian
analysis is performed. This analysis is schematically illustrated in Figure 2,
which highlights the linking between the full input image and individual
output pixels via rows of the network’s Jacobian matrix. The analysis is
performed as follows:

  1. 1.

1000 random points (pixels) are selected from the network’s prediction, herein
referred to as probe points;

  2. 2.

the network’s Jacobian matrix is evaluated with respect to the probe points,
in other words, 1000 rows of the Jacobian matrix is computed;

  3. 3.

the computed Jacobian rows are each reshaped to the dimensions of the input
image and a window of $31\times 31$ is extracted around the probe points
(where this window covers all neighbouring pixels that provide a significant
contribution towards the predicted central value);

  4. 4.

finally, the absolute values of the 1000 extracted windows are averaged to
give a single matrix that represents the average contribution of the
neighbouring pixels to any predicted pixels.

Figure 2: Probing of the neural network’s Jacobian Matrix for the creation of
the Jacobian Map that illustrates the contributions from input pixels onto a
predicted pixel in the output.

The final averaged matrix is herein referred to as the Jacobian Map. For
comparative purposes, we normalise the map such that the contribution from all
pixels sums to one. It is important to note that the Jacobian Map represents
the total contribution of the pixels without any separation between the pixels
that aid the noise recreation and those that are vital for the signal
recreation. However, we hypothesise that the signal contributions that the
blind-spot network will learn will average out across the probe points, while
the contributions from noise oriented in a particular manner will amplify, as
schematically illustrated in Figure 3.

Figure 3: Separate contributions from signal and noise towards the full
averaged Jacobian Map. The three annotated points on the single and noise
panels correspond to the succeeding Jacobian Maps, paired via their colouring
scheme. Red dashed lines on the signal Jacobian Maps indicate the orientation
of pixels with a high contribution.

After the computation of the Jacobian Map, two approaches can be taken to
determine the shape of the mask for the subsequent blind-mask network - both
of which are investigated in this study. First, the mask can be manually
created using the results of the contribution analysis alongside prior
knowledge of the pixels that strongly contribute to the signal recreation and
the expected noise. This prior knowledge can easily be obtained by training a
blind-spot network and computing the Jacobian Map for a clean synthetic
dataset with the same frequency content and survey design of the field data of
interest - similar to the top row of Figure 3. We term this approach “XAI-
guided mask design” to emphasise the additional information considered in the
mask design process. The alternative approach is to set a contribution limit
and all pixels that contribute above such limit in the Jacobian Map become
part of the mask. As this incorporates no additional information, we term this
“XAI-driven mask design”. The former approach has the benefit that it reduces
the removal of pixels that are heavy contributors to the signal prediction,
however, its drawback is that it is a manual approach. The latter approach may
result in a mask that is larger than necessary and hence removes pixels that
provide useful signal information, however, it is an automated approach that
only requires the setting of the contribution limit.

After the mask is selected, the blind-mask training procedure is performed, to
teach the network to utilise only unmasked pixels in the prediction of the
central pixel. Once the network is trained, the entire raw, noisy data can be
passed through the network, without any requirement of pre-processing or
masking. The final result is the denoised product arising from the self-
supervised, blind-mask scheme.

##  3 Numerical Analysis

To validate our proposed procedure, the approach is first tested on a
realistic synthetic dataset contaminated by different types of noise. The
dataset has been modelled using a realistic synthetic velocity model that
mimics the subsurface structure of the Volve field and a ocean-bottom
configuration (see [24] for details). The noise types under investigation
include White Gaussian Noise (WGN), time-correlated noise, both isotropic and
anisotropic coloured Gaussian noise (CGN), and pseudo-rig noise. The isotropic
CGN has been obtained by filtering a WGN realization with a filter length of 3
in both time and space directions, whilst the anisotropic CGN used a filter
length of 4 and 2 in time and space, respectively. The time-correlated noise
is generated by applying a trace-wise, band-pass filter to a WGN realisation,
and the pseudo-rig noise is generated by creating linear events mimicking
waves travelling at a constant speed (i.e., through the sea) convolved with a
non-compact wavelet extending for 10 seconds.

In all the following examples, the data undergoes log1p amplitude scaling,
similar to that of [17]. The benefit of such a scaling method, as opposed to
the conventional amplitude gain correction (AGC), is that it can be easily
reversed, allowing the amplitudes to return to their original values after the
denoising task is complete. In order to obtain an adequate number of training
samples, the shots are sliced into overlapping patches of size $64\times 64$.
To ensure patches contain an adequate percentage of signal (as opposed to pre-
direct arrival empty contributions), patches whose mean energy is in the
$25$th percentile or below are excluded from the training dataset. A 4-layer
UNet is trained, with an initial filter size of 32 that doubles per layer.
From all of the selected patches, $4096$ are used for training and $1024$ are
retained for validation purposes, utilising a batch size of $128$. For both
the blind-spot and blind-mask procedures, $2\%$ of the pixels within each
input patch are identified as active pixels.

The application of blind-spot denoising is provided in Figure 4, where the
first column shows the raw, noisy data and the second column shows the data
after blind-spot denoising. The other two columns show the difference between
the denoised data with the noisy and clean data, respectively. Whilst the
final column illustrates the averaged Jacobian matrix for 1000 randomly
selected pixels in the denoised product. As expected, the network successfully
suppresses the WGN in the first example, due to the noise’s incoherency
between pixels. The Jacobian associated with this result therefore highlights
which neighbouring pixels provide a strong contribution towards replicating
the signal in the central pixel. Note that due to the powerful denoising
capabilities of the blind-spot approach for data contaminated with WGN, this
scenario will not be considered any further in this study.

The following four scenarios all contain varying levels of correlated noise
and their associated blind-spot products are almost identical to the raw,
noisy data. The near-perfect recreation of the noisy data is unsurprising
given the coherent nature of the noise components. For the time-correlated
noise (second row), the Jacobian analysis clearly indicates a higher
contribution from pixels along the same trace as the central pixel. For both
CGN examples, the pixels within the filter length of the noise are shown to be
high contributors, evident by the rectangular shape of the pixels with a high
value in the Jacobian Map. Finally, the Jacobian analysis of the pseudo-rig
noise highlights the linear, slightly rotated nature of the noise with
stronger contributions of pixels to the top left and bottom right of a central
pixel.

To conclude, minimal signal leakage is observed in every example of Figure 4,
as determined by the lack of coherency in the differences between the
predictions and their associated clean shots. This highlights the network’s
suitability for recreating the complex seismic wavefronts present in such a
realistic synthetic dataset.

###  3.1 XAI-guided mask design

The XAI-guided methodology allows us to incorporate prior knowledge of the
acquisition environment, alongside a modicum of trial-and-error, in order to
identify the optimal denoising mask. This allows us to keep highly
contributing pixels out of the masked area if we believe them to be likely to
provide only signal information. The time-correlated noise is a prime example
of this: with the knowledge that the noise is independent between traces and
therefore the highly contributing pixels in nearby traces likely provide
signal information, then a noise mask can be built only masking highly
contributing pixels along the trace direction, as highlighted in the top, left
plot of Figure 5. The noise masks for the remaining noise scenarios are shown
in the left column of Figure 5 with their associated denoised products
displayed in the middle column of the same figure. The two final columns
illustrate the difference between the original noisy data and the clean data
used for benchmarking. A conventional noise metric, Peak Signal-to-Noise Ratio
(PSNR), is provided in the titles for the original noisy data and the denoised
products.

Figure 4: Self-supervised, blind-spot application to synthetic data
contaminated with WGN, time-correlated noise, isotropic CGN, anisotropic CGN
and pseudo-rig noise, from top to bottom respectively. From left to right, the
noisy training data, the ‘denoised’ product, and its difference with the
original noisy and clean shots, and finally the averaged Jacobian matrix.
Figure 5: XAI-Guided blind-mask implementation on synthetic data contaminated
with time-correlated noise, isotropic CGN, anisotropic CGN and pseudo-rig
noise, from top to bottom respectively. The first column represents the XAI-
guided noise masks where the red cross marks the active pixel location. The
second column shows the noisy data and the third the denoised product. The
fourth and fifth columns show the difference between the denoised product and
the original noisy data and the clean shot, respectively.

For all of the considered scenarios, the manually selected masks have resulted
in suppression of almost all the noise present in the shot gathers,
significantly increasing the PSNR values. However, the noise masks cover a
significantly larger area than the previous single blind-spots and therefore
some useful signal information is also being masked. This is evident for the
bottom two noise masks which are substantial in size and result in a slightly
higher level of signal leakage in comparison to the blind-spot approach, when
considering the difference between the denoised product and the clean data.

###  3.2 XAI-driven mask design

Unlike the XAI-guided approach, the only variable to be optimised in the XAI-
driven approach is the contribution threshold (from the Jacobian analysis) at
which a pixel does, or does not, become part of the blind-mask. To ensure
consistency across examples, the total contribution of all pixels in the
Jacobian map is always provided as a normalised sum; in other words, the sum
of the contribution of all pixels within the input image sums to one. The
spatial extent of the noise coherency dictates the masks’ contribution
threshold (due to the number of significantly contributing pixels). In the
following analysis, the cut-off percentage has been chosen based on a sweeping
analysis.

The left column of Figure 6 illustrates the automatically selected noise masks
and their selected thresholds. In general, these masks cover a larger area
than their XAI-guided counterparts and their shapes are less obviously rooted
in the noise characteristics. For example, the time-correlated noise mask now
includes the masking of nearby traces, whilst the pseudo-rig noise mask
includes pixels outwith the general masked area. The remainder of Figure 6
portrays the results of the blind-mask denoising procedure. As in the XAI-
guided approach, the XAI-driven procedure has managed to suppress almost all
the noise present in each of the four coherent noise scenarios, resulting in
significant gains in the PSNR values. Interestingly, despite the masks
typically covering a larger area than the XAI-guided approach, a very similar
level of signal leakage is visually observed.

As before, the suppression of the pseudo-rig noise poses the biggest challenge
to the blind-mask network. As the interfering waves propagate through the same
medium as those from the actual survey of interest (i.e., our seismic signal),
directionality becomes the only discriminative factor. It is therefore
particularly challenging to find a threshold that adequately masks the noise
contribution whilst retaining enough information for the network to recreate
the desired seismic signal. Whilst the noise has been efficiently suppressed,
the high-coherency of the noise and the shared propagation properties between
the noise and the signal, has led to higher signal leakage than would be
desired. This feature is particularly noticeable in the early arrivals on the
right-side of the denoised product (and its corresponding difference plots),
where the noise and signal propagate at similar velocities. We conclude that
this scenario may be pushing the limits on what can be adequately tackled via
self-supervised, blind-mask networks.

Figure 6: XAI-Driven blind-mask implementation on synthetic data contaminated
with time-correlated noise, isotropic CGN, anisotropic CGN and pseudo-rig
noise, from top to bottom respectively. The first column represents the
automatically selected, XAI-driven noise masks, where the red cross marks the
active pixel location and the CCO value represents the contribution cut-off
threshold for inclusion in the mask. The second column shows the noisy data
and the third the denoised product. The fourth and fifth columns show the
difference between the denoised product and the original noisy data and the
clean shot, respectively.

##  4 Field Results

Following on from the previous examples, which utilised a synthetic dataset
inspired from the Volve subsurface, in this section we apply the same
methodology to shot gathers from a survey collected over the Volve oil field,
previously released by Equinor 111https://www.equinor.com/energy/volve-data-
sharing. Figure 7 displays five shot gathers (top row); we observe how such
gathers are contaminated with a small volume of trace-wise noise. As with the
synthetic examples, the data undergoes log1p amplitude scaling prior to being
sliced into overlapping patches of size $64\times 64$, with patches whose mean
energy is below the $15$th percentile being excluded from the training
dataset. As before, a 4-layer UNet is trained, with an initial filter size of
32 that doubles per layer. All remaining parameters (e.g., training sample
size, number of active pixels, batch size, etc.) are kept consistent with the
earlier synthetic studies. To avoid utilising synthetics to aid the mask
design, or performing extensive analysis of the field noise, only the
automatic, XAI-Driven approach is implemented for the field data, with any
pixel contributing more than 2% to the Jacobian Map being included into the
blind-mask.

Figure 7: Common shot gathers from the Volve field dataset (top row) and their
respective denoised products (middle row) following the XAI-driven masking
procedure. The final row illustrates the difference between the raw and
denoised shots.

Figure 8 details the steps of the mask selection and application with panel b
illustrating the denoising performance of the blind-spot network and its
associated pixel contributions as defined by the Jacobian Map (panel c). The
mask (panel f) is defined as anywhere where the pixel contribution is greater
than $2\%$ of the total pixel contributions. Finally, panels d and g
illustrate the blind-mask denoising performance and the associated residual.
Similarly, the middle row of Figure 7 illustrates the blind-mask denoising
performance on the five shot gathers selected earlier, with their associated
residuals (i.e., suppressed noise) displayed in the final row. Whilst the
original shots are relatively clean, the computed linear mask has successfully
suppressed the trace-wise noise that appeared at different traces across all
of the shot gathers. Considering the residual plots, very little coherent
energy is observed indicating that minimal signal damage has occurred.

Figure 8: A single shot example of the XAI-driven denoising procedure for the
Volve field dataset. (a) the raw data, (b) after blind-spot denoising, and (d)
after blind-mask denoising. (c) is the averaged Jacobian matrix and (f) the
respective XAI-Driven noise mask. (e) and (g) are the residuals for the blind-
spot and blind-mask procedures, respectively.

###  4.1 Application to post-stack data

Until now, our proposed methodology has been applied on seismic data prior to
migration (i.e., shot-gathers); in other words, denoising is applied to data
at the early stages of the seismic processing chain. This is typically where
the majority of geophysicists focus their denoising efforts - removing noise
before it has a chance to affect down-the-line procedures. However, the
proposed methodology is not restricted to any specific stage of processing. To
illustrate this, the same methodology is applied to a post-stack time-migrated
image from a land field dataset. Similar to the previous examples, a UNet is
used; however due to the reduced density of signals in post-stack data, the
network is made up of only 2 levels prior to the bottle neck. As only one 2D
image is available, the training data is created by dividing it into $7470$
patches of size $32\times 32$, with an overlap of $4$ samples in the $x-$axis
and $6$ samples in the depth-axis.

Figure 9(b) illustrates the results of the blind-spot denoising with 2% of the
patch pixels being selected as active pixels and after training the network
over 250 epochs. The accompanying Jacobian Map, averaged over all the pixels
on inference, is provided in Figure 9(c). The subsequent denoising mask is
selected by identifying any pixels that contribute more than $2\%$, resulting
in the mask portrayed in Figure 9(f). A significant volume of noise is
removed, as evident in the residual Figure 9(g), however an observable level
of signal leakage is now present with the linear interfaces in the upper
section visible within the residual plot.

Figure 9: Application of the XAI-driven denoising procedure on a 2D slice of
land data after post-stack time migration. (a) the raw data, (b) after blind-
spot denoising, and (d) after blind-mask denoising. (c) is the averaged
Jacobian matrix and (f) the respective XAI-driven noise mask. (e) and (g) are
the residuals for the blind-spot and blind-mask procedures, respectively.

##  5 Discussion

###  5.1 Denoising

In the original implementation of N2V, [13] selected between 0.5 and $2\%$ of
the pixels within a training patch to become active pixels in the training
loss. For the seismic denoising scenario, where ‘random’ noise still exhibits
loose coherency with neighbouring points, [14] illustrated how the assumption
of independent and identically distributed noise could be loosened by
significantly increasing the number of active pixels (which will provide a
highly corrupted input to the network), alongside adopting an early stopping
strategy in the network’s training process. In this study, instead of trying
to loosen the requirement of noise being random (or corrupted to appear as
so), we explicitly design a mask that ensures neighbouring pixels cannot
contribute to a central pixel’s predicted value. As such, we follow the
original approach of [13] where only a small percentage of the patch is
corrupted during the training pre-processing and we allow the training process
to run to convergence. One of the benefits of this approach is that the
network is being allowed as much time as is needed to fully learn to replicate
the seismic signal, without much concern being afforded to the network
potentially learning to replicate the coherent noise elements of the data.
Figure 10 provides a comparison for the post-stack field data, using the
initial proposal of [14] where $1/3$ of the total number of pixels are
selected as active pixels and the network is trained for only 15 epochs. By
selecting an appropriate blinding mask, as opposed to excessive corruption of
the training input, the signal leakage has almost been completely eliminated,
whilst the overall volume of noise suppressed has actually increased. Both
examples have the exact same network architecture, training samples and
network hyperparameters.

Figure 10: Comparison of the denoising performance of the blind-spot procedure
of (b) [14] and (c) this study, and (d,e) their associated residuals,
respectively. (a) Original, noisy post-stack field dataset.

Unlike the alternative blind-spot implementation of [25] \- and utilised for
seismic deblending by [18] \- where the receptive field of the neural network
is explicitly restricted to remove the central pixel, the pre-processing
approach of [13] assumes that the replaced central pixel value is never useful
and therefore the network learns implicitly not to use it. Similarly, the
extension to a blind-mask assumes that the network learns not to use the
masked pixels. Throughout, the Jacobian has been used to analyse the
contribution of the neighbouring pixels on a central pixel’s prediction.
Therefore, it can be similarly used to evaluate the efficiency of the masking
procedure. Table 1 details the influence of the active pixel and associated
mask area as a percentage of the total influence the pixels within the
receptive field have on a central pixel’s predicted value. In all scenarios
the active pixel does provide a small contribution towards it’s own predicted
value; thereby, illustrating the network hasn’t fully learnt to ignore the
influence of the active pixel, despite it being corrupted during the training
process. For the masked areas, considering the total contribution, the masked
pixels’ contributions are significantly large, with many contributing more
than $10\%$ towards the predicted value. However, the average contribution per
masked pixel remains below $1\%$ per masked pixel - lower than some of the
contributions of the blind-spot in the standard blind-spot denoising. Whilst
this analysis highlights that the network has not fully learnt to ignore the
masked area as had been hoped, the influence of these pixels does not have a
meaningful value due to their true amplitudes being corrupted (i.e., replaced)
during the training. Therefore, it can be concluded that these pixels’
contributions only serve to introduce noise into the prediction procedure.
Given the promising denoising results, future work will investigate the level
to which the mask efficiency is required in order to achieve the best possible
denoised product.

Table 1: Contribution of masked pixels for the blind-spot and blind-mask denoising schemes, given as a percentage of the total overall pixel contributions within the Jacobian Map.  |  | CGN |
---|---|---|---
| Time Corr. | Iso. | Aniso. | Rig noise
| Active Pixel | 0.8 | 1.2 | 0.9 | 0.5
Blind-spot |  [HTML]C0C0C0Masked Area |  [HTML]C0C0C0 |  [HTML]C0C0C0 |  [HTML]C0C0C0 |  [HTML]C0C0C0
| Active Pixel | 1.6 | 0.7 | 0.7 | 0.7
|  | Total | 13.2 | 6.0 | 11.2 | 14.8
XAI-Guided Blind-mask | Masked Area | Avg./Pixel | 0.4 | 0.7 | 0.7 | 0.6
| Active Pixel | 1.0 | 0.8 | 0.5 | 0.5
|  | Total | 6.0 | 8.8 | 12.9 | 16.2
XAI-Driven Blind-mask | Masked Area | Avg./Pixel | 0.9 | 0.7 | 0.5 | 0.5

###  5.2 Jacobian Analysis

One of the challenges with utilising the Jacobian matrix for self-supervised,
blind-spot denoising is that there is no clear delineation between pixels that
contribute towards the noise component of the central pixel versus those that
contribute towards the signal component. In our synthetic analysis, we
proposed the use of an XAI-guided procedure where the researcher could utilize
their prior knowledge of the pixel contributions as a complement to the output
of the Jacobian matrix analysis. Such prior knowledge could be knowledge of
the noise statistics helping the researcher to identify any highly
contributing pixels that may be providing only signal information.
Alternatively, it could be knowledge from a prior synthetic investigation
which identified the pixels that are likely to contribute signal information.
Whilst the XAI-guided approach allows the incorporation of prior knowledge as
a means to reduce the masking of pixels that provide useful information on the
seismic signal, this introduces some manual labour in the process. To the
contrary, the XAI-driven methodology is more heavy-handed with respect to the
selection of the mask, however it is a fully automated procedure that requires
no prior knowledge on the monitoring conditions or the properties of the
seismic signals. As the mask design in the XAI-guided approach is subjective,
a direct comparison between XAI-guided and XAI-driven is not possible.
However, in the numerical experiment section of this paper we show that both
are suitable to the synthetic case. Due to a lack of prior knowledge, only the
XAI-driven approach was implemented for the field data studies.

##  6 Conclusions

Self-supervised, blind-spot networks have been shown to be powerful denoisers
for seismic data contaminated by random noise; however, they require an
additional masking step for the suppression of coherent noise. The mask design
is often non-trivial, requiring extensive knowledge of the noise
characteristics and how the network will discriminate between them and the
signal that we wish to preserve. In this study, we have proposed the use of
bias-free networks, which allow the denoising process to be interpreted as a
linear combination of the input pixels, where the weights can be identified by
computing the Jacobian matrix of the network itself; ultimately allowing the
identification of pixels with a strong influence on any predicted value. This
knowledge is utilised to automatically design noise masks to hide said pixels
during the denoising networks training process. The resulting blind-mask
networks are shown to successfully suppress a range of coherent noise types on
realistic synthetic seismic shot gathers. Furthermore, the methodology is
applied to two field datasets at different ends of the processing pipeline -
common shot gathers and post-stack images. The successful suppression of noise
in the field datasets highlights not only the robust, denoising capability of
the proposed procedure, but also its general applicability to a variety
seismic denoising applications. To conclude, the proposed procedure offers an
automated, deep learning-based workflow that operates directly on the raw
field data with no requirement of prior knowledge on the site conditions.

##  7 Acknowledgements

The authors thank the KAUST Deep Imaging Group and Seismic Wave Analysis Group
for insightful discussions, and in particular, Professor Tariq Alkhalifah, Dr
Nick Luiken, and Ms Sixiu Liu. For computer time, this research used the
resources of the Supercomputing Laboratory at King Abdullah University of
Science & Technology (KAUST) in Thuwal, Saudi Arabia. CB kindly acknowledges
funding from Saudi Aramco for aspects of this work.

##  8 Data Availability Statement

All data used in this study are openly available within the public domain. The
Volve synthetic waveform data is available on the Zenodo data platform, at
http://doi.org/10.5281/zenodo.6572286. The Volve field data is available
directly from Equinor ASA at https://www.equinor.com/energy/volve-data-
sharing. The post-stack image is available in the Madagascar data repository,
at https://github.com/ahay/src. The code to generate the synthetic noise
instances is available from the authors by request.

## References

  * [1] Chongyi Li, Chunle Guo, Linghao Han, Jun Jiang, Ming-Ming Cheng, Jinwei Gu, and Chen Change Loy.  Low-light image and video enhancement using deep learning: A survey.  IEEE transactions on pattern analysis and machine intelligence, 44(12):9396–9416, 2021.
  * [2] Kenji Suzuki.  Overview of deep learning in medical imaging.  Radiological physics and technology, 10(3):257–273, 2017.
  * [3] Ellappan Mohan, Arunachalam Rajesh, Gurram Sunitha, Reddy Madhavi Konduru, Janagaraj Avanija, and Loganathan Ganesh Babu.  A deep neural network learning-based speckle noise removal technique for enhancing the quality of synthetic-aperture radar images.  Concurrency and Computation: Practice and Experience, 33(13):e6239, 2021.
  * [4] Yuxing Zhao, Yue Li, Xintong Dong, and Baojun Yang.  Low-frequency noise suppression method based on improved dncnn in desert seismic data.  IEEE Geoscience and Remote Sensing Letters, 16(5):811–815, 2018\.
  * [5] Feng Qian, Wei Guo, Zhangbo Liu, Hongtao Yu, Gulan Zhang, and Guangmin Hu.  Unsupervised erratic seismic noise attenuation with robust deep convolutional autoencoders.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–16, 2022\.
  * [6] Harpreet Kaur, Sergey Fomel, and Nam Pham.  Seismic ground-roll noise attenuation using deep learning.  Geophysical Prospecting, 68(7):2064–2077, 2020.
  * [7] Xintong Dong, Hongzhou Wang, Tie Zhong, and Yue Li.  Desert mixed seismic noise suppression by using multiple forward models and a supervised deep-learning method.  Exploration Geophysics, 52(4):431–445, 2021.
  * [8] Tie Zhong, Ming Cheng, Shaoping Lu, Xintong Dong, and Yue Li.  Rcen: a deep-learning-based background noise suppression method for das-vsp records.  IEEE Geoscience and Remote Sensing Letters, 19:1–5, 2021.
  * [9] Claire Birnie and Tariq Alkhalifah.  Transfer learning for self-supervised, blind-spot seismic denoising.  Frontiers in Earth Science, 10:1053279, 2022.
  * [10] Omar M Saad and Yangkang Chen.  Deep denoising autoencoder for seismic random noise attenuation.  Geophysics, 85(4):V367–V376, 2020.
  * [11] Rosita Hamidi, Abdul Halim Latif, and Wei Yi Lee.  Autoencoder neural network application for coherent noise attenuation in high frequency shallow marine seismic data.  In Offshore Technology Conference Asia, page D012S001R086. OTC, 2020\.
  * [12] Magdalena Markovic, Reza Malehmir, and Alireza Malehmir.  Diffraction pattern recognition using deep semantic segmentation.  Near Surface Geophysics, 20(5):507–518, 2022.
  * [13] Alexander Krull, Tim-Oliver Buchholz, and Florian Jug.  Noise2void-learning denoising from single noisy images.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 2129–2137, 2019.
  * [14] Claire Birnie, Matteo Ravasi, Sixiu Liu, and Tariq Alkhalifah.  The potential of self-supervised networks for random noise suppression in seismic data.  Artificial Intelligence in Geosciences, 2:47–59, 2021.
  * [15] Coleman Broaddus, Alexander Krull, Martin Weigert, Uwe Schmidt, and Gene Myers.  Removing structured noise with self-supervised blind-spot networks.  In 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI), pages 159–163. IEEE, 2020.
  * [16] Sixiu Liu, Claire Birnie, and Tariq Alkhalifah.  Coherent noise suppression via a self-supervised blind-trace deep learning scheme.  arXiv preprint arXiv:2206.00301, 2022.
  * [17] Shirui Wang, Wenyi Hu, Pengyu Yuan, Xuqing Wu, Qunshan Zhang, Prashanth Nadukandi, German Ocampo Botero, and Jiefu Chen.  Seismic deblending by self-supervised deep learning with a blind-trace network.  In SEG/AAPG/SEPM First International Meeting for Applied Geoscience & Energy. OnePetro, 2021.
  * [18] Nick Luiken, Matteo Ravasi, and Claire E Birnie.  A hybrid approach to seismic deblending: when physics meets self-supervision.  arXiv preprint arXiv:2205.15395, 2022.
  * [19] L Mosser and B Alaei.  Deep learning-based unsupervised denoising for coherent noise attenuation in seismic data.  In Second EAGE Subsurface Intelligence Workshop, volume 2022, pages 1–5. EAGE Publications BV, 2022.
  * [20] Ryan Benkert, Oluwaseun Joseph Aribido, and Ghassan AlRegib.  Explainable seismic neural networks using learning statistics.  In First International Meeting for Applied Geoscience & Energy, pages 1425–1429. Society of Exploration Geophysicists, 2021.
  * [21] Kyubo Noh, Dowan Kim, and Joongmoo Byun.  Intrinsically explainable deep learning for supervised seismic facies classification using intrinsic method.  IEEE Transactions on Geoscience and Remote Sensing, 2023.
  * [22] Danilo Chamorro, Jiahua Zhao, Claire Birnie, Myrna Staring, Fliedner Moritz, and Matteo Ravasi.  Deep learning-based extraction of surface wave dispersion curves from seismic shot gathers.  arXiv preprint arXiv:2305.13990, 2023.
  * [23] Antoni Buades, Bartomeu Coll, and J-M Morel.  A non-local algorithm for image denoising.  In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 2, pages 60–65. Ieee, 2005.
  * [24] Matteo Ravasi, Tamin Selvan, and Nick Luiken.  Stochastic multi-dimensional deconvolution.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–14, 2022\.
  * [25] Samuli Laine, Tero Karras, Jaakko Lehtinen, and Timo Aila.  High-quality self-supervised deep image denoising.  Advances in Neural Information Processing Systems, 32, 2019.


11institutetext: Space Science Division, Naval Research Laboratory, Code
7674L, Washington DC 20375

# On the Detectability of 57Fe Axion-Photon Mode Conversion in the Sun

J. Martin Laming 11

###### Abstract

Aims. The purpose of this paper is to assess the feasibility of axion
detection by X-ray spectroscopy of the sun.

Methods. We review the theory of axion-photon mode conversion with special
attention to axions emitted in the 14.4 keV M1 decay of 57Fe at the solar
center. These then mode convert to photons in the outer layers of the solar
envelope, and may in principle be detected subsequently as X-rays.

Results. For axion masses above about $10^{-4}$ eV, resonant mode conversion
at a layer where the axion mass matches the local electron plasma frequency is
necessary. For axion masses above about $10^{-2}$ eV, this mode conversion
occurs too deep in the solar atmosphere for the resulting photon to escape the
solar surface and be detected before Compton scattering obscures the line. At
the (detectable) axion masses below this, the flux of mode converted photons
predicted by axion models appears to be too low for detection to be feasible
with current instrumentation. Nonresonant mode conversion for axion masses
below $10^{-4}$ eV is also plausible, but with still lower predicted fluxes,
since the axion coupling constant is related to it mass.

Conclusions. Prospects for meaningful constraints on massive axion parameters
from X-ray observations of this transition from the Sun do not appear to be
promising. However parameters for massless counterparts (e.g. the “arion”) may
still result from such observations. It may mode convert in the outer layers
of the solar atmosphere, but is not restricted by this to have a small
coupling constant.

###### Key Words.:

Elementary Particles – Sun: X-rays, gamma rays – Sun: Particle Emission

##  1 Introduction

The basic form of the Lagrangian for the theory of strong interactions known
as Quantum Chromodynamics (QCD) contains a pseudo-scalar term involving the
gluon field strength tensor that violates symmetry under the actions of charge
and parity reversal (CP violation). One practical consequence of this is the
prediction of an electric dipole moment for the neutron in excess of the
experimental limits by about ten orders of magnitude. The most generally
accepted solution to this problem involves the introduction of a new global
chiral symmetry (Peccei & Quinn (1977a, b)). Weinberg (1978) and Wilczek
(1978) pointed out that such a symmetry, being spontaneously broken,
necessitates the existence of a new boson, the axion, with mass in the range
$10^{-12}$ eV to 1 MeV, depending on the scale at which the global symmetry is
broken. The axion mass is related to the symmetry breaking scale, $f_{PQ}$, by
(e.g. Moriyama (1995))

$m=0.62\left(10^{7}{\rm GeV}/f_{PQ}\right)\quad{\rm eV.}$ (1)

Limits on the axion mass are placed by various astrophysical constraints.
There are two most often discussed classes of axion; “hadronic” axions that do
not couple to leptons and quarks at tree level (Kim (1979); Shifman et al.
(1980)) and DFSZ axions that do (Dine et al. (1981); Zhitnitskii (1980)). Both
types of axion are constrained to have masses less than about
$10^{-3}-10^{-2}$ eV based on analyses of the neutrino burst from SN 1987A
(Burrows et al. (1989); Keil et al. (1997)), while hadronic axions may also
have a mass of a few eV, bounded from below by treatments of axion trapping in
SN 1987A (Burrows et al. (1990)) and from above by studies of red giant
evolution (Haxton & Lee (1991)). Studies of white dwarf cooling and pulsation
also indicate DFSZ axion masses of order $10^{-3}$ eV (Isern et al. (2010)).
Cosmological arguments (e.g. Kolb & Turner (1994)) give a lower limit on the
axion mass in the range $10^{-6}-10^{-5}$ eV, with a mass of $\sim 10^{-5}$ eV
supplying the measured cold dark matter density (Erken et al. (2011)).

The sun is the only astrophysical object besides SN 1987A from which neutrinos
have been directly detected, and so it is natural to consider it also as a
potential source of axions. While thermal processes at the solar center
produce axions by Compton and Primakoff photoproduction processes, illustrated
in Figure 1, (Chanda et al. (1988)), considerable attention has also been
devoted to a transition at 14.4 keV in the 57Fe nucleus, which lies in a
spectral region with essentially zero background from the nonflaring solar
atmosphere. Being an M1 transition, it has a branching ratio for axion
emission (DFSZ or hadronic) instead of photon radiation (Haxton & Lee (1991)),
and has a sufficiently low threshold to allow thermal excitation. An estimate
of the axion flux in this line has been given by Andriamonje et al. (2009) as
$F_{A}=4.56\times 10^{23}\left(g_{aN}^{eff}\right)^{2}$ cm-2s-1, where
$g_{aN}^{eff}$ is the effective axion-nucleon coupling. Moriyama (1995) gives
$F_{A}=5.3\times 10^{8}\left(m/1{\rm eV}\right)^{2}$ cm-2s-1 assuming the
naive quark model (NQM) where $g_{aN}^{eff}=3.53\times 10^{-8}\left(m/1{\rm
eV}\right)$. Andriamonje et al. (2009) give $g_{aN}^{eff}=1.48\times
10^{-8}\left(m/1{\rm eV}\right)$ (in the calculation of the 57Fe M1 axion
branching ratio), which would give a flux of $F_{A}=10^{8}\left(m/1{\rm
eV}\right)^{2}$ cm-2s-1 but also discuss higher possible values for
$g_{aN}^{eff}$.

Figure 1: Feynman diagrams illustrating the two Compton scattering diagrams
(left) and the Primakoff process (right).

Andriamonje et al. (2009) describe the CERN Axion Solar Telescope (CAST),
designed to detect solar 57Fe axions by observing the 14.4 keV photons
produced when the axions mode convert inside the detector. The parameters of
CAST; a 9T magnetic field extending over a length of $s=9.26$ m, give a
product $Bs\simeq 8\times 10^{7}$ G cm, which is smaller than the product of a
typical solar magnetic field and density scale height, and so the question of
whether 14.4 keV photons might be detectable in the solar spectrum as a
signature of in situ axion mode conversion becomes interesting to consider.

##  2 Axion-Photon Interactions

The axion-photon interaction Lagrangian is

$L=-{g_{a\gamma\gamma}\over
4}aF^{\mu\nu}\tilde{F}_{\mu\nu}=g_{a\gamma\gamma}a{\bf E}\cdot{\bf B}$ (2)

where $g_{a\gamma\gamma}$ is the axion-photon coupling constant, $a$ is the
axion field, $F^{\mu\nu}$ is the electromagnetic field tensor and
$\tilde{F}_{\mu\nu}$ is its dual. The product
$F^{\mu\nu}\tilde{F}_{\mu\nu}=-4{\bf E}\cdot{\bf B}$, which leads to the final
form. Here, ${\bf E}$ is interpreted as the photon electric field and ${\bf
B}$ the ambient magnetic field. Thus axions only interact with photons when
propagating across the ambient magnetic field, and then only with the mode
polarized along the magnetic field, known as the ordinary mode.

The effective Lagrangian for an axion-photon system is then

$L={1\over 2}{\bf E}\cdot{\bf\epsilon}\cdot{\bf E}-{1\over 2}{\bf B}\cdot{\bf
B}+{1\over 2}\left(\partial_{\mu}a\right)^{2}-{1\over
2}m^{2}a^{2}+g_{a\gamma\gamma}a{\bf E}\cdot{\bf B}$ (3)

where $m$ is the axion mass and ${\bf\epsilon}=1-\omega_{p}^{2}/\omega^{2}$ is
the dielectric tensor in terms of the electron plasma frequency $\omega_{p}$.
We are working in the limit of small magnetic fields. By way of contrast,
Yoshimura (1988) includes the magnetic field terms in ${\bf\epsilon}$. The
equations of motion are

$\displaystyle{\bf\epsilon}\cdot{\bf E}-{k^{2}\over\omega^{2}}{\bf
E}+g_{a\gamma\gamma}a{\bf B}=0$ (4) $\displaystyle{g_{a\gamma\gamma}{\bf
E}\cdot{\bf
B}\over\omega^{2}}+\left(1-{m^{2}\over\omega^{2}}-{k^{2}\over\omega^{2}}\right)a=0,$
(5)

which for nontrivial solutions require

$\displaystyle n^{2}$ $\displaystyle={k^{2}\over\omega^{2}}$ (6)
$\displaystyle={1\over
2}\left(\epsilon+1-{m^{2}\over\omega^{2}}\right)\pm{1\over
2}\sqrt{\left(\epsilon-1+{m^{2}\over\omega^{2}}\right)^{2}+4{g_{a\gamma\gamma}^{2}B^{2}\over\omega^{2}}}.$
(7)

These two solutions represent axion and photon propagation across the magnetic
field. We rewrite the photon-axion mass matrix from equation 4 in symmetric
form as

$M=\left(\begin{array}[]{cc}{m^{2}-\omega_{p}^{2}\over
2\omega}&g_{a\gamma\gamma}B\\\ g_{a\gamma\gamma}B&-{m^{2}-\omega_{p}^{2}\over
2\omega}\\\ \end{array}\right).$ (8)

Pre and post multiplying by a rotation matrix,

$\displaystyle\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\\\
\end{array}\right)M\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\\\ \end{array}\right)=$ (13)
$\displaystyle\left(\begin{array}[]{cc}{m^{2}-\omega_{p}^{2}\over 2\omega}\cos
2\theta-g_{a\gamma\gamma}B\sin 2\theta&{m^{2}-\omega_{p}^{2}\over 2\omega}\sin
2\theta+g_{a\gamma\gamma}B\cos 2\theta\\\ {m^{2}-\omega_{p}^{2}\over
2\omega}\sin 2\theta+g_{a\gamma\gamma}B\cos
2\theta&-{m^{2}-\omega_{p}^{2}\over 2\omega}\cos
2\theta+g_{a\gamma\gamma}B\sin 2\theta\\\ \end{array}\right),$ (16)

which is diagonal for a mixing angle $\theta$ given by $\tan
2\theta=2g_{a\gamma\gamma}B\omega/\left(\omega_{p}^{2}-m^{2}\right)$.

The evolution of the axion-photon mixed state with time $t$ is described by

$\displaystyle\left|a\left(t\right)\right>=e^{-iE_{a}t}\left|a^{\prime}\left(0\right)\right>\cos\theta+e^{-iE_{\gamma}t}\left|\gamma^{\prime}\left(0\right)\right>\sin\theta$
(17)
$\displaystyle=e^{-iE_{a}t}\left|a\left(0\right)\right>\cos^{2}\theta+e^{-iE_{\gamma}t}\left|a\left(0\right)\right>\sin^{2}\theta$
(18)
$\displaystyle+\cos\theta\sin\theta\left(e^{-iE_{a}t}-e^{-iE_{\gamma}t}\right)\left|\gamma\left(0\right)\right>.$
(19)

The probability of axion to photon conversion is then

$\displaystyle
P\left(a\rightarrow\gamma\right)=\left|\left<\gamma|a\left(t\right)\right>\right|^{2}$
$\displaystyle=\cos^{2}\theta\sin^{2}\theta\left|e^{-iE_{a}t}-e^{-iE_{\gamma}t}\right|^{2}$
(20)
$\displaystyle=\sin^{2}2\theta\sin^{2}{\left(E_{a}-E_{\gamma}\right)t\over 2}$
(21)

which upon substituting for the mixing angle becomes

$P\left(a\rightarrow\gamma\right)={g_{a\gamma\gamma}^{2}B^{2}\over
g_{a\gamma\gamma}^{2}B^{2}+\left(\omega_{p}^{2}-m^{2}\right)^{2}/4\omega^{2}}\sin^{2}{\left(E_{a}-E_{\gamma}\right)t\over
2}.$ (22)

We put $\left(E_{a}-E_{\gamma}\right)t=\omega s\left(n_{a}-n_{\gamma}\right)$
with $n_{a,\gamma}^{2}$ given by the two solutions in equation 5, and $s$ the
distance traveled through the medium. Then
$E_{a}-E_{\gamma}\simeq\sqrt{\left(m^{2}-\omega_{p}^{2}\right)^{2}/\omega^{2}+4g_{a\gamma\gamma}^{2}B^{2}}/2$
and

$\displaystyle P\left(a\rightarrow\gamma\right)$
$\displaystyle={g_{a\gamma\gamma}^{2}B^{2}\over
g_{a\gamma\gamma}^{2}B^{2}+\left(\omega_{p}^{2}-m^{2}\right)^{2}/4\omega^{2}}\times$
(23) $\displaystyle\sin^{2}{{s\over
2}\sqrt{\left(m^{2}-\omega_{p}^{2}\right)^{2}/4\omega^{2}+g_{a\gamma\gamma}^{2}B^{2}}}$
(24) $\displaystyle\rightarrow s^{2}g_{a\gamma\gamma}^{2}B^{2}/4$ (25)

when
$\left(s/2\right)\sqrt{g_{a\gamma\gamma}^{2}B^{2}+\left(\omega_{p}^{2}-m^{2}\right)^{2}/4\omega^{2}}\rightarrow
0$.

The forgoing estimate assumed that the electron density in constant. When
crossing a resonance, the rapid variation in $\tan 2\theta$ with
$\omega_{p}^{2}$ first greater than $m^{2}$ becoming equal and then less than
$m^{2}$, as the axion moves from high density to low density, means that the
axion-photon oscillation cannot keep up with the changing Hamiltonian. In this
limit, the equations of motion (4) admit solutions in terms of parabolic
cylinder functions, yielding an axion transmission probability of
$\exp\left(-\pi g_{a\gamma\gamma}^{2}B^{2}kh/\omega_{p}^{2}\right)$, (Cairns &
Lashmore-Davies (1983); Cally (2006)). The corresponding mode conversion
probability is then

$P\left(a\rightarrow\gamma\right)=1-\exp\left(-\pi
g_{a\gamma\gamma}^{2}B^{2}kh/\omega_{p}^{2}\right)\simeq\pi
g_{a\gamma\gamma}^{2}B^{2}kh/\omega_{p}^{2},$ (27)

and as will be seen below, is generally larger than the nonresonant conversion
probability given in equation 11.

##  3 In Real Numbers …

We have used natural units, where $\hbar=c=1$, above. To put variables in
units more familiar to solar physicists, we use Andriamonje et al. (2009)

$g_{a\gamma\gamma}=8.36\times 10^{-4}/f_{PQ}\quad{\rm GeV}^{-1},$ (28)

concentrating on DFSZ axions for the time being, where $f_{PQ}$ is the
symmetry breaking scale in GeV, and hence from equation 1

$m=0.62\left(10^{7}{\rm
GeV}/f_{PQ}\right)=0.74\left(g_{a\gamma\gamma}/10^{-10}{\rm
GeV}^{-1}\right){\rm eV}.$ (29)

With these definitions

$\displaystyle g_{a\gamma\gamma}B$
$\displaystyle=10^{-16}B\left(g_{a\gamma\gamma}/10^{-10}{\rm GeV}^{-1}\right)$
(30) $\displaystyle=1.3\times 10^{-16}B\left(m/1{\rm eV}\right){\rm cm}^{-1}$
(31)

where $B$ is given in Gauss, and the conversion probability is

$P\left(a\rightarrow\gamma\right)=4.5\times 10^{-33}B^{2}\left({\rm
G}\right)\left(m/1{\rm eV}\right)^{2}s^{2}.$ (32)

At a photon energy of 14.4 keV,

$\displaystyle\omega_{p}^{2}/\omega$ $\displaystyle=4.8\times
10^{-21}n_{e}{\rm cm}^{-1}$ (33) $\displaystyle m^{2}/\omega$
$\displaystyle=3.53\left(m/1{\rm eV}\right)^{2}{\rm cm}^{-1}.$ (34)

For $m>10^{-4}$ eV and $s>10^{8}$ cm (a typical solar atmosphere scale size)
the small angle approximation for the sine in equation 11 is not valid. Both
terms in equation 17 are $>>g_{a\gamma\gamma}^{2}B^{2}$ and so unless
$\omega_{p}\simeq m$, the conversion probability is significantly reduced from
the value in equation 16. The electron density at which resonant mode
conversion can occur is (see also Yoshimura (1988))

$n_{e}\simeq 7.35\times 10^{20}\left(m/1{\rm eV}\right)^{2}.$ (36)

The electron scattering optical depth between this layer and the solar surface
is

$\tau\simeq n_{e}\sigma_{T}h\simeq 3\times 10^{-17}n_{e}$ (37)

where $\sigma_{T}$ is the Thomson cross section and $h\simeq 500$ km is the
atmosphere density scale height. Figure 2 shows emergent line profiles at the
top of the solar atmosphere after traveling through electron scattering
optical depths of 0.1, 0.3, 1.0, 3.0, and 5.0, calculated following the
treatment of Chandrasekhar (1960) and Münch (1948) for plane parallel
atmospheres. As can be seen, the detectable line disappears once $\tau\sim 5$,
meaning a maximum density at which detectable mode conversion will occur is
$\sim 10^{17}$ cm-3, with a maximum detectable axion mass of $1.2\times
10^{-2}$ eV. At higher opacities, a broad Compton scattered feature may still
be seen. Figure 3 shows this feature at opacities of 5, 15, 50, and 150. It is
visible as a 1-2 keV wide feature out to an opacity of between 15 and 50. This
would increase the maximum detectable axion mass to about $4-7\times 10^{-2}$
eV, but would probably require more photons detected due to the broader
feature.

Figure 2: Line profiles of the emergent 14.4 keV 57Fe line, originally emitted
as an axion, and mode converted at various optical depths in the solar
atmosphere. The electron scattering opacities through which the various
profiles have propagated are noted. The dashed curve gives the original
unscattered profile. At an opacity of $\simeq 5$, the central narrow line has
essentially disappeared.

Figure 3: Same as Figure 1, but for electron scattering opacities of 5, 15,
50, and 150. The $\tau=5$ curve is truncated for numerical reasons.

Equation 12 gives the mode conversion probability as $P=7.5\times
10^{-25}B^{2}$, independent of axion mass. For $B\sim 10^{3}$G, $P\sim
10^{-18}$, and the photon flux at earth is $10^{-10}\left(m/{\rm
1eV}\right)^{2}$ cm-2s-1, assuming the axion flux calculated by Andriamonje et
al. (2009). For axions with mass up to $10^{-2}$ eV (the high end of the mass
range allowed by SN 1987A observations, and close to the limit allowed by
Compton scattering in the solar atmosphere), this leads to a count rate of
order $10^{-11}$ s-1 in a instrument of effective area $10^{3}$ cm2 (e.g.
NUSTAR; Harrison et al. (2010, 2013)); too low for detection. Significantly
higher conversion probabilities would be expected in pulsar strength magnetic
fields, (Yoshimura (1988); Perna et al. (2012)).

For $s\sim 10^{8}$ cm, the argument of the sine in equation 11 becomes small
for axion masses at the low end of the allowed range, $m<10^{-4}$ eV. In this
case nonresonant conversion can occur with probability
$s^{2}g_{a\gamma\gamma}^{2}B^{2}/4\simeq 10^{-19}$ for $B\simeq 10^{3}$ G and
$m\simeq 10^{-4}$ eV. This yields a flux at earth of mode converted photons of
$5\times 10^{4}\left(g_{aN}^{eff}\right)^{2}$. Taking $g_{aN}^{eff}=1.48\times
10^{-8}\left(m/1{\rm eV}\right)$ as before gives a negligible mode converted
photon flux at earth of $\sim 10^{-19}$ cm-2s-1 with $m=10^{-4}$ eV. The lower
axion mass leads to a significantly lower flux, through its coupling to
$g_{a\gamma\gamma}$. In hadronic axion models, the parameters
$g_{a\gamma\gamma}$ and $g_{aN}^{eff}$ are less tightly constrained. Yoshimura
(1988) suggests $g_{a\gamma\gamma}$ may increase by a factor of ten, which
would lead to a commensurate increase in the count rates estimated above,
though still insufficient to change the conclusions above.

Figure 4: The $g_{a\gamma\gamma}-g_{aN}^{eff}$ parameter space. The dark
region at the top and right indicates the region ruled out by CAST
(Andriamonje et al. 2009). The horizontal portion labeled (luminosity) comes
from a constraint on the total Fe57 axion luminosity as a fraction of the
solar luminosity, and the vertical section labeled (Primakoff) comes from CAST
constraints on axions produced by the Primakoff process, both from Andriamonje
et al. (2009). The dashed contour below this indicates the where this limit
moves with the nondetection of a count rate of $10^{-5}$ s-1 by NUSTAR (a flux
of $10^{-8}$ photons cm-2s-1) for resonant axion-photon mode conversion. The
shaded contours indicate the “forbidden” mass regions ($>10^{-2}$ eV) towards
the upper right for DFSZ or hadronic axions corresponding to the values of
$g_{a\gamma\gamma}$ and $g_{aN}^{eff}$ (assuming the naive quark model). The
lower dashed contour indicates the limit achieved by the same NUSTAR count
rate for nonresonant axion-photon mode conversion in the outer layers of the
solar atmosphere. The rectangular boxes indicate the “allowed” mass limits
($<10^{-4}$ eV) for this case. These “forbidden” or “allowed” mass regions do
not apply in the case of massless arions; only the CAST limits apply. Clearly,
substantial parameters space is available to solar observations detecting
nonresonant mode conversion of arions.

These potential constraints are illustrated in Fig. 4., showing the
$g_{a\gamma\gamma}-g_{aN}^{eff}$ parameter space. The top and right hand axes
show the corresponding axion masses in the cases of DFSZ axions (equation 1),
and those coming from the naive quark model from $g_{aN}^{eff}$ respectively.
The dark region at the top and right indicates the region ruled out by CAST
(Andriamonje et al. 2009), and other constraints discussed therein. The dashed
contour below this indicates the where this limit moves with the nondetection
of a count rate of $10^{-5}$ s-1 by NUSTAR (a flux of $10^{-8}$ photons
cm-2s-1, and about 10 times the background count rate in this rate estimated
from Table 2 in Harrison et al. 2013) for resonant axion-photon mode
conversion. The shaded contours indicate the allowed mass region ($<10^{-2}$
eV)in the lower left for DFSZ or hadronic axions. The lower dashed contour
indicates the limit achieved by the same NUSTAR count rate for nonresonant
axion-photon mode conversion in the outer layers of the solar atmosphere. The
rectangular boxes indicate the mass limits ($<10^{-4}$ eV) for this case. As
can be seen, in either resonant or nonresonant cases, the expected flux at the
appropriate axion mass is well below plausible detection limits.

An exception might arise in the case of pseudo-Goldstone bosons arising from
spontaneously broken family symmetry, which are naturally massless (Wilczek
(1982)), but may otherwise have similar properties to the axion (Anselm
(1988)). Berezhiani & Khlopov (1990a), Berezhiani & Khlopov (1990b), and
Berezhiani & Khlopov (1991) give further discussion of the origin of such a
particle, dubbed the “arion”, in spontaneously broken family symmetry. Being
massless, it may be detected following mode conversion anywhere in the solar
atmosphere where the electron scattering opacity $\tau<1$, which roughly
corresponds to the nonresonant case in Figure 4. In this case
$g_{a\gamma\gamma}$ and $g_{aN}^{eff}$ are disconnected from the mass, and the
forbidden regions of Figure 4. for DFSZ and Hadronic axions with mass above
$10^{-2}$ eV no longer apply (the CAST forbidden region is still applicable).
In such a case, NUSTAR might provide limits parameters competitive with those
derived by other methods. For DFSZ axions, studies of the Primakoff process in
massive stars limit $g_{a\gamma\gamma}<0.8\times 10^{-10}$ GeV-1 (Friedland et
al. (2012)), a limit that would be accessible to solar observations. A similar
limit comes from CAST observations of Primakoff solar axions (Andriamonje et
al. (2007)).

As mentioned above, axions may also be produced by thermal processes at the
solar center (Chanda et al. (1988)). Hadronic axions may also be produced by
the Primakoff process. Compton scattering and electron-nucleus bremsstrahlung
can also produce DFSZ axions. The Primakoff process alone produces around 10 -
100 times more axions (Andriamonje et al. (2009)) than the 57Fe decay, while
the other two process increase the axion flux by approximately $10^{3}$
(Moriyama (1995); Derbin et al. (2011)). The count rate so produced according
to axion models would be $\sim 10^{-8}$ s-1 at a mass of $10^{-2}$ eV, still
too low for observational feasibility, but possibly adequate to detect the
arion. However these axions, and consequently the photons produced by their
mode conversion, dominantly have energies in the 1 - 3 keV range (Derbin et
al. (2011)). At these energies, the dominant opacity is photoelectric
absorption rather than Compton scattering, and the depth of solar atmosphere
from which mode converted photons may be detected is smaller. In this spectral
band, the sun also emits a thermal bremsstrahlung spectrum, as well as
spectral lines from active regions and flares, in contrast to the 14.4 keV
region, where no solar emission is expected outside of flare. The fundamental
problem would be that if a definitive detection is not expected, to be of
value, the observation needs to be designed so that a null result produces a
further restriction of the observationally allowed parameter space. Hudson et
al. (2012) discuss progress and limitations in this approach. Davoudiasl &
Huber (2006) suggest observing axion mode conversion in the magnetic field of
the earth, in such a manner that the earth itself blocks out the thermal
continuum from the sun. While this estimate uses the observationally allowed
axion parameter space rather than a specific model to arrive at a count rate,
it does have the advantage that a null result would be more readily
interpreted.

##  4 Conclusions

The detection of photons from axion mode conversion within the solar
atmosphere is practically limited to axion masses at the higher end of the
mass window, $\sim 10^{-5}-10^{-2}$ eV, (bounded by cosmological arguments
from below, and by stellar luminosities and SN 1987A from above). Detection of
solar axions is also necessarily limited to masses below $10^{-2}$ eV due to
the Compton scattering expected between the mode conversion layer (deeper in
the solar atmosphere for higher mass axions) and the solar surface. This
restriction means that current instrumentation is unable to provide meaningful
constraints on theoretical axion models, because axion coupling constants
scale as the axion mass, $m$, and so the mode converted signal scales as
$m^{4}$, leading to an intrinsically weak signal for low mass axions. A
possible exception to this might a massless counterpart to the axion, the
arion, where $g_{a\gamma\gamma}$ and $g_{aN}^{eff}$ are now decoupled from the
mass, and consequently the “forbidden” regions in Figure 4. do not apply.

###### Acknowledgements.

This work has been supported by basic research funds of the Office of Naval
Research.

## References

  * Andriamonje et al. (2007) Andriamonje, S., et al. (CAST Collaboration) 2007, J. Cosmology and Astroparticle Phys., 4, 10
  * Andriamonje et al. (2009) Andriamonje, S., Aune, S., Autiero, D., et al. 2009, J. Cosmology and Astroparticle Phys., 12, 2
  * Anselm (1988) Anselm, A. A. 1988, Phys. Rev. D, 37, 2001
  * Berezhiani & Khlopov (1990a) Berezhiani, Z. G, & Khlopov, M. Yu., 1990a, Sov. J. Nucl. Phys. 51, 739
  * Berezhiani & Khlopov (1990b) Berezhiani, Z. G, & Khlopov, M. Yu., 1990b, Sov. J. Nucl. Phys. 51, 935
  * Berezhiani & Khlopov (1991) Berezhiani, Z. G, & Khlopov, M. Yu., 1991, Z. Phys. C 49, 73
  * Burrows et al. (1989) Burrows, A., Turner, M. S., & Brinkmann, R. P. 1989, Phys. Rev. D, 39, 1020
  * Burrows et al. (1990) Burrows, A., Ressell, M. T., & Turner, M. S. 1990, Phys. Rev. D, 42, 3297
  * Cairns & Lashmore-Davies (1983) Cairns, R. A., & Lashmore-Davies, C. N. 1983, Phys. Fluids, 26, 1268
  * Cally (2006) Cally, P. S. 2006, Phil. Trans. R. Soc. A, 364, 333
  * Chanda et al. (1988) Chanda, R., Nieves, J. F., & Pal, P. B. 1988, Phys. Rev. D, 37, 2714
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Radiative Transfer, (Dover: New York)
  * Davoudiasl & Huber (2006) Davoudiasl, H., & Huber, P. 2006, Phys. Rev. Lett., 97, 141302
  * Derbin et al. (2011) Derbin, A. V., Kayunov, A. S., Muratova, V. V., Semenov, D. A., & Unzhakov, E. V. 2011, Phys. Rev. D, 83, 3505
  * Dine et al. (1981) Dine, M., Fischler, W., & Srednicki, M. 1981, Phys. Lett., 104B, 199
  * Erken et al. (2011) Erken, O., Sikivie, P., Tam, H., & Yang, Q. 2011, Phys. Rev. Lett., 108, 061304
  * Friedland et al. (2012) Friedland, A., Giannotti, M., & Wise, M. 2012, arXiv:1210.1271
  * Harrison et al. (2010) Harrison, F. A., Boggs, S., Christensen, F., et al. 2010, Proc. SPIE 7732, 77320S-1
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F., et al. 2013, ApJ, submitted, arXiv.1301:7307
  * Haxton & Lee (1991) Haxton, W. C., & Lee, K. Y. 1991, Phys. Rev. Lett., 66, 2557
  * Hudson et al. (2012) Hudson, H. S., Acton, L. W., DeLuca, E., Hannad, I. G., Reardon, K., & Van Bibber, K. 2012, in Proc. 4th Hinode Science Meeting ”Unsolved Problems and Recent Insights”, ed. L. R. Bellot Rubio, F. Reale, & M. Carlsson, arXiv:1201.4607
  * Isern et al. (2010) Isern, J., García-Berro, E., Althaus, L. G., & Córsico, A. H., 2010, A&A, 512, A86
  * Keil et al. (1997) Keil, W., Janka, H.-T., Schramm, D. N., Sigl, G., Turner, M. S., & Ellis, J. 1997, Phys. Rev. D, 56, 2419
  * Kolb & Turner (1994) Kolb, E. W., & Turner, M. S. 1994, The Early Universe, (Addison Wesley: Boston)
  * Kim (1979) Kim, J. E. 1979, Phys. Rev. Lett., 43, 103
  * Moriyama (1995) Moriyama, S. 1995, Phys. Rev. Lett., 75, 3222
  * Münch (1948) Münch,G. 1948, ApJ, 108, 116
  * Peccei & Quinn (1977a) Peccei, R. D., & Quinn, H. R. 1977a, Phys. Rev. Lett., 38, 1440
  * Peccei & Quinn (1977b) Peccei, R. D., & Quinn, H. R. 1977b, Phys. Rev. L, 38, 1440
  * Perna et al. (2012) Perna, R., Ho, W. C. G., Verde, L., van Adelsberg, M., & Jimenez, R. 2012, ApJ, 748, 116
  * Shifman et al. (1980) Shifman, M. A., Vainshtein, A. I., & Zakharov, V. I. 1980, Nucl. Phys. B, 166, 493
  * Weinberg (1978) Weinberg, S. 1978, Phys. Rev. Lett., 40, 223
  * Wilczek (1978) Wilczek, F. 1978, Phys. Rev. Lett., 40, 279
  * Wilczek (1982) Wilczek, F. 1982, Phys. Rev. Lett., 49, 1549
  * Yoshimura (1988) Yoshimura, M. 1988, Phys. Rev. D, 37, 2039
  * Zhitnitskii (1980) Zhitnitskii, A. R. 1980, Sov. J. Nucl. Phys., 31, 260