# UncertaintyPlayground: A Fast and Simplified Python Library for Uncertainty
Estimation

Ilia Azizi
University of Lausanne
Lausanne, Switzerland
<EMAIL_ADDRESS>

###### Abstract

This paper introduces UncertaintyPlayground, a Python library built on PyTorch
and GPyTorch for uncertainty estimation in supervised learning tasks. The
library offers fast training for Gaussian and multi-modal outcome
distributions through Sparse and Variational Gaussian Process Regressions
(SVGPRs) for normally distributed outcomes and Mixed Density Networks (MDN)
for mixed distributions. In addition to model training with various
hyperparameters, UncertaintyPlayground can visualize the prediction intervals
of one or more instances. Due to using tensor operations, the library can be
trained both on CPU and GPU and offers various PyTorch-specific techniques for
speed optimization. The library contains unit tests for each module and
ensures multi-platform continuous integration with GitHub Workflows (online
integration) and Tox (local integration). Finally, the code is documented with
Google-style docstrings and offers a documentation website created with MkDocs
and MkDocStrings.

_Keywords_ Uncertainty Estimation $\cdot$ Python Library $\cdot$ Gaussian
Processes $\cdot$ Mixed Density Network

##  1 Introduction

Uncertainty estimation is a critical aspect of machine learning, providing
insights into the reliability of predictions. While efficient libraries such
as PyTorch [1] and GPyTorch [2] offer great flexibility, they are designed to
be relatively low-level and, deliberately, do not offer process abstraction.
This paper introduces
UncertaintyPlayground111https://unco3892.github.io/UncertaintyPlayground/, a
Python library that facilitates uncertainty estimation in supervised learning
tasks, offering a user-friendly yet performant interface to powerful,
parallelized uncertainty estimation techniques.

The advent of Neural Networks (NNs) and Graphics Processing Units (GPUs) has
revolutionized uncertainty estimation. NNs, with their ability to model
complex, non-linear relationships, provide a flexible framework for
uncertainty estimation. This is further enhanced by Mixed Density Networks
(MDNs), which allow for the modeling of complex, multi-modal distributions.
The parallel processing capabilities of GPUs expedite the training of these
large, complex models, making it feasible to implement computationally
intensive techniques such as Variational Inference and Kullback–Leibler (KL)
divergence. UncertaintyPlayground leverages these advancements, offering
Gaussian Processes (GPs) and MDNs for regression tasks. The library abstracts
away the complexities of these techniques, allowing users to focus on their
data and tasks. Designed to be user-friendly yet performant, it serves as a
valuable tool for researchers and practitioners in machine learning, as well
as less technical users and experts in statistics, economics, and
bioinformatics.

The rest of the paper is organized as follows: Section 2 briefly covers the
research question, necessary theoretical foundations, existing literature, and
similar libraries for uncertainty estimation. Section 3 discusses the
library’s algorithm and workflow, including modules for model definition,
training, and plotting prediction distributions. This section also discusses
the parallelization capacity of the package . Section 4 describes the unit
tests necessary for ensuring the build and outlines code maintenance
strategies, including a local and an online workflow for continuous
integration. Section 5 discusses the process of documenting the code and
creating a website for the code documentation. Section 6, through an example,
demonstrates the output of UncertaintyPlayground for the user. Section 7
provides additional ideas for further developing the library, and Section 8
concludes the paper.

##  2 Research Question and Relevant Literature

###  2.1 Research Question

This project’s primary research question is how to estimate uncertainty in
supervised learning tasks in an easy-to-use manner. The relevant literature
includes works on GPs [3] and MDNs [4]. The GP implementation is based on
Sparse and Variational Gaussian Process techniques (SVGP) [5]. We briefly
highlight SVGP and why it was chosen for this library.

Neural Networks provide a flexible framework for modeling complex patterns in
data. In particular, SVGPRs and MDNs, which are neural network-based
implementations of Gaussian Processes and Mixture Models, offer a powerful
approach to modeling uncertainty. SVGPRs allow for efficient and scalable
Gaussian Process Regression, while MDNs enable the modeling of complex, multi-
modal distributions. UncertaintyPlayground aims to provide these neural
network-based techniques and tackle issues with these models’ algorithmic
complexity through methods such as Variational Inference and KL divergence.

###  2.2 Relevant Literature

####  2.2.1 Gaussian Processes

Gaussian Processes (GPs) have been recognized as a fundamental tool in machine
learning, providing function priors for various tasks. In cases with Gaussian
likelihoods, inference can be performed in a closed-form manner. However, for
non-Gaussian likelihoods, posterior and marginal likelihood approximations are
required [6, 7]. Traditional GP regressions have a significant drawback: the
computational cost of the exact calculation of the posterior and marginal
likelihood scales as $O(N^{3})$ in time and $O(N^{2})$ in memory, where $N$ is
the number of training examples.

Sparse and Variational Gaussian Process Regression (SVGPR) is a particular
implementation of GPs that allows for efficient and scalable Gaussian Process
Regression [5, 8]. SVGPR addresses this issue by choosing $M$ inducing
variables that are $M\ll N$ to summarize the entire posterior. This reduces
the computational cost to $O(NM^{2}+M^{3})$ in time and $O(NM+M^{2})$ in
memory [9]. SVGPR works by minimizing the Kullback-Leibler (KL) divergence
between the approximate and the true posterior, which allows us to learn the
model parameters via gradient descent. This is significantly lower than the
$O(N^{3})$ time complexity of standard Gaussian Process Regression, making
SVGPR a more scalable alternative for large datasets [8].

####  2.2.2 Mixed Density Networks

Mixed Density Networks (MDNs) represent another powerful approach in machine
learning, particularly for modeling complex, multi-modal distributions. MDNs
combine the flexibility of neural networks with the robustness of mixture
models, enabling them to capture intricate patterns in data [4].

The computational complexity of MDNs is primarily dependent on the number of
mixture components and the dimensionality of the data. Unlike GPs, the
complexity of MDNs does not increase cubically with the number of data points.
However, as the number of mixture components or the dimensionality of the data
increases, the computational cost of training an MDN can become substantial.
Despite this, the parallel processing capabilities of modern hardware, such as
Graphics Processing Units (GPUs), can significantly expedite the training of
these models, making MDNs a feasible approach for large-scale, high-
dimensional tasks.

##  3 Code Base

The methodology involves the implementation of SVGP and MDN using the GPyTorch
[2] and PyTorch [1] libraries. As mentioned, the GP implementation is based on
SVGP, hence SVGPR, while MDN is implemented purely in PyTorch. Other
dependencies for this project are Numpy [10], Scikit-learn [11], Matplotlib
[12], and Seaborn [13]. Details regarding the version of these packages can be
found in the README.md and requirements.txt files.

The entire code base (version 0.1) has been depicted in Figure 1. The workflow
of the package goes through models, trainers (supported via supplementary sub-
modules in utils), and then the predictions are plotted via predplot. After
the development of the library, several unit tests were designed to ensure its
functionality and placed in the tests directory. Finally, the documentation
and some examples can be found in docs and examples folders.

.1 UncertaintyPlayground/. .2 README.md. .2 setup.py. .2 LICENSE. .2
MANIFEST.in. .2 tox.ini. .2 mkdocs.yml. .2 docs/. .3 bib.md. .3 example.md. .3
gen_ref_pages.py. .3 index.md. .3 README.md. .2 examples/. .3
compare_models_example.py. .3 mdn_example.py. .3 svgp_example.py. .2
uncertaintyplayground/. .3 requirements.txt. .3 models/. .4 mdn_model.py. .4
svgp_model.py. .3 trainers/. .4 base_trainer.py. .4 mdn_trainer.py. .4
svgp_trainer.py. .3 predplot/. .4 grid_predplot.py. .4 mdn_predplot.py. .4
svgp_predplot.py. .3 utils/. .4 early_stopping.py. .4 generate_data.py. .3
tests/. .4 test_early_stopping.py. .4 test_generate_data.py. .4
test_mdn_model.py. .4 test_mdn_predplot.py. .4 test_mdn_trainer.py. .4
test_svgp_model.py. .4 test_svgp_predplot.py. .4 test_svgp_trainer.py.

Figure 1: Project structure of UncertaintyPlayground v0.1

###  3.1 Models

As mentioned already, in the current version of the package, two models are
implemented: an MDN and an SVGP, defined in the modules mdn_model and
svgp_model, respectively.

The MDN model is designed to predict multi-modal distributions. It is
particularly useful in scenarios where the data does not adhere to a simple
Gaussian distribution. This model includes a neural network with one hidden
layer comprised of 10 neurons and a Tanh (hyperbolic tangent) activation
function. This network generates parameters for a mixture of Gaussians, which
include mixture weights (z_pi), means (z_mu), and standard deviations
(z_sigma).

The model provides three options for prediction methods: ‘max_weight_mean’,
‘max_weight_sample’, and ‘average_sample’. These approaches represent
different strategies for choosing components from the Gaussian mixture and
generating predictions:

  * •

Max Weight Mean: This method selects the component with the highest weight and
uses its mean as the prediction. It is a deterministic strategy and tends to
produce the most probable prediction.

  * •

Max Weight Sample: This method selects a component from the mixture based on
the weights, then samples from the selected Gaussian component to yield a
prediction. This is a stochastic strategy, introducing variability into the
predictions.

  * •

Average Sample: This method generates multiple samples from the mixture, each
time selecting a component based on the weights and sampling from the selected
Gaussian component. The final prediction is the average of these samples,
providing a balance between the deterministic ’max_weight_mean’ and the
stochastic ’max_weight_sample.’

The prediction method selection depends on the particular use case and the
trade-off between variability and predictability in the predictions. In all
cases, a loss function, mdn_loss, computes the target variable’s negative log-
likelihood given the mixture’s predicted parameters. The forward pass of the
model computes these parameters based on the input, and a sample method is
included to generate samples from the output distribution based on the input.

The SVGP model is a Sparse Variational Gaussian Process model designed for
large-scale regression tasks where traditional Gaussian Process models are
computationally infeasible. It employs a subset of data points, referred to as
inducing points, to approximate the full Gaussian process. This approximation
reduces the computational complexity from cubic, making the uncertainty
estimation more scalable to large datasets. The SVGP model inherits from the
ApproximateGP class provided by the gpytorch library. The model is initialized
with a variational distribution strategy, defined using the provided inducing
points and data type. The model is designed to learn the inducing locations
during the training process. The SVGP model includes a constant mean module
(mean_module) and a scaled radial basis function (RBF) kernel (covar_module).
The forward pass of the model computes a multivariate normal distribution with
a mean and covariance provided by the mean and covariance modules.

While both models compute the output distribution given input during the
forward pass, they differ significantly in their representation of the output
distribution. The MDN uses a mixture of Gaussians to represent complex, multi-
modal distributions. In contrast, SVGP employs Gaussian processes, more
specifically, a sparse approximation, suitable for modeling functions with
smooth, continuous outputs in large-scale settings.

###  3.2 Utilities

The utils module in the project comprises two smaller sub-modules:
early_stopping and generate_data. These sub-modules offer vital functionality
to control the model training process and generate synthetic data for testing,
respectively.

The sub-module early_stopping encompasses the EarlyStopping class, which is
designed to halt the training procedure of a model when a specified
performance metric ceases to improve over several consecutive epochs. The
extent of tolerance for non-improvement is decided by the patience parameter.
Furthermore, a comparison function, compare_fn, compares the metric values. By
default, the function evaluates with a logical ’less than’ comparison,
implying that lower metric values are superior since we minimize the loss
functions in both situations. In case of an improvement in the validation
metric, the model’s state is saved. This feature is crucial as it aids in the
prevention of overfitting and reduces unnecessary computation.

The generate_data sub-module provides a function, generate_multi_modal_data,
which generates multi-modal data helpful in testing the models. This function
accepts the total number of data samples to be generated, along with a list of
dictionaries specifying the modes of the distribution. Each dictionary
represents a mode defined by its mean, standard deviation, and weight. Here,
the weight dictates the proportion of total samples that will be drawn from
the specific mode. The function returns a NumPy array of the generated data
samples. Such functionality is pivotal for generating synthetic data that
mimics the complex, multi-modal distributions the models aim to handle.

###  3.3 Trainers

The training aspect of the project is catered by a class BaseTrainer and the
two model-specific trainers classes: MDNTrainer and SVGPTrainer. The base
trainer provides the general functionality for training machine learning
models, which is then specialized for training the MDN and SVGP models in
their respective trainers.

The BaseTrainer class provides a basic structure for the training procedure,
defining common functionalities that are universally required for model
training, such as setting up the training, test, and validation datasets via
PyTorch DataLoaders, defining the loss function and optimizer, and setting up
the training loop just to name a few arguments. This base class is designed to
be extended by more specific trainer classes to train particular models, and
the training process itself is implemented within the child classes. One
important capability of this class is providing sample weights to have a
weighted loss function. By default, all the weights are set to 1.

One child class of BaseTrainer is the MDNTrainer. This class is specialized
for training the MDN model, as defined in the mdn_model module. The trainer
initializes the MDN model with a specified number of hidden neurons and
Gaussian components and employs a specified optimizer for training. During
training, the trainer iterates over the data batches, using the MDN’s loss
function to compute the training loss and update the model parameters.
Additionally, the trainer calculates the validation loss after each epoch to
track the model’s performance on unseen data. It utilizes the EarlyStopping
mechanism mentioned earlier to halt training if the validation loss does not
improve over a specified number of epochs, which helps prevent overfitting.

The MDNTrainer also includes a method for making predictions with uncertainty,
predict_with_uncertainty. This method feeds new data into the trained MDN
model and returns the parameters of the predicted Gaussian mixture, including
the mixture weights, means, and standard deviations, as well as samples from
the predicted distribution.

The SVGPTrainer follows a similar structure and role as the MDNTrainer,
providing functionality for training the SVGP model defined in the svgp_model
module. During the initialization of the SparseGPTrainer, the class takes in
the number of inducing points as an argument and initializes the SVGP model
with these inducing points. It also sets up a Gaussian likelihood for the
model. The training method in this class contains a training loop that uses
the VariationalELBO loss function, which is related to the KL divergence
mentioned earlier. Similar to MDNTrainer, the method also includes early
stopping functionality, checking the validation loss after each epoch and
stopping the training if the validation loss does not improve for a set number
of epochs.

The function predict_with_uncertainty varies between SVGPTrainer and
MDNTrainer due to the nature of the models they handle. The trainer in
SVGPTrainer also contains a predict_with_uncertainty method to make
predictions and the associated uncertainty using the trained model. In
SVGPTrainer, this method passes the input tensor through the trained SVGP
model and its Gaussian likelihood, producing a Gaussian predictive
distribution for each data point. The mean and variance of these distributions
are then interpreted as the predicted output and its associated uncertainty.
In contrast, MDNTrainer works with MDNs, yielding a Gaussian Mixture Model
(GMM) as output. This results in multiple Gaussian components per data point
(three per component), each contributing to the final prediction and its
uncertainty. Consequently, despite both methods aiming to estimate prediction
uncertainty, the structure of the output distributions and the procedure of
extracting uncertainty information differ considerably.

###  3.4 Prediction Plots

The predplot module is dedicated to visualizing the models’ predictions and
the corresponding uncertainty. This module includes three sub-modules;
svgp_predplot, mdn_predplot and grid_predplot. As the names suggest, the first
two sub-modules are for plotting the predictions with MDN and SVGP models for
a single instance, while grid_predplot takes as input any of the two functions
and plots the same function for two or more instances. In all cases,
Matplotlib is the main library for visualizations.

The svgp_predplot, mdn_predplot sub-modules each contain a
compare_distributions_ (svgpr or mdn) function which takes a trained MDN
model, an input instance, and optionally the actual outcome(s), and plots the
predicted distribution for that given instance, and compare the final scalar
prediction against the actual outcome value (given that the latter is
provided). The predicted distribution is shown as a Kernel Density Estimate
(KDE) with the Seaborn library, and the predicted value and actual value(s)
are represented as vertical lines on the plot.

The two functions were separated for several reasons. For the MDN model, the
dedicated function plots the learned GMM for each test point, indicating (and
printing) the predicted means, variances, and weights of the Gaussian
components. For the SVGP model, the function plots the predicted mean and
variance at each test point. As a result, there are many more possibilities
for the MDN techniques regarding how the different Gaussian distributions can
be represented (for instance, separately or KDE as it currently stands).

The grid_predplot sub-modules contains a function and a class:
plot_results_grid and DisablePlotDisplay. The plot_results_grid method
generates a grid of scatter plots for given data sets. Each subplot presents
the ground truth data as a scatter plot while superimposing the model’s
prediction and the associated uncertainty as a shaded area. This method is
generic and can be used for both models. Additionally, this sub-module
contains a DisablePlotDisplay class, a special utility class included in the
codebase to control the display of plots during the execution of the script.
This class is especially useful when running the code in an environment where
the graphical display of plots is not possible or not desirable, such as with
our automated tests for the plots.

###  3.5 Parallelization

The UncertaintyPlayground codebase effectively leverages parallelization at
various stages to speed up the computations. PyTorch’s native support for GPU
acceleration is used to parallelize computations involved in model training.
This functionality is implemented within the BaseTrainer class, where the
device attribute determines whether computations are to be performed on a CPU
or GPU.

In addition to GPU acceleration, the software utilizes parallelization when
loading data, significantly reducing the data-loading time, especially for
large datasets. The data-loading step in PyTorch can be easily parallelized
across multiple CPUs using the num_workers parameter of the DataLoader class.
By setting num_workers to be greater than zero, the data loading tasks are
divided among multiple subprocesses. This approach allows to load data in
parallel and ensures that the GPU is not idle. In contrast, the data is being
loaded, thus maximizing the GPU utilization and overall performance.

While currently, the project does not explicitly support multi-GPU
configurations for model training; this could be added in the future using
PyTorch’s DataParallel wrapper for distributed computation. This could further
speed up model training, especially for large models and datasets. However,
this enhancement would require careful management of memory and
synchronization of model parameters across the different GPUs, which would be
a complex task warranting further investigation.

##  4 Code Maintenance

The codebase is shared and maintained using GitHub, which facilitates
collaborative development. The repository includes a .gitignore file,
containing specific files and directories that Git should ignore, such as the
package build files, logs, and local configurations that are not intended to
be part of the shared codebase. Additionally, we ignore the files that are
generated when running local integration with Tox and Pyenv.

Unit tests are implemented using Python’s built-in unittest framework. There
are eight unit tests in total, ensuring the correctness and robustness of the
codebase. The use of unit tests also contributes to the modularity of the
code, as each component can be tested independently.

Continuous Integration (CI) was designed to work both locally (offline) and
online. The CI system automatically tests the code with multiple Python
versions (3.8, 3.9, and 3.10), ensuring that the code works as expected across
different platforms and Python versions. Local integration is implemented with
Tox (tox.ini) and Pyenv, and online integration is set up through GitHub
Workflows (.github/workflows/ci_cd.yml). The dual implementation of CI is due
to the advantages and disadvantages of these techniques. On the one hand,
GitHub Workflow integration tests compatibility on Mac, Windows, and Linux OS,
while Tox is specific to the local OS. On the other hand, GitHub workflows
have limited testing calls for private repositories (lifted for public
projects), while Tox can be used free of charge.

##  5 Code Documentation

The code documentation is based on Google-style docstrings, which are used
throughout the code to describe the purpose and behavior of classes, methods,
and functions. The MkDocs tool, a fast and simple static site generator,
produces the project documentation from these docstrings and additional
markdown files. The docstrings themselves are converted via the MkDocStrings
software extension. The markdown files used for MkDocs are located in the docs
directory, and the mkdocs.yml configuration file is placed at the project’s
root directory. This approach was specifically chosen to ensure that the
documentation always stays in sync with the code.

##  6 Usage

In this section, we discuss the output for the users and some capabilities of
the package. We shed light on the capability of this data with a real data set
commonly used in ML, namely the California housing data. This dataset provides
a realistic scenario for the models. However, the outcome variable, the
average house value, is not necessarily multi-modal. In the examples folder of
the package and the ‘Usage’ section of the documentation, you find a better
example for MDN where simulated data which contain multi-modal distributions
can be used to test specific aspects of the MDN models.

First, we load California housing data and convert it to the desired floating-
point format as shown in Figure 2. Then, we initialize and train an SVGP model
with 100 inducing points as depicted in Figure 3 and obtain the output for the
training. It can be observed that the model starts learning and then stops at
the 30th epoch. Please note here that we do not set a seed for Numpy and
PyTorch, but when applying this library, it is good practice to set the two
seeds.

⬇

>>> from sklearn.datasets import fetch_california_housing

>>> from sklearn.model_selection import train_test_split

>>> import numpy as np

>>> california = fetch_california_housing()

>>> X = np.array(California.data, dtype=np.float32)

>>> y = np.array(California.target, dtype=np.float32)

>>> X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
random_state=42)

Figure 2: Loading the California Housing dataset

⬇

>>> from uncertaintyplayground.trainers.svgp_trainer import SparseGPTrainer

>>> california_trainer_svgp = SparseGPTrainer(

X=X_train,

y=y_train,

num_inducing_points=100,

num_epochs=30,

batch_size=512,

lr=0.1,

patience=3)

>>> california_trainer_svgp.train()

Epoch 1/30, Weighted Loss: 1.731, Val MSE: 1.380385, Val R2: -0.000

Epoch 2/30, Weighted Loss: 1.676, Val MSE: 1.381376, Val R2: -0.001

...

Epoch 30/30, Weighted Loss: 1.514, Val MSE: 0.940440, Val R2: 0.319

Figure 3: Initializing and Training SVGPR

We can apply the same kind of pipeline for the MDN model as illustrated in
Figure 4. Aside from the function arguments, the outputs differ from those of
Figure 3. The first difference is that other outputs are not given aside from
the loss metric. It may be argued that inference with multiple Gaussian
distributions at every epoch may be undesirable. However, this feature will be
added in the next versions of the package. The second difference is that the
training stopped earlier than expected. This is due to the patience argument,
which triggers the early stopping at the default value of 10, meaning that
validation loss did not improve after the 11th epoch for another ten epochs.
Hence the training was stopped, and the best model was returned.

⬇

>>> from uncertaintyplayground.trainers.mdn_trainer import MDNTrainer

>>> california_trainer_mdn = MDNTrainer(

X=X_train,

y=y_train,

num_epochs=100,

lr=0.001,

dense1_units=50,

n_gaussians=10)

>>> california_trainer_mdn.train()

Epoch 1/100, Training Loss: 1.450, Validation Loss: 1.522

Epoch 2/100, Training Loss: 1.382, Validation Loss: 1.458

...

Epoch 11/100, Training Loss: 1.381, Validation Loss: 1.261

...

Epoch 21/100, Training Loss: 1.412, Validation Loss: 1.293

Early stopping after 21 epochs

Figure 4: Initializing and Training MDN

We must note that we do not show an example of predict_with_uncertainty in
this paper for brevity. These functions are also applied for generating
prediction plots and will be discussed briefly. Please note that this function
for both models produces the learned parameters of the model. More information
can be found on this topic in the documentation.

After the models have been trained, we visualize the prediction intervals. For
the SVGP model, use the function compare_distributions_svgpr and the function
plot_results_grid.

Next, we visualize the predictions of the MDN model using the function
compare_distributions_mdn and the function plot_results_grid. These workflows
can be found in Figure 5. The four produced plots are shown in Figure 6,
Figure 7 and Figure 8 and Figure 9.

⬇

>>> from uncertaintyplayground.predplot.mdn_predplot import
compare_distributions_mdn

>>> from uncertaintyplayground.predplot.svgp_predplot import
compare_distributions_svgpr

>>> from uncertaintyplayground.predplot.grid_predplot import plot_results_grid

>>> compare_distributions_svgpr(

trainer=california_trainer_svgp,

x_instance=X_test[900, :],

y_actual=y_test[900])

>>> plot_results_grid(

trainer=california_trainer_svgp,

compare_func=compare_distributions_svgpr,

X_test=X_test,

Y_test=y_test,

indices=[900, 500],

ncols=2)

>>> compare_distributions_mdn(

trainer=california_trainer_mdn,

x_instance=X_test[900, :],

y_actual=y_test[900])

>>> plot_results_grid(

trainer=california_trainer_mdn,

compare_func=compare_distributions_mdn,

X_test=X_test,

Y_test=y_test,

indices=[900, 500],

ncols=1)

Figure 5: Plotting results from both SVGPR and MDN models Figure 6: SVGPR
compare_distributions_svgpr for instance number 900 Figure 7: SVGPR
plot_results_grid on two instances Figure 8: MDN compare_distributions_mdn for
instance number 900 Figure 9: MDN plot_results_grid on two instances

The plots and their predictive distributions illustrated how to diagnose the
model’s performance and better understand the underlying data structure. We do
not discuss the quality of the model in this report. Nonetheless, it is
helpful to see that the MDN model may find skewed data that SVGPR does not
capture. Whether this model is better than SVGPR needs proper metrics and
validation and is beyond the scope of this report.

##  7 Further Development

The UncertaintyPlayground project, as a dynamic platform for uncertainty
estimation, has numerous avenues for future enhancements and expansions. The
opportunities for improvements are distributed across different elements of
the package:

  1. 1.

Flexible Model Architectures: For the MDN model, it would be beneficial to
make the neural network architecture modular, allowing for the incorporation
of different layers and activation functions. For the SVGP model, adding the
option to use different kernel functions could extend the model’s flexibility.

  2. 2.

Improved Noise Modelling: Introducing the capability to use different types of
noise in the MDN model could significantly improve the quality of uncertainty
estimates.

  3. 3.

Classification Capabilities: Both the MDN and SVGP models could be extended to
support binary and multi-class classification. This would involve modifying
the likelihood function and the performance metric. This would, however,
require extensive theoretical and empirical validation as MDNs and SVGPRs are
not traditionally used for classification tasks [14].

  4. 4.

Hardware Utilization: The package could benefit from implementing multi-GPU
support, which would allow for more efficient training of large models on
large datasets. Optimizing the parallel data loading process for maximized CPU
utilization could significantly improve overall performance.

  5. 5.

Improved Code Documentation: The addition of type hints to the docstrings
would offer better clarity and type checking, enhancing the readability and
maintainability of the codebase. Additionally, the documentation can benefit
from more examples.

  6. 6.

Benchmarking Performance: In the further iterations of this package, the
performance, both in terms of speed and accuracy of prediction, can be
measured against other models. For instance, one can compare our approach with
a traditional GPR for larger and smaller datasets since, as already discussed,
GPR has an algorithmic complexity of $O(N^{3})$ and does not scale well beyond
a few hundred observations.

##  8 Conclusion

This paper comprehensively reviews a newly developed Python package,
UncertainityPlayground, a library built for simplified uncertainty estimation
with PyTorch and GPyTorch. The two machine learning algorithms, MDN and SVGP
models, are implemented for fast and easy uncertainty estimation of continuous
outcomes. The results demonstrate that both models can model complex data
distributions and provide meaningful uncertainty estimates. The built-in
plotting functionalities allow users to study the inferred distribution of
given outcomes. The codebase is well-maintained, well-documented, and designed
for future development.

## References

  * [1] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., S
  * [2] Gardner, J.R., Pleiss, G., Bindel, D., Weinberger, K.Q. and Wilson, A.G., 2018. GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems.
  * [3] Rasmussen, C.E. and Williams, C.K.I., 2005. Gaussian Processes for Machine Learning. The MIT Press. doi:10.7551/mitpress/3206.001.0001
  * [4] Bishop, C.M., 1994. Mixture Density Networks. Aston University.
  * [5] M. Titsias.  Variational learning of inducing variables in sparse Gaussian processes.  In D. van Dyk and M. Welling, editors, _Proceedings of the Twelfth International Workshop on Artificial Intelligence and Statistics_ , volume 5, pages 567–574, Clearwater Beach, FL, 16-18 April 2009. JMLR W&CP 5\.
  * [6] M. Kuss and C. E. Rasmussen.  Assessing approximate inference for binary Gaussian process classification.  _The Journal of Machine Learning Research_ , 6:1679–1704, 2005.
  * [7] H. Nickisch and C. E. Rasmussen.  Approximations for binary Gaussian process classification.  _Journal of Machine Learning Research_ , 9:2035–2078, 2008\.
  * [8] J. Hensman, N. Fusi, and N. D. Lawrence.  Gaussian processes for big data.  In A. Nicholson and P. Smyth, editors, _Uncertainty in Artificial Intelligence_ , volume 29. AUAI Press, 2013.
  * [9] David Burt, Carl Edward Rasmussen, and Mark Van Der Wilk.  Rates of Convergence for Sparse Variational Gaussian Process Regression.  In _Proceedings of the 36th International Conference on Machine Learning_ , pages 862–871. PMLR, 2019.
  * [10] Stéfan van der Walt, S. Chris Colbert and Gaël Varoquaux.  The NumPy Array: A Structure for Efficient Numerical Computation.  _Computing in Science & Engineering_, 13(2):22-30, 2011.
  * [11] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12:2825-2830, 2011.
  * [12] John D. Hunter.  Matplotlib: A 2D Graphics Environment.  _Computing in Science & Engineering_, 9(3):90-95, 2007.
  * [13] Michael Waskom.  Seaborn: statistical data visualization.  https://seaborn.pydata.org/, 2022.
  * [14] James Hensman, Alex Matthews, and Zoubin Ghahramani.  Scalable Variational Gaussian Process Classification.  _arXiv preprint arXiv:1411.2005_ , 2015.


# An Analysis of X-Ray Hardness Ratios Between Asynchronous and Non-
Asynchronous Polars

Eric Masington  Student  Department of Physics and Astronomy, Texas Tech
University, Lubbock, TX 79409, USA  Thomas J. Maccarone  Department of Physics
and Astronomy, Texas Tech University, Lubbock, TX 79409, USA  Liliana Rivera
Sandoval  Department of Physics and Astronomy, Texas Tech University, Lubbock,
TX 79409, USA  Department of Physics, University of Alberta, Edmonton, Alberta
T6G 2R3, CA  Craig Heinke  Department of Physics, University of Alberta,
Edmonton, Alberta T6G 2R3, CA  Arash Bahramian  Department of Physics and
Astronomy, Curtin University, Perth, Western Australia 6845, AU  Aarran W.
Shaw  Department of Physics, University of Nevada, Reno, NV 89557, USA

###### Abstract

The subclass of magnetic Cataclysmic Variables (CV), known as asynchronous
polars, are still relatively poorly understood. An asynchronous polar is a
polar in which the spin period of the white dwarf is either shorter or longer
than the binary orbital period (typically within a few percent). The
asynchronous polars have been disproportionately detected in soft gamma-ray
observations, leading us to consider the possibility that they have
intrinsically harder X-ray spectra. We compared standard and asynchronous
polars in order to examine the relationship between a CV’s synchronization
status and its spectral shape. Using the entire sample of asynchronous polars,
we find that the asynchronous polars may, indeed, have harder spectra, but
that the result is not statistically significant.

(stars:) novae, cataclysmic variables — X-rays: binaries — catalogs

##  1 Introduction

One of the first results on accreting white dwarfs with the International
Gamma-Ray Astrophysical Laboratory (INTEGRAL) was that the asynchronous polars
represented a disproportionate fraction of its detected cataclysmic variables
(Barlow et al., 2016). Asynchronous polars are magnetically accreting white
dwarfs with deviations between the orbital and spin period (unlike standard
polars) and with streamlike accretion rather than accretion disks (unlike the
intermediate polars). The difference between the spin and orbital period in
the asynchronous polars is typically about 5% or less. It is unclear whether
INTEGRAL preferentially detected these objects because the X-ray and
$\gamma-$ray spectra of the asynchronous polars are different from those of
the standard polars, or merely because they tend to be more luminous, hence at
higher fluxes within the well-understood samples. Here, we test whether the
spectral indices of these sources in the soft $\gamma-$ray band alone, and
between X-ray and $\gamma$-ray, are systematically harder for the asynchronous
polars than for the standard polars.

The asynchronous polars are often suggested to have been driven out of
synchronization by classical novae, which can affect both the orbital and spin
periods of cataclysmic variables, motivated by the association between one of
the asynchronous polars, V1500 Cyg with a classical nova in 1975(Campbell &
Schwope, 1999). Searches for additional nova shells around other asynchronous
polars have not yielded any new evidence for the nova hypothesis (Pagnotta &
Zurek, 2016), but it remains a viable one, as nova shells may have lifetimes
shorter than the synchronization timescales of the asynchronous polars.
Because the sample sizes of the asynchronous polars are quite small, and
rather long, well-sampled light curves are needed to identify that there are
two separate, but similar periods in the light curves, it is worth exploring
new methods that might work to find new members of the class, and the INTEGRAL
discoveries of these objects suggest that gamma-ray surveys might be an
interesting approach. With this in mind, we undertake an exploration of
whether the high energy spectra of asynchronous polars are fundamentally
different from those of the standard polars.

##  2 Data used

We obtained a set of cataclysmic variables from the Ritter and Kolb catalog,
update 7.24 (Ritter, Kolb, 2003), hard X-ray data from the 2018 Swift-BAT
105-month All-Sky Survey catalogue (Oh et al., 2018) and soft X-ray data from
the ROSAT All-Sky Survey (Boller et al., 2016). Matching all three catalogues
with a maximum separation of 3 arc minutes (typical for Swift-BAT for faint
sources) yields 10 objects, including 4 asynchronous polars and 6 standard
polars. We define an asynchronous polar to be a system with a spin period
within 5% of the orbital period. The asynchronous polars are BY Cam, CD Ind,
V1432 Aql, and the recently identified IGR J19552+0044 111We consider a 5
sigma upper limit for V1500 Cyg and Rx J0838.7-2827 based on the survey depth,
but do not consider upper limits for the much larger class of synchronous
polars. (Tovmassian et al., 2017). V1500 Cyg and RX J0838.7-2827 were added to
the analysis to provide for a more complete sample. The synchronous polars
included are AM Her, Swift J231930.4+261517, 1RXS 145341.1-552146, IW Eri,
V2301 Oph, and V834 Cen.

Figure 1: Shown here is the ratio between the flux (Swift-BAT 15-150 keV band)
and count rate (ROSAT 0.1-2.0 keV band) plotted against luminosity (calculated
using Gaia distances and Swift flux values). As can be seen, no significant
relationship can be established between whether or not a polar is asynchronous
and the hardness of its flux ratio. It must be considered that the ROSAT and
Swift data are taken non-simultaneously, so because the ROSAT count numbers
are generally too small for spectral fitting, they cannot be reliably
converted into fluxes. Thus, this means that only very strong trends could
have been detected using this combination of data. Such trends are not
present, but the method would not be particularly sensitive to subtle
systematic variations. Also shown are the gamma-ray spectral indices. We can
see that the asynchronous polars seem to be harder (lower spectral indices),
but more data is needed to confirm this. (Gaia Collaboration et al., 2018)

##  3 Analysis and conclusions

First, we looked at the ratio of hard X-rays from Swift-BAT (15-150 keV), to
soft X-rays from ROSAT (0.1-2.0 keV). This comparison is done between a count
rate for ROSAT and a flux for Swift BAT because the standard ROSAT data
include only a count rate, and the standard BAT data include only a flux. The
fluxes do, thus, show some model dependence, but since the spectra are all
steep power laws, with photon index greater than 2.0, in all cases, the BAT
flux is dominated by counts near the lower end of the band, and this
comparison is nearly equivalent to a count rate-to-count rate comparison.
However, the ROSAT and Swift data are taken non-simultaneously, and because
the ROSAT count numbers are generally too small for spectral fitting, they
cannot be reliably converted into fluxes. This thus means that only very
strong trends could have been detected using this combination of data. Such
trends are not present, but the method would not be particularly sensitive to
subtle systematic variations. No trend is found in this ratio between the two
classes of polars. Also, there exists much uncertainty in the soft X-ray flux
values for V1500 Cyg and RX J0838.7-2827, given that they are not in the Swift
catalog. Because the plot shows strong scatter between the Swift and ROSAT
data, and the ROSAT data in most cases are insufficient for detailed spectral
analysis, we simply leave this plot as a ratio of a flux to a count rate.

Next, we consider the spectra within the $\gamma-$ray band alone. Within the
Swift band, the mean spectral index is 2.73 for the asynchronous polars and
3.34 for the standard polars. We apply the Anderson-Darling test to the
distributions of spectral indices. This is a cumulative statistic, similar to
the Kolomogorov-Smirnov test, but with greater diagnostic power in cases where
the differences are strongest near the edges of the distributions, and at
least equal power in all cases. The Anderson-Darling (AD) test statistic here
is 2.28 (computed using https://www.real-statistics.com/non-parametric-
tests/goodness-of-fit-tests/two-sample-anderson-darling-test/). For this
sample size, the critical value of the AD test statistic is 3.38 for a 99%
confidence level detection of a difference. The asynchronous polars do show a
different $\gamma$-ray spectral index at the 95% confidence level. Since this
is a marginally significant difference, we expect that a larger sample of
objects would have a reasonable probability of establishing a difference.
Unfortunately, doing so will require finding new asynchronous polars, as we
have already investigated the properties of the whole sample, with only V1500
Cyg and RX J0838.7-2827 undetected in the Swift-BAT data. If more asynchronous
polars can be found from optical or soft X-ray searches, NuSTAR would easily
be capable of measuring $\gamma$-ray spectral indices for objects much fainter
than the BAT survey can, so searches for more asynchronous polars would be
well-motivated.

## References

  * Barlow et al. (2016) Barlow, E. J., Knigge, C., Bird, A. J., et al. 2006, MNRAS, 372, 224
  * Boller et al. (2016) Boller T., Freyberg, M. J., Trumper, J., Haberl, F., Voges, W., Nandra, K. 2016, A&A, 588, A103
  * Campbell & Schwope (1999) Campbell, C. G. & Schwope, A. D. 1999, A&A, 343, 132
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Oh et al. (2018) Oh, K., Koss, M., Markwardt, C. B., Schawinski K., Baumgartner W. H., Barthelmy, S. D., Cenko, S. B., Gehrels, N., Mushotzky, R., Petulante A., Ricci, C., Lien, A., & Trakhtenbrot, B. 2018, ApJS, 235, 1
  * Pagnotta & Zurek (2016) Pagnotta, A. & Zurek, D. 2016, MNRAS, 458, 1833. doi:10.1093/mnras/stw424
  * Ritter, Kolb (2003) Ritter H., Kolb U. 2003, A&A, 404, 301
  * Tovmassian et al. (2017) Tovmassian, G., González-Buitrago, D., Thorstensen, J., Kotze, E., Breytenbach, H., Schwope, A., Bernardini, F., Zharikov, S. V., Hernandez, M. S., Buckley, D. A. H., de Miguel, E., Hambsch, F. J., Myers, G., Goff, W., Cejudo, D., Starkey, D., Campbell, T., Ulowetz, J., Stein, W., Nelson, P., Reichart, D. E., Haislip, J. B., Ivarsen, K. M., LaCluyze, A. P., Moore, J. P., Miroshnichenko, A. S. 2017, A&A, 608, A36


# Electron Preacceleration at Weak Quasi-perpendicular Intracluster Shocks:
Effects of Preexisting Nonthermal Electrons

Ji-Hoon Ha Department of Physics, College of Natural Sciences UNIST, Ulsan
44919, Korea  Dongsu Ryu Department of Physics, College of Natural Sciences
UNIST, Ulsan 44919, Korea  Hyesung Kang Department of Earth Sciences, Pusan
National University, Busan 46241, Korea  Sunjung Kim Department of Physics,
College of Natural Sciences UNIST, Ulsan 44919, Korea Ji-Hoon Ha, Dongsu Ryu
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Radio relics in the outskirts of galaxy clusters imply the diffusive shock
acceleration (DSA) of electrons at merger-driven shocks with Mach number
$M_{s}\lesssim 3-4$ in the intracluster medium (ICM). Recent studies have
suggested that electron preacceleration and injection, prerequisite steps for
DSA, could occur at supercritical shocks with $M_{s}\gtrsim 2.3$ in the ICM,
thanks to the generation of multiscale waves by microinstabilities such as the
Alfvén ion cyclotron (AIC) instability, the electron firehose instability
(EFI), and the whistler instability (WI). On the other hand, some relics are
observed to have subcritical shocks with $M_{s}\lesssim 2.3$, leaving DSA at
such weak shocks as an outstanding problem. Reacceleration of preexisting
nonthermal electrons has been contemplated as one of possible solutions for
that puzzle. To explore this idea, we perform Particle-in-Cell (PIC)
simulations for weak quasi-perpendicular shocks in high-$\beta$ ($\beta=P_{\rm
gas}/P_{B}$) plasmas with power-law suprathermal electrons in addition to
Maxwellian thermal electrons. We find that suprathermal electrons enhance the
excitation of electron-scale waves via the EFI and WI. However, they do not
affect the ion reflection and the ensuing generation of ion-scale waves via
the AIC instability. The presence of ion-scale waves is the key for the
preacceleration of electrons up to the injection momentum, thus the shock
criticality condition for electron injection to DSA is preserved. Based on the
results, we conclude that preexisting nonthermal electrons in the preshock
region alone would not resolve the issue of electron preacceleration at
subcritical ICM shocks.

acceleration of particles – cosmic rays – galaxies: clusters: general –
methods: numerical – shock waves

††journal: The Astrophysical Journal

##  1 Introduction

Structure formation shocks are induced by supersonic flow motions during the
hierarchical clustering of the large-scale structures of the universe
including galaxy clusters (e.g., Miniati et al., 2000; Ryu et al., 2003;
Pfrommer et al., 2006; Skillman et al., 2008; Vazza et al., 2009; Hong et al.,
2014; Schaal & Springel, 2015). In particular, major mergers of clusters drive
shocks with low sonic Mach numbers, $M_{s}\lesssim 3-4$ (e.g., Ha et al.,
2018), in the hot intracluster medium (ICM) of high plasma beta, $\beta\sim
100$ (e.g., Ryu et al., 2008; Porter et al., 2015). Here, $\beta=P_{\rm
gas}/P_{B}$, and $P_{\rm gas}$ and $P_{B}$ are the gas and magnetic pressures,
respectively. Such merger-driven shocks are energetic enough to be detected in
radio and X-ray observations (e.g., Markevitch & Vikhlinin, 2007; van Weeren
et al., 2010; Brüggen et al., 2012; Brunetti & Jones, 2014).

In the so-called radio relics, the diffuse radio emission has been interpreted
as synchrotron radiation from the cosmic-ray (CR) electrons produced by
diffusive shock acceleration (DSA; e.g., Bell, 1978; Blandford & Ostriker,
1978; Drury, 1983) at merger-driven shocks (see, e.g., van Weeren et al.,
2019, for a review). Those shocks in the ICM are collisionless as in other
astrophysical environments. With the shock obliquity angle between the shock
normal and the background magnetic field, $\theta_{Bn}$, collisionless shocks
are classified into “quasi-perpendicular” ($Q_{\perp}$, hereafter) if
$\theta_{Bn}\gtrsim 45^{\circ}$ and “quasi-parallel” ($Q_{\parallel}$,
hereafter) if $\theta_{Bn}\lesssim 45^{\circ}$. It is known from in-situ
observations of Earth’s bow shocks that electrons are energized mostly at
$Q_{\perp}$ shocks (e.g., Gosling et al., 1989). So the acceleration of CR
electrons is expected to be effective primarily at the $Q_{\perp}$-population
of merger shocks. On the other hand, using a high-resolution, long-
integration, one-dimensional Particle-in-Cell (PIC) simulation, Park et al.
(2015) showed that, in “strong” $Q_{\parallel}$-shocks ($M_{s}\approx 40$ or
so), electrons could be accelerated via the Fermi acceleration mediated by
ion-induced non-resonant Bell waves.

One of the key issues for electron acceleration at collisionless shocks is the
so-called injection to DSA. This is because thermal electrons have the
gyroradius much smaller than the shock thickness, and hence need to be
“preaccelerated” to the “injection momentum”, $p_{\rm inj}$, in order to
diffuse across the shock transition and participate in the first-order Fermi
process of DSA (e.g., Kang et al., 2002; Caprioli et al., 2015). Here, the
injection momentum is $p_{\rm inj}\gtrsim{\rm a~{}few}\times p_{i,{\rm th}}$,
where $p_{i,{\rm th}}=\sqrt{2m_{i}k_{B}T_{i2}}$, $m_{i}$ is the ion mass,
$k_{B}$ is the Boltzmann constant, and $T_{i2}$ is the downstream ion
temperature. Hereafter, the subscripts 1 and 2 denote the quantities in the
shock upstream and downstream regions, respectively.

To study the electron preacceleration at ICM shocks, PIC simulations have been
performed for $Q_{\perp}$-shocks with $M_{s}\lesssim$ several in plasmas of
$\beta\approx 20-100$ (e.g., Matsukiyo et al., 2011; Matsukiyo & Matsumoto,
2015; Guo et al., 2014a, b; Kang et al., 2019; Niemiec et al., 2019; Ha et
al., 2021; Kobzar et al., 2021). In those simulations, two mechanisms
involving wave-particle interactions have been identified for electron
preacceleration: (1) multiple cycles of the shock drift acceleration (SDA)
mediated by upstream self-excited waves, that is, the Fermi-like acceleration
in the preshock region, and (2) the stochastic SDA (SSDA) mediated by waves of
multiscales, ranging from electron to ion scales, in the shock transition
zone. While the SDA of electrons, even after multiple cycles, would not be
sufficient for electron injection, the SSDA could energize electrons all the
way up to $p_{\rm inj}$.

The plasma kinetic waves involved in these preacceleration processes are
produced by a variety of microinstabilities. In collisionless
$Q_{\perp}$-shocks, the instabilities are governed by the ions and electrons
reflected or compressed at the shock ramp (see, e.g., Kim et al., 2020, 2021,
for more details). In the shock upstream, the reflected electrons induce the
electron temperature anisotropy, $T_{e\parallel}/T_{e\perp}>1$, which in turn
excites the electron firehose instability (EFI). Scattering off the EFI-
induced waves moderates a Fermi-like acceleration between the shock ramp and
the upstream region (e.g. Guo et al., 2014a, b; Kang et al., 2019). In the
shock transition zone, on the other hand, the ion temperature anisotropy is
induced by the shock-reflected ions that are advected downstream, while the
electron temperature anisotropy is induced by magnetic field compression at
the shock ramp (e.g. Guo et al., 2017; Katou & Amano, 2019). Then, those
ion/electron temperature anisotropies could trigger instabilities of the
following three types: (1) the Alfvén ion cyclotron (AIC) instability due to
$T_{i\perp}/T_{i\parallel}>1$, (2) the whistler instability (WI) due to
$T_{e\perp}/T_{e\parallel}>1$, and (3) the ion and electron mirror
instabilities due to $T_{i\perp}/T_{i\parallel}>1$ and
$T_{e\perp}/T_{e\parallel}>1$, respectively.

The SSDA in the shock transition zone, which would be necessary for electron
injection to DSA, relies on the extended confinement of electrons within the
zone by scattering off multiscale waves and the resulting gradient-drift along
the motional electric field, and hence requires the effective excitation of
the AIC instability (Niemiec et al., 2019; Kobzar et al., 2021). Ha et al.
(2021, Paper I, hereafter) showed that the reflection of incoming ions and the
ensuing excitation of the AIC instability is expected to be effective only at
$Q_{\perp}$-shocks with $M_{s}\gtrsim M_{\rm crit}^{*}\approx 2.3$ in the ICM,
implying that DSA would work only at those supercritical shocks.

The properties of merger-driven shocks in radio relics have been studied
through observations. For instance, $M_{s}$ of the so-called “radio relic
shocks” has been inferred either from the radio spectral index in radio
observations (e.g., van Weeren et al., 2010, 2019) or using the
discontinuities in temperature or surface brightness in X-ray observations
(e.g., Markevitch et al., 2002; Markevitch & Vikhlinin, 2007). In some radio
relics, the Mach numbers inferred from radio observations is larger than those
estimated from X-ray observations, i.e., $M_{\rm radio}\gtrsim M_{\rm X}$
(e.g., Akamatsu & Kawahara, 2013; van Weeren et al., 2019; Wittor et al.,
2021). It has been argued that the surfaces of merger shocks in the turbulent
ICM are not uniform and consist of parts with different $M_{s}$, and radio and
X-ray observations preferentially pick up the parts of high and low Mach
numbers, respectively (e.g., Hong et al., 2015; Botteon et al., 2020; Wittor
et al., 2021). This implies that there are intrinsic and inevitable
uncertainties in $M_{s}$ estimated from observations of merger shocks.
Nevertheless, some radio relics seem to have shocks with $M_{s}$ close to or
even smaller than $M_{\rm crit}^{*}\approx 2.3$, and examples include well-
studied radio relics, such as those in ZwCL 0008.8+5215 (e.g., Kierdorf et
al., 2017) and Abell 521 (e.g., Macario et al., 2013). Hence, to explain the
origin of radio relics, additional ingredients, other than the preacceleration
and injection due to self-excited waves, would need to be considered.

As a possible solution for radio relics with low $M_{s}$ shocks, the so-called
“reacceleration model” was suggested, where preexisting relativistic fossil
electrons are injected to DSA and further accelerated (e.g., Kang et al.,
2012; Pinzke et al., 2013; Kang, 2016a, b; Kang & Ryu, 2016). Those fossil
electrons could enhance the acceleration efficiency, especially at weak shocks
with $M_{s}\lesssim 3$ (see, e.g., Kang & Ryu, 2011). The presence of
relativistic electron populations in the ICM is expected, possibly being
supplied by radio jets of AGNs or produced through previous episodes of
shock/turbulence acceleration (e.g., Kang & Ryu, 2015, 2016, and references
therein). As a matter of fact, some radio relics are observed to be associated
with nearby radio galaxies (e.g., van Weeren et al., 2017), and low-frequency
observations indicate the possibility of preexisting nonthermal electrons
(e.g., Mandal et al., 2020). However, the reacceleration scenario at weak
shocks would operate, only if preexisting nonthermal electrons with $p\gtrsim
p_{\rm inj}$ could be scattered diffusively across the shock transition. It
requires either the triggering of microinstabilities and the self-generation
of scattering waves or the presence of preexisting kinetic waves in the
turbulent ICM.

In this paper, we revisit our previous work of PIC simulations for weak
$Q_{\perp}$ ICM shocks, described in Paper I, by including power-law
nonthermal electrons in addition to Maxwellian thermal electrons in the
upstream plasma. Due to the limitation of computational capacities, the power-
law component extends only up to the Lorenz factor of $\gamma\approx$ several,
representing a suprathermal population, rather than the bonafide nonthermal
population that continues to $\gamma\gg 1$. We study the effects of the
preexisting suprathermal electrons on kinetic plasma processes in high-$\beta$
$Q_{\perp}$-shocks. In particular, we examine whether the suprathermal
electrons could enhance the microinstabilities that generate multiscale waves,
leading to the preacceleration of electrons up to $p_{\rm inj}$, even in
subcritical shocks with $M_{s}\lesssim 2.3$.

The paper is organized as follows. In the next section, simulations are
described. The results are presented in Section 3. Summary and discussion
follow in Section 4.

Table 1: Model Parameters for Simulations  Model Name | $M_{s}$ | $M_{A}$ | $u_{0}/c$ | $\theta_{Bn}$ | $\beta$ | $m_{i}/m_{e}$ | $\alpha$ | $n_{p}/n_{0}$ | $\gamma_{\rm cut}$ | $L_{x}[\lambda_{se}]$ | $L_{y}[\lambda_{se}]$ | $t_{\rm end}[\Omega_{ci}^{-1}]$
---|---|---|---|---|---|---|---|---|---|---|---|---
M2.0 | 2.0 | 8.2 | 0.038 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 5 | 2000 | 200 | 22
M2.3 | 2.3 | 9.4 | 0.046 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 5 | 2000 | 200 | 22
M2.5 | 2.5 | 10.2 | 0.053 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 5 | 2000 | 200 | 22
M3.0 | 3.0 | 12.3 | 0.068 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 5 | 2000 | 200 | 22
M2.0-$\alpha$4 | 2.0 | 8.2 | 0.038 | $63^{\circ}$ | 20 | 50 | 4 | 0.01 | 5 | 2000 | 200 | 22
M2.0-$\alpha$4.5 | 2.0 | 8.2 | 0.038 | $63^{\circ}$ | 20 | 50 | 4.5 | 0.01 | 5 | 2000 | 200 | 22
M2.0-np0.1 | 2.0 | 8.2 | 0.038 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.1 | 5 | 2000 | 200 | 22
M2.0-$\gamma$cut10 | 2.0 | 8.2 | 0.038 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 10 | 2000 | 200 | 22
M2.0-m25 | 2.0 | 8.2 | 0.053 | $63^{\circ}$ | 20 | 25 | 4.2 | 0.01 | 5 | 2000 | 140 | 22
M2.0-$\beta$50-m25 | 2.0 | 12.9 | 0.053 | $63^{\circ}$ | 50 | 25 | 4.2 | 0.01 | 5 | 2000 | 225 | 22
M3.0-$\alpha$4 | 3.0 | 12.3 | 0.068 | $63^{\circ}$ | 20 | 50 | 4 | 0.01 | 5 | 2000 | 200 | 22
M3.0-$\alpha$4.5 | 3.0 | 12.3 | 0.068 | $63^{\circ}$ | 20 | 50 | 4.5 | 0.01 | 5 | 2000 | 200 | 22
M3.0-np0.1 | 3.0 | 12.3 | 0.068 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.1 | 5 | 2000 | 200 | 22
M3.0-$\gamma$cut10 | 3.0 | 12.3 | 0.068 | $63^{\circ}$ | 20 | 50 | 4.2 | 0.01 | 10 | 2000 | 200 | 22
M3.0-m25 | 3.0 | 12.3 | 0.096 | $63^{\circ}$ | 20 | 25 | 4.2 | 0.01 | 5 | 2000 | 140 | 22
M3.0-$\beta$50-m25 | 3.0 | 19.4 | 0.096 | $63^{\circ}$ | 50 | 25 | 4.2 | 0.01 | 5 | 2000 | 225 | 22

##  2 Numerics

Two-dimensional (2D) PIC simulations have been performed using TRISTAN - MP, a
parallel electromagnetic PIC code (Buneman, 1993; Spitkovsky, 2005). Here, we
brief simulations, while further details can be found in Paper I. In the
downstream rest frame, the background plasma composed with ions and electrons
moves with the velocity $-u_{0}\hat{x}$, and a collisionless shock propagating
to the $+x$-direction is generated via interaction with reflected particles.
The simulation domain is in the $x$-$y$ plane, and a uniform background
magnetic field, ${\bf{B_{0}}}$, lies in the plane. We adopt the plasma
parameters relevant for the ICM, the ion and electron temperatures
$k_{B}T=k_{B}T_{i}=k_{B}T_{e}=0.0168~{}m_{e}c^{2}=8.6$ keV (or
$T_{i}=T_{e}=10^{8}$ K), the ion and electron number densities
$n_{0}=n_{i}=n_{e}=10^{-4}{\rm cm}^{-3}$, and the plasma beta $\beta=16\pi
n_{0}k_{B}T/B_{0}^{2}=20-50$. For given $\beta$ and $k_{B}T$, the sonic and
Alfvén Mach numbers are calculated as $M_{s}=u_{\rm sh}/\sqrt{2\Gamma
k_{B}T_{i}/m_{i}}$ and $M_{A}=u_{\rm sh}/(B_{0}/\sqrt{4\pi
n_{0}m_{i}})=M_{s}\sqrt{\Gamma\beta/2}$, respectively. Here, $\Gamma=5/3$ is
the adiabatic index. The shock velocity, $u_{\rm sh}=ru_{0}/(r-1)$, is
estimated using the shock compression ratio,
$r=(\Gamma+1)/(\Gamma-1+2/M_{s}^{2})$. The sonic Mach number of simulated
shocks is in the range of $M_{s}\sim 2-3$, covering subcritical and
supercritical shocks, and $\theta_{Bn}=63^{\circ}$ is chosen to emulate
$Q_{\perp}$-shocks.

The spatial resolution is $\Delta x=\Delta y=0.1\lambda_{se}$, and the time
step is $\Delta t=0.045\omega_{pe}^{-1}$, where $\omega_{pe}=\sqrt{4\pi
e^{2}n_{0}/m_{e}}$ and $\lambda_{se}=c/\omega_{pe}$ are the electron plasma
frequency and the electron skin depth, respectively. Reduced mass ratios,
$m_{i}/m_{e}=25-50$, are employed to alleviate the demanding computational
requirement. We note that $\Delta x$, $\Delta y$ and $\Delta t$ are expressed
in units of the plasma skin depth, $\lambda_{se}$, and the electron plasma
oscillation period, $\omega_{pe}^{-1}$. On the other hand, the shock changes
on the time scale of the gyroperiod of incoming ions,
$\Omega_{ci}^{-1}=m_{i}c/eB_{0}\propto m_{i}\sqrt{\beta}$, and the shock
transition develops on the length scale of their gyroradius, $r_{L,i}\equiv
u_{0}/\Omega_{ci}\approx
57\lambda_{se}\cdot(M_{s}/3)\sqrt{\beta/20}\sqrt{(m_{i}/m_{e})/50}$. Hence, we
use $\Omega_{ci}$ and $r_{L,i}$, when we describe the evolution and structure
of the shock.

Figure 1: Energy spectra of electrons in the far upstream region
$10\leq(x-x_{\rm sh})/r_{L,i}\leq 11$ (solid lines) and in the initial
incoming plasmas (dashed lines) for two models M3.0 and M3.0-np0.1 (see Table
1). Here, $x_{\rm sh}$ is the shock position, and $r_{L,i}\approx
57\lambda_{se}$. The energy spectrum, $dN/d\gamma$, is related to the momentum
spectrum, $f(p)\propto p^{-\alpha}$, as $4\pi f(p)p^{2}dp/dE\propto
dN/d\gamma\propto(\gamma-1)^{-s}$, where $s=\alpha-2$.

The main difference between this work and Paper I lies in the inclusion of
preexisting power-law electrons (PPEs, hereafter). In this work, we add a
nonthermal electron population of “isotropic” power-law distribution,
$f(p)\propto p^{-\alpha}$, to the background thermal plasma, reducing the
amount of thermal electrons accordingly. The power-law slope we consider
covers the range of $\alpha=4-4.5$, and the number fraction of PPEs ranges
$n_{p}/n_{0}=0.01-0.1$. For the fiducial model, $\alpha=4.2$ is adopted,
assuming that the preexisting nonthermal electrons are produced by
relativistic shocks such as those in radio jets (e.g., Kirk et al., 2000;
Achterberg et al., 2001). The minimum and maximum (or cutoff) Lorenz factors
of the nonthermal component are $\gamma_{\rm min}=1.15$ (corresponding to the
velocity $\sim 0.5c$) and $\gamma_{\rm cut}=5-10$, respectively. Figure 1
shows the energy spectra of electrons measured in the far upstream region of
$10\leq(x-x_{\rm sh})/r_{L,i}\leq 11$ at $\Omega_{ci}t\sim 3$, and compares
them with those of the initial incoming plasmas. Here, $x_{\rm sh}$ denotes
the position of the shock. We see a smoothing in the upstream energy spectra,
resulting in the continuous merging of the power-law component to the
Maxwellian component. Although it is not further discussed here, the whistler
mode is responsible for it. Even with the smoothing, the power-law
distribution is pretty well maintained.

In PIC simulations, the limited number of particles induces numerical
dissipations (e.g., Melzani et al., 2013), and hence a sufficiently large
number of particles is required to maintain the power-law distribution of
electrons for a sufficiently long time. In this work, 120 particles per cell
(60 per species) are placed, which are much larger than 32 particles per cell
(16 per species) used in Paper I. Even with this number of particles,
$\gamma_{\rm cut}=5-10$ is the maximum value that we can accommodate. We point
that this value of $\gamma_{\rm cut}$ is set by the computational limitation,
rather than the physical argument. With $\gamma_{\rm cut}=5-10$, the PPE
component should be regarded in effect as a suprathermal population.

The model parameters of our simulations are given in Table 1. Models with
different $M_{s}$ are named with the combination of the letter ‘M’ and the
sonic Mach number (for instance, the M3.0 model has $M_{s}=3$). The four
models in the top group are the fiducial models, which have $\beta=20$,
$m_{i}/m_{e}=50$, $\alpha=4.2$, $n_{p}/n_{0}=0.01$, and $\gamma_{\rm cut}=5$.
Models with parameters different from those of the fiducial models have the
names that are appended by a character for the specific parameter and its
value; for example, the M3.0-np0.1 model has $n_{p}/n_{0}=0.1$, while the
M3.0-$\alpha 4$ has the power-law slope $\alpha=4$. For the box size, $L_{x}$
and $L_{y}$, the end time of simulations, $t_{\rm end}$, and $\beta$, we adopt
the values somewhat smaller than those of Paper I to compensate the longer
computational time due to the larger number of particles. Yet, the adopted
values of these parameters should be pertinent enough to capture the main
results of this work. In the next section, the fiducial models in the current
work will be compared with the models with $\beta=20$ in Paper I (without
PPEs).

Figure 2: Fraction of reflected electrons (panel (a)) and average energy
gained by reflected electrons in a single cycle of SDA (panel (b)), based on
the relativistic SDA description (Guo et al., 2014a; Kang et al., 2019).
Electron temperature anisotropy, $T_{e\parallel}/T_{e\perp}$, estimated in the
immediate upstream region, $0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$ (panel (c)),
and $T_{e\perp}/T_{e\parallel}$ estimated in the immediate downstream region,
$-1\leq(x-x_{\rm sh})/r_{L,i}\leq 0$ (panel (d)). Instability parameters for
the EFI, $I_{\rm EFI}$ (panel (e)), and for the WI, $I_{\rm WI}$ (panel (f)).
The red solid lines draw the fiducial models with PPEs, while the black solid
lines draw the corresponding models without PPEs from Paper I, at
$\Omega_{ci}t\sim 22$. The red dashed lines in panels (a) and (b) show the
quantities estimated only with PPEs, that is, the fraction of reflected PPEs
and the average energy gained by reflected PPEs in a single cycle of SDA.

Figure 3: Ion number density (top panels), $n_{i}$, normalized to the upstream
ion number density, $n_{0}$, and $z$-magnetic field (bottom panels), $B_{z}$,
normalized to the upstream magnetic field, $B_{0}$, in the region of
$-60\leq(x-x_{\rm sh})/\lambda_{se}\leq 140$ around the shock at
$\Omega_{ci}t\sim 22$. Here, $x_{\rm sh}$ is the shock position. The fiducial
M2.0 and M3.0 models with $\alpha=4.2$ and $n_{p}/n_{0}=0.01$ (panels (a),
(e), (c), (g)) are compared to the corresponding models without PPEs from
Paper I (panels (b), (f), (d), (h)).

Figure 4: Magnetic field power spectra,
$P_{B_{z}}(k_{y})\propto(k_{y}/2\pi)(\delta B_{z}(k_{y})^{2}/B_{0}^{2})$, in
the immediate upstream region, $0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$ (top
panels), and the immediate downstream region, $-1\leq(x-x_{\rm
sh})/r_{L,i}\leq 0$ (bottom panels), at $\Omega_{ci}t\sim 22$, for the models
with PPEs (left panels) and without PPEs (right panels). Here, $x_{\rm sh}$ is
the shock position, and $r_{L,i}\approx 57\lambda_{se}\cdot(M_{s}/3)$. The
magenta and gray boxes denote the ranges of the gyroradius of electrons with
$p_{\rm inj}$ and $p_{\rm spt}\equiv\sqrt{m_{e}/m_{i}}~{}p_{\rm inj}$,
respectively, in the upstream and downstream regions.

##  3 Results

###  3.1 Particle Reflection and Temperature Anisotropies

According to the description of relativistic SDA theory for oblique shocks
(e.g., Guo et al., 2014a; Kang et al., 2019), the amount of reflected
electrons is expected to increase, if PPEs with higher velocities are added to
the Maxwellian distribution. We first examine its consequences and effects on
SDA and instabilities in Figure 2. In Figure 2(a,b), the fraction of reflected
electrons and the average energy gained by those reflected electrons in a
single SDA cycle are shown. Indeed, the fraction of reflected PPEs (red dashed
line) is larger than that of reflected thermal electrons (black solid line),
so is the fraction of total reflected electrons in the models with PPEs (red
solid line). The average energy gained by reflected PPEs in a single SDA cycle
is larger than that of thermal electrons as well, because electrons with
larger $\gamma$ drift longer distances along the shock surface. Hence, with
PPEs, the electron temperature anisotropies, $T_{e\parallel}/T_{e\perp}$,
estimated in the shock upstream, $0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$, increase
(Figure 2(c)), and then based on the instability condition for the EFI,
$I_{\rm EFI}\approx 1-(T_{\rm e\parallel}/T_{\rm
e\perp})^{-1}-{1.27}{\beta_{\rm e\parallel}^{-0.95}}>0$ (see, e.g., Kang et
al., 2019), the growth of the EFI is expected to be enhanced in the upstream
region (Figure 2(e)). On the other hand, the SDA-accelerated electrons and
PPEs, which are advected downstream, would lead to a slight increase of the
electron temperature anisotropy, $T_{e\perp}/T_{e\parallel}$, in the shock
transition zone, $-1\leq(x-x_{\rm sh})/r_{L,i}\leq 0$, (Figure 2(d)). Then,
according to the instability condition for the WI, $I_{\rm WI}\approx T_{\rm
e\perp}/T_{\rm e\parallel}-1-{0.21}{\beta_{\rm e\parallel}^{-0.6}}>0$ (see,
e.g., Guo et al., 2017), the growth of the WI would be slightly enhanced in
the transition zone (Figure 2(f)).

The instability condition for the AIC, $I_{\rm AIC}\approx(T_{\rm
i\perp}/T_{\rm i\parallel})-1-{1.6}{\beta_{\rm i\parallel}^{-0.72}}>0$ (Paper
I), depends on the temperature anisotropy of ions, which is induced by the
gyrating reflected ions in the plane perpendicular to the shock-compressed
magnetic field (Guo et al., 2017). Obviously, the presence of PPEs would not
affect the ion reflection, so we expect that it would not increase the ion
temperature anisotropy, nor enhance the growth of the AIC in the shock
transition zone.

###  3.2 Generation of Waves by Microinstabilities

Figure 3 shows the distributions of the ion number density, $n_{i}$ (top
panels), and the self-excited magnetic field, $B_{z}$ (bottom panels), around
the shock for the M2.0 and M3.0 models with and without PPEs. Oblique waves
with $\lambda\sim 20-30\lambda_{se}$ are generated in the preshock region of
the M2.0 model with PPEs (Figure 3(e)), while such waves are absent in the
M2.0 model without PPEs (Figure 3(f)). They are expected to be induced by the
EFI, and present also in the upstream region of the M3.0 models both with and
without PPEs (Figure 3(g, h)). Kang et al. (2019) demonstrated that the
fraction of reflected, backstreaming electrons, the energy gain via SDA, and
the temperature anisotropy, $T_{e\parallel}/T_{e\perp}$, are large enough to
trigger the EFI only in supercritical shocks with $M_{s}\gtrsim 2.3$. If PPEs
are added to the incoming plasma, on the other hand, the fraction of reflected
electrons and $T_{e\parallel}/T_{e\perp}$ increase as shown in Figure 2(a,c),
and thus the EFI is excited even in the subcritical shock with $M_{s}=2$
(Figure 3(e)). Moreover, the comparison of the models with and without PPEs
reveals that electron-scale waves along the overshoot, excited by the WI, are
somewhat enhanced by PPEs as well.

In Paper I, it was shown that in supercritical shocks, the AIC instability due
to the ion temperature anisotropy could excite ion-scale waves with
$\lambda_{\rm ripple}\sim 70\lambda_{se}$, leading to ripples propagating
along the shock surface (Figure 3 (c,d)). As a result, multiscale plasma
waves, ranging from electron to ion scales, appear in the supercritical M3.0
models (Figure 3 (g,h)). On the contrary, the emergence of such ripples and
ion-scale waves are not observed in the subcritical M2.0 models, regardless of
PPEs (Figure 3 (a,b)).

The presence of plasma waves on relevant kinetic scales is essential for
electron preacceleration and subsequent injection to DSA. This is because for
electron energization up to $p_{\rm inj}$, wave-particle interactions mediated
by waves with wavelengths as long as the gyroradius of electrons with $p_{\rm
inj}$ are required. The injection momentum ranges $p_{\rm inj}\approx
3\sqrt{2m_{i}k_{B}T_{2}}\approx 5.6-7.5m_{e}c$ for the fiducial models with
$M_{s}=2-3$, $T_{1}=10^{8}$ K, and $m_{i}/m_{e}=50$. Hence, the gyroradius of
electrons with $\gamma_{\rm inj}$ is $r_{\rm inj,1}\equiv\gamma_{\rm
inj}c/\Omega_{ce,1}\approx 92-123\lambda_{se}$ and $r_{\rm
inj,2}\equiv\gamma_{\rm inj}c/\Omega_{ce,2}\approx 40-41\lambda_{se}$ in the
preshock and postshock regions, respectively. Here, $\Omega_{ce,1}$ and
$\Omega_{ce,2}$ are the gyrofrequencies of electrons. Figure 3 demonstrates
that in subcritical shocks, such ion-scale waves with $\lambda\sim r_{\rm
inj,1}$ or $r_{\rm inj,2}$ are not excited by the addition of PPEs (Figure
3(a,e)).

Figure 5: Electron energy spectra in the immediate upstream region,
$0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$, at $\Omega_{ci}t\sim 22$, for the models
with PPEs (panel (a)) and without PPEs (panel (b)). The dark purple dashed
lines draw either the far upstream spectrum shown in Figure 1 (panel (a)) or
the Maxwellian spectrum (panel (b)). The magenta and gray boxes denote the
ranges of $\gamma_{\rm inj}$ and $\gamma_{\rm spt}$, respectively. Panel (c)
shows the injection fraction of electrons defined in Equation (1), estimated
for the fiducial models with PPEs in panel (a). Here, $x_{\rm sh}$ is the
shock position, and $r_{L,i}\approx 57\lambda_{se}\cdot(M_{s}/3)$.

For the further characterization of waves, the power spectra of $B_{z}$,
$P_{B_{z}}$, in the immediate upstream (top panels) and downstream (bottom
panels) regions of the shock are shown in Figure 4, for all the fiducial
models with PPEs in Table 1 (left panels) and the corresponding models without
PPEs (right panels). The ranges of $r_{\rm inj}$ (magenta) and the gyroradius,
$r_{\rm spt}$, for the suprathermal momentum, $p_{\rm
spt}\equiv\sqrt{m_{e}/m_{i}}~{}p_{\rm inj}$, (gray) are plotted for reference.
Note that $p_{\rm spt}$ roughly marks the momentum above which electrons
change from thermal to suprathermal distributions (Kang et al., 2019). The
comparison of the subcritical shocks (green and blue lines) with and without
PPEs shows that PPEs enhance waves mainly with $\lambda\sim
20-30\lambda_{se}\ll r_{\rm inj,1},r_{\rm inj,2}$, which cannot mediate the
energization of electrons up to $p_{\rm inj}$. The downstream waves with
$r_{\rm spt,2}<\lambda<r_{\rm inj,2}$ appear partly due to the WI at the shock
transition and partly via the advection of upstream waves.

On the other hand, waves with $\lambda\gtrsim r_{\rm inj,2}$ are produced only
in the supercritical shocks (black and red lines), regardless of the presence
of PPEs. An interesting point is that in the marginally critical shock of
$M_{s}=2.3$, waves with $\lambda\gtrsim r_{\rm inj,2}$ are enhanced slightly
with PPEs, possibly owing to the modification of magnetic field fluctuations
by the EFI mode. Again, Figure 4 confirms that the multiscale plasma waves
that can resonate with electrons of $p_{\rm spt}\lesssim p\lesssim p_{\rm
inj}$ are produced only at supercritical shocks with $M_{s}\gtrsim 2.3$.

###  3.3 Energy Spectrum and Injection Fraction

The consequences of wave-particle interactions should be manifested in the
electron energy spectra. Figure 5 shows the spectra measured in the shock
upstream for all the models with PPEs (Figure 5 (a)) and without PPEs (Figure
5 (b)). The spectra in Figure 5(b) are basically the same as those in Paper I.
In the M2.0 model without PPEs (green line in Figure 5(b)), the EFI is not
triggered, and so only a single SDA cycle occurs. In the M3.0 model without
PPEs (black line in Figure 5(b)), on the other hand, thermal electrons are
energized up to $\gamma\sim 2$ via the Fermi-like acceleration and SSDA.

In the spectra for the models with PPEs in Figure 5(a), two points are noted.
(1) The suprathermal population grows in all the models, even in the M2.0
model. Especially, the population with $\gamma\lesssim 1.5$ represents the
incoming thermal electrons energized via the Fermi-like acceleration mediated
by EFI-driven waves. (2) In the supercritical models, the spectra stretch
beyond $\gamma_{\rm inj}$ through the SSDA mediated by multiscale waves.
However, in the M2.0 model, the high-energy end of the spectrum does not
change, because PPEs do not aid the excitation of AIC-driven ion-scale waves.
These results demonstrate that the presence of PPEs alone does not modify the
critical Mach number, $M_{\rm crit}\approx 2.3$, for electron preacceleration
and subsequent injection to DSA.

In PIC simulations, the maximum energy of electrons is limited by
computational constraints; hence, with a larger simulation box and a longer
integration time, the spectra would extend to higher energies, as noted in
Paper I. Previously, Trotta & Burgess (2019) performed hybrid simulations,
implemented with test-particle electrons of initially kappa distributions. In
their work, electron-scale waves, induced by the EFI and whistler modes, are
not present due to the lack of electron dynamics. Yet, they showed that in
supercritical shocks, electrons could be energized well above $p_{\rm inj}$
through interactions with multiscale waves, which are accompanied by the shock
surface rippling triggered by the ion dynamics.

To quantify electron preacceleration, we estimate the injection fraction,
$\xi$, defined as

$\xi\equiv\int_{\gamma_{\rm
inj}}^{\infty}\frac{dN}{d\gamma}d\gamma~{}\Bigg{/}\int_{1}^{\infty}\frac{dN}{d\gamma}d\gamma,$
(1)

in the shock upstream, $0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$. Figure 5(c) shows
$\xi(M_{s})$ for the fiducial models in Table 1. As expected, $\xi$ increases
sharply at $M_{s}\gtrsim 2.3$, marking the critical Mach number. We point that
$\xi$ shown in Figure 5(c) is not necessarily the injection fraction of
electrons at $Q_{\perp}$-shocks in the ICM. Our PIC simulations are limited
only to reproduce processes on kinetic scales. With a much larger spatial
domain and a much longer dynamical time, both thermal and suprathermal
electrons are expected to be preaccelerated and injected to DSA in
supercritical shocks, hence the injection fraction would be substantially
larger than $\xi$ in Figure 5(c). While the estimation of realistic injection
fractions is beyond the scope of this paper, we conjecture that $\xi$ should
be small at subcritical shocks unless there are additional processes to boost
it.

Figure 6: Magnetic field power spectra,
$P_{B_{z}}(k_{y})\propto(k_{y}/2\pi)(\delta B_{z}(k_{y})^{2}/B_{0}^{2})$, in
the immediate downstream region, $-1\leq(x-x_{\rm sh})/r_{L,i}\leq 0$, at
$\Omega_{ci}t\sim 22$, for the M2.0 (top) and M3.0 (bottom) models with
different parameters (see Table 1). Here, $x_{\rm sh}$ is the shock position,
and $r_{L,i}\approx
57\lambda_{se}\cdot(M_{s}/3)\sqrt{\beta/20}\sqrt{(m_{i}/m_{e})/50}$. The
magenta and gray boxes denote the ranges of the gyroradius of electrons with
$p_{\rm inj}$ and $p_{\rm spt}$, respectively, in the downstream region.

Figure 7: Electron energy spectra in the immediate upstream region,
$0\leq(x-x_{\rm sh})/r_{L,i}\leq 1$, at $\Omega_{ci}t\sim 22$, for the M2.0
(upper panels) and M3.0 (lower panels) models with different parameters (see
Table 1). The dashed lines draw the electron energy spectra in the far
upstream region, $10\leq(x-x_{\rm sh})/r_{L,i}\leq 11$, at $\Omega_{ci}t\sim
3$, for comparison. Here, $x_{\rm sh}$ is the shock position, and
$r_{L,i}\approx
57\lambda_{se}\cdot(M_{s}/3)\sqrt{\beta/20}\sqrt{(m_{i}/m_{e})/50}$. Some of
the blue solid and dashed lines marked by ($\times 0.1$) are shifted
vertically for clarity.

###  3.4 Dependence on Parameters

We next examine how our findings depend on the parameters that specify PPEs,
such as $n_{p}/n_{0}$, $\alpha$, and $\gamma_{\rm max}$, as well as the
simulation parameters, such as $m_{i}/m_{e}$ and $\beta$. Figure 6 shows the
magnetic field power spectra, $P_{B_{z}}$, in the immediate downstream, for
the M2.0 models in the second group (Figure 6(a)) and the M3.0 models in the
third group (Figure 6(b)) of Table 1. The comparison of all the lines in
Figure 6(b) with the black line in Figure 4(c) indicates that the generation
of waves does not strongly depend on these parameters in the case of the
supercritical shocks with $M_{s}=3$. In all the M3.0 models, $P_{B_{z}}$ for
$r_{\rm spt,2}\lesssim\lambda\lesssim r_{\rm inj,2}$ is substantial. By
contrast, there are some differences in the case of the subcritical M2.0
models. For instance, in the M2.0-np0.1 model (green line in Figure 6(a)),
with $n_{p}/n_{0}$ ten times larger than that of the fiducial model (green
line in Figure 4(a)), the excitation of waves is enhanced. In the
M2.0-$\beta$50-m25 model (orange line in Figure 6(a)), with a weaker
background magnetic field, waves with longer wavelengths appear. As noted
above, waves in the M2.0 models are mostly electron-scale waves, induced by
the EFI and whistler modes.

Finally, Figure 7 compares the electron energy spectra in the immediate
upstream for the M2.0 models in the second group (upper panels) and the M3.0
models in the third group (lower panels) of Table 1. The energy spectra extend
up to $\gamma-1\sim 10$ in all the M3.0 models, independent of the parameters.
This is because the ripples along the shock surface have $\lambda_{\rm
ripple}\sim 70\lambda_{se}$ (see Figure 3(c,d)), which marks the longest
wavelength of AIC-driven waves. Hence, the maximum energy in the spectra is
given as $\gamma_{\rm max}\approx\lambda_{\rm ripple}\Omega_{ce,2}/c\approx
11$, in our simulations. Again, the main results described in the previous
subsection remain valid, independent of the model parameters.

##  4 Summary and Discussion

Electrons are expected to be preaccelerated and injected to the Fermi-I, DSA
process in supercritical $Q_{\perp}$-shocks with $M_{s}\gtrsim 2.3$ in the
high-$\beta$ ICM (Kang et al., 2019, and Paper I). Here, we explore if the
presence of preexisting nonthermal electrons could facilitate the generation
of kinetic waves via microinstabilities and assist the electron injection to
DSA, especially in the subcritical regime. We have performed 2D PIC
simulations for the following parameters: $M_{s}=2-3$, $\theta_{\rm
Bn}=63^{\circ}$, and $\beta=20$ (see Table 1). Preexisting nonthermal
electrons are represented by a power-law component (PPEs) that extends up to
$\gamma_{\rm cut}=5-10$.

The main findings are summarized as follows. (1) The presence of PPEs enhances
the EFI owing to the increased fraction of the electrons reflected at the
shock ramp. As a result, waves with $\lambda\sim 20-30\lambda_{se}$, generated
by the EFI, appear even in the upstream region of the subcritical shocks in
our M2.0 models. At the same time, the suprathermal population increases
through the preacceleration via the Fermi-like acceleration mediated by the
EFI waves in the M2.0 models, as well as in other models with supercritical
shocks. (2) Ion-scale waves with $\lambda\gtrsim r_{L,i}\approx
57\lambda_{se}\cdot(M_{s}/3)$ are excited only in the supercritical shocks of
the M2.3 - M3.0 models, as the consequence of the AIC instability, regardless
of the presence of PPEs. In those models, the energy spectrum of electrons
continues beyond $p_{\rm inj}$ via the SSDA mediated by multiscale waves
ranging from electron to ion scales. On the other hand, ion-scale waves are
not produced in the subcritical M2.0 models, since PPEs do not facilitate the
excitation of ion-scale waves via the AIC instability. Hence, in the M2.0
models, the high-energy end of the electron energy spectrum does not change
even with PPEs.

We further comment on the potential effects of the three ingredients that are
not fully accounted for in our simulations: (1) realistic preexisting
nonthermal electrons, (2) preexisting nonthermal ions, and (3) preexisting
kinetic turbulence in the ICM.

As noted above, while the preexisting power-law component of electrons extends
only up to $\gamma\lesssim\gamma_{\rm inj}$ in our simulations due to limited
computational capacities, the nonthermal population in the ICM is expected to
stretch to higher energies of $\gamma\sim 10^{2}-10^{3}$, considering the
physical condition there (e.g., Sarazin, 1999; Mandal et al., 2020). Hence, in
reality, some of the electrons reflected at the shock ramp may have energies
higher than those reproduced in our simulations. On the other hand, the
properties of the EFI induced by the reflected electrons, such as the
wavelength and growth rate of the fastest growing mode, are not very sensitive
to the characteristics of the electron beam that induces the instability (see,
e.g., Kim et al., 2020); as long as the bulk of the beam energy resides in the
low-energy part of the power-law distribution, most of the EFI-induced waves
would be still on electron scales. Although the details should be further
investigated, we suspect that the power-law nonthermal electrons that extend
to higher energies would not resolve the issue of the generation of ion-scale
waves in subcritical shocks.

In the ICM, not just relativistic electrons, but also relativistic CR ions
(mostly protons) are expected to be present. The ratio of CR to thermal ions,
which is constrained by non-detection of gamma-rays from nearby clusters, for
instance, in Fermi-LAT observations, seems to be fairly small, less than
$10^{-4}-10^{-3}$, averaged over all the cluster volume (e.g., Ryu et al.,
2019; Ha et al., 2020; Wittor et al., 2020). Yet, we may postulate a scenario
that about the same numbers of nonthermal electrons and ions are present in
the upstream of merger shocks. However, the presence of nonthermal power-law
ions would not change the picture of electron preacceleration described in
this paper. In weak $Q_{\perp}$-shocks, unlike suprathermal electrons, higher
energy ions are more likely to pass through the shock ramp, rather than to be
reflected at the shock potential barrier; hence, the fraction of reflected
ions would not increase and the AIC instability would not be enhanced by
adding nonthermal power-law ions. Thus, we conclude that the nonthermal
population of CR ions in the upstream plasma would not change the results of
this work. Although not shown here, it has been confirmed in PIC simulations
we have additionally performed.

A more promising possibility would be the presence of broadband magnetic
fluctuations, possibly produced by the turbulence in the ICM. Preexisting
nonthermal electrons with energies up to $\gamma\sim 10^{2}-10^{3}$ could be
directly injected into DSA via resonant scatterings off those preexisting
turbulent waves. As a matter of fact, previous studies using test-particle
simulations for electrons showed that electrons could be accelerated via SDA
and DSA, regardless of the shock obliquity angle, by interacting with
preexisting large-amplitude magnetic fluctuations in interplanetary shocks
(e.g., Giacalone, 2005; Guo & Giacalone, 2015). While such acceleration
requires efficient wave-particle interactions in the ICM environment, the
details of involved processes should be investigated through further
simulations.

Additionally, we point that cosmological hydrodynamic simulations showed that
radio relics, formed in the turbulent ICM, normally consist of shock surfaces
with varying $M_{\rm s}$ and $\theta_{\rm Bn}$ (e.g., Hong et al., 2015; Roh
et al., 2019; Wittor et al., 2021). As suggested by several previous papers
including Paper I, it would be interesting to consider a scenario, in which
nonthermal electrons, accelerated in the locally supercritical portions of the
shock surface, are injected to DSA at subcritical portions. Furthermore, the
variation of $\theta_{Bn}$ in the shock surface could contribute to the
generation of the large-scale modulation in the surface. In the interplanetary
shocks, for instance, the surface ripplings on scales larger than the AIC-
induced waves have been detected (e.g., Kajdič et al., 2019), which are
thought to be generated due to the upstream magnetic field fluctuations
produced by the backstreaming ions at $Q_{\parallel}$-portions.

Finally, we thus conjecture that there is still the possibility of electron
DSA in radio relics with subcritical shocks, when additional processes and/or
ingredients other than those considered in this work are included. The
investigation of such elements is beyond the scope of this article, so we
leave it for future work. For now, we conclude that the presence of nonthermal
electrons in the ICM “alone” would not resolve the issue of electron
preacceleration and injection into DSA, and hence could not explain the
production of CR electrons in radio relics with subcritical shocks.

This work was supported by the National Research Foundation (NRF) of Korea
through grants 2016R1A5A1013277, 2020R1A2C2102800, 2020R1F1A1048189, and
2020R1C1C1012112. J.-H. H. was also supported by the Global PhD Fellowship of
the NRF through grant 2017H1A2A1042370. Some of simulations were performed
using the high performance computing resources of the UNIST Supercomputing
Center.

## References

  * Achterberg et al. (2001) Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W. 2001, MNRAS, 328, 393, doi: 10.1046/j.1365-8711.2001.04851.x
  * Akamatsu & Kawahara (2013) Akamatsu, H., & Kawahara, H. 2013, PASJ, 65, 16, doi: 10.1093/pasj/65.1.16
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29, doi: 10.1086/182658
  * Botteon et al. (2020) Botteon, A., Brunetti, G., Ryu, D., & Roh, S. 2020, A&A, 634, A64, doi: 10.1051/0004-6361/201936216
  * Brüggen et al. (2012) Brüggen, M., Bykov, A., Ryu, D., & Röttgering, H. 2012, Space Sci. Rev., 166, 187, doi: 10.1007/s11214-011-9785-9
  * Brunetti & Jones (2014) Brunetti, G., & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007, doi: 10.1142/S0218271814300079
  * Buneman (1993) Buneman, O. 1993, in Computer Space Plasma Physics: Simulation Techniques and Software, ed. H. Matsumoto & Y. Omura, 67
  * Caprioli et al. (2015) Caprioli, D., Pop, A.-R., & Spitkovsky, A. 2015, ApJ, 798, L28, doi: 10.1088/2041-8205/798/2/L28
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973, doi: 10.1088/0034-4885/46/8/002
  * Giacalone (2005) Giacalone, J. 2005, ApJ, 624, 765, doi: 10.1086/429265
  * Gosling et al. (1989) Gosling, J. T., Thomsen, M. F., Bame, S. J., & Russell, C. T. 1989, J. Geophys. Res., 94, 10011, doi: 10.1029/JA094iA08p10011
  * Guo & Giacalone (2015) Guo, F., & Giacalone, J. 2015, ApJ, 802, 97, doi: 10.1088/0004-637X/802/2/97
  * Guo et al. (2014a) Guo, X., Sironi, L., & Narayan, R. 2014a, ApJ, 794, 153, doi: 10.1088/0004-637X/794/2/153
  * Guo et al. (2014b) —. 2014b, ApJ, 797, 47, doi: 10.1088/0004-637X/797/1/47
  * Guo et al. (2017) —. 2017, ApJ, 851, 134, doi: 10.3847/1538-4357/aa9b82
  * Ha et al. (2021) Ha, J.-H., Kim, S., Ryu, D., & Kang, H. 2021, ApJ, 915, 18, doi: 10.3847/1538-4357/abfb68
  * Ha et al. (2018) Ha, J.-H., Ryu, D., & Kang, H. 2018, ApJ, 857, 26, doi: 10.3847/1538-4357/aab4a2
  * Ha et al. (2020) —. 2020, ApJ, 892, 86, doi: 10.3847/1538-4357/ab7c5b
  * Hong et al. (2015) Hong, S. E., Kang, H., & Ryu, D. 2015, ApJ, 812, 49, doi: 10.1088/0004-637X/812/1/49
  * Hong et al. (2014) Hong, S. E., Ryu, D., Kang, H., & Cen, R. 2014, ApJ, 785, 133, doi: 10.1088/0004-637X/785/2/133
  * Kajdič et al. (2019) Kajdič, P., Preisser, L., Blanco-Cano, X., Burgess, D., & Trotta, D. 2019, ApJ, 874, L13, doi: 10.3847/2041-8213/ab0e84
  * Kang (2016a) Kang, H. 2016a, Journal of Korean Astronomical Society, 49, 83, doi: 10.5303/JKAS.2016.49.3.83
  * Kang (2016b) —. 2016b, Journal of Korean Astronomical Society, 49, 145, doi: 10.5303/JKAS.2016.49.4.145
  * Kang et al. (2002) Kang, H., Jones, T. W., & Gieseler, U. D. J. 2002, ApJ, 579, 337, doi: 10.1086/342724
  * Kang & Ryu (2011) Kang, H., & Ryu, D. 2011, ApJ, 734, 18, doi: 10.1088/0004-637X/734/1/18
  * Kang & Ryu (2015) —. 2015, ApJ, 809, 186, doi: 10.1088/0004-637X/809/2/186
  * Kang & Ryu (2016) —. 2016, ApJ, 823, 13, doi: 10.3847/0004-637X/823/1/13
  * Kang et al. (2019) Kang, H., Ryu, D., & Ha, J.-H. 2019, The Astrophysical Journal, 876, 79, doi: 10.3847/1538-4357/ab16d1
  * Kang et al. (2012) Kang, H., Ryu, D., & Jones, T. W. 2012, ApJ, 756, 97, doi: 10.1088/0004-637X/756/1/97
  * Katou & Amano (2019) Katou, T., & Amano, T. 2019, ApJ, 874, 119, doi: 10.3847/1538-4357/ab0d8a
  * Kierdorf et al. (2017) Kierdorf, M., Beck, R., Hoeft, M., et al. 2017, A&A, 600, A18, doi: 10.1051/0004-6361/201629570
  * Kim et al. (2020) Kim, S., Ha, J.-H., Ryu, D., & Kang, H. 2020, ApJ, 892, 85, doi: 10.3847/1538-4357/ab7cd9
  * Kim et al. (2021) —. 2021, ApJ, 913, 35, doi: 10.3847/1538-4357/abf1e1
  * Kirk et al. (2000) Kirk, J. G., Guthmann, A. W., Gallant, Y. A., & Achterberg, A. 2000, ApJ, 542, 235, doi: 10.1086/309533
  * Kobzar et al. (2021) Kobzar, O., Niemiec, J., Amano, T., et al. 2021, ApJ, 919, 97, doi: 10.3847/1538-4357/ac1107
  * Macario et al. (2013) Macario, G., Venturi, T., Intema, H. T., et al. 2013, A&A, 551, A141, doi: 10.1051/0004-6361/201220667
  * Mandal et al. (2020) Mandal, S., Intema, H. T., van Weeren, R. J., et al. 2020, A&A, 634, A4, doi: 10.1051/0004-6361/201936560
  * Markevitch et al. (2002) Markevitch, M., Gonzalez, A. H., David, L., et al. 2002, ApJ, 567, L27, doi: 10.1086/339619
  * Markevitch & Vikhlinin (2007) Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1, doi: 10.1016/j.physrep.2007.01.001
  * Matsukiyo & Matsumoto (2015) Matsukiyo, S., & Matsumoto, Y. 2015, in Journal of Physics Conference Series, Vol. 642, Journal of Physics Conference Series, 012017, doi: 10.1088/1742-6596/642/1/012017
  * Matsukiyo et al. (2011) Matsukiyo, S., Ohira, Y., Yamazaki, R., & Umeda, T. 2011, ApJ, 742, 47, doi: 10.1088/0004-637X/742/1/47
  * Melzani et al. (2013) Melzani, M., Winisdoerffer, C., Walder, R., et al. 2013, A&A, 558, A133, doi: 10.1051/0004-6361/201321557
  * Miniati et al. (2000) Miniati, F., Ryu, D., Kang, H., et al. 2000, ApJ, 542, 608, doi: 10.1086/317027
  * Niemiec et al. (2019) Niemiec, J., Kobzar, O., Amano, T., et al. 2019, in International Cosmic Ray Conference, Vol. 36, 36th International Cosmic Ray Conference (ICRC2019), 368
  * Park et al. (2015) Park, J., Caprioli, D., & Spitkovsky, A. 2015, Phys. Rev. Lett., 114, 085003, doi: 10.1103/PhysRevLett.114.085003
  * Pfrommer et al. (2006) Pfrommer, C., Springel, V., Enßlin, T. A., & Jubelgas, M. 2006, MNRAS, 367, 113, doi: 10.1111/j.1365-2966.2005.09953.x
  * Pinzke et al. (2013) Pinzke, A., Oh, S. P., & Pfrommer, C. 2013, MNRAS, 435, 1061, doi: 10.1093/mnras/stt1308
  * Porter et al. (2015) Porter, D. H., Jones, T. W., & Ryu, D. 2015, ApJ, 810, 93, doi: 10.1088/0004-637X/810/2/93
  * Roh et al. (2019) Roh, S., Ryu, D., Kang, H., Ha, S., & Jang, H. 2019, ApJ, 883, 138, doi: 10.3847/1538-4357/ab3aff
  * Ryu et al. (2008) Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Science, 320, 909, doi: 10.1126/science.1154923
  * Ryu et al. (2019) Ryu, D., Kang, H., & Ha, J.-H. 2019, ApJ, 883, 60, doi: 10.3847/1538-4357/ab3a3a
  * Ryu et al. (2003) Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, ApJ, 593, 599, doi: 10.1086/376723
  * Sarazin (1999) Sarazin, C. L. 1999, in Diffuse Thermal and Relativistic Plasma in Galaxy Clusters, ed. H. Boehringer, L. Feretti, & P. Schuecker, 185.  https://arxiv.org/abs/astro-ph/9906272
  * Schaal & Springel (2015) Schaal, K., & Springel, V. 2015, MNRAS, 446, 3992, doi: 10.1093/mnras/stu2386
  * Skillman et al. (2008) Skillman, S. W., O’Shea, B. W., Hallman, E. J., Burns, J. O., & Norman, M. L. 2008, ApJ, 689, 1063, doi: 10.1086/592496
  * Spitkovsky (2005) Spitkovsky, A. 2005, in AIP Conf. Proc., Astrophysical Sources of High Energy Particles and Radiation, ed. T. Bulik, B. Rudak, & G. Madejski (New York: AIP), 801, 345, doi: 10.1063/1.2141897
  * Trotta & Burgess (2019) Trotta, D., & Burgess, D. 2019, MNRAS, 482, 1154, doi: 10.1093/mnras/sty2756
  * van Weeren et al. (2019) van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, Space Sci. Rev., 215, 16, doi: 10.1007/s11214-019-0584-z
  * van Weeren et al. (2010) van Weeren, R. J., Röttgering, H. J. A., Brüggen, M., & Hoeft, M. 2010, Science, 330, 347, doi: 10.1126/science.1194293
  * van Weeren et al. (2017) van Weeren, R. J., Andrade-Santos, F., Dawson, W. A., et al. 2017, Nature Astronomy, 1, 0005, doi: 10.1038/s41550-016-0005
  * Vazza et al. (2009) Vazza, F., Brunetti, G., & Gheller, C. 2009, MNRAS, 395, 1333, doi: 10.1111/j.1365-2966.2009.14691.x
  * Wittor et al. (2021) Wittor, D., Ettori, S., Vazza, F., et al. 2021, MNRAS, 506, 396, doi: 10.1093/mnras/stab1735
  * Wittor et al. (2020) Wittor, D., Vazza, F., Ryu, D., & Kang, H. 2020, MNRAS, 495, L112, doi: 10.1093/mnrasl/slaa066