11institutetext:  Paris-Saclay University, CentraleSupélec, 91190, Gif-sur-
Yvette, France 22institutetext: Paris-Saclay University, Gustave Roussy,
Inserm, 94800, Villejuif, France 33institutetext: TheraPanacea, 75014, Paris,
France

# Sparse convolutional context-aware multiple instance learning for whole
slide image classification

Marvin Lerousseau 1122 Maria Vakalopoulou 11 Nikos Paragios 33 Eric Deutsch 22

###### Abstract

Whole slide microscopic slides display many cues about the underlying tissue
guiding diagnostic and the choice of therapy for many diseases. However, their
enormous size often in gigapixels hampers the use of traditional neural
network architectures. To tackle this issue, multiple instance learning (MIL)
classifies bags of patches instead of whole slide images. Most MIL strategies
consider that patches are independent and identically distributed. Our
approach presents a paradigm shift through the integration of spatial
information of patches with a sparse-input convolutional-based MIL strategy.
The formulated framework is generic, flexible, scalable and is the first to
introduce contextual dependencies between decisions taken at the patch level.
It achieved state-of-the-art performance in pan-cancer subtype classification.
The code of this work will be made available at censored for review.

###### Keywords:

Multiple instance learning Whole slide image classification Weakly supervised
learning.

##  1 Introduction

Histopathology, or the microscopic examination of tissues, is essential in
diagnostic medicine. The advent of digital scanners is leading to an
increasing digitization of microscopic slides into whole slide images (WSI)
[22]. These images encompass a great deal of information about the phenotypes
of cells and reach billions of pixels, preventing the use of traditional
convolutional neural networks (CNN) due to memory limitations. Downscaling is
not viable since most of the cues about a tumor are displayed at the cell-
level (or nuclei-level), which are lost with downsampling. Besides, WSI cannot
be classified from a contiguous extracted region since pertinent information
may not be evenly spread in those images due to tumor heterogeneity [12].
Therefore, the scientific community has predominantly adopted the strategy of
classifying patches extracted throughout WSI using the multiple instance
learning [9] paradigm.

In this work, we propose a novel type of framework built on top of multiple
instance learning that exploits the spatial structure of WSI with the use of
sparse-input convolutional neural networks. In particular, our methodology is
designed to harvest the spatial context of patches as additional information
for improving WSI classification. On another hand, the relative localization
of patches should enhance the consistency of features from neighboring patches
which often display the same tissue phenotypes. Our contributions are twofold:
(i) we propose a novel sparse embedding map to fully process huge WSI keeping
the spatial relations of the different instances, (ii) large scale experiments
with other state of the art MIL schemes prove the high potentials of our
formulation. Our framework is open-sourced at censored for review.

Figure 1: Visual representation of the proposed approach. First, a set of
patches are randomly sampled throughout a WSI, and are then concurrently and
independently forwarded into a shared patch embedder $f_{\theta_{1}}$. Then, a
sparse map is built by placing each resulting embedding at the location of its
associated patch. This map is forwarded into a sparse-input CNN producing a
bag embedding which is finally classified into probabilities.

##  2 Background

###  2.1 Multiple instance learning

MIL deals with problems where labels are associated with groups of instances,
or bags (whole slide images in our case), whereas instance (tiles or patches)
labels are unknown or non-existent. In the most general case, MIL models can
be decomposed into 3 parts: an instance embedder $f_{\theta_{1}}$, a pooling
operator $g_{\theta_{2}}$ and a classifier $h_{\theta_{3}}$ such that a
prediction for a bag $(x_{1},\dots,x_{K})$ is obtained with
$h_{\theta_{3}}\Big{(}g_{\theta_{2}}\big{(}f_{\theta_{1}}(x_{1}),\dots,f_{\theta_{1}}(x_{K})\big{)}\Big{)}$.
Instance-level MIL is a particular case of MIL where embeddings are one-
dimensional i.e. apparent to probabilities. Most of the MIL community efforts
have revolved around the pooling operator $g_{\theta_{2}}$ which can be placed
in two categories: permutation invariant operators, and others.

####  2.1.1 Permutation invariant pooling operators

In some applications, such as the drug discovery problem [3], instances do not
exhibit dependency, ordering, or spatial information among each other, i.e.
they are independent and identically distributed (i.i.d.). Equivalently, for
any permutation $\sigma$, the output of an i.i.d. pooling operator is the same
for $(x_{1},\dots,x_{K})$ and $(\sigma(x_{1}),\dots,\sigma(x_{K}))$, i.e. they
are permutation invariant, such as:

$\begin{array}[]{ll}\textnormal{mean}&\colon(x_{1},...,x_{K})\mapsto\dfrac{(x_{1}+...+x_{K})}{K}\\\
\textnormal{max}&\colon(x_{1},...,x_{K})\mapsto\max(x_{1},...,x_{K})\\\
\textnormal{log-sum-
exp~{}\cite[cite]{[\@@bibref{}{ramon2000multi}{}{}]}}&\colon(x_{1},...,x_{K})\mapsto\dfrac{1}{M}\log\big{(}\sum\exp(M\times
x_{k})\big{)}\\\
\textnormal{attention~{}\cite[cite]{[\@@bibref{}{ilse2018attention}{}{}]}}&\colon(x_{1},...,x_{K})\mapsto{\sum_{k}}\dfrac{\exp\big{(}w^{\top}\tanh(Vx_{k}^{\top})\big{)}}{\sum\exp\big{(}w^{\top}\tanh(Vx_{j}^{\top})\big{)}}\cdot
x_{k}\\\ \textnormal{gated-
attention~{}\cite[cite]{[\@@bibref{}{ilse2018attention}{}{}]}}&\colon(x_{1},...,x_{K})\mapsto{\sum_{k}}\dfrac{\exp\big{(}w^{\top}(\tanh(Vx_{k}^{\top})\odot\mathop{\mathrm{sigm}}(Ux_{k}^{\top}))\big{)}}{\sum\exp\big{(}w^{\top}(\tanh(Vx_{j}^{\top})\odot\mathop{\mathrm{sigm}}(Ux_{j}^{\top})\big{)}}\cdot
x_{k}\\\ \end{array}$

where $M\in\mathbb{R}^{+}$, $V\in\mathbb{R}^{L\times\dim(\mathbb{H})}$,
$U\in\mathbb{R}^{L\times\dim(\mathbb{H})}$, $w\in\mathbb{R}^{L\times 1}$,
$L\in\mathbb{N}^{*}$ are parameters, $\odot$ is the elementwise
multiplication, and $\mathop{\mathrm{sigm}}$ is the elementwise sigmoid
function.

####  2.1.2 Permutation noninvariant pooling operators

The i.i.d. assumption does not hold for applications where there is inherent
structural information about instances, such as document classification from
sentences [1], or WSI classification from tiles. Zhou et al. [24] have notably
achieved state-of-the-art performance over i.i.d. MIL poolings by building a
graph from instance embeddings, and then performing classification with kernel
methods. With recent improvements of graph neural networks [20], further
iterations have been proposed by the community [18, 21, 23].

###  2.2 Sparse-input convolutional neural network

Convolutional layers are optimized for dense inputs such as natural images or
3D medical volumes but are not efficient for sparse data in terms of both
memory consumption and FLOPs. Sparse data can be defined as variables
containing few active cells, i.e. most of their volume is empty \- examples
are cloud points from LiDAR or tiles extracted from WSI. In particular and for
clarity, the sparsity in sparse-input CNN is related to the input, and not to
the weights sparsity which is another matter. Several convolutional
implementations designed for sparse input data have been proposed [4, 5, 14,
6] revolving around the idea that convolutions should only be performed on
active data points, i.e. by ignoring empty input regions, therefore decreasing
the number of computations. Additionally, sparse data have their own tensor
structure which considers only active data points, therefore using virtually
no memory for non-active tensor cells.

##  3 Methods

We propose to use a sparse-input convolutional neural network to learn the MIL
operator in whole slide image classification applications; we coin our
framework SparseConvMIL (Figure 1). More specifically, our pooling operator
first builds a sparse map of instance embeddings, which is then processed into
a bag embedding with a sparse-CNN.

Let us consider one input WSI $X\in[0,1]^{w\times h\times 3}$, where $w$ and
$h$ are typically in dozens of thousands. $b$ locations $l_{X}=(i_{k},j_{k})$
are randomly sampled from $X$. A tile $(x_{i_{k},j_{k}})$ is extracted from
each location. Providing an instance embedder $f_{\theta_{1}}$, all tiles are
independently and concurrently inferred, producing the set of $b$ embeddings
$f_{\theta_{1}}(x_{i_{k},j_{k}})\in\mathbb{H}$. The sparse embedding map
$S^{X}$ is built by corresponding each tile embedding with their location:

$S^{X}=\Big{\\{}\big{[}(i_{k},j_{k}),f_{\theta_{1}}(x_{i_{k},j_{k}})\big{]};1\leq
k\leq b\Big{\\}}\subset(\mathbb{N}\times\mathbb{N}\times\mathbb{H})^{b}$

Alternatively, $S^{X}$ can be characterized by:

$\forall i\in\llbracket 1,w\rrbracket,\forall j\in\llbracket
1,h\rrbracket:S^{X}_{i,j}=\begin{cases}\O&\textnormal{if }(i,j)\notin l_{X}\\\
f_{\theta_{1}}(x_{i,j})&\mathop{\mathrm{otherwise}}\end{cases}$

Once $S^{X}$ is reconstructed, it is forwarded into a sparse-input CNN which
can be made of convolutional layers, pooling layers, and global pooling
layers. For a convolution layer
$U\in\mathbb{R}^{2f+1}\times\mathbb{R}^{2f+1}\times\mathbb{R}^{o}$ with a
filter of half-size $f$, stride $s$, and $o$ output channels, the output
$U\circledast S^{X}$ is:

$\forall i\in\left\llbracket
1,\left\lfloor\frac{w}{s}\right\rfloor\right\rrbracket,\forall
j\in\left\llbracket 1,\left\lfloor\frac{h}{s}\right\rfloor\right\rrbracket:\\\
(U\circledast
S^{X})_{i,j}=\sum\limits_{m=-f}^{f}\sum\limits_{n=-f}^{f}U_{m+f,n+f}\cdot
S^{X}_{si+m,sj+n}$

where a term $S^{X}_{si+m,sj+n}=\text{\O}$ has no impact in the sum. In
particular, for a position $(i,j)$, if no embeddings are present in a
neighboring region of half-width $f$, i.e. when the set
$\\{S_{X;m,n};|m-i|\leq f,|n-j|\leq f\\}$ is empty, the output of the
convolution is Ø which implies that the output $U\circledast S^{X}$ is also
sparse.

Similarly, for a function $p$ such as the functions max or average, a pooling
layer with filter size $f$ and stride $s$ produces the output $p(S^{X})$ from
an input sparse map $S^{X}$ such that:

$\forall i\in\left\llbracket
1,\left\lfloor\frac{w}{s}\right\rfloor\right\rrbracket,\forall
j\in\left\llbracket
1,\left\lfloor\frac{h}{s}\right\rfloor\right\rrbracket:p(S^{X})_{i,j}=p\big{(}\\{S^{X}_{si+m,sj+n};m\leq
f,n\leq f\\}\big{)}$

In particular, adaptive global pooling layers with output size $o$ are defined
as pooling layers of stride $1$ and of filter size
$\left\lfloor\frac{w}{o}\right\rfloor$ and
$\left\lfloor\frac{h}{o}\right\rfloor$. Noteworthily, locations with no data
do not interfere with the pooling function.

Eventually, after a succession of sparse-input layers, a bag embedding is
obtained. Bag probabilities are then computed with any type of classifier,
including multi-layer perceptrons.

We introduce an additional parameter called the downsampling factor. In
theory, a sparse map thus built has the same size as its associated WSI.
However in practice, few tiles can be extracted from a WSI, implying that most
convolutional operations would involve at most one activate size which would
cancel out the effect of convolutions. To address this issue, sparse maps can
be downsampled, i.e. all tile locations can be divided by a shared factor -
the downsampling factor. We evaluate its impact in our experiments.

Finally, SparseConvMIL benefits from additional data augmentation over
permutation noninvariant pooling strategies, precisely because it treats
instances as non i.i.d. Spatial augmentations (e.g. flipping, rotations, local
shuffling or elastic deformations) performed on tiles locations, or
equivalently on the sparse map, can help reduce the burden of overfitting by
increasing the input data to the pooling CNN. Besides, these augmentations can
be performed after tile embedding inference, implying that multiple sparse map
spatial augmentations can be done with a low additional memory footprint.

##  4 Experimental validation

We have used two well studied data examples to demonstrate the extreme
potentials of our method. The first is CRCHistoPhenotype for classifying
whether power fields contain epithelial cells, and the second is The Cancer
Genome Atlas for pan-cancer subtype classification.

###  4.1 Classical MIL dataset

####  4.1.1 Dataset

The CRCHistoPhenotype [17] dataset consists of 100 haematoxylin and eosin-
stained (H&E) $500\times 500$ pixels histology images of colorectal
adenocarcinomas. A total of 22,444 nuclei were annotated with the position of
their center and their type amid epithelial, inflammatory, fibroblast, or
miscellaneous. In this context, bags are made of $27\times 27$ pixel patches
extracted around each nuclei center, and are considered positive iff they
contain at least one epithelial cell, resulting in 51 positive and 49 negative
bags.

####  4.1.2 Implementation details

All methods share the same training parameters which are detailed thereafter.
The instance embedding model $f_{\theta_{1}}$ (supplementary data Table 1)
proposed in [17] and used in [8] is employed. Binary cross-entropy computes
error signals which update the models parameters using Adam optimizer [10]
with default momentum values, learning rate of $1e^{-4}$, weight decay of
$5e^{-4}$, batch size of $1$ for 100 epochs. For data augmentation, tiles are
randomly rotated and flipped. Besides, each tile is first color deconvoluted
in the HE space using [15], then reconstructed into RGB space after
multiplication of both H and E stains with two independent random gaussian
variables with mean 1 and standard deviation of 3. Finally each tile is
standard scaled with mean and standard deviation extracted from the training
set. Sparse maps of the proposed approach are downsampled by $2$, i.e. there
are of width 250 pixels. All methods share the same data split with 55, 20 and
25 samples for resp. the training, validation and testing sets. The validation
set is used to select the snapshot with least validation error for inference
on the testing set.

Attention-based approaches use a two-layer neural network for attention with
$128$ hidden neurons as used in [8], resulting in 66304 pooling parameters and
131968 for the gated version. The proposed SparseConvMIL is implemented with
two 12-channels convolutional filters with filter size 3, resulting in a
module with 56628 parameters, or less than both previous. The implementation
of other methods are detailed in supplementary Table 2. The training/testing
process was performed 5 times for each method to derive confidence intervals.

Method |  Accuracy |  Precision |  Recall |  F1-score |  AUC
---|---|---|---|---|---
Instance+max |  0.842$\pm$0.021 |  0.866$\pm$0.017 |  0.816$\pm$0.031 |  0.839$\pm$0.023 |  0.914$\pm$0.010
Instance+mean |  0.772$\pm$0.012 |  0.821$\pm$0.011 |  0.710$\pm$0.031 |  0.759$\pm$0.017 |  0.866$\pm$0.008
max |  0.824$\pm$0.015 |  0.884$\pm$0.014 |  0.753$\pm$0.020 |  0.813$\pm$0.017 |  0.918$\pm$0.010
mean |  0.860$\pm$0.014 |  0.911$\pm$0.011 |  0.804$\pm$0.027 |  0.853$\pm$0.016 |  0.940$\pm$0.010
Attention [8] |  0.904$\pm$0.011 |  0.953$\pm$0.014 |  0.855$\pm$0.017 |  0.901$\pm$0.011 |  0.968$\pm$0.009
Gated-Attention [8] |  0.898$\pm$0.020 |  0.944$\pm$0.016 |  0.851$\pm$0.035 |  0.893$\pm$0.022 |  0.968$\pm$0.010
Proposed |  0.944$\pm$0.019 |  0.929$\pm$0.021 |  0.944$\pm$0.019 |  0.932$\pm$0.024 |  0.958$\pm$0.008
Table 1: Results on CRCHistoPhenotype in mean $\pm$ standard deviation of 5
runs. Otherwise specified, methods are embedding-level MIL.

####  4.1.3 Results

Results are reported in Table 1. The proposed approach achieves the best
performance in terms of balanced accuracy and f1-score. Although the proposed
obtained a slightly lower precision than the attention-based methods, it
achieved a significantly higher recall. This is preferable for clinical
practice since carcinomas arise from epithelial cells [19], which implies that
no epithelial tissue should be discarded for exhaustive cancer screening.

###  4.2 Large-scale whole slide image dataset

####  4.2.1 Dataset

9,200 WSI were downloaded from The Cancer Genome Atlas from all available
solid cancer subtypes, resulting in 5.57Tb of data. These WSI were randomly
selected among the slides with some neoplastic tissue, i.e. by not considering
normal slides which provide no cue about a cancer subtype given their location
since they have no tumor material. The number WSI per subtype and per location
is reported in supplementary Table 3. All WSI were tiled into 512 pixels width
tiles with 128 pixels overlap on both sides at $10\times$ magnification using
code from [11]. The cohort was split on a patient-basis in 3597, 2001 and 3602
slides for respectively training, validation and testing sets.

####  4.2.2 Implementation details

All benchmarked methods use a ResNet34 architecture [7] pre-trained on
Imagenet [2] as the instance embedding $f_{\theta_{1}}$, and use one fully
connected layer for classifier (mapping each bag embedding to a probability).
For embedding-based methods, the last fully connected layer of ResNet34 was
removed resulting in $512$ output channels per tile instead of $1$
probability. During training, 200 randomly cropped $128\times 128$ pixel tiles
were randomly sampled within each WSI. Data augmentation was the same as in
the first experiment (subsection 4.1). Adam [10] was used with cross-entropy
loss, learning rate of $1e^{-4}$ and batch size of 10 WSI (or 2000 taking into
account the number of tiles per WSI). Due to significant imbalance in the
class distribution, oversampling was employed during training with frequencies
equal to the inverse of the classes counts. For each method, the training
process lasted approximately 1 week on 2 Nvidia V100. SparseConvMIL used two
convolutional layers with varying numbers of channels. Downsampling was 128.
Attention-based also used varying hidden neurons as specified.

Method | #Params | Accuracy | Precision | F1-score |  AUC |  CE
---|---|---|---|---|---|---
Random | N/A | 0.031 | 0.031 | 0.031 |  0.500 |  3.506
Instance+max† | 16k | 0.417 | 0.365 | 0.360 |  0.879 |  2.027
Instance+mean† | 16k | 0.463 | 0.417 | 0.414 |  0.905 |  1.783
Instance+LSE† | 16k | 0.451 | 0.406 | 0.403 |  0.898 |  1.819
max† | 70k | 0.441 | 0.434 | 0.403 |  0.913 |  1.821
mean† | 70k | 0.488 | 0.463 | 0.456 |  0.917 |  1.604
Attention-128 [8] | 82k | 0.481 | 0.448 | 0.449 |  0.913 |  1.619
Attention-512 [8] | 279k | 0.487 | 0.451 | 0.453 |  0.912 |  1.616
Attention-2048 [8] | 1055k | 0.472 | 0.452 | 0.452 |  0.909 |  1.621
Gated-Attention-128 [8] | 147k | 0.492 | 0.452 | 0.456 |  0.916 |  1.613
Gated-Attention-512 [8] | 542k | 0.487 | 0.447 | 0.450 |  0.911 |  1.629
Gated-Attention-2048 [8] | 2105k | 0.483 | 0.457 | 0.459 |  0.911 |  1.604
Graph-CNN [18] | 789k | 0.464 | 0.439 | 0.436 |  0.907 |  1.673
Proposed-c32,c32 | 157k | 0.523 | 0.508 | 0.504 |  0.935 |  1.386
Proposed-c128,c128 | 742k | 0.568 | 0.568 | 0.553 |  0.944 |  1.267
Proposed-c256,c256 | 1772k | 0.567 | 0.571 | 0.560 |  0.948 |  1.221
Table 2: Results of the 32 classes classification on the TCGA dataset. Params
are the number of pooling parameters (in thousand). Accuracy is balanced. CE
stands for cross-entropy. Random is the random performance. Otherwise
specified, methods are embedding-level MIL. † denotes parameters-unscalable
pooling methods.

#### Results

Table 2 reports results for several metrics computed by averaging one-vs-all
metrics for each class. The proposed approach achieved superior results in all
metrics and for all configurations, even the small one with 32 channels
convolutions. Furthermore, SparseConvMIL seemed to scale better than
attention-based MIL. Indeed as inferred in Table 2, increasing the parameters
of SparseConvMIL improved significantly performance. On the contrary,
attention-based MIL stagnated with additional parameters.

(a) varying the downsampling factor

(b) varying the number of sampled tiles

Figure 2: Performance of SparseConvMIL-c32,c32 by varying the input sparsity

We aimed at understanding the limitations of the proposed approach. 2(a) plots
the performance of SparseConvMIL-c32,c32 with various downsampling. In
particular, the method performed significantly worse for low downsampling,
which we conjecture is due to the fact that tiles are too far apart to exploit
spatial context with convolution. Furthermore, its performance decreased for
high downsamplings, e.g. 256 or 512, which is probably due to uncoalesced
sparse maps, where there may be duplicate coordinates for several tiles
provoking a loss of input tiles. 2(b) plots the performance by varying the
number of tiles sampled per WSI for downsampling of 128. The performance of
SparseConvMIL increased with respect to the number of sampled tiles, which is
not surprising since more tiles provide more information about WSI.

##  5 Conclusion

In this paper, we proposed a flexible and promising sparse-input convolutional
MIL approach for classifying gigapixel images. We showed on one small-scale
dataset that our method performs on par or better than other approaches for
multiple metrics. We also demonstrated that the proposed approach achieves
significant better performances for all metrics on very large-scale pan-cancer
whole slide images classification from the public dataset The Cancer Genome
Atlas. The limitations of our method have been measured with respect to the
level of sparsity in input data.

There are several ways to iterate over this work. First, CNN visualization
techniques, e.g. GradCAM [16], can be employed to decipher the key instances
and improve interpretability of the framework. Besides, more sophisticated
sparse-input CNN architectures could further improve performance, e.g. with
more layers and residual connections [7].

## References

  * [1] Angelidis, S., Lapata, M.: Multiple instance learning networks for fine-grained sentiment analysis. Transactions of the Association for Computational Linguistics 6, 17–31 (2018)
  * [2] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [3] Dietterich, T.G., Lathrop, R.H., Lozano-Pérez, T.: Solving the multiple instance problem with axis-parallel rectangles. Artificial intelligence 89(1-2), 31–71 (1997)
  * [4] Engelcke, M., Rao, D., Wang, D.Z., Tong, C.H., Posner, I.: Vote3deep: Fast object detection in 3d point clouds using efficient convolutional neural networks. In: 2017 IEEE International Conference on Robotics and Automation (ICRA). pp. 1355–1361. IEEE (2017)
  * [5] Graham, B.: Sparse 3d convolutional neural networks. arXiv preprint arXiv:1505.02890 (2015)
  * [6] Graham, B., van der Maaten, L.: Submanifold sparse convolutional networks. arXiv preprint arXiv:1706.01307 (2017)
  * [7] He, K., Zhang, X., Ren, S., Sun, J.: Identity mappings in deep residual networks. In: European conference on computer vision. pp. 630–645. Springer (2016)
  * [8] Ilse, M., Tomczak, J., Welling, M.: Attention-based deep multiple instance learning. In: International conference on machine learning. pp. 2127–2136. PMLR (2018)
  * [9] Keeler, J.D., Rumelhart, D.E., Leow, W.K.: Integrated segmentation and recognition of hand-printed numerals. Microelectronics and Computer Technology Corporation (1991)
  * [10] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [11] Lerousseau, M., Vakalopoulou, M., Classe, M., Adam, J., Battistella, E., Carré, A., Estienne, T., Henry, T., Deutsch, E., Paragios, N.: Weakly supervised multiple instance learning histopathological tumor segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 470–479. Springer (2020)
  * [12] Marusyk, A., Polyak, K.: Tumor heterogeneity: causes and consequences. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer 1805(1), 105–117 (2010)
  * [13] Ramon, J., De Raedt, L.: Multi instance neural networks. In: Proceedings of the ICML-2000 workshop on attribute-value and relational learning. pp. 53–60 (2000)
  * [14] Riegler, G., Osman Ulusoy, A., Geiger, A.: Octnet: Learning deep 3d representations at high resolutions. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 3577–3586 (2017)
  * [15] Ruifrok, A.C., Johnston, D.A., et al.: Quantification of histochemical staining by color deconvolution. Analytical and quantitative cytology and histology 23(4), 291–299 (2001)
  * [16] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-cam: Visual explanations from deep networks via gradient-based localization. In: Proceedings of the IEEE international conference on computer vision. pp. 618–626 (2017)
  * [17] Sirinukunwattana, K., Raza, S.E.A., Tsang, Y.W., Snead, D.R., Cree, I.A., Rajpoot, N.M.: Locality sensitive deep learning for detection and classification of nuclei in routine colon cancer histology images. IEEE transactions on medical imaging 35(5), 1196–1206 (2016)
  * [18] Tu, M., Huang, J., He, X., Zhou, B.: Multiple instance learning with graph neural networks. arXiv preprint arXiv:1906.04881 (2019)
  * [19] Warburg, O.: The metabolism of carcinoma cells. The Journal of Cancer Research 9(1), 148–163 (1925)
  * [20] Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., Philip, S.Y.: A comprehensive survey on graph neural networks. IEEE transactions on neural networks and learning systems (2020)
  * [21] Yi, Y., Lin, M.: Human action recognition with graph-based multiple-instance learning. Pattern Recognition 53, 148–162 (2016)
  * [22] Zarella, M.D., Bowman, D., Aeffner, F., Farahani, N., Xthona, A., Absar, S.F., Parwani, A., Bui, M., Hartman, D.J.: A practical guide to whole slide imaging: a white paper from the digital pathology association. Archives of pathology & laboratory medicine 143(2), 222–234 (2019)
  * [23] Zhao, Y., Yang, F., Fang, Y., Liu, H., Zhou, N., Zhang, J., Sun, J., Yang, S., Menze, B., Fan, X., et al.: Predicting lymph node metastasis using histopathological images based on multiple instance learning with deep graph convolution. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4837–4846 (2020)
  * [24] Zhou, Z.H., Sun, Y.Y., Li, Y.F.: Multi-instance learning by treating instances as non-iid samples. In: Proceedings of the 26th annual international conference on machine learning. pp. 1249–1256 (2009)


Astronomy Reports, 2023, Vol. , No.

# Schwinger Pair Production and Vacuum Birefringence around High Magnetized
Neutron Stars111Paper presented at the Fifth Zeldovich meeting, an
international conference in honor of Ya. B. Zeldovich held in Yerevan, Armenia
on June 12–16, 2023. Published by the recommendation of the special editors:
R. Ruffini, N. Sahakyan and G. V. Vereshchagin.

Chul Min Kim<EMAIL_ADDRESS>Advanced Photonics Research Institute, Gwangju
Institute of Science and Technology, Gwangju 61005, Korea  Center for
Relativistic Laser Science, Institute for Basic Science, Gwangju 61005, Korea
Sang Pyo Kim<EMAIL_ADDRESS>Department of Physics, Kunsan National
University, Gunsan 54150, Korea  Asia Pacific Center for Theoretical Physics,
Pohang 37673, Korea  ELI Beamlines, ELI-ERIC, Dolní Břežany, Czech Republic

###### Abstract

Received: ; Revised: ; Accepted: .

Highly magnetized neutron stars have magnetic fields of order of the critical
field and can lead to measurable QED effects. We consider the Goldreich-Julian
pulsar model with supercritical magnetic fields, induced subcritical electric
fields, and a period of milliseconds. We then study the strong field physics,
such as Schwinger pair production and the vacuum birefringence including the
wrench effect, whose X-ray polarimetry will be observed in future space
missions.

##  I Introduction

Ultra-strong electromagnetic (EM) fields abound in astrophysics. Neutron stars
have strong magnetic fields that are stable for long duration and extend to
astrophysical scale (for review, see Shapiro and Teukolsky (1983) and also
Meszaros (1992)), whose strength goes beyond the present laboratory
achievements, $I=1.1\times 10^{23}\,{\rm W/cm^{2}}$ Yoon et al. (2021).
Chandrasekhar and Fermi Chandrasekhar and Fermi (1953) derived the upper bound
for the magnetic fields from the virial theorem of magnetohydrostatic
equilibrium that the magnetic energy of compact object can not exceed the
gravitational binding energy, $B\leq 10^{18}\bigl{(}M_{\rm
N}/1.4M_{\odot}\bigr{)}\bigl{(}10\,{\rm km}/R_{\rm N}\bigr{)}^{2}\,{\rm G}$.
The dynamo can further amplify the magnetic field beyond $10^{15}\,{\rm G}$,
though the flux conservation during gravitational collapse of stars gives the
magnetic fields about $10^{12}\,{\rm G}$ Thompson and Duncan (1993). Recently,
magnetars, highly magnetized neutron stars, have been observed Vasisht and
Gotthelf (1997), and those with supercritical magnetic fields and their
astrophysical properties are listed in McGill Catalog Olausen and Kaspi (2014)
and Kaspi and Beloborodov (2017).

It has been known since the QED action was obtained by Heisenberg-Euler and
Schwinger (HES) Heisenberg and Euler (1936); Schwinger (1951) that strong EM
fields make the Dirac vacuum polarized as media, and that strong electric
fields spontaneously create pairs of electrons and positrons. The HES action
results in nonlinear electrodynamics that modifies the linear Maxwell theory.
The criteria for necessity of nonlinear electrodynamics are given as follows
(in the Lorentz-Heaviside units with $c=\hbar=1$): (i) in plasma and matter
systems, the magnetic or electric length $(1/\sqrt{eB},1/\sqrt{eE})$ for
charge $e$ is comparable to or shorter than the characteristic length scale of
physical systems, (ii) in the Minkowski spacetime, the characteristic length
is the Compton wavelength $(1/\sqrt{eB_{c}}=1/\sqrt{eE_{c}}=1/m)$ and QED
loops contribute to the EM theory. The critical electric and magnetic fields
$E_{c}=m^{2}/e=B_{c}$ for a charge $e$ with mass $m$ have the electrostatic
potential energy over a Compton wavelength or the Landau levels equal to the
rest mass: $E_{c}=1.3\times 10^{16}\,{\rm V/cm}$ and $B_{c}=4.4\times
10^{13}\,{\rm G}$ or the intensity $I_{c}=2.3\times 10^{29}\,{\rm W/cm^{2}}$.
When electric or magnetic fields become comparable to or stronger than the
critical field, the QED loop corrections should be included for accuracy.

QED loop corrections result in nonlinear electrodynamics in which the vacuum
acquires field-dependent permittivity and permeability, in particular, the
magneto-electric response. In the presence of strong electric fields, the
vacuum decays due to spontaneous creation of electron-positron pairs. The
nonlinear theory that respects the Lorentz- and gauge-invariance should be a
functional of the Maxwell scalar ${\cal F}$ and the pseudo-scalar ${\cal G}$:

$\displaystyle{\cal
F}=\frac{1}{4}F_{\mu\nu}F^{\mu\nu}=\frac{1}{2}(\vec{B}^{2}-\vec{E}^{2}),\quad{\cal
G}=\frac{1}{4}F_{\mu\nu}F^{*\mu\nu}=-\vec{E}\cdot\vec{B}.$ (1)

In this paper we will consider the dipole pulsar model by Goldreich and Julian
(GJ) Goldreich and Julian (1969), as shown in Fig. 1(a), whose ${\cal F}$ and
${\cal G}$ are

$\displaystyle{\cal F}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}}{8}\Bigl{(}\frac{R_{N}}{r}\Bigr{)}^{6}\Bigl{[}1+3\cos^{2}\theta_{n}-\Bigl{(}\frac{R_{N}}{r}\Bigr{)}^{2}\Bigl{(}\frac{R_{N}\Omega}{c}\Bigr{)}^{2}\bigl{(}1-2\cos^{2}\theta_{n}+5\cos^{4}\theta_{n}\bigr{)}\Bigr{]},$
$\displaystyle{\cal G}$ $\displaystyle=$ $\displaystyle
B_{0}^{2}\Bigl{(}\frac{R_{N}}{r}\Bigr{)}^{7}\Bigl{(}\frac{R_{N}\Omega}{c}\Bigr{)}\cos^{3}\theta_{n}.$
(2)

Here, $R_{N}$ and $\Omega$ are the radius and angular frequency of a neutron
star, and $\theta_{n}$ is the polar angle from the north pole (for simplicity,
the magnetic polar axis assumed to be along the rotational axis, but see Ref.
Turolla et al. (2015) for realistic configurations of magnetic fields). Note
that the maximal surface velocity of neutron star (in units of $c$) can be
written as

$\displaystyle\frac{R_{N}\Omega}{c}=\frac{\pi}{15}\times\Bigl{(}\frac{R_{N}}{R_{10}}\Bigr{)}\Bigl{(}\frac{P_{-3}}{P}\Bigr{)},$
(3)

where $R_{10}=10\,\mathrm{km}$ and $P_{-3}=10^{-3}\,\mathrm{sec}$ refers to
millisecond neutron stars. Most neutron stars have
$(R_{N}/R_{10})(P_{-3}/P)\leq 1$, so the maximal surface velocity does not
violate the special relativity. On the equatorial plane $(\theta_{n}=\pi/2)$,
${\cal G}=0$, but
$E^{2}(r)/B^{2}(r)=\bigl{(}{R_{N}}/{r}\bigr{)}^{2}\bigl{(}{R_{N}\Omega}/{c}\bigr{)}^{2}$.
In general, ${\cal G}/{\cal
F}<2\bigl{(}{R_{N}}/{r}\bigr{)}\bigl{(}{R_{N}\Omega}/{c}\bigr{)}$, so we may
assume ${\cal G}\ll{\cal F}$ away from the neutron stars.

(a)

(b)

(c)

Figure 1:  Electromagnetic fields around a pulsar. (a) EM fields of the
Goldreich-Julian dipole model for pulsars. The red arrows denote the magnetic
field (in units of $B_{0}$), and the blue ones denote the electric field (in
units of $B_{0}R_{N}\Omega/c$). $B_{0}$ refers to the field strength at the
north pole, $R_{N}$ to the pulsar’s radius, and $\Omega$ to the angular
frequency of rotation; (b) $\mathcal{F}(r,\theta_{n})$ and (c)
$\mathcal{G}(r,\theta_{n})$ in units of $B_{0}^{2}$ for
$R_{N}\Omega/c=\pi/15/1000$, e.g., a neutron star with a radius of 10 km and a
period of 1 s.

The purpose of this paper is to further explore strong QED effects: the vacuum
polarization and Schwinger pair production in supercritical magnetic fields
and subcritical electric fields, motivated by highly magnetized neutron stars,
in particular, magnetars, which have recently been intensively studied in
Refs. Kim and Kim (2021, 2023). To do so, we find a new analytical expression
for HES action by keeping supercritical magnetic fields in Refs. Dittrich
(1976); Dittrich and Reuter (1985); Kim (2011); Kim and Lee (2019) but by
expanding the part for the subcritical fields. We then study the vacuum
birefringence for photons in EM fields of neutron stars. This study can be
used in space missions that have recently been proposed to probe the QED
regime, such as Compton Telescope Wadiasingh et al. (2019) and the enhanced
X-ray Timing and Polarimetry (eXTP) Santangelo et al. (2019). Physical
processes in supercritical magnetic fields combined with induced electric
fields differ from those in weak-field regime, for instance the curvature
synchrotron radiation and QED cascade of pairs Meszaros (1992); Battesti and
Rizzo (2013); Harding and Lai (2006), not to mention the electron-positron
pair production in astrophysics Ruffini et al. (2010).

##  II QED Action in a Supercritical Magnetic Field and a Subcritical Electric
Field

The vacuum polarization and Schwinger pair production in a constant EM field
has been described by Heisenberg, Euler, and Schwinger’s exact one-loop QED
action

$\displaystyle{\mathcal{L}^{(1)}}(a,b)=-\frac{1}{8\pi^{2}}\int_{0}^{\infty}ds\frac{e^{-m^{2}s}}{s^{3}}\Bigl{[}(eas)\coth(eas)(ebs)\cot(ebs)-1-\frac{(es)^{2}}{3}(a^{2}-b^{2})\Bigr{]},$
(4)

where $a$ and $b$ are the gauge and Lorentz invariants

$\displaystyle a=\sqrt{\sqrt{{\cal F}^{2}+{\cal G}^{2}}+{\cal F}},\quad
b=\sqrt{\sqrt{{\cal F}^{2}+{\cal G}^{2}}-{\cal F}}.$ (5)

The QED action respects the CP invariance, ${\cal L}^{(1)}(\pm a,\pm b)={\cal
L}^{(1)}(a,b)$, and the EM duality, ${\cal L}^{(1)}(a,b)={\cal
L}^{(1)}(ib,-ia)$ Kim (2011); Kim and Lee (2019).

Given an EM field configuration, $a$ and $b$ characterize in a gauge- and
Lorentz-invariant way the magnetic field strength and the electric field
strength, respectively. When ${\mathcal{G}}=0$, called electromagnetically
wrenchless, a Lorentz frame can be found to remove the weaker between the
electric field and the magnetic field: $a=B,b=0$ when $\mathcal{F}>0$ and
$a=0,b=E$ when $\mathcal{F}<0$. When ${\mathcal{G}}\neq 0$, called
electromagnetically wrenched, a Lorentz frame can be found to make the
magnetic field parallel to the electric field: $a=B,b=E$. Since $a$ and $b$
determine the effective action and thereby the resultant phenomena, we may
regard $a$ the intrinsic magnetic field strength and $b$ the intrinsic
electric field strength of a given field configuration, irrespective of
reference frames. Hence, the configuration of a supercritical magnetic field
and a subcritical electric field is unequivocally represented by the condition
$a>E_{c}>b$.

In a supercritical magnetic field with or without a subcritical electric
field, (4) can be used for a numerical evaluation. Hence, it will be helpful
to have an analytical expression for (4) in terms of known special functions.
To do so, we expand $(ebs)\cot(ebs)$ in (4) using the power series 1.411 in
Gradshteĭn and Zwillinger (2015)

$x\cot x=1-\sum_{k=1}^{\infty}\frac{2^{2k}|B_{2k}|}{(2k)!}x^{2k},$ (6)

where $B_{2k}$ are the Bernoulli numbers. The leading term in the expansion of
${\cal L}^{(1)}(a,b)$, i.e., ${\cal L}^{(1)}_{(0)}(a,0)$, was obtained in a
closed form by Dittrich through dimensional regularization Dittrich (1976);
Dittrich and Reuter (1985) and also by Kim and Lee using the gamma function
regularization of the in-out formalism Kim and Lee (2019):

$\displaystyle{\cal
L}^{(1)}_{(0)}(a,0)=\frac{m^{4}}{8\pi^{2}\bar{a}^{2}}\left[\zeta^{\prime}(-1,\bar{a})-\frac{1}{12}+\frac{\bar{a}^{2}}{4}-\left(\frac{1}{12}-\frac{\bar{a}}{2}+\frac{\bar{a}^{2}}{2}\right)\ln\bar{a}\right].$
(7)

Here, $\bar{a}={m^{2}}/{2ea}$,
$\zeta^{\prime}(-1,\bar{a})=d\zeta(s,\bar{a})/ds|_{s=-1}$, and
$\zeta(s,\bar{a})$ is the Hurwitz zeta function. The next leading term,
proportional to $(eb)^{2}$, is

$\displaystyle{\cal
L}^{(1)}_{(1)}(a,b)=\frac{(eb)^{2}}{24\pi^{2}}\left[-\frac{1}{2\bar{a}}+\ln{\bar{a}}-\psi(\bar{a})\right],$
(8)

where $\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function. The
higher order terms are given by

$\displaystyle{\cal L}^{(1)}_{(\geq
2)}(a,b)=\frac{(ea)(eb)}{4\pi^{2}}\sum_{k\geq
2}^{\infty}\frac{|B_{2k}|}{(2k)(2k-1)}\left[2\Bigl{(}\frac{b}{a}\Bigr{)}^{2k-1}\zeta\bigl{(}2k-1,\bar{a}\bigr{)}-\Bigr{(}\frac{2eb}{m^{2}}\Bigl{)}^{2k-1}\right].$
(9)

Here, we have used the integral formula 3.551 in Gradshteĭn and Zwillinger
(2015) and the definition of $\Gamma$-function:

$\int_{0}^{\infty}\mathrm{d}x\ x^{\mu-1}e^{-\beta x}\coth
x=\Gamma(\mu)\left[2^{1-\mu}\zeta\left(\mu,\frac{\beta}{2}\right)-\beta^{-\mu}\right]\quad\left[\mathrm{Re}(\mu)>1,\
\mathrm{Re}(\beta)>0\right].$ (10)

Note that the QED action was obtained in powers of ${\cal G}$ by Heyl and
Hernquist Heyl and Hernquist (1997a, b).

The importance of the correction terms ${\cal L}^{(1)}_{(1)}$ and ${\cal
L}^{(1)}_{(2)}$ with respect to ${\cal L}^{(1)}_{(0)}$ is exemplified in Fig.
2 for $b=0.1E_{c},~{}0.01E_{c}$. We may take the ratio of 0.1 as the threshold
of small corrections. When $b=0.1E_{c}$, the first-order correction ${\cal
L}^{(1)}_{(1)}$ is small compared to ${\cal L}^{(1)}_{(0)}$ for $a\geq
0.66B_{c}$, whereas for $b=0.01E_{c}$, it is small for $a\geq 0.068B_{c}$.
These results empirically imply that $a\geq 7b$ is necessary to keep the
correction terms less than 10% of the leading term when $b\leq 0.1E_{c}$,
which is consistent with the assumption of a supercritical magnetic field and
a subcritical electric field. When ${\cal L}^{(1)}_{(1)}$ can be considered
small, the correction ${\cal L}^{(1)}_{(2)}$ can be neglected as it is smaller
than ${\cal L}^{(1)}_{(1)}$ by two orders, as shown in Fig. 2. Thus, including
only the subleading correction ${\cal L}^{(1)}_{(1)}$ would be sufficient for
an accurate calculation in the valid parameter range. At a higher value of $b$
closer to $E_{c}$, the constraint on $a$ becomes more demanding: $a>8b$ is
necessary to keep the correction small when $b=0.5E_{c}$. In such a case,
however, the pair production due to the imaginary part of ${\cal
L}^{(1)}(a,b)$ becomes significant and thus should be taken into account (see
the next section).

Figure 2: ${\cal L}^{(1)}_{(1)}(a,b)/{\cal L}^{(1)}_{(0)}(a,b)$ (denoted by
$r_{1}$) and ${\cal L}^{(1)}_{(2)}(a,b)/{\cal L}^{(1)}_{(0)}(a,b)$ (denoted by
$r_{2}$) for $b=0.1E_{c},~{}0.01E_{c}$. The ratio of 0.1, marked by the red
thick horizontal line, may be considered as the threshold of small
corrections.

Remarkably, Ritus obtained the QED action in supercritical magnetic field up
to the two-loop level Ritus (1976)

$\displaystyle{\cal L}^{(1)}+{\cal
L}^{(2)}=\frac{(eB)^{2}}{24\pi^{2}}\Bigl{[}\Bigl{(}1+\frac{3\alpha}{4\pi}\Bigr{)}\ln\Bigl{(}\frac{eB}{\pi\gamma
m^{2}}\Bigr{)}+\frac{6}{\pi^{2}}\zeta^{\prime}(2)+\cdots\Bigr{]}.$ (11)

Thus the two-loop action is smaller by a factor of $\alpha$ than the one-loop
action, which implies the one-loop action can be used even for supercritical
magnetic fields up to this small deviation. In reality, an electric field can
hardly exceed the critical strength because of the electrical short induced by
Schwinger pair production, near the critical field. The dominant one-particle
irreducible diagrams of all loops to the HES action by Karbstein Karbstein
(2019) replace $\alpha$ by $\alpha/(1-(\alpha/3\pi)\ln(eB/m^{2}))$, which is
illustrative in astrophysics applications.

##  III Schwinger Pair Production

The Schwinger pair production and vacuum decay originates from the poles of
HES action (4). The pair production rate per unit Compton four-volume is

$\displaystyle
R(r,\theta_{n};B_{0}/B_{c})=\frac{\tilde{a}\tilde{b}m^{4}}{8\pi^{2}}\coth\Bigl{(}\pi\frac{\tilde{a}}{\tilde{b}}\Bigr{)}e^{-\pi/\tilde{b}},$
(12)

where $\tilde{a}=a/B_{c}$ and $\tilde{b}=b/B_{c}$ are dimensionless
parameters. The vacuum decay rate is given by the vacuum persistence amplitude

$\displaystyle 2\,\mathrm{Im}{\cal
L}^{(1)}=\frac{\tilde{a}\tilde{b}m^{4}}{4\pi^{2}}\sum_{n=1}^{\infty}\coth\Bigl{(}\pi\frac{\tilde{a}}{\tilde{b}}\Bigr{)}e^{-\pi
n/\tilde{b}},$ (13)

and in the case of $\tilde{b}=0$ the sum can be expressed in terms of the
polylogarithm function
$(\tilde{a}m^{2})^{2}/(4\pi^{2})L_{2}(e^{-\pi/\tilde{b}})$.

We consider only the EM field of the GJ model with (I), in which ${\cal F}$
$({\cal G})$ is symmetric (antisymmetric) with respect to the equatorial
plane. Furthermore, the parameters $\tilde{a}$ and $\tilde{b}$ are
proportional to $B_{0}/B_{c}$ and depend only the distance $r$ from the
surface and the polar angle $\theta_{n}$. In this paper, we study a
theoretical model for pulsars with the millisecond period and supercritical
magnetic fields. The pair production is dominant in regions near the north and
south poles because ${\cal G}<{\cal F}$ on the surface and ${\cal G}\ll{\cal
F}$ away from the surface, so $\tilde{b}\simeq|{\cal
G}/B_{c}^{2}|/\sqrt{2{\cal F}/B_{c}^{2}}$.

As shown in Figs. 1(b) and 1(c), ${\cal F}$ is maximal along polar axis and
decreases away from the neutron star surface, and ${\cal G}$ vanishes on the
equatorial plane and rapidly increases near the polar axis. This means that
Schwinger pair production is significant close to the polar axis, as shown in
Fig. 3 (left panel). Because of the electric field direction at poles in Fig.
1(a), near the north pole, positrons are emitted outward and electrons fall on
the neutron star, and vice versa near the south pole. The emission forms
“jets” of different charges along polar axis. The neutron star remains charge
neutral because the pair production rate is symmetric with respect to the
equatorial plane. The pair production rate increase by two order when the
magnetic field strength increases from $B_{0}/B_{c}=40$ to $B_{0}/B_{c}=50$
and the rate rapidly decreases away from the surface of neutron star as shown
in Fig. 3 (right panel).

The mean number of electrons or positrons must be multiplied by a huge factor
from the astrophysical scale in units of the Compton four-volume
$m^{4}=1/\lambda_{\rm C}^{3}T_{\rm C}=4.5\times 10^{57}/{\rm
m^{3}\,sec}=4.5\times 10^{66}/{\rm km^{3}\,sec}$ at the north or south pole
region. Roughly, one expects the number of electrons or positron on a pole per
second ${\cal N}/{\mathrm{sec}}=4.5\times 10^{66}\times N$, so ${\cal
N}/{\mathrm{sec}}\sim 10^{20}$ even when $N=10^{-44}$ which is produced by
$B_{0}/B_{c}=6.08$.

Figure 3: The pair production rate per unit Compton four-volume against
$r/R_{N}$ and $\cos\theta_{n}$ when $B_{0}/B_{c}=50$ [left] and at poles
against $r/R_{N}$ and $B_{0}/B_{c}$ [right].

##  IV Vacuum birefringence for low-energy photons

The one-loop effective action $\mathcal{L}^{(1)}$ renders the vacuum under
strong EM fields into an EM medium, of which total Lagrangian is given as

$\mathcal{L}=-\mathcal{F}+\mathcal{L}^{(1)}.$ (14)

Then, the polarization ($\vec{P}$) and magnetization ($\vec{M}$), a response
of the medium, can be obtained by differentiating $\mathcal{L}$ with respect
to $\vec{E}$ and $\vec{B}$ Greiner and Reinhardt (2009):

$\vec{D}=\vec{E}+\vec{P}=\frac{\partial\mathcal{L}}{\partial\vec{E}},\quad\vec{H}=\vec{B}-\vec{M}=-\frac{\partial\mathcal{L}}{\partial\vec{B}},$
(15)

which, for $\mathcal{L}=\mathcal{L}(\mathcal{F},\mathcal{G})$, are reduced to

$\vec{D}=-\mathcal{L}_{\mathcal{F}}\vec{E}-\mathcal{L}_{\mathcal{G}}\vec{B},\quad\vec{H}=-\mathcal{L}_{\mathcal{F}}\vec{B}+\mathcal{L}_{\mathcal{G}}\vec{E}.$
(16)

The subscripts denote the differentiation with respect to the variable. These
expressions imply that the vacuum is a nonlinear medium as $\mathcal{L}$ is a
complicated function of $\vec{E}$ and $\vec{B}$.

We are interested in the propagation of a photon in such a nonlinear vacuum.
When the photon energy is much smaller than $m$, the photon cannot produce
electron-positron pairs by interacting with the strong background field. Then,
the photon can be treated as a weak classical perturbation to the background
field. Splitting all the fields in (16) into a strong background and a weak
perturbation, $\vec{A}=\vec{A}_{0}+\delta\vec{A}$, we can obtain the
constitutive relation for the photon field:

$\delta\vec{D}=\tensor{\epsilon}_{E}\cdot\delta\vec{E}+\tensor{\epsilon}_{B}\cdot\delta\vec{B},\quad\delta\vec{H}=\tensor{\bar{\mu}}_{B}\cdot\delta\vec{B}+\tensor{\bar{\mu}}_{E}\cdot\delta\vec{E}.$
(17)

It is assumed that the constitutive relations (16) hold separately for the
background quantities and perturbation quantities. The tensor coefficients
$\tensor{\epsilon}_{E}$, $\tensor{\epsilon}_{B}$, $\tensor{\bar{\mu}}_{B}$,
and $\tensor{\bar{\mu}}_{E}$ are given in terms of $\vec{E}_{0}$,
$\vec{B}_{0}$, and $\mathcal{L}_{\mathcal{F}}$, $\mathcal{L}_{\mathcal{G}}$,
$\mathcal{L}_{\mathcal{FF}}$, $\mathcal{L}_{\mathcal{FG}}$,
$\mathcal{L}_{\mathcal{GG}}$, evaluated for $\vec{E}=\vec{E}_{0}$ and
$\vec{B}=\vec{B}_{0}$ Kim and Kim (2022). Note that $\tensor{\epsilon}_{B}$
and $\tensor{\bar{\mu}}_{E}$ give rise to the magneto-electric effect, which
is rarely observed in usual optical media but has recently been observed in
multiferroic materials Eerenstein et al. (2006); Fiebig (2005).

Then, the Maxwell equations for the photon field are obtained as

$\frac{\partial\delta\vec{B}}{\partial
t}=-\nabla\times\delta\vec{E},\quad\frac{\partial\delta\vec{D}}{\partial
t}=\nabla\times\delta\vec{H}.$ (18)

To find the photon propagation modes, we assume that all the perturbation
quantities vary harmonically as $\exp(i\omega(n\hat{n}\cdot\vec{x}-t))$, where
$\hat{n}$ is the direction of photon propagation, $\omega$ the photon
frequency, and $n$ the refractive index to be determined. Then, the Maxwell
equations are reduced to

$\delta\vec{B}=\vec{n}\times\delta\vec{E},\quad\delta\vec{D}=-\vec{n}\times\delta\vec{H},$
(19)

where $\vec{n}=n\hat{n}$. Combining (17) and (19) yields a homogeneous matrix-
vector equation, whose solution is the polarization vector $\delta\vec{E}$ and
the refractive index $n$:

$\tensor{\Lambda}\cdot\delta\vec{E}=\tensor{\epsilon}_{E}\cdot\delta\vec{E}+\tensor{\epsilon}_{B}\cdot\vec{n}\times\delta\vec{E}+\vec{n}\times(\tensor{\bar{\mu}}_{B}\cdot\vec{n}\times\delta\vec{E}+\tensor{\bar{\mu}}_{E}\cdot\delta\vec{E})=0$
(20)

Albeit complicated, the equation can be solved exactly through a systematic
procedure Kim and Kim (2022).

For example, we consider a weak, purely magnetic background field, for which
the refractive indices and polarization vectors are obtained in simple forms.
In such a case,
$\mathcal{L}=-\mathcal{F}+\eta_{1}\mathcal{F}^{2}+\eta_{2}\mathcal{G}^{2}$ ,
where $\eta_{1}/4=\eta_{2}/7=e^{4}/(360\pi^{2}m^{4})$ Adler (1971). Assuming
$\vec{B}_{0}=B_{0}\hat{z}$ and $\hat{n}=(\sin\theta,0,\cos\theta)$ without
loss of generality, we obtain the modes as follows Kim and Kim (2022):

$\begin{split}n_{1}=\sqrt{\frac{1-\eta_{1}B_{0}^{2}}{1-\eta_{1}B_{0}^{2}-2\eta_{1}B_{0}^{2}\sin^{2}\theta}}\;\mathrm{with}\;\delta\vec{E}_{1}&=(0,1,0)\\\
n_{2}=\sqrt{\frac{1-\eta_{1}B_{0}^{2}+2\eta_{2}B_{0}^{2}}{1-\eta_{1}B_{0}^{2}+2\eta_{2}B_{0}^{2}\cos^{2}\theta}}\;\mathrm{with}\;\delta\vec{E}_{2}&=\left(-\frac{(1-\eta_{1}B_{0}^{2}+2\eta_{2}B_{0}^{2})}{1-\eta_{1}B_{0}^{2}}\cot\theta,0,1\right).\end{split}$
(21)

As $n_{1}\neq n_{2}$, birefringence occurs. Note that
$\hat{n}\cdot\delta\vec{E}_{1}=\delta\vec{E}_{1}\cdot\delta\vec{E}_{2}=0$ but
$\hat{n}\cdot\delta\vec{E}_{2}\neq 0$.

When the background field is strong,
$\mathcal{L}^{(1)}=\mathcal{L}^{(1)}_{(0)}+\mathcal{L}^{(1)}_{(1)}+\cdots$ in
Sec. II should be used to obtain the modes. For the same configuration of
$\vec{B}_{0}$ and $\hat{n}$ as in the weak field case, the modes are given as
follows Kim and Kim (2023):

$\begin{split}n_{1}=\sqrt{\frac{1}{1+\mu\sin^{2}\theta}}\;\mathrm{with}\;\delta\vec{E}_{1}&=(0,1,0)\\\
n_{2}=\sqrt{\frac{1+\epsilon}{1+\epsilon\cos^{2}\theta}}\;\mathrm{with}\;\delta\vec{E}_{2}&=\left(-(1+\epsilon)\cot\theta,0,1\right),\end{split}$
(22)

where

$\mu=-\frac{B_{0}^{2}\mathcal{L}^{(1)}_{\mathcal{FF}}}{1-\mathcal{L}^{(1)}_{\mathcal{F}}},\quad\epsilon=\frac{B_{0}^{2}\mathcal{L}^{(1)}_{\mathcal{GG}}}{1-\mathcal{L}^{(1)}_{\mathcal{F}}}.$
(23)

As the background field is purely magnetic, one may be tempted to drop the
terms $\mathcal{L}^{(1)}_{\geq 1}$, which is null without an electric field,
from the beginning. However, the fields in $\mathcal{L}$ are the combinations
of the background field and the photon field, which definitely has an electric
part due to its propagating nature. Thus, the expansion presented in Sec. II
is crucial even for the mode analysis with a purely magnetic background field.
If an additional electric field is added to have a component parallel to the
strong magnetic field, the birefringence effect, $|n_{1}-n_{2}|$, increases,
and the polarization vectors rotate Kim and Kim (2023).

The birefringence experienced by low-energy photons around highly magnetized
neutron stars is a characteristic phenomenon of the nonlinear quantum vacuum.
If observed in such environments, it will not only reveal the nature of the
quantum vacuum under strong fields but also estimate the magnetic field around
compact astronomical objects. The formalism developed in this paper and a
systematic mode analysis will provide an accurate framework for such goals.

##  V Conclusion

Ultra-strong magnetic fields of neutron stars induce relatively strong
electric fields due to their rapid rotation from milliseconds to seconds.
Hence the strong QED plays an important role in understanding physics, such as
vacuum polarization effects, electron-positron pair production, photon-photon
interactions, interactions of charges in strong EM fields, particularly, the
acceleration of charges and radiations and plasma QED etc. In this paper we
have studied strong field QED physics in supercritical magnetic fields and
induced subcritical electric fields, and particularly, two of the prominent
features: the Schwinger pair production and the vacuum birefringence.

A theoretical model studied in this paper with a millisecond period and
supercritical magnetic fields leads to the Schwinger pair production, the
dominant emission mechanism of electron-positron pairs due to the
astrophysical scale volume $\mathrm{km}^{3}$ in Compton units around the north
and south poles. This model provides us with the most extreme environments for
strong field QED. However, the confirmed magnetars with supercritical magnetic
fields have the period of second or longer and millisecond neutron stars have
magnetic fields weaker by order than the critical field, so the Schwinger
effect is exponentially suppressed by a huge factor even for the astrophysical
volume. On the other hand, the vacuum birefringence may give rise to
observable polarization effects, which can be measured through X-ray
polarimetry: space telescopes like Compton Telescope and eXTP are planned for
that purpose. The precision observations of X-ray polarimetry will enable one
to reconstruct EM field configurations of neutron stars and charged and/or
magnetized black holes.

###### Acknowledgements.

S.P.K. would like to thank Sergei Bulanov, Tae Moon Jung, Remo Ruffini, Rashid
Shaisultanov, Gregory Vereshchagin for useful discussions. This work was
supported by Institute for Basic Science (IBS) under IBS-R012-D1. The work of
C.M.K. was also supported by Ultrashort Quantum Beam Facility operation
program (140011) through APRI, GIST and GIST Research Institute (GRI) grant
funded by GIST in 2023. The work of S.P.K. was in part supported by National
Research Foundation of Korea (NRF) funded by the Ministry of Education
(2019R1I1A3A01063183).

## References

  * Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky, _Black holes, white dwarfs, and neutron stars : the physics of compact objects_ (Wiley, New York, 1983).
  * Meszaros (1992) P. Meszaros, _High-energy radiation from magnetized neutron stars_ (The University Of Chicago Press, Chicago, 1992).
  * Yoon et al. (2021) J. W. Yoon, Y. G. Kim, I. W. Choi, J. H. Sung, H. W. Lee, S. K. Lee, and C. H. Nam, Optica 8, 630 (2021).
  * Chandrasekhar and Fermi (1953) S. Chandrasekhar and E. Fermi, Astrophys. J. 118, 116 (1953).
  * Thompson and Duncan (1993) C. Thompson and R. C. Duncan, Astrophys. J. 408, 194 (1993).
  * Vasisht and Gotthelf (1997) G. Vasisht and E. Gotthelf, Astrophys. J. Lett. 486, L129 (1997).
  * Olausen and Kaspi (2014) S. A. Olausen and V. M. Kaspi, Astrophys. J. Suppl. Ser. 212, 6 (2014).
  * Kaspi and Beloborodov (2017) V. M. Kaspi and A. M. Beloborodov, Annu. Rev. Astron. Astrophys. 55, 261 (2017).
  * Heisenberg and Euler (1936) W. Heisenberg and H. Euler, Zeitschr. Phys 98, 714 (1936).
  * Schwinger (1951) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Goldreich and Julian (1969) P. Goldreich and W. H. Julian, Astrophys. J. 157, 869 (1969).
  * Turolla et al. (2015) R. Turolla, S. Zane, and A. L. Watts, Rep. Prog. Phys. 78, 116901 (2015).
  * Kim and Kim (2021) C. M. Kim and S. P. Kim, in _17th Italian-Korean Symposium on Relativistic Astrophysics_ (2021), eprint 2112.02460.
  * Kim and Kim (2023) C. M. Kim and S. P. Kim, Eur. Phys. J. C 83, 104 (2023).
  * Dittrich (1976) W. Dittrich, J. Phys. A: Math Gen. 9, 1171 (1976).
  * Dittrich and Reuter (1985) W. Dittrich and M. Reuter, _Effective Lagrangians in Quantum Electrodynamics_ (Springer, Berlin, Heidelberg, 1985).
  * Kim (2011) S. P. Kim, Phys. Rev. D 84, 065004 (2011).
  * Kim and Lee (2019) S. P. Kim and H. K. Lee, J. Korean Phys. Soc. 74, 930 (2019).
  * Wadiasingh et al. (2019) Z. Wadiasingh, G. Younes, M. G. Baring, A. K. Harding, P. L. Gonthier, K. Hu, A. v. d. Horst, S. Zane, C. Kouveliotou, A. M. Beloborodov, et al., Bulletin of the AAS 51 (2019).
  * Santangelo et al. (2019) A. Santangelo, S. Zane, H. Feng, R. Xu, V. Doroshenko, E. Bozzo, I. Caiazzo, F. C. Zelati, P. Esposito, D. Gonzalez-Caniulef, et al., SCIENCE CHINA Physics, Mechanics & Astronomy 62, 1 (2019).
  * Battesti and Rizzo (2013) R. Battesti and C. Rizzo, Rep. Prog. Phys. 76, 016401 (2013).
  * Harding and Lai (2006) A. K. Harding and D. Lai, Rep. Prog. Phys. 69, 2631 (2006).
  * Ruffini et al. (2010) R. Ruffini, G. Vereshchagin, and S.-S. Xue, Phys. Rep. 487, 1 (2010).
  * Gradshteĭn and Zwillinger (2015) I. S. Gradshteĭn and D. Zwillinger, _Table of integrals, series, and products_ (Elsevier, Amsterdam, 2015).
  * Heyl and Hernquist (1997a) J. S. Heyl and L. Hernquist, Phys. Rev. D 55, 2449 (1997a).
  * Heyl and Hernquist (1997b) J. S. Heyl and L. Hernquist, J. Phys. A: Math Gen. 30, 6485 (1997b).
  * Ritus (1976) V. I. Ritus, Sov. Phys. JETP 42, 774 (1976).
  * Karbstein (2019) F. Karbstein, Phys. Rev. Lett. 122, 211602 (2019).
  * Greiner and Reinhardt (2009) W. Greiner and J. Reinhardt, _Quantum Electrodynamics_ (Springer, 2009).
  * Kim and Kim (2022) C. M. Kim and S. P. Kim (2022), eprint 2210.12890.
  * Eerenstein et al. (2006) W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 442, 759 (2006).
  * Fiebig (2005) M. Fiebig, J. Phys. D: Appl. Phys. 38, R123 (2005).
  * Adler (1971) S. L. Adler, Ann. Phys. 67, 599 (1971).