11institutetext: Department of Biomedicine and Prevention
University of Rome Tor Vergata (IT) 11email<EMAIL_ADDRESS>
22institutetext: Martinos Center For Biomedical Imaging
MGH and Harvard Medical School (USA)

# Through their eyes: multi-subject Brain Decoding with simple alignment
techniques

Matteo Ferrante 11 Tommaso Boccato 11 Nicola Toschi 1133

###### Abstract

To-date, brain decoding literature has focused on single-subject studies, i.e.
reconstructing stimuli presented to a subject under fMRI acquisition from the
fMRI activity of the same subject. The objective of this study is to introduce
a generalization technique that enables the decoding of a subject’s brain
based on fMRI activity of another subject, i.e. cross-subject brain decoding.
To this end, we also explore cross-subject data alignment techniques. Data
alignment is the attempt to register different subjects in a common anatomical
or functional space for further and more general analysis.

We utilized the Natural Scenes Dataset, a comprehensive 7T fMRI experiment
focused on vision of natural images. The dataset contains fMRI data from
multiple subjects exposed to 9841 images, where 982 images have been viewed by
all subjects. Our method involved training a decoding model on one subject’s
data, aligning new data from other subjects to this space, and testing the
decoding on the second subject based on information aligned to first subject.
We also compared different techniques for fMRI data alignment, specifically
ridge regression, hyper alignment, and anatomical alignment.

We found that cross-subject brain decoding is possible, even with a small
subset of the dataset, specifically, using the common data, which are around
$10\%$ of the total data, namely 982 images, with performances in decoding
compararble to the ones achieved by single subject decoding. Cross-subject
decoding is still feasible using half or a quarter of this number of images
with slightly lower performances. Ridge regression emerged as the best method
for functional alignment in fine-grained information decoding, outperforming
all other techniques.

By aligning multiple subjects, we achieved high-quality brain decoding and a
potential reduction in scan time by $90\%$. This substantial decrease in scan
time could open up unprecedented opportunities for more efficient experiment
execution and further advancements in the field, which commonly requires
prohibitive (20 hours) scan time per subject.

###### Keywords:

Brain Decoding Neuroscience Cross subject decoding

##  1 Introduction

Figure 1: Scheme for cross-subject decoding: The procedure involves the
following steps: In the first step, Subj01 (top row) is selected as the target
subject. A decoding model is trained to reconstruct seen images based on the
brain activity of Subj01 on the training set of Subj01 images (8859 images).
These images showed uniquely to this subject. Next, we decode the brain
activity of a second subject, Subj02, who was exposed to a share of the same
stimuli that Subj01 was exposed to (982 images). Using the shared images we
can compare the brain activity related to the same stimuli across different
subjects. We used this shared information to align the functional activity of
Subj02 with that of the target subject. Once the alignment transformation is
learned, we can align the complete dataset (including unseen data) and we can
utilize the pretrained decoder to reconstruct images from Subj02 without
training a decoding model specifically for Subj02.

Deep learning has revolutionized numerous fields, including neuroscience. The
application of deep learning techniques in neuroscience has led to significant
advancements in understanding brain function and decoding the intricate
workings of the human mind. Brain decoding, in particular, has emerged as a
crucial area where deep learning plays a pivotal role.

Brain decoding involves the extraction of meaningful information from recorded
brain activity, allowing researchers to infer mental states, perceptual
experiences, or cognitive processes. For example, deep learning algorithms
have been used to decode brain activity and predict whether an individual is
looking at a face or an object based on their neural responses [23, 2].

The potential applications of brain decoding are vast, from understanding
various aspects of brain function, such as information processing strategies,
decision-making, memory formation, and consolidation, to potential uses in
neurofeedback, neuroaesthetics, or neuromarketing strategies [10]. Moreover,
successful brain decoding could lead to novel strategies for diagnosing and
treating neurological or neuropsychiatric conditions, and potentially
contribute to the development of radically new algorithmic learning
strategies. However, these promising endeavors are not without challenges.
Noninvasive data, for instance, have lower temporal or spatial resolution than
neural firing, which may limit the granularity of information that can be
retrieved. Furthermore, physiological noise and signal/image artifacts can
affect both fMRI and EEG data, which can only be imperfectly removed after the
data are acquired. Nevertheless, several brain decoding studies have achieved
impressive results [23, 2].

One of the key challenges in brain decoding is the subject-specific nature of
all models developed thus far. This means that the models are tailored to
individual subjects, which can lead to significant variability in the results,
given that the amount of data collection per subject could be limited by
external factors like time and acquisition costs. Moreover, intrinsic inter-
individual variability poses further challenges and every model has to be
built from scratch for each new subject. This variability is a consequence of
the unique functional and anatomical structure of each individual’s brain, and
implies the need to acquire an entire dataset and train an individual model
for each subject. This technology’s use is limited by a bottleneck requiring
extensive data collection—typically thousands of stimulus images—to function
properly. This complexity stems from the unique brain anatomy, information
processing methods, and functional responses each individual possesses,
complicating the training of a universally applicable brain activity decoding
model. Despite individual differences in brain anatomy and function, common
structures enable reliable neuroscience analysis using template matching
techniques, such as anatomical and functional alignment. Anatomical alignment
transforms individual brain images to match a standard ’average’ brain
template or ’atlas’, aligning the size, shape, and orientation of brain
images. This facilitates meaningful cross-comparisons of brain images,
although it is more effective for larger, well-defined structures and may lack
precision for smaller, variable regions. Additionally, it does not account for
functional differences across brains. Thus, anatomical alignment is often
supplemented with functional alignment, which synchronizes brain activity
patterns across individuals, aiding the comparison and analysis of functional
data. This method is vital as activity locations can differ among individuals.
Numerous functional alignment methods exist, each with distinct applications
and limitations.

Figure 2: Example results: The first column, ”Stimulus”, presents the stimuli
from the fMRI experiment. The ”Subj01” column (in gray) displays the decoded
activity from Subj01, providing an upper performance baseline using a subject-
specific decoder model. All other columns (in red) show results from
functional alignment using Ridge Regression with $100\%$ common data (952
images), meaning subjects were functionally aligned to Subj01 and decoded
using Subj01’s trained decoder. To ensure robust visual comparisons, none of
the displayed images were used in learning alignment transformation,
demonstrating functional alignment on unseen data not used for decoder
training or alignment function learning.

In this study, we explore and contrast three methods for cross-subject brain
decoding of visual stimuli, using the established, cutting-edge Brain-Diffuser
decoding procedure. This approach mitigates variability from new decoding
procedures, enabling straightforward quantitative and qualitative comparisons
of cross-subject decoding outcomes.

We train a visual stimuli decoding model for one subject (Subj01) and employ
anatomical alignment, hyper alignment, and functional alignment with ridge
regression for three others, decoding their activity using the pretrained
model. Our goal is to demonstrate the feasibility of fine-grained cross-
subject decoding for visual stimuli reconstruction, potentially reducing scan
times significantly by only acquiring data necessary for alignment, thereby
reaching state-of-the-art in image reconstruction. Figure 1 outlines our
proposed pipeline, while Figure 2 provides examples of cross-subject decoding
results.

##  2 Related Work

Figure 3: The Brain-Diffuser pipeline, the state-of-the-art decoder for brain
activity used in this study, begins with brain activity from viewing an image
stimulus. A model is trained to estimate the latent representation of the
VDVAE autoencoder as well as the text and visual embeddings of the CLIP model,
using linear models. These estimated vectors and an initial guess image
obtained by decoding the autoencoder latents, are fed into Versatile
Diffusion—a latent diffusion model—to reconstruct the final image.

In the burgeoning field of deep learning-based brain decoding, a range of
models has been utilized to scrutinize preprocessed fMRI time series as input,
particularly focusing on visual stimuli decoding. This involves reconstructing
images that could have triggered specific fMRI patterns—termed brain activity.
Here, we review major works in this domain.

Some methods have employed variational autoencoders with a generative
adversarial component (VAE-GAN) to encode latent human face representations,
estimating these encoded representations from fMRI activity using a linear
model [22]. Sparse linear regression has also been utilized on preprocessed
fMRI data to predict features from the early convolutional layers of a pre-
trained convolutional neural network (CNN) for natural images [17].

Unsupervised and adversarial strategies have been used for image
reconstruction, including dual VAEGAN and unsupervised methods for decoding
fMRI stimuli, utilizing multiple encoder and decoder approaches [20, 19, 13].
Pretrained architectures like BigBiGAN and IC-GAN have optimized latent
spaces, significantly enhancing high-fidelity image reconstruction from fMRI
patterns [9, 6].

Recently, diffusion models have become prominent in the decoding pipeline,
providing superior image generation performance [21, 8]. These models often
incorporate semantic-based strategies [11] and multi-step decoding strategies
[18, 12].

Regarding alignment techniques, there are several approaches [3].
Hyperalignment [15, 16] aligns functional brain activity across individuals in
a high-dimensional space, enhancing the precision of cross-subject brain
activity predictions, but requires extensive high-quality data and complex
computational resources. The Shared Response Model (SRM) [7] aligns brain
activity by identifying a common response pattern across subjects, ideal for
shared experiences, but assumes uniform responses, which individual perception
and cognition differences may contradict. Independent Component Analysis (ICA)
[5], separates multivariate signals into additive subcomponents, identifying
common brain activity patterns, but requires statistical independence of
subcomponents, which may not always apply to brain data.

While functional alignment methods provide powerful tools for comparing brain
activity, their limitations and assumptions require careful result
interpretation. These methods align and compare functional brain data,
complementing, not replacing, anatomical alignment. Various other methods,
each with its pros and cons, have been proposed.

In this paper, we compare anatomical alignment, hyperalignment-based
functional alignment, and ridge regression-based alignment methods for cross-
subject brain decoding.

##  3 Material and Methods

Figure 4: More example results. Format and conventions as in Figure 2

In this section, we describe the proposed method and the data we used. The
data are publicly available and can be requested at
https://naturalscenesdataset.org/. All experiments and models were trained on
a server equipped with four NVIDIA A100 GPU cards (80GB RAM each connected
through NVLINK) and 2 TB of System RAM.

The study utilizes the Natural Scenes Dataset (NSD) [1], a vast fMRI data set
from eight subjects exposed to images from the COCO21 dataset. We focused on
four subjects, forming a unique training dataset of 8,859 images and 24,980
fMRI trials, and a common dataset of 982 images and 2,770 trials. To reduce
spatial dimensionality, we applied a mask to the fMRI signal (resolution of
1.8mm isotropic) using the NSDGeneral ROI, targeting various visual areas.
This strategic ROI selection enhanced the signal-to-noise ratio and simplified
data complexity, enabling exploration of both low-level and high-level visual
features. Temporal dimensionality was reduced using precomputed betas from a
general linear model (GLM) with a fitted hemodynamic response function (HRF)
and a denoising process as detailed in the NSD paper. Data from Subj01,
Subj02, Subj05, Subj07, warped into the Montreal Neurological Institute common
space (MNI) and downsampled at 2mm, represented the brain activity of each
subject and helped decrease computational time and cost. We used the common
dataset as alignment, keeping out 30 images for visual comparison, so there
are 8859 unique images for each subject. We only used them for training the
decoding model for Subj01. Then there are 952 common images across all
subjects that was used to functionally align them to the activity of Subj01,
and 30 common images kept out for visual comparison on images neither used in
the training or in the alignment procedure. These 30 images were chosen
because they’re used as visual qualitative evaluation of decoding results in
other papers and could help the reader to compare results across different
methods. Decoding metrics are evaluated on the 952 images which correspond to
our test set for each one of the subjects, since these images are never seen
by the decoder model, so the evaluation is still fair and on unseen images.
When we refer to $100\%$ of common data we are pointing to these 952 images.

###  3.1 Decoding Pipeline: Brain-Diffuser

The ”Brain-Diffuser” [18] model is a two-stage framework for reconstructing
natural scenes from fMRI signals. Initially, a Very Deep Variational
Autoencoder (VDVAE) provides an ”initial guess” of the reconstruction,
focusing on low-level details. This guess is refined using high-level semantic
features from CLIP-Text and CLIP-Vision models, employing a latent diffusion
model (Versatile Diffusion) for final image generation. The model, represented
in Fig. 3, takes fMRI signals as input and generates reconstructed images,
capturing low-level properties and overall layout. As a state-of-the-art
procedure, Brain-Diffuser was trained using data from Subj01 in the MNI space
(cross-subject decoding requires a of a common space). Further details about
the decoding model are available in the original paper.

Figure 5: More example results. Format and conventions as in Figure 2

###  3.2 Alignment

This study investigates three alignment strategies to evaluate cross-subject
fine-grained brain decoding’s feasibility: anatomical alignment, functional
alignment via hyper alignment, and functional alignment through ridge
regression.

Anatomical alignment, our baseline, relies solely on anatomical details,
transforming functional aspects using pre-computed structural image
transformations. On the other hand, functional alignment necessitates a more
comprehensive approach. Consider the scenario where the brain activity of a
source subject $\mathbf{S}$ needs to align with a target subject $\mathbf{T}$.
These activities, responses to numerous stimuli, are matrices of shape ($\\#$
stimuli, $\\#$ voxels). Given that subjects encounter several common stimuli
(i.e., they view identical images in the fMRI scanner), we can divide the
datasets into $\mathbf{T}{common},\mathbf{T}{different}$ and
$\mathbf{S}{common},\mathbf{S}{different}$. Our goal is to leverage the
$common$ dataset portion to learn a mapping from $S$ to $T$, aligning the
entire $\mathbf{S}$ dataset with the $\mathbf{T}$ functional space. The NSD
experiment’s structure, with separate training and test sets (the latter
containing identical images for each subject), provides a common stimuli set
for alignment purposes.

####  3.2.1 Anatomical Alignment

Anatomical alignment, a common neuroscience method, aligns to a standard
template, here, the MNI space, facilitating anatomical structure comparison.
This alignment typically involves a linear coregistration of anatomical images
between native and common spaces, followed by a nonlinear warping to match
common brain structures. Several software options like FSL and ANTs can
perform this task. The NSDData authors [1] elaborate on this process in their
released code, providing betas (i.e. coefficients obtained by theressing the
stimuls waveform against the fMRI data) for all subjects in the MNI common
space (1mm). We downsampled these to 2mm to approximate the resolution used in
the original Brain-Diffuser decoding paper (1.8mm) and to reduce spatial
dimensionality.

####  3.2.2 HyperAlignment

HyperAlignment [15, 16], a functional data alignment technique, models
functional data as high-dimensional points, with each voxel representing a
dimension with betas ranging in $\mathcal{R}$. This method, based on
Procrustes Analysis [14], presents a high-dimensional model of the
representational space in the human ventral temporal (VT) cortex, wherein
dimensions are response-tuning functions common across individuals.

To perform the Procrustes analysis for functional brain alignment, we aim to
find a rotation matrix $\mathbf{R}$ and a scale factor $c$ such that the
difference $|c\mathbf{S}\mathbf{R}-\mathbf{T}|^{2}$ is minimized.

This is achieved by computing the matrix product
$\mathbf{P}=\mathbf{S}^{T}_{common}\mathbf{T}_{common}$, Performing the
singular value decomposition of $\mathbf{P}$ to obtain left and right
eigenvector matrices $\mathbf{U}$ and $\mathbf{V}$, Computing
$\mathbf{R}=\mathbf{U}\mathbf{V}^{T}$ and the scaling factor
$c=\frac{trace(\mathbf{T_{common}}^{T}(\mathbf{S_{common}}\mathbf{R}))}{trace(\mathbf{S}^{T}_{common}\mathbf{S}_{common})}$.
Finally, we can apply the matrix $\mathbf{R}$ and the scaling $c$ to both
common and non-common source data to align them with the target subject. We
computed these values for Subj02, Subj05, and Subj07 as source subjects, using
Subj01 as the target, to align all subjects to the functional space of the
first one. For detailed mathematical proofs and other insights, please refer
to the original articles [16, 15, 14].

Figure 6: More example results. Format and conventions as in Figure 2

####  3.2.3 Ridge Regression

Our third approach embraces a simple assumption: even in different subjects,
all functional data contain the information for the same stimuli, albeit
possibly spread across different voxels. This suggests that one subject’s
activity (source) might be expressed as a linear combination of the activity
of another subject (target) for the same stimuli. By deriving a linear
combination for each voxel of the target from all possible voxels of the
source, we can create a linear map from the source to the target, facilitating
functional alignment. The target subject activity can be expressed as
$\mathbf{t}i=\sum{j}\mathbf{w}{j}\mathbf{s}{j}$, where $\mathbf{t}i$ is the
$i$-th activity of the target voxel for each common dataset stimulus. Here,
$\mathbf{t}i$ represents the $i$-th column of $\mathbf{T}{common}$, expressed
as a linear combination of all $\mathbf{S}{common}$ columns. The challenge
lies in finding the vector of $\mathbf{w}$ values. When extended to all the
target subject voxels, the $w$ vector morphs into a square matrix
$\mathbf{W}$, each column of which contains weights to estimate one target
subject voxel from a combination of source values. The objective can be
redefined as minimizing
$|\mathbf{S}{common}\mathbf{W}^{T}-\mathbf{T}{common}|^{2}$.

We employed Ridge Regression from [4] to determine the W matrix, conducting a
5-fold cross-validation to select the optimal hyper-parameter $\alpha$ from
the values $[0,1,10,1e2,1e3,1e4]$. Our findings indicated that $\alpha=1000$
yields superior performance, hence we adopted it as our final regularization
parameter. We computed these values to align all subjects to the initial
functional space For the sources Subj02, Subj05, and Subj07, and using Subj01
as the target.

###  3.3 Evaluation

Our research seeks to evaluate visual stimuli’s detailed brain decoding
feasibility, scrutinizing the alignment methods and shared data ratio at play.
We examined how the alignment performance fluctuates when the shared data
makes up $10\%$, $25\%$, $50\%$, and $100\%$ of the total common data (952
images).

Our shared dataset, or ”test dataset,” comprises 982 images, all viewed by
every subject. In order to allow visual comparison, we excluded 30 images from
the original Brain-Diffuser paper. Thus, these excluded images neither
influenced the training of the decoding pipeline nor the alignment process.
The remaining 952 images serve as the shared dataset. We computed
transformations for each alignment method (anatomical, hyperalignment, ridge
regression) and shared dataset proportion, applying the linear transformation
to the complete dataset. This procedure aligns the images with Subject 01’s
functional space. We then used the pre-trained Brain-Diffuser pipeline for
decoding the aligned fMRI activity and reconstructing the images. We assessed
our image reconstruction process through both basic and advanced metrics,
including PixCorr, SSIM, and 2-way accuracy in AlexNet, Inception, and CLIP
latent spaces. This comprehensive evaluation approach allows us to benchmark
our results against other brain decoding studies. However, the goal here is
not merely comparison, but rather the examination of performance in relation
to the shared data fraction and alignment method, given a fixed decoding
pipeline, trained solely on Subj01 as a reference target.

##  4 Results

Figure 7: A: This figure compares Ridge Regression functional alignment over
different fractions of common data. The ”Stimulus” column presents the images
used during the experiment, with the remaining columns illustrating the
decoded, aligned activity of Subj02 on Subj01. B: Here, we compare various
alignment modalities. Again, the ”Stimulus” column displays the images used
during the experiment, while the other columns visualize the decoded activity
of Subj02 aligned to Subj01 using different methods.

Figures 4, 5, 2, 6 provide a comparison between stimuli and decoded images
from Subj01 (on which the decoding model is trained). These figures also
display the aligned activity of all other subjects using Ridge Regression.
Figure 7 compares fractions of common data used for Ridge Regression-based
alignment and different alignment methods. Lastly, Figure 8 illustrates each
quantitative metric, computed and averaged over the entire test set for each
aligned subject (2,5,7). Metrics are expressed as a fraction of the entire
dataset, which contains approximately 10k images per subject (8859 unique +
982 common across subjects). As common images are necessary, the maximum
amount of images that can be included in the alignment process is 952 (termed
the ”test set”, except 30 images left out for visualization purposes),
representing around 10% of the dataset. This represents the maximum data that
can be incorporated into the procedure, and we experimented with half, a
quarter, and a tenth of this data.

Anatomical and Hyperalignment methods fail to yield satisfactory results,
demonstrating just above chance performance levels for 2-way classification
accuracy and poor performance for low-level metrics such as SSIM and PixCorr.
However, Ridge Regression exhibits an increasing performance based on the
volume of data used for alignment mapping function learning. This method
reaches performance levels comparable with the within-subject decoder in both
low-level and high-level metrics, using all the common data (approximately 10%
of the entire dataset).

Our findings are encapsulated in the following key points:

Functional alignment’s critical role in fine-grained brain decoding: Our
research emphasizes the pivotal role of functional alignment in fine-grained
brain decoding. This process, which interprets neural signals to reconstruct
perceived images or thoughts, greatly benefits from precise functional
alignment of brain activity. Accurate alignment ensures that neural signals
are matched correctly to their corresponding brain regions, thus enhancing the
decoding accuracy.

Anatomical method’s inefficacy: As corroborated by previous studies [16], our
research found that anatomical methods for brain decoding are ineffective.
Relying on the physical structure of the brain for alignment and decoding does
not deliver the requisite precision for fine-grained decoding tasks. This
could be attributed to inherent anatomical variability across different
individuals, which may not necessarily align with functional differences. The
specialized areas in the brain with functional selectivity can sometimes yield
performance above chance levels. However, in most cases, decoded images do not
correlate with the stimulus, undermining the reliability of this method for
cross-subject brain decoding.

Overfitting tendency of complex techniques: We noted that more sophisticated
techniques, like hyperscanning for brain decoding, tend to overfit the data.
This results in poor generalization to unseen data, with metrics measuring
n-way accuracy reaching only chance levels. While these techniques might seem
to offer superior decoding accuracy initially, their lack of generalizability
limits their practical utility. Of course, room for improvement exists,
perhaps through the incorporation of regularization techniques.

Ridge Regression’s efficacy: Our results demonstrate that Ridge Regression-
based methods for brain decoding can achieve above-chance performance levels
with as little as 1% of the entire dataset. Furthermore, these methods near
baseline performance levels with around 10% of the dataset. This crucial
finding implies that reliable brain decoding results can be achieved while
significantly reducing scan time. This efficiency could be instrumental in
making brain decoding research and applications more feasible and cost-
effective.

These results contribute to our comprehension of the challenges and potential
remedies in brain decoding and emphasize the need for additional research to
refine these techniques and augment their effectiveness.

Figure 8: This figure illustrates the performance of various methods evaluated
using different metrics. Blue lines represent metrics from the target
subject’s decoded images, derived from their test set brain activity. Green
lines denote the mean and standard deviation (std) of performance on test sets
from other subjects, aligned using hyperalignment. Gray lines present results
achieved using anatomical alignment, while orange lines display outcomes using
Ridge Regression. Remarkably, the Ridge Regression approach yields positive
results even when using a tiny fraction of the entire dataset. Furthermore, as
this fraction approaches roughly $10\%$ of the total set, resulting in 952
images the performance becomes comparable with those obtained by the within-
subject model.

##  5 Discussion

Our study underscores the intricacies and potential of cross-subject fine-
grained brain decoding, a field promising to enhance our understanding of the
human brain and cognition.

We identified the criticality of functional alignment for successfully
executing brain decoding. This alignment, which maps neural signals to their
corresponding brain regions, is vital for accurately decoding neural activity
from other individuals using a model pre-trained on separate subject data.
This insight holds promise for constructing large studies with a high-accuracy
decoding pipeline, subsequently requiring only alignment data acquisition for
new subjects. This approach negates the need for an entire experimental
reproduction each time, streamlining the process.

Our research reveals the limitations of anatomical methods for brain decoding,
which rely on the physical brain structure for alignment and decoding. These
methods underperformed due to inherent brain anatomical variability across
individuals, which may not align with functional differences. Thus, this study
emphasizes the need for functional, not merely anatomical, considerations in
decoding studies.

Excitingly, our results suggest significant reductions in scan-time are
possible. Ridge Regression-based methods were found to provide reliable brain
decoding results with just a fraction of the entire dataset, implying
practical implications for brain decoding research feasibility and cost-
effectiveness.

Our study also highlights the qualitative similarities in decoded images
across subjects. While these images largely match high-level semantic content,
intra-subject differences appear minimized. This observation prompts us to
consider whether the decoding procedure is fully captured. Given that high-
level concepts are generally aligned, we propose a possible brain activity
decomposition into brain activity $=$ concept \+ individual perception. Such a
model might only capture the concept while treating differences as noise,
offering new research directions to explore fine-grained inter-subject
differences.

Apart from speculation, our findings suggest that despite individual brain
structure and function differences, common neural activity patterns can be
decoded across individuals. This opens up intriguing opportunities for
developing generalized brain decoding models applicable across different
individuals. However, our study also reveals the limitations of current
functional alignment methods. Future research might consider exploring more
complex models, such as neural networks, known for their ability to capture
intricate, non-linear relationships, potentially improving functional
alignment.

Looking ahead, several promising future work directions emerge, such as
training models across subjects and machines, which could lead to more robust
and generalizable brain decoding models. The development of new techniques and
methodologies could potentially address current limitations of brain decoding,
heralding more accurate and efficient decoding of brain activity. As we delve
into fine-grained brain decoding, addressing potential privacy concerns and
ethical implications is paramount. Current research suggests that while
certain brain activity aspects can be decoded across subjects, the process is
not yet a comprehensive or intrusive ’mind-reading’ tool. A key finding
highlights the disruptive role of attention mechanisms, suggesting that brain
decoding is only possible with actively participating, aware subjects.

While our methods currently prevent involuntary or covert ’mind reading’, as
the field advances, maintaining strong ethical frameworks for brain decoding
research becomes even more critical. Informed consent, strict data privacy
protocols, and potential societal implications consideration remain key.
Decoding brain activity raises broader ethical questions, such as its
potential use to enhance communication for individuals with speech or motor
impairments or its potential misuse for coercive or manipulative purposes.
These critical questions must be confronted by the scientific community and
society as we continue to explore brain decoding potential.

##  6 Conclusions

In this study, we present a method for conducting brain decoding of visual
stimuli across different subjects. We delved into the role of functional
alignment in decoding neural signals with accuracy, uncovering the limitations
inherent in anatomical methods and the pitfalls of complex decoding techniques
that are susceptible to overfitting.

Our study saw a significant breakthrough with the successful implementation of
Ridge Regression-based methods. This technique showed remarkable efficiency in
decoding neural activity, utilizing only a fraction of the dataset. This
discovery indicates the possibility of substantial reductions in scan time -
with potential reductions nearing 90$\%$. Such an advancement could
revolutionize the feasibility and cost-effectiveness of brain decoding
research and its applications.

In our research, we achieved successful cross-subject brain decoding by
training the decoding pipeline on one subject and effectively decoding neural
activity across multiple individuals. This significant finding points to the
existence of shared neural activity patterns, paving the way for the
development of generalized brain decoding models. We also revealed a
hierarchical structure in the brain’s processing and representation of
information, separating brain decoding into concept and perception components.
Despite these advancements, we recognize the limitations of current functional
alignment methods and advocate for exploring future research directions, such
as training models across subjects and machines.

In conclusion, our study illuminates the challenges and potential solutions
within the sphere of fine-grained cross-subject brain decoding, particularly
in relation to visual stimuli.

## References

  * [1] Allen, E.J., St-Yves, G., Wu, Y., Breedlove, J.L., Prince, J.S., Dowdle, L.T., Nau, M., Caron, B., Pestilli, F., Charest, I., Hutchinson, J.B., Naselaris, T., Kay, K.: A massive 7t fmri dataset to bridge cognitive neuroscience and artificial intelligence. Nature Neuroscience 25(1), 116–126 (Jan 2022). https://doi.org/10.1038/s41593-021-00962-x, https://doi.org/10.1038/s41593-021-00962-x
  * [2] Awangga, R.M., Mengko, T.L.R., Utama, N.P.: A literature review of brain decoding research. IOP Conference Series: Materials Science and Engineering 830(3), 032049 (Apr 2020). https://doi.org/10.1088/1757-899X/830/3/032049, https://iopscience.iop.org/article/10.1088/1757-899X/830/3/032049
  * [3] Bazeille, T., DuPre, E., Richard, H., Poline, J.B., Thirion, B.: An empirical evaluation of functional alignment using inter-subject decoding. NeuroImage 245, 118683 (2021). https://doi.org/https://doi.org/10.1016/j.neuroimage.2021.118683, https://www.sciencedirect.com/science/article/pii/S1053811921009563
  * [4] Buitinck, L., Louppe, G., Blondel, M., Pedregosa, F., Mueller, A., Grisel, O., Niculae, V., Prettenhofer, P., Gramfort, A., Grobler, J., Layton, R., VanderPlas, J., Joly, A., Holt, B., Varoquaux, G.: API design for machine learning software: experiences from the scikit-learn project. In: ECML PKDD Workshop: Languages for Data Mining and Machine Learning. pp. 108–122 (2013)
  * [5] Calhoun, V.D., Liu, J., Adalı, T.: A review of group ica for fmri data and ica for joint inference of imaging, genetic, and erp data. NeuroImage 45(1, Supplement 1), S163–S172 (2009). https://doi.org/https://doi.org/10.1016/j.neuroimage.2008.10.057, https://www.sciencedirect.com/science/article/pii/S1053811908012032, mathematics in Brain Imaging
  * [6] Casanova, A., Careil, M., Verbeek, J., Drozdzal, M., Romero-Soriano, A.: Instance-conditioned gan (2021). https://doi.org/10.48550/ARXIV.2109.05070, https://arxiv.org/abs/2109.05070
  * [7] Chen, P.H.C., Chen, J., Yeshurun, Y., Hasson, U., Haxby, J., Ramadge, P.J.: A reduced-dimension fmri shared response model. In: Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 28. Curran Associates, Inc. (2015)
  * [8] Chen, Z., Qing, J., Xiang, T., Yue, W.L., Zhou, J.H.: Seeing beyond the brain: Conditional diffusion model with sparse masked modeling for vision decoding (2022)
  * [9] Donahue, J., Simonyan, K.: Large scale adversarial representation learning (2019). https://doi.org/10.48550/ARXIV.1907.02544, https://arxiv.org/abs/1907.02544
  * [10] Du, B., Cheng, X., Duan, Y., Ning, H.: fmri brain decoding and its applications in brain and computer interface: A survey. Brain Sciences 12(2) (2022). https://doi.org/10.3390/brainsci12020228, https://www.mdpi.com/2076-3425/12/2/228
  * [11] Ferrante, M., Boccato, T., Toschi, N.: Semantic brain decoding: from fmri to conceptually similar image reconstruction of visual stimuli (2023)
  * [12] Ferrante, M., Ozcelik, F., Boccato, T., VanRullen, R., Toschi, N.: Brain captioning: Decoding human brain activity into images and text (2023)
  * [13] Gaziv, G., Beliy, R., Granot, N., Hoogi, A., Strappini, F., Golan, T., Irani, M.: Self-supervised Natural Image Reconstruction and Large-scale Semantic Classification from Brain Activity. NeuroImage 254, 119121 (Jul 2022). https://doi.org/10.1016/j.neuroimage.2022.119121, https://linkinghub.elsevier.com/retrieve/pii/S105381192200249X
  * [14] Gower, J.C.: Generalized procrustes analysis. Psychometrika 40(1), 33–51 (Mar 1975). https://doi.org/10.1007/BF02291478, https://doi.org/10.1007/BF02291478
  * [15] Haxby, J.V., Guntupalli, J.S., Connolly, A.C., Halchenko, Y.O., Conroy, B.R., Gobbini, M.I., Hanke, M., Ramadge, P.J.: A common, high-dimensional model of the representational space in human ventral temporal cortex. Neuron 72(2), 404–416 (Oct 2011)
  * [16] Haxby, J.V., Guntupalli, J.S., Nastase, S.A., Feilong, M.: Hyperalignment: Modeling shared information encoded in idiosyncratic cortical topographies. eLife 9, e56601 (jun 2020). https://doi.org/10.7554/eLife.56601, https://doi.org/10.7554/eLife.56601
  * [17] Horikawa, T., Kamitani, Y.: Generic decoding of seen and imagined objects using hierarchical visual features. Nature Communications 8(1), 15037 (Aug 2017). https://doi.org/10.1038/ncomms15037, http://www.nature.com/articles/ncomms15037
  * [18] Ozcelik, F., VanRullen, R.: Brain-diffuser: Natural scene reconstruction from fmri signals using generative latent diffusion (2023)
  * [19] Ren, Z., Li, J., Xue, X., Li, X., Yang, F., Jiao, Z., Gao, X.: Reconstructing Perceived Images from Brain Activity by Visually-guided Cognitive Representation and Adversarial Learning (Oct 2019), http://arxiv.org/abs/1906.12181, arXiv:1906.12181 [cs]
  * [20] Shen, G., Dwivedi, K., Majima, K., Horikawa, T., Kamitani, Y.: End-to-end deep image reconstruction from human brain activity. Front. Comput. Neurosci. 13, 21 (Apr 2019)
  * [21] Takagi, Y., Nishimoto, S.: High-resolution image reconstruction with latent diffusion models from human brain activity. bioRxiv (2023). https://doi.org/10.1101/2022.11.18.517004, https://www.biorxiv.org/content/early/2023/03/11/2022.11.18.517004
  * [22] VanRullen, R., Reddy, L.: Reconstructing faces from fmri patterns using deep generative neural networks. Communications Biology 2(1), 193 (May 2019). https://doi.org/10.1038/s42003-019-0438-y, https://doi.org/10.1038/s42003-019-0438-y
  * [23] Zafar, R., Malik, A.S., Kamel, N., Dass, S.C., Abdullah, J.M., Reza, F., Abdul Karim, A.H.: Decoding of visual information from human brain activity: A review of fMRI and EEG studies. Journal of Integrative Neuroscience 14(02), 155–168 (Jun 2015). https://doi.org/10.1142/S0219635215500089, http://www.worldscientific.com/doi/abs/10.1142/S0219635215500089


# A Quarter Century of Guitar Nebula/Filament Evolution

Martijn de Vries  Department of Physics/KIPAC, Stanford University, Stanford,
CA 94305-4060, USA  Roger W. Romani  Department of Physics/KIPAC, Stanford
University, Stanford, CA 94305-4060, USA  Oleg Kargaltsev  Department of
Physics, George Washington University, Washington, DC 20052, USA  George
Pavlov  Department of Astronomy & Astrophysics, Pennsylvania State University,
University Park, PA 16802, USA  Bettina Posselt  Department of Astronomy &
Astrophysics, Pennsylvania State University, University Park, PA 16802, USA
Department of Physics, University of Oxford, OX1 3PU Oxford, UK  Patrick Slane
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA  Niccolo’ Bucciantini  INAF - Osservatorio Astrofisico di Arcetri,
L.go Fermi 5, 50125, Firenze, Italy  C.-Y. Ng  Department of Physics, The
University of Hong Kong, Pokfulam Road, Hong Kong  Noel Klingler  Department
of Astronomy and Astrophysics, The Pennsylvania State University, 525 Davey
Laboratory, University Park, PA 16802, USA

###### Abstract

We have collected a new deep Chandra X-ray Observatory (CXO) exposure of PSR
B2224+65 and the ‘Guitar Nebula’, mapping the complex X-ray structure. This is
accompanied by a new HST H$\alpha$ image of the head of the Guitar. Comparing
the HST and CXO structures in 4 epochs over 25 years, we constrain the
evolution of the TeV particles that light up the filament. Cross-field
diffusion appears to be enhanced, likely by the injected particles, behind the
filament’s sharp leading edge, explaining the filament width and its evolving
surface brightness profile.

stars: neutron — pulsars: individual (PSR B2224+65)

††facilities: HST(ACS/WFC), Gemini(GMOS-N), CXO

##  1 Introduction

A handful of fast-moving pulsars have been seen to have narrow X-ray
structures (‘filaments’) extending from the pulsar point source at large angle
to the proper motion axis, itself often marked by a Pulsar Wind Nebula (PWN)
trail. The first discovered, and arguably most spectacular, is that associated
with PSR B2224+65. This pulsar has a very large proper motion $\mu=194.1\pm
0.2$ mas/y, which at its VLBI-measured $0.83^{+0.17}_{-0.10}$ kpc parallax
distance (Deller et al., 2019) implies a highly supersonic $v_{\rm PSR}=765$
km/s. This leads to the formation of its remarkable H$\alpha$ bow shock, ‘The
Guitar Nebula’ (Cordes et al., 1993), which has a very small $\sim
0.1^{\prime\prime}$ standoff angle, and a long trail forming the neck and body
of the Guitar. Projecting from near the pulsar at $\sim 115^{\circ}$ to its
proper motion is an X-ray filament, with a variable width of $\sim
20^{\prime\prime}$ and length $\sim 2.5^{\prime}$. It has a sharp leading
edge, in the direction of the pulsar motion with surface brightness trailing
off behind.

In the picture sketched by Bandiera (2008) and explored numerically by Barkov
et al. (2019) and Olmi & Bucciantini (2019), pulsar filaments are created by
multi-TeV pulsar $e^{\pm}$ leaking out near the bow shock apex and escaping to
external ISM field lines. Electrons whose gyroradius $r_{c}$ approaches or
exceeds the standoff distance $r_{0}=[{\dot{E}}/(4\pi\mu m_{p}n_{0}cv_{\rm
PSR}^{2})]^{1/2}$ may escape to the filaments. Small $r_{0}$ requires some
combination of low pulsar $\dot{E}$, high velocity $v_{\rm PSR}$, and high
ambient ISM density $n_{0}$. The presence of an H$\alpha$ bow shock in B2224
(implying high ISM density), as well as its large transverse velocity seem to
support this picture.

In de Vries & Romani (2022), we argued that the morphology of the filament
associated with PSR J2030+4415 can be explained by a variable particle
injection rate: the H$\alpha$ morphology suggests a temporary decrease in the
standoff distance, subsequently enabling enhanced particle injection into the
ISM over a period of approximately a decade. In order to connect the
properties of the X-ray filament to the level of particle injection, multi-
epoch H$\alpha$ bow shock observations are crucial as they allow for direct
measurement of $r_{0}$ and provide important clues to fluctuations in
$\dot{E}$ or $n_{0}$ at earlier times. The Guitar, which has been observed
several times by both Chandra and HST over a period of more than 25 years,
thus provides us with a unique opportunity to witness the time evolution of a
pulsar filament.

The filament of B2224 has been the subject of three previous Chandra ACIS
exposures in 2000, 2006 and 2012. Since there were significant changes between
these epochs, we have collected a new, deeper ACIS exposure to provide a
fiducial structure and context for the earlier observations. The evolution and
X-ray spectral parameters give information about the multi-TeV $e^{\pm}$
injection which lights up the filament. Concurrent with the new epoch of X-ray
observations, we have obtained a new HST ACS/WFC H$\alpha$ image of the apex
of the nebula.

Figure 1: Overview of the Guitar/filament complex in the most recent HST and
Chandra epochs. Red: 2020 HST ACS/WFC H$\alpha$, green: PanSTARRS2 $r$, blue:
2021 CXO 1-5 keV X-rays. The pulsar point source lies at the tip of the
H$\alpha$ nebula. Filament X-ray emission extends primarily to the right (NW).
This has a sharp leading edge and extension behind. The Guitar ‘body’ shows
faintly in the H$\alpha$ in the lower half of the image.
Figure 2: GMOS-N IFU H$\alpha$ velocity channel images, covering
$5.5^{\prime\prime}\times 6.1^{\prime\prime}$. The combined IFU image from
2016 is shown in the middle panel. This is flanked by velocity slices (central
velocity, in km/s, at upper right in each frame), with the H$\alpha$ limb as a
red outline, for comparison. The oval region near the apex provides the line
spectrum shown in Figure 3.

##  2 Observations and Data reduction

###  2.1 Chandra Observations

Table 1: Overview of B2224 CXO observations used in this paper. The ‘Aim’ column indicates aimpoint on the ACIS-I or ACIS-S array. ‘Exp’ lists the exposure times in kiloseconds, after filtering out periods of high background. Date | Obs | Aim | Exp | Date | Obs | Aim | Exp
---|---|---|---|---|---|---|---
|  |  | [ks] |  |  |  | [ks]
2000-10-21 | 755 | S | 48.8 | 2021-04-21 | 24433 | I | 25.7
2006-08-28 | 6691 | S | 10.0 | 2021-04-23 | 24431 | I | 25.7
2006-10-06 | 7400 | S | 36.6 | 2021-04-25 | 24429 | I | 24.5
2012-07-28 | 14467 | S | 14.6 | 2021-07-04 | 24428 | I | 29.7
2012-07-29 | 14353 | S | 34.6 | 2021-07-25 | 24427 | I | 24.7
2012-08-01 | 13771 | S | 49.2 | 2021-10-09 | 24430 | I | 29.7
2021-02-19 | 24437 | I | 24.7 | 2021-10-20 | 23537 | I | 57.2
2021-03-14 | 24434 | I | 29.5 | 2021-11-14 | 24432 | I | 29.6
2021-03-15 | 24435 | I | 14.6 | 2022-02-21 | 24426 | I | 20.9
2021-03-16 | 24992 | I | 14.9 | 2022-02-24 | 26336 | I | 18.0
2021-04-03 | 24436 | I | 24.3 |  |  |  |

The Guitar, at its high Northern declination, requires short dwell times for
the thermal health of CXO. Thus we obtained fifteen $15-60\,$ks visits between
2021 February 19 and 2022 February 24 to collect 393.8 ks of exposure. For all
observations we used the I3 chip of ACIS-I array, with the aimpoint positioned
so that the full filament was covered at any roll angle. An overview of the
new and archival observations is given in Table 1. In addition to the
2021-2022 epoch of 393.8 ks, we re-analyze the three archival epochs: 2000.89
(48.8ks total, Ep1), 2006.80 (46.5 ks total, Ep2), and 2012.66 (98.4 ks total,
Ep3).

All data were reprocessed with the standard CIAO reprocessing tools, using
CIAO 4.12 and CALDB 4.9.1. We performed a relative astrometric correction by
using wavdetect on the aimpoint chip of each observation (S3 for the archival
epochs, I3 for the latest epoch) and then using reproject_aspect to minimize
the point source offsets between observations. ObsID 23537 was used as the
reference observation, because it has the longest exposure time in the new
epoch. The pulsar was excluded from the point source list, because its large
proper motion makes it unsuitable for astrometry. For most ObsIDs, at least 4
sources could be used for registration. However the degraded soft X-ray
response of ACIS and the short exposure time in the later observations (most
notably 24435, 24436, 24992, and 26336) combined with the lack of large
numbers of bright field stars means that for these obervations only 2-3
matching sources could be found. The bright source directly NW of the PSR
unfortunately decreased in brightness over the course of 2021, making this
source difficult to use for astrometry in the later observations of the 2021
epoch as well. We estimated the error on the frame registration from the RMS
residuals after source matching, and summed the errors of individual frames
weighted by exposure time to find the error on the relative astrometric
correction in each epoch. We find $1\sigma$ errors of 0.08, 0.12, 0.11, and
0.16 arcsec for the 2000, 2006, 2012, and 2021-2022 epochs respectively. After
the astrometric correction, we combined the event files for each epoch and
generated $1-5$ keV exposure maps using a power law with $\Gamma=1.5$
(appropriate to the PWN emission, see §3.3) as an input energy weighting, with
the tool merge_obs.

In our X-ray analysis, we compare the exposure-corrected data from the S and I
chips, which have significantly different particle background levels. We
therefore estimated the particle background contribution by retrieving the
ACIS ‘stowed’ background maps from the calibration database, and scaling them
by the number of 9.5-12 keV counts of the observation in question (no bright
sources were present to contribute significant photon counts to this band).
The scaled backgrounds were subtracted from the data before making the
exposure-corrected image for each epoch.

###  2.2 New HST Apex Image

To probe the current state of the bow shock we obtained a new HST ACS
H$\alpha$ image of the nebula apex (the ‘head’ of the Guitar) under Program
16426. This is a challenging observation, since, although the head is the
brightest portion of the nebula, its surface brightness is still low. And
while the ACS/WFC has the highest H$\alpha$ sensitivity of the present HST
cameras, at low light levels it suffers severely from a degraded Charge
Transfer Efficiency (CTE). Happily at $DEC=+65^{\circ}$, the Guitar lies far
enough North to be occasionally in the Continuous Viewing Zone and CVZ
observation were kindly granted for the three awarded orbits by HST, allowing
much longer exposures, and greater photoelectron count per pixel at readout,
than would otherwise be available. In the end we were able to schedule five
2910 s F658N exposures (and two 338 s F625W frames for continuum monitor and
subtraction) on 2020 November 8 (MJD 59161). With few exposures we also reduce
the total read-noise cost of the observation but slightly decrease the CR
rejection efficiency. The second mitigation is to place the Guitar apex near
the WFC1-CTE position, where the number of row transfers is minimal,
decreasing CTE degradation. The cost is that WFC field of view is cropped
closely around the apex, decreasing the number of upstream field star
detections for registration and context. In addition, the geometrical
distortions are largest near the array corners, such as the CTE1 position.
Happily, at the observation epoch the default roll angle placed the body of
the Guitar farther on to the WFC array. Its surface brightness is very low so
that the only useful measurements of the body at HST resolution are the limb
position in some of the brighter areas. Nevertheless it is gratifying to
detect the whole structure (Figure 1).

In the end these mitigation measures were successful and we have obtained the
best-ever image of the Guitar head. Since the pulsar had moved
$2.74^{\prime\prime}$ since the last HST exposure (§3.1), there are quite
substantial changes. All the HST data used in this paper can be found in
MAST:http://dx.doi.org/10.17909/ytdx-cf49 (catalog 10.17909/ytdx-cf49).

###  2.3 Archival IFU apex data

Inspecting the Gemini archive, we found unpublished GMOS-N IFU observations of
the Guitar nebula taken on 2008 June 10 and July 3 (Program GN-2008A-Q-3, van
Kerkwijk, PI). The observations included $4\times 3863$ s exposures with the
IFU-2 mask, the B1200 grating and RG610 filter, covering the head of the
nebula with a $0.2^{\prime\prime}$ fiber grid. A few 300 s direct acquisition
images using the G0310 H$\alpha$ filter and associated calibration files were
also obtained. Conditions were good, meeting the 20% best seeing criterion.

To improve S/N on this faint nebula, the GMOS-N detector was binned $2\times
2$ during this observation. As it happens, the $2\times$ spatial binning left
the fiber traces poorly resolved. This is generally not recommended and,
indeed, meant that the Gemini pipeline software failed to trace the spectra
and extract the data. To complete reductions, we therefore had to mark the
1000 fiber spectra positions by hand near the position of the H$\alpha$
feature in each spectral exposure flat field image and force a low-order trace
centered on these fiber centroids. Importing these traces to the arcs we were
able to use line features (again initially identified by hand) to establish a
good wavelength solution. The traces applied to the target integrations could
then be used to extract and calibrate the 1-D spectra. These were spatially
flattened using sky lines and assembled into position-velocity cubes. The
final weighted combination of these data gave a data cube with a spatial scale
of $0.1^{\prime\prime}$/pixel and a wavelength scale of 0.4729Å/image plane
(21.6 km/s/image plane). The spectral resolution of 1.23Å(56 km/s) was
confirmed by measuring sky lines and the $0.6^{\prime\prime}$ FWHM spatial
resolution of the pre-images is maintained in the data cubes, as indicted by
the width of the narrow H$\alpha$ limb $\sim 2.5^{\prime\prime}$ behind the
apex.

Figure 3: Velocity profile just behind the nebula apex (oval in Figure 2). A
reflected version of the profile (red) shows that it is centered at $\sim-12$
km/s. The high velocity $|v|=150-300$ km/s wings (blue points, shifted from
the red profile to match original apex profile wings) show a slight additional
blue shift.

We find that the nebula is best detected in low velocity channels near the
nebula limb, as expected from neutral excitation in the post shock gas (and by
projection effects). Larger velocities are found principally right behind the
nebula apex, where charge exchange allows accelerated post-shock protons to
obtain electrons and emit H$\alpha$. The line spectrum from a region just
behind the apex is shown in Figure 3. The similarity of the channel maps red
and blue of the central velocity supports the inference from prior bow shock
image fits that the pulsar space velocity lies close to the plane of the sky.
Focusing on the velocity extrema we see an offset of $-32-(-12)=-20$ /km/s for
the wing components peaking at $\sim\pm 250$ km/s, giving
$\theta_{v}\approx{\rm arctan}(-20{\rm km/s/250km/s})\approx 5^{\circ}$ out of
the plane of the sky.

Assuming approximate axial symmetry, we can check these velocities by
examining the transverse expansion of the head boundary in the HST images. By
comparing the 2020 and 2006 HST images we can see that across the bulk of the
head region, starting $5^{\prime\prime}$ back from the 2020 apex, the lateral
expansion is $\approx 0.3^{\prime\prime}$ over 14y or $\approx\pm 85{\rm
km/s}$ at $d=0.83$ kpc, in good agreement with the brightest emission in the
channel maps. It is a bit more difficult to measure the transverse expansion
corresponding to the elliptical apex region marked in Figure 2, but comparing
the 2008 GMOS-N pre-image with the 2006 HST frame, we see transverse expansion
of $\sim 0.1^{\prime\prime}$ over 1.7y or $\sim\pm 230{\rm km/s}$, in
reasonable agreement with the wing component speeds in Figure 3.

##  3 Multi-Epoch Comparison

Over the past 25 years we now have 4 epochs of observation each with HST in
H$\alpha$ and CXO in X-rays (Table 1). Other ground-based H$\alpha$ images
exist, but the very small angular scale of the bow shock stand-off requires
HST resolution for serious study of the evolving morphology. Comparison of the
shock between these epochs gives important clues to the nature of the
filament.

Figure 4: 1994-2021 Guitar/Filament Evolution. Left panels show four
$22^{\prime\prime}\times 21^{\prime\prime}$ cutouts from the HST H$\alpha$
images (1994, 2001, 2006, 2020). The top panel has the pulsar positions at the
four HST and four CXO epochs marked on lines representing 50y
($9.7^{\prime\prime}$) of proper motion. The lower panels have the outline of
the 2020 HST H$\alpha$ limb marked, for comparison. The right panels show the
filament evolution over the four CXO epochs ($4^{\prime}\times 1^{\prime}$
cutouts). The filament leading edge from the 2021 epoch is marked by the blue
line, while the red arrow in the top and bottom panels shows the 50y proper
motion, and a stationary background source is marked for reference.

###  3.1 Optical Evolution

Table 2: The standoff distances (in mas) measured for each HST epoch, for thin
and wide shock models; shown are the median values from the posterior
distributions, with the 14th and 86th percentiles as the errors

. The estimates of Chatterjee & Cordes (2004) and Ocker et al. (2021) are also
shown.  Epoch $r_{0,\rm{thin}}$ $r_{0,\rm{wide}}$ $r_{0,\rm{O2021}}$
$r_{0,\rm{CC2004}}$ 1994 $86_{-9}^{+12}$ $94_{-7}^{+11}$ $77\pm 4$ $120\pm 40$
2001 $116_{-15}^{+16}$ $112_{-14}^{+14}$ $110\pm 10$ $100\pm 40$ 2006
$97_{-4}^{+4}$ $93_{-3}^{+3}$ $94\pm 6$ 2020 $96_{-2}^{+4}$ $92_{-2}^{+3}$

Comparing the optical images at the left hand side of Figure 4 with the line
marking the limb of the 2020 image we see that as the pulsar advances, the
perpendicular expansion is rapid at the apex but slows by a few arcsec behind.
This is also visible in the IFU data cube. The general structure of the Guitar
head is best seen in our new high S/N 2020 image; it is roughly symmetric,
with indentations, especially a ‘pinch’ $\sim 7^{\prime\prime}$ behind the
pulsar, and higher limb brightness regions, e.g. $\sim 3^{\prime\prime}$ and
$\sim 6^{\prime\prime}$ behind the pulsar. Thus the geometry of the bow shock
apex and its expansion rate must vary. The most extreme illustrations of this
are, of course, the apparently closed bubble of the Guitar head and the double
cavity of the Guitar body itself. The approximate bi-lateral symmetry of the
overall nebula indicates either that the central pulsar wind varies or that
the perturbations producing these structures have a coherence scale
substantially larger than the width of the nebula. However, there is also
significant right-left asymmetry, which implicates instabilities in the shock
flow or variations in the external medium on the few arcsec scale of the head
width.

The spectrum of such perturbations have recently been explored by Ocker et al.
(2021), who, following Chatterjee & Cordes (2004) discuss apparent changes in
the bow shock standoff distance in the previous three HST epochs. These
$\theta_{0}$ were, however, estimated by marking the apparent bow shock limb
by hand and then fitting to these marked points. This, of necessity,
introduces substantial subjectivity. We therefore have sought to fit Wilkin
(2000) apex models directly to the HST images. This model computes the locus
of the contact discontinuity, which for a ‘thin’ shock marks the H$\alpha$
front. In practice post-shock pressure widens the structure; the H$\alpha$
emission standoff should be $\sim 1.3r_{0}$ at the apex and this factor should
grow downstream. A simple approximation increases the transverse scale by
$1.25\times$ (Brownsberger & Romani, 2014) for a ‘wide’ shock model. For the
first three epochs we were able to register the frames to Gaia stars to
determine the position of the pulsar in the frame with an $1\sigma$
uncertainty of 0.07, 0.07, and 0.06 pixels respectively. For 2020, however,
WFC corner distortions defeated such registration, so we have let the pulsar
position adjust over a 1 pixel ($\sim 50$mas) range. We used the affine
invariant Markov Chain Monte Carlo (MCMC) algorithm of Goodman & Weare (2010),
implemented through the python package emcee, to sample the likelihood
function and obtain posterior distributions for the standoff distance in each
epoch. The MCMC routine was run using 50 walkers and 5000 steps for each
walker. In order to run MCMC efficiently, we first performed a maximum
likelihood analysis, and started the walkers with initial parameters close to
the best-fit parameters. Upon visual inspection of the chains, we further
excluded the first 500 steps of each walker to ‘burn in’ the chains, ensuring
convergence. We calculate the integrated autocorrelation time $\tau_{f}$ to be
$\sim 50$ steps, meaning that the walker length of 4500 steps should be
sufficient.111In the emcee documentation on autocorrelation analysis
(https://emcee.readthedocs.io/en/stable/tutorials/autocorr/), it is suggested
that each walker should have a length of $>\,50\tau_{f}$ steps, so that enough
independent samples can be obtained to yield accurate results..A visual
comparison of the data and the Wilkin thin shock model for each epoch is shown
in Figure 5. The posterior distributions for $r_{0}$ are plotted in Figure 6,
with values listed in Table 2. For the 2020 epoch, the position uncertainty
from the posterior (the 14th and 86th percentiles of the distribution) is
$\approx 0.15$ pixel ($\sim 8$ mas) in both the $x$ and $y$ directions.

Figure 5: The Wilkin-model bow shock fits for the 4 different H$\alpha$
epochs. All images are shown in native resolution of $0.05\arcsec$/pixel. The
left column shows a cut-out of the bow shock apex for the 4 epochs, using
pixels up to $\approx 1^{\prime\prime}$ behind the apex. The 1994 and 2001
epochs have been lightly smoothed for visualization. The middle column shows
the Wilkin model, and the rightmost column shows the residual (data-
model)2/$\sigma^{2}$. The red cross shows the location of the pulsar - a free
parameter for the 2020 epoch, and referenced to Gaia using a set of reference
stars for the other three epochs. Figure 6: The posterior distributions for
the standoff distance $r_{0}$ in each epoch (left to right, top to bottom:
1994, 2001, 2006, 2020. Shown for reference are the estimated standoff
distances of Ocker et al. (2021) and Chatterjee & Cordes (2004)

The Ocker et al. (2021) estimates generally lie within the $r_{0}$ uncertainty
ranges, but have nominal errors much smaller than we find for a direct fit,
especially for the 1994 and 2001 data. Alas our more realistic errors mean
that direct evidence for $r_{0}$ variation is poor. The bulk of our
uncertainty range suggests that $r_{0}$ was larger in 2001, but even this
result is weak. Additional images of the quality of our new ACS/WFC exposure
are needed to probe stand-off variation at the required $\sim 5$mas level.

Nevertheless the head limb shape does suggest that the standoff was small when
pulsar was at the position of the head’s closed base, $\sim
15.5^{\prime\prime}$ behind the present apex (i.e. in $\sim 1940$). The
transition into the head bubble may be similar to the ‘break-through’ inferred
for the PSR J2030+4415 H$\alpha$ nebula and filament (de Vries & Romani,
2022). The ‘pinch’ $\sim 7^{\prime\prime}$ back (i.e. in 1985), and the
increased limb brightness $\sim 3.2^{\prime\prime}$ behind the apex (in $\sim
2004$) suggest weaker compression events.

Figure 7: Regions used for spectral analysis of the filament in the 2021
Chandra epoch. We divided the filament into three main sections: the inner
section where the leading edge is sharpest (0–0.7′), the middle section where
the filament appears to become more diffuse (0.7–1.3′) and the outer, most
diffuse section where the sharp leading edge has largely disappeared
(1.3–2.4′). Additionally, the inner and middle sections are divided by the
green dashed line into the ’Leading’ and ’Trailing’ regions. We also identify
a ‘Remnant’ region of bright emission around $12\arcsec$ behind the leading
edge in the middle section. The red regions show the contours of the Guitar
nebula head (solid line) and body (dashed line) in the 2020 HST H$\alpha$
image.

###  3.2 X-ray Evolution

Figure 8: Light-curves of the leading edge ($1^{\prime\prime}$ width) of the
inner filament (0–0.7′ segment), showing the flux in each epoch and the flux
of that same region of the sky in following epochs. The y-axis indicates the
1-5 keV photon surface brightness. Figure 9: Light-curves of the inner
counter-filament and inner and middle filament sections. Data points have been
slightly offset from each other for legibility. The y-axis indicates the 1-5
keV photon flux per arcmin length of filament. The lengths of the inner
counter-filament, and inner and middle filament sections are 0.25′, 0.65′, and
0.68′ respectively and fluxes have been integrated across the width of the
main filament.

In Figure 7 we define several regions useful in describing the filament’s
spectrum and its evolution. The $e^{\pm}$ injection site shifts with the
steady pulsar motion, and in Figure 4 it is apparent that the filament leading
edge marches along with the pulsar, as also noted by Wang (2021). Our deeper
2021 exposure provides a much better view of the counter-filament (CF) than
earlier epochs. It extends at least $20^{\prime\prime}$ and likely
$40^{\prime\prime}$ from the pulsar. Interestingly it does not line up well
with the filament leading edge, instead intersecting the proper motion axis
some $1.5^{\prime\prime}$ behind the pulsar position. Both it and the filament
have substantial curvature near the bow shock. This is likely a field line
‘draping’ effect or field distortion from supra-thermal particles as most
clearly seen in the ‘lighthouse’ PWN filament and counter-filament (Pavan et
al., 2016).

We checked to see if the PWN PSR trail is detected in our deep 2021 image.
Using the ‘head region’ of the Guitar (see Figure 7) as an aperture and
subtracting similar flanking regions as background, we find an excess of
$11\pm 4.5$ counts in the 0.7–5 keV, range, a marginal $2.5\sigma$ detection.
This gives a filament/trail flux ratio $>100$, the largest among known
filaments.

In Figure 8 we measure the surface brightness at the filament leading edge in
each epoch and compare the flux in the same aperture in subsequent epochs. In
general the region corresponding to the edge shows an initial rapid decrease
in surface brightness in the following epoch (see also Figure 4). We infer a
rapid change in the electron population as the pulsar moves ahead to the next
set of field lines, due to cooling, advection or diffusion. The subsequent
brightness decrease, if any, is much smaller. Fortuitously the leading edge
was much brighter than usual during the original 2000 epoch, which helped in
the filament’s discovery. This may be related to enhanced injection around
this epoch. For example if $r_{0}$ decreases, then more pulsar/PWN shock
particles have gyroradii exceeding $r_{0}$, so escape to the filament might
increase and the filament surface brightness may temporarily increase.
Averaged over the full width the fluxes per unit length seem quite constant
(Figure 9); although the inner counter-filament appears more prominently in
the 2021 image, its flux per unit length remains consistent with that of the
filament, within errors.

Although the statistics are limited in the early images, there appear to be
changes behind the leading edge. In Figure 4 the most notable changes are in
the ‘Middle’ section of the filament where the emission spreads behind the
leading edge as a shifting ridge. We quantify this trend in Figure 10, where
fits to Gaussian distributions transverse to the filament show a progressive
shift and broadening of the maximum. Note that the integral flux is consistent
with constant across the four epochs. A fit to such regions in the ‘Inner’
zone gives similar evolution with nearly identical parameters, but lower
statistical significance. We attempt to interpret these results in the
conclusions.

Figure 10: Spread of the emission behind the leading edge in the filament
middle section across four epochs. The data are well described by a steady
shift of the peak behind the leading edge, a steady increase in the width, and
a constant integrated flux. The y-axis indicates the 1–5 keV photon surface
brightness. Labels indicate the best-fit parameter value (with $1\sigma$
errors ): $x_{2000}$ and $\sigma_{2000}$ are the peak position and standard
deviation of the Gaussian component in 2000 respectively; $\dot{\sigma}$
indicates the increase in $\sigma$ over time; and $\mu$ indicates the shift of
the peak away from the leading edge over time.

###  3.3 Spectral fits

We have extracted spectra for the several regions of Figure 7 using the
standard CIAO tools. To each of the spectra, we have fit a power law
multiplied by Galactic absorption, which we have set at $2.7\times 10^{21}$
cm-2. The results of the spectral analysis are shown in Table 3. There are no
significant differences in spectral index between the regions. The weak
evidence for spectral softening with distance from the pulsar would require
much deeper observation for a serious test.

Additionally, we have estimated the magnetic field strength under the
assumption of equipartition. For an optically thin region filled with
relativistic electrons and magnetic field emitting synchrotron radiation

$B=46\left[\frac{J_{\rm-20}(E_{1},E_{2})\sigma}{\phi}\frac{C_{1.5-\Gamma}(E_{m},E_{M})}{C_{2-\Gamma}(E_{1},E_{2})}\right]^{2/7}\mu
G$ (1)

where

$C_{q}(x_{1},x_{2})=\frac{x_{2}^{q}-x_{1}^{q}}{q}.$ (2)

$J_{\rm-20}(E_{1},E_{2})=4\pi f_{\rm-20}(E_{1},E_{2})d^{2}/V$ is the observed
emissivity (in $10^{-20}$ erg s-1 cm-3, between $E_{1}$ keV and $E_{2}$ keV),
$\sigma=w_{B}/w_{e}$ is the magnetization parameter, $\phi$ the filling
factor, and $E_{m}$ and $E_{M}$ the minimum and maximum energies, in keV, of
the synchrotron spectrum with photon $\Gamma$. We assume that the structures
are cylindrical, with diameter set to the observed region width. We list the
derived equipartition fields in Table 3 for $\sigma=\phi=1$,
$E_{m}=0.01\,\rm{keV}$ and $E_{M}=10\,\rm{keV}$.

Table 3: Spectral fit results (with $1\sigma$ errors) for the filament in the
2021 Chandra epoch (see Figure 7 for the regions). The ‘Leading’ and
‘Trailing’ regions are composed of the front and back halves, respectively of
the combined inner and middle regions. Beq was computed assuming cylindrical
volumes for each region.

Region | Counts | $\Gamma$ |  $f_{-15}$b |  $\chi^{2}$/DoF | Beq
---|---|---|---|---|---
|  |  |  |  | [$\mu G$]
Inner | $214\pm 17$ | $1.31\pm 0.16$ | 9.9 | 29.3/27 | 13
Middle | $209\pm 17$ | $1.37\pm 0.17$ | 10.2 | 24.3/24 | 14
Outer | $489\pm 32$ | $1.58\pm 0.15$ | 24.1 | 53.2/48 | 8
CF | $86\pm 11$ | $1.71\pm 0.30$ | 3.5 | 23.8/24 | 17
Leading | $273\pm 19$ | $1.39\pm 0.14$ | 13.6 | 22.7/33 | 19
Trailing | $154\pm 16$ | $1.60\pm 0.20$ | 7.1 | 30.5/30 | 17
Remnant | $174\pm 19$ | $1.40\pm 0.27$ | 7.2 | 38.0/33 | 9

a $N_{H}$ fixed at $2.7\times 10^{21}{\rm cm^{-2}}$.

b $0.5-7\,$keV unabsorbed fluxes in units of $10^{-15}{\rm
erg\,cm^{-2}s^{-1}}$.

##  4 Discussion and Conclusions

The shape of the filament is complex and the epoch-to-epoch changes are
subtle. We seek to explain these through a combination of variable particle
injection at the moving pulsar, particle flow along field lines, particle
diffusion across field lines and possible cooling. In practice cooling is
likely not important on the scale of the observed filament since standard
synchrotron theory gives a cooling time of

$\tau\approx 7.6\times 10^{4}E_{\rm keV}^{-1/2}B_{\mu G}^{-3/2}{\rm y}.$ (3)

With an observed photon energy of $\sim 2$ keV and $B_{\mu G}\sim 15$
estimated in §2.3, we get a cooling time $\tau\approx 930$ y, so over our four
epochs we expect no significant cooling. Since the pulsar moves $\sim
3^{\prime}$ (twice the size of the Guitar body) in time $\tau$, cooling
predicts a fading on this scale. Accordingly, the smaller scale morphology
changes must be due to variable injection, advection and diffusion.

In the original Bandiera (2008) picture the $r_{c}$ relevant for escape was
that of the shocked pulsar wind. Since the mean field in that wind increases
as $r_{0}$ decreases, $r_{c}/r_{0}$ is essentially constant, and does not
control the particle escape; in this picture most bow-shock pulsars should
produce filaments and they should do so at all epochs independent of the bow
shock size. This does not appear to be the case, since filaments are rare and
preferentially associated with pulsars with small $r_{0}$. Instead we argue
that energetic $e^{\pm}$ are produced via reconnection throughout the shocked
pulsar wind and that $r_{c}$ beyond the contact discontinuity, in the shocked
ISM and external medium, controls escaping particle motion. The near-apex
external field is modified by the draping effect to have a characteristic
curvature radius $r_{0}$ and thus $r_{c}/r_{0}$ in this medium can control
which particles move far enough in a gyroradius to encounter different
external field orientations, and escape.

The curvature of the filament leading edge implies that the ambient field
lines are not completely straight, although the similarity of the edge from
epoch to epoch suggests that they are locally approximately parallel. The
leading edge is quite sharp. The $e^{\pm}$ gyroradius $r_{c}$ in the local
field subtends an angle of

$\theta=r_{c}/d\approx 26^{\prime\prime}E_{\rm keV}^{1/2}B_{\mu
G}^{-3/2}d_{\rm kpc}^{-1}$ (4)

for particles producing a peak photon energy $E_{\rm keV}$. For a leading edge
field of $20\mu G$ (Table 3) we get $\theta\approx 0.5^{\prime\prime}$ for the
filament. This is comparable to the CXO resolution (but substantially larger
than $r_{0}$). The filament leading edge stays sharp for the inner and middle
zones, spreading primarily in the outer zone. This implies that the cross
field diffusion coefficient ahead of the leading edge in the ambient ISM is
small.

With an estimate of the flow speed $v_{e^{\pm}}$ along the leading edge, one
could use the broadening with distance to get an estimate of this forward
diffusion coefficient. Noting that this edge is actually the front reached by
particles moving rapidly along a set of field lines, we see that the far
filament represents earlier injection, onto field lines behind that connect to
the pulsar at its current position. Thus the filament front follows an angle
$\theta_{f}\sim v_{\rm PSR}{\rm cos}\Psi/v_{e^{\pm}}$ behind the ISM field
lines, with the field lines themselves at an angle $\Psi\sim 25^{\circ}$ to
the proper motion. If the filament and counter-filament propagation speeds are
equal, we can account for $\Psi$ by comparing the PAs of the two sides; these
should differ by $2\theta_{f}$. In practice this measurement is difficult
since the counter-filament is short and the section closest to the bow shock
suffers PWN-induced distortion. Very roughly, we estimate $\theta_{f}\lesssim
2^{\circ}$, and thus $v_{e^{\pm}}\approx v_{\rm PSR}{\rm
cos}\Psi/\theta_{f}\gtrsim c/13$.

With small cross-field diffusion, we would expect particles confined to their
injection field line and the filament would present an approximately uniform
band, shifted increasingly farther from the Guitar axis, since particles on
field lines to the rear would have more time to propagate away. This band
would have brighter ridges marking times (field lines) of enhanced particle
injection and a smooth fading on arcmin scales behind the leading edge due to
synchrotron cooling. This is not what we see. Instead the emission behind the
leading edge is patchy and seems to evolve on times short compared to the
cooling times.

This may be understood if cross-field advection and diffusion increase behind
the leading edge. From Figure 10 we estimate the 2021 surface brightness peak
as having position $x_{2021}\approx 12.2^{\prime\prime}$ behind the 2021
leading edge, with bulk motion of $\mu_{\rm ridge}\sim 0.15^{\prime\prime}{\rm
y^{-1}}$ and spread of $\sigma(t)\sim[3.0+0.25(t-2000)]^{\prime\prime}$. We
can attribute these increased rates to turbulence induced behind the leading
edge by the injected particles; this leads to increased scattering and easier
cross-field propagation. It then becomes interesting to trace the origin of
the ridge that moves through Figure 10. With a coordinate increasing normal to
and behind the leading edge, we can write the pulsar position at year $t$ as
$x_{p}=\mu_{PSR}{\rm cos}\Psi(2021-t)$. Similarly $x_{\rm
ridge}=x_{2021}+\mu_{ridge}(t-2021)$, with the proper motions in arcsec/y.
Finally the propagation time between the pulsar and the middle zone $l\sim
75^{\prime\prime}d\sim 1$ lt-y away is $t_{\parallel}\approx l/v_{e^{\pm}}\sim
c/v_{e^{\pm}}$ years. Thus the date for the enhanced injection of the
$e^{\pm}$ that we see in 2021 as a ridge moving through the filament is

$t_{\rm inj}\approx 2021-x_{\rm r,2021}/(\mu_{\rm PSR}{\rm cos}\Psi+\mu_{\rm
ridge})+t_{\parallel}.$ (5)

From our fit to the ridge evolution we get $t_{\rm inj}\approx
1993_{-16}^{+7}+t_{\parallel}$, so to identify the moving ridge with particles
injected when the pulsar was at the ‘pinch’ in the Guitar head,
$7^{\prime\prime}$ behind the apex, in 1985, we would want $\mu_{\rm ridge}$
low in the fit range and $t_{\parallel}=l/v_{e^{\pm}}<8$ y. Note that $\sigma$
decreases to 0 at $1990_{-10}^{+5}$, so consistent with $\sim 1985$, as well.

It is likely a coincidence that the back-propagation of the ridge brings it
nearly parallel with the base of the Guitar head in the 2021 epoch. Although
we don’t see strong emission at this position in our earlier epochs, those
images are shallow, and it is possible that future deep observations will show
that this ridge is a permanent feature fixed in space. In that case it would
be compatible with the simpler hypothesis that it is the fossil of strong
injection at the point the pulsar broke into the head region in $\sim 1940$.
More generally the lack of such ‘fossil’ X-ray emission parallel with the
Guitar body suggests that when the pulsar was blowing the bubbles
corresponding to the body structure, the standoff $r_{0}$ was large and that
little or no $e^{\pm}$ escape occurred. Thus the Guitar may have had an X-ray
filament only since the very compact head region was formed.

Although forward propagation of the pulsar-generated cosmic rays is severely
limited at the leading edge, it seems much freer behind, rearranging the
injected particles long before they cool. While this means that the filament
surface brightness profile is not a simple historical record of injection
history, it does offer the opportunity to probe the diffusion of multi-TeV
$e^{\pm}$ through the ISM and, more importantly, their effect, via induced MHD
waves, on the local particle propagation. For B2224, injection seems to be
effective over the $\sim 20^{\prime\prime}$ region of the Guitar head and neck
where the bow shock standoff was evidently small, leading to a wide filament.
In contrast, the filament of PSR J2030+4415 stays narrow since the injected
period was short and the pulsar covered little distance in this time. This
picture of variable injection may certainly be tested by finding more example
filaments, and connecting them with bow shock properties. Numerical
simulations can also be useful in determining whether external field-
controlled escape is viable or some other peculiarity of small $r_{0}$ bow
shocks, such as enhanced local turbulence or asymmetric reconnection to the
external fields, needs to be invoked.

The effect of injected $e^{\pm}$ may be especially important in connection
with the recently observed TeV halos around nearby young pulsars (Abeysekara
et al., 2017). Interesting experiments to probe these effects would be enabled
by identifying enhanced filament injection events (plausibly via H$\alpha$ bow
shock monitoring) followed by a decade of sensitive X-ray images to observe,
via synchrotron emission, the TeV $e^{\pm}$ pulse propagating and spreading.
Such a campaign would be expensive in observation time, but would yield a rich
harvest of information of cosmic ray and magnetic field dynamics, spread out
before the observer in evolving filament images.

We wish to thank the observatory staff who helped in planning the exposures
described in this paper, especially Jean Connelly of the CfA for help with CXO
and Alison Vick and Ray Lucas of STScI for help with the ACS/WFC. We also wish
to thank Marten van Kerkwijk for advice on the GMOS-N IFU data set. MdV and
RWR were supported in part by NASA grant G08-19050A, through the Smithsonian
Astrophysical Observatory. GGP was supported by NASA grant G08-19050B. OK was
supported by NASA grant GO8-19050C and ADAP grant 80NSSC19K0576. Support for
this work was provided by the National Aeronautics and Space Administration
through Chandra Award Number GO8-19050 issued by the Chandra X-ray Observatory
Center, which is operated by the Smithsonian Astrophysical Observatory for and
on behalf of the National Aeronautics Space Administration under contract
NAS8-03060.

## References

  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911, doi: 10.1126/science.aan4880
  * Bandiera (2008) Bandiera, R. 2008, A&A, 490, L3, doi: 10.1051/0004-6361:200810666
  * Barkov et al. (2019) Barkov, M. V., Lyutikov, M., Klingler, N., & Bordas, P. 2019, MNRAS, 485, 2041, doi: 10.1093/mnras/stz521
  * Brownsberger & Romani (2014) Brownsberger, S., & Romani, R. W. 2014, ApJ, 784, 154, doi: 10.1088/0004-637X/784/2/154
  * Chatterjee & Cordes (2004) Chatterjee, S., & Cordes, J. M. 2004, ApJ, 600, L51, doi: 10.1086/381498
  * Cordes et al. (1993) Cordes, J. M., Romani, R. W., & Lundgren, S. C. 1993, Nature, 362, 133, doi: 10.1038/362133a0
  * de Vries & Romani (2022) de Vries, M., & Romani, R. W. 2022, ApJ, 928, 39, doi: 10.3847/1538-4357/ac5739
  * Deller et al. (2019) Deller, A. T., Goss, W. M., Brisken, W. F., et al. 2019, ApJ, 875, 100, doi: 10.3847/1538-4357/ab11c7
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65, doi: 10.2140/camcos.2010.5.65
  * Ocker et al. (2021) Ocker, S. K., Cordes, J. M., Chatterjee, S., & Dolch, T. 2021, ApJ, 922, 233, doi: 10.3847/1538-4357/ac2b28
  * Olmi & Bucciantini (2019) Olmi, B., & Bucciantini, N. 2019, MNRAS, 490, 3608, doi: 10.1093/mnras/stz2819
  * Pavan et al. (2016) Pavan, L., Pühlhofer, G., Bordas, P., et al. 2016, A&A, 591, A91, doi: 10.1051/0004-6361/201527703
  * Wang (2021) Wang, Q. D. 2021, Research Notes of the American Astronomical Society, 5, 5, doi: 10.3847/2515-5172/abd854
  * Wilkin (2000) Wilkin, F. P. 2000, ApJ, 532, 400, doi: 10.1086/308576