11institutetext: Department of Electrical and Computer Engineering
The University of Iowa, Iowa City, USA
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Globally Optimal Surface Segmentation using Deep Learning with Learnable
Smoothness Priors

Leixin Zhou  Xiaodong Wu

###### Abstract

Automated surface segmentation is important and challenging in many medical
image analysis applications. Recent deep learning based methods have been
developed for various object segmentation tasks. Most of them are a
classification based approach, e.g. U-net, which predicts the probability of
being target object or background for each voxel. One problem of those methods
is lacking of topology guarantee for segmented objects, and usually post
processing is needed to infer the boundary surface of the object. In this
paper, a novel model based on convolutional neural network (CNN) followed by a
learnable surface smoothing block is proposed to tackle the surface
segmentation problem with end-to-end training. To the best of our knowledge,
this is the first study to learn smoothness priors end-to-end with CNN for
direct surface segmentation with global optimality. Experiments carried out on
Spectral Domain Optical Coherence Tomography (SD-OCT) retinal layer
segmentation and Intravascular Ultrasound (IVUS) vessel wall segmentation
demonstrated very promising results.

##  1 Introduction

Automated segmentation of objects or equivalently boundary surfaces plays a
very import role in quantitative image analysis. In several years, deep
learning based method for semantic segmentation has become very popular in
computer vision and medical imaging. The fully convolutional networks (FCN)
[11] , and then U-net [15] for medical image segmentation have been proposed.
All these methods model the segmentation problem as a pixel-wise or voxel-wise
classification problem. Popular loss functions for training these networks
include cross entropy (CE), weighted cross entropy (WCE) for imbalanced
classes, and dice similarity coefficient (DSC). As there is no explicit
constraint among the labeling of different pixels, usually post processing,
such as morphological operation, is required to get reasonable prediction.

It is well known that prior knowledge, like shape priors, can help
significantly improve segmentation performance, especially in medical scenario
the target objects often lack strong boundaries, and/or present similar
intensity profiles with multiple crowed nearby tissues. Some studies have
investigated along this line. For example, Milletari et al. [12] proposed to
integrate into CNN the statistical shape priors, which were obtained from
principal components analysis (PCA). Although the robustness of segmentation
was improved, the mean DSC did not increase. Ravishankar et al. [14] utilized
a shape regularization network to project the preliminary prediction into
shape space, and a shape loss was added to make the prediction more fit into
the desired shape space. Chen et al. [5] proposed an end-to-end approach with
CNN and Conditional Random Fields (CRFs) for segmentation. However, the
inference of CRFs was approximate and not globally optimal.

Inspired by the graph search method [10], Shah et al. [17] first modeled the
terrain-like surfaces segmentation as direct surface identification using
regression. The network consists of an FCN followed by fully-connected layers.
The network is very light weighted and no post processing is required.
Surprisingly the results are very promising [17]. It is well known that U-net
outperforms FCN because U-net has an additional expansive path such that
features of high resolution can be learned and then better prediction accuracy
can be achieved. To improve segmentation accuracy further, it is natural to
think to replace the FCN with a U-net. However, it is not reasonable to
concatenate a U-net with fully-connected layers, as the invariance of feature
maps in the original resolution is supposed to be much less than that in the
low resolution, such that there would be much more chance that the fully-
connected layers heavily overfit to the training data. Another drawback of
Shah et al’s work [17] is that the surface smoothness is implicitly learned
within the whole network (mainly within fully-connected layers) as a black
box. It is hard to decode after training.

To resolve problems mentioned above, we propose to explicitly model the
surface segmentation problem as a quadratic programming, which can be solved
with guaranteed global optimality. The whole network can be trained end-to-
end. Our contributions are in three folds: 1) The first time to parameterize
the output of the U-net as Gaussians (the mean represents surface position
prediction from the U-net and the standard deviation encodes the prediction
confidence), which converts the description from the discrete to the
continuous space; 2) The parameter that controls smoothness can be learned and
has a clear meaning; 3) The method works in the continuous space and enables
sub-pixel segmentation.

##  2 Method

Figure 1: Surface segmentation definition.

We first define the surface segmentation problem. A $3$-D image can be viewed
as a $3$-D tensor $\mathcal{I}$. A terrain-like surface in $\mathcal{I}$ is
oriented and shown in Fig. 1. Let $N_{1}$, $N_{2}$ and $N_{3}$ denote the
image sizes in three dimensions, respectively. Let all column index be a set
$\Omega=\\{(1,1),(1,2),\cdots,(N_{1},N_{3})\\}$. The surface $\mathbf{x}$ is
defined by $x_{i}\in[1,N_{2}]$, $\forall i\in\Omega$. Thus any feasible
surface in $\mathcal{I}$ intersects with each column exactly once. Generally
in surface segmentation [10], the problem is formulated as minimizing the
energy function $E(\mathbf{x}):\mathbb{R}^{N_{1}\times
N_{3}}\rightarrow\mathbb{R}$

$\displaystyle E(\mathbf{x})=E_{u}(\mathbf{x})+E_{p}(\mathbf{x}),$ (1)

where the unary term $E_{u}$ is the energy when considering each column
independently, and the pairwise energy term $E_{p}$ penalizes discontinuity of
surface position among adjacent columns. The design of $E_{u}$ and $E_{p}$
will be detailed in Section 2.1.

###  2.1 Proposed Inference Pipeline

One should note that the proposed method can be applied in both $2$-D and
$3$-D. For the purpose of proof-of-concept and clear explanation, the detailed
description of the proposed method and all experiments are done in $2$-D. The
inference pipeline of the proposed method is demonstrated in Fig. 2(a). The
trained U-net takes in the original image
$\mathcal{I}\in\mathbb{R}^{N_{2}\times N_{1}}$, and outputs the discrete
probability map $P\in\mathbb{R}^{N_{2}\times N_{1}}$. Ideally, for each image
column, the probability for the target surface position is high, and it is
gradually reduced on the positions away from the one on the surface, as
demonstrated in Fig. 2(b). We thus propose a block to convert the discrete
probability map $P$ to a Gaussian parameterization
$G\equiv(\mathbf{\gamma},\mathbf{\sigma})\in\mathbb{R}^{N_{1}\times 2}$, where
$\gamma_{i}$ specifies the mean surface position on each column $i$ and
$\sigma_{i}$ is the corresponding standard deviation. The Gaussian
parameterization $G$ is then fed into the trained smoothing block (SB) , which
incorporate the learned surface smoothness priors to infer the optimal target
surface. Next, we detail the novel blocks in our deep optimal surface
segmentation neural network.

(a)

(b)

Figure 2: (a) Inference pipeline of the proposed method. (b) Target surface
position (left) and its relaxed surface possibility map (right) for one
column.

D2C block The D2C block is basically designed to convert the discrete
probability map of each column $P_{i}\in\mathbb{R}^{N_{2}}$ $\forall
i\in\\{1,2,\cdots,N_{1}\\}$, which is output from U-net, to a continuous
representation $G_{i}\equiv(\gamma_{i},\sigma_{i})\in\mathbb{R}^{2}$ (Gaussian
is utilized in our design). This enables optimizing directly on the surface
position and sub-pixel accuracy prediction. The proposed conversion is
realized by fitting a continuous Gaussian function to the discrete probability
map $P_{i}$, which can be thought as discrete samples of a continuous Gaussian
probability density function. Recall that one dimensional Gaussian function
has the formula $f(j)=A\exp(\frac{-(j-\gamma)^{2}}{2\sigma^{2}}),$ and then
$\ln(f(j))=\ln(A)+\frac{-(j-\gamma)^{2}}{2\sigma^{2}}=\ln(A)-\frac{\gamma^{2}}{2\sigma^{2}}+\frac{2\gamma
j}{2\sigma^{2}}-\frac{j^{2}}{2\sigma^{2}}=a+bj+cj^{2},$ where
$a=\ln(A)-\frac{\gamma^{2}}{2\sigma^{2}}$, $b=\frac{\gamma}{\sigma^{2}}$, and
$c=-\frac{1}{2\sigma^{2}}$. In our setting, for each column, we have $N_{2}$
samples of $(j,f(j))$. We can define an error function [8], namely
$\varepsilon=\sum_{j=1}^{N_{2}}f(j)^{2}(\ln(f(j))-(a+bj+cj^{2}))^{2}.$ Then
minimizing the weighted mean square error (MSE) $\varepsilon$, one can get the
estimates of $a,b,c$ by solving a set of three linear equations and then
$A,\gamma,\sigma$. The problem is very similar to least square problem [8]. In
our implementation, a linear transform is utilized to normalize the
probability map for each column to the range [0,1], then we can ignore the
magnitude $A$. As the computation of each column is independent, it is
straightforward to be extended to $3$-D.

Smoothing Block (SB) To integrate the surface segmentation model (Eq. 1) with
smoothness priors, we define the energy function
$E(\mathbf{x}):\mathbb{R}^{N_{1}}\rightarrow\mathbb{R}$ as

$\displaystyle E(\mathbf{x})$
$\displaystyle=\sum_{i\in\Omega}\psi_{u}(x_{i})+\sum_{\begin{subarray}{c}i,j\in\Omega,(i,j)\in\mathcal{N}\end{subarray}}\psi_{p}(x_{i},x_{j}),$
(2)

where $\psi_{u}(x_{i})$ and $\psi_{p}(x_{i},x_{j})$ are defined as
$\psi_{u}(x_{i})=\frac{(x_{i}-\gamma_{i})^{2}}{2\sigma_{i}^{2}},~{}~{}\psi_{p}(x_{i},x_{j})=w_{comp}\cdot(x_{i}-x_{j})^{2},$
and $\mathcal{N}$ is the set of neighbor columns.

Figure 3: The proposed SB architecture. Only one parameter $w_{comp}$ needs to
be trained.

For simplicity, the nearest neighbor pairs, i.e.
$\mathcal{N}=\\{(1,2),(2,3),\cdots,(N_{1}-1,N_{1})\\}$, are considered as the
set of neighbor columns. The whole energy in Eq. 2 can be reformulated as the
standard quadratic form
$E(\mathbf{x})=\frac{1}{2}\mathbf{x}^{T}H\mathbf{x}+\mathbf{c}^{T}\mathbf{x}+\text{const}.$
It can be proved that the Hessian matrix $H$ is positive definite by using
Gershgorin circle theorem and then the energy function is convex. The gradient
is $\nabla=H\mathbf{x}+\mathbf{c}.$ Let the gradient to be zero, we have the
global optimal solution $\mathbf{x}^{*}=-H^{-1}\mathbf{c}.$ Another advantage
of the proposed energy formulation is the optimal solution can be calculated
in one step, i.e. we do not need to make use of a recurrent neural network
(RNN) to implement SB. It is also straightforward to implement the smoothing
block in 3D.

###  2.2 Training Strategy

####  2.2.1 U-net Pre-training

A common U-net architectures is utilized to generate the discrete probability
map for the input image ($N_{2}\times N_{1}$). The detailed architecture is
demonstrated in Fig. 10. In the proposed method, the softmax layer works on
each column, not on each pixel. The rational is that we assume the target
surface intersects with each column by exactly once, and so the probabilities
are normalized within each column. Also we assume the U-net should output a
Gaussian shaped probability map for each column, which mimics the Bayesian
learning for each column and shares merits with knowledge distillation [9] and
distillation defense [13]. To encourage the U-net outputs reasonable
probability maps, the Kullback–Leibler divergence (KLD) loss is utilized for
the U-net pre-training. KLD is a measure of how one probability distribution
is different from a reference probability distribution. It is equivalent to
the Cross Entropy when the reference is a Kronecker delta function. We propose
to relax the delta function to a Gaussian distribution such that the proposed
D2C block can work properly. We pick $\sigma$ to be around 0.1 times of the
column length and our method is insensitive to $\sigma$. Then we will have one
Gaussian distribution ground truth for each column, denoted as
$T^{\sigma}_{i}\in\mathbb{R}^{N_{2}}$. One illustration is demonstrated in
Fig. 2(b). We denote the output from the U-net as $P\in\mathbb{R}^{N_{2}\times
N_{1}}$ and the loss for the pre-training is formulated as
$loss_{pre}(P,T^{\sigma})=-\sum_{i\in\Omega}D_{KL}(T^{\sigma}_{i}||P_{i})).$

####  2.2.2 Fine Tuning

During the fine tunning phase, the mean square surface position error (MSE) is
utilized as the loss function, which is formulated as
$loss_{fine}(\mathbf{x},\mathbf{t})=\sum_{i\in\Omega}(x_{i}-t_{i})^{2},$ where
$\mathbf{t}\in\mathbb{R}^{N_{1}}$ denotes the ground truth surface positions.
The fine tuning of the whole network proceeds in an alternation fashion (Fig.
4).

Figure 4: Two phases of fine tuning of the proposed network. The fine tuning
alternates between training U-net for $EP_{\text{U-net}}$ epochs and training
SB for $EP_{\text{SB}}$ epochs. The shaded blocks are kept fixed during the
respective training phase.

The validation data is utilized to train the SB, and the training data is for
U-net training. As SB only has one parameter-$w_{comp}$ to train, the
overfitting chance of it is very low. Also the U-net is not trained on
validation data, the learned $w_{comp}$ should be more representative in the
wild. Otherwise if fine tuning U-net+SB simultaneously on the training data,
the learned $w_{comp}$ is generally smaller than necessary, as the pre-trained
U-net generally has fit the training data well and SB does not play an
important role.

##  3 Experiments

###  3.1 SD-OCT Retinal Layer Segmentation

Figure 5: Sample SD-OCT data. Ground truth (left) and predictions (right).
Red: IRPE; Blue: OBM.

The proposed method was applied to retinal layer segmentation in SD-OCT
images, which were obtained from the public dataset [6]. Since the manual
tracings were only available for a region centered at the fovea, subvolumes of
size $400\times 60\times 512$ were extracted around the fovea. The dataset was
randomly divided into 3 sets: 1) Training set - 263 volumes (79 normal, 184
with age-related macular degeneration (AMD)); 2)Validation set - 59 volumes
(18 normal, 41 AMD); 3) Testing set - 58 volumes (18 normal, 40 AMD). The
surfaces considered are S2-Inner Aspect of Retinal Pigment Epithelium Drusen
Complex (IRPE) and S3-Outer Aspect of Bruch Membrane (OBM) as shown in Fig. 5,
which are very challenging to segment.

Pre-processing and Augmentation The intensity of each slice was normalized to
the range [-1, 1]. No additional pre-processing methods were utilized. For the
purpose of pre-training the U-Net, the standard deviation of the Gaussian
model of the surface position on each column was set $\sigma=50$. We augmented
the training data by applying mirroring along the horizontal dimension, and 2
random translations along the axial dimension on the original image and the
mirrored image.

Hyperparameters All training utilized Adam optimizer. For U-net pre-training,
the learning rate was $10^{-4}$, and the training ran for 2000 epochs. For
fine-tuning, the learning rate of U-net ($LR_{\text{U-net}}$) was $10^{-5}$,
the learning rate of SB ($LR_{\text{SB}})$ was $10^{-2}$, the total epochs for
the two phases were: $EP_{\text{U-net}}=10$, and $EP_{\text{SB}}=10$. The
initial smoothness parameter $w_{comp}^{init}$ was set to $10^{-5}$.

Surface | W/O, normal | W/, normal | W/O, AMD | W/, AMD |  [17] normal |  [17] AMD
---|---|---|---|---|---|---
IRPE | 2.31$\pm$0.54 | 2.21$\pm$0.58 | 3.80$\pm$2.01 | 3.63$\pm$1.71 | 3.84$\pm$0.58 | 6.07$\pm$1.84
OBM | 2.97$\pm$0.35 | 2.93$\pm$0.37 | 6.09$\pm$3.10 | 5.83$\pm$2.89 | 4.97$\pm$1.01 | 5.85$\pm$1.80
Table 1: Unsigned mean surface positioning errors (UMSP) for the results on
the SD-OCT test dataset. The unit is in $\mu m$. W/O: without SB; W/: with SB.

Results Unsigned mean surface positioning error (UMSP) [7] is utilized for
evaluation of segmentation accuracy. The quantitative results are summarized
in Table 1. We compare to another deep learning based method [17], which is
the state-of-the-art on this dataset. As for the IRPE surface, the proposed
method without SB outperforms [17]. With SB plugged in, the performance can be
improved further: for the normal cases, the UMSP can be improved by 39%; for
the AMD cases, we can achieve a similar improvement (40%). The segmentation of
the OBM surface is more challenging than that of the IRPE surface. One can
notice that the learned SB can improve the segmentation accuracy consistently.
As for the OBM surface, compared to [17], the proposed method achieves 41%
improvement on the normal cases and performs comparably on the AMD cases. An
interesting observation is that the learned $w_{comp}$ are around 0.05 and 1.0
for IRPE and OBM, respectively, which coincides with the human observation
that generally less context information exists for OBM. Sample segmentation
results are illustrated in Fig. 5, Fig. 6 and Fig. 7.

Figure 6: Sample segmentations of the OBM surface. The four columns illustrate
original images, ground truth, predictions without SB (red: Gaussian means,
blue: std bars) and that with SB (red: prediction, blue: ground truth). For
illustrations, the standard deviation bars are downsampled. The errors shown
are in the unit of pixel.

Figure 7: Sample segmentations of the IRPE surface. The four columns
illustrate original images, ground truth, predictions without SB (red:
Gaussian means, blue: std bars) and that with SB (red: prediction, blue:
ground truth). For illustrations, the standard deviation bars are downsampled.
The errors shown are in the unit of pixel.

###  3.2 IVUS Vessel Wall Segmentation

The data used for this experiment was obtained from the standardized
evaluation of IVUS image segmentation database [2]. In this experiment, the
dataset B was used. This dataset consists of 435 images with a size of
$384\times 384$, as well as the respective expert manual tracings of lumen and
media surfaces. It comprises two groups - a training set (109 slices) and a
testing set (326 slices).

Figure 8: Sample IVUS data. Ground truth (left) and predictions (right) in
Cartesian (top) and Polar system (bottom). Red: lumen; Blue: media.

The experiment with the proposed method was conducted in conformance with the
directives provided for the IVUS challenge. In our experiment, we randomly
split the 109 training slices into 100 slices for training and 9 slices for
validation. Pre-processing and Augmentation Each slice was transformed to be
represented in the polar coordinate system with a size of $256\times 128$, as
illustrated in Fig. 8. The intensity of each slice was normalized to have a
zero mean and a unit standard deviation. The Gaussian truth was generated for
each column using $\sigma=15$, as it has a shorter column (128 vs 512 in the
SD-OCT data). As the number of training data was limited, we augmented the
data on the fly by random combinations of various operations including
mirroring, circulation shifting along the polar dimension, adding Gaussian
noises (mean=0, std=0.1), adding Salt and Pepper noises (5%), and cropping
(90% of the original size) and then resizing ($256\times 128$). All training
parameters were the same as those for the SD-OCT segmentations, except that
the number of epochs for the pre-training was 4000.

Methods | Lumen | Media
---|---|---
| JM | PAD | HD | JM | PAD | HD
P3[3, 4] | 0.88$\pm$0.05 | 0.06$\pm$0.05 | 0.34$\pm$0.14 | 0.91$\pm$0.04 | 0.05$\pm$0.04 | 0.31$\pm$0.12
VGG-U-net[1] | 0.80(-) | - | - | 0.81(-) | - | -
[16] | 0.86$\pm$0.04 | 0.09$\pm$0.03 | 0.37$\pm$0.14 | 0.90$\pm$0.03 | 0.07$\pm$0.03 | 0.43$\pm$0.12
W/O SB | 0.87$\pm$0.11 | 0.08$\pm$0.12 | 0.43$\pm$1.90 | 0.88$\pm$0.10 | 0.08$\pm$0.10 | 0.41$\pm$0.31
W/ SB | 0.88$\pm$0.06 | 0.07$\pm$0.07 | 0.28$\pm$0.18 | 0.89$\pm$0.07 | 0.07$\pm$0.07 | 0.37$\pm$0.27
Table 2: Segmentation results comparison on the IVUS dataset.

Results Jaccard Measure (JM), Percentage of Area Difference (PAD) and
Hausdroff Distance (HD) are utilized to evaluate segmentation accuracy. The
results are summarized in Table 2. P3 [3, 4] is the state-of-the-art method
for this IVUS dataset. It is an expectation maximization based method and is
semi automated. VGG-U-net denotes a deep learning based method [1]. The state-
of-the-art fully automation method is a graph search based method working in
irregularly sampled space [16], the unary energy of which is learned with
random forests. From Table 2, one can find that the proposed method
outperforms the sub-voxel graph based method [16] for the lumen surface and
the performance on the media surface is comparable. Compared to the state-of-
the-art semi-automate method P3, the performance of the proposed is comparable
for the lumen surface and is marginally inferior on the media surface. The
proposed method outperforms the VGG-U-net by a big margin. Sample segmentation
results are illustrated in Fig. 8 and Fig. 9.

Figure 9: Sample segmentation results on IVUS data. Red: lumen ground truth;
Blue: media ground truth; Green: lumen prediction; Purple: media prediction.
The four columns correspond to the original image, the ground truth, the
prediction without SB, the prediction with SB.

As to the inference computation time, the proposed method needs more overhead
than the VGG-U-net method (0.21 vs 0.09 sec/slice). The overhead is mainly
from the Hessian matrix computation, as well as two separate runs of the
program for two surfaces. While compared to P3 [3, 4] (8.6 sec/slice) and the
graph based method [16] (187.4 sec/slice), the proposed method is highly
efficient.

##  4 Conclusion

In this paper, a novel segmentation model based on a convolutional neural
network (CNN) and a learnable surface smoothing block is proposed to tackle
the surface segmentation problem with end-to-end training. To the best of our
knowledge, this is the first study for surface segmentation which can achieve
guaranteed globally optimal solutions using deep learning. Experiments on SD-
OCT retinal layer segmentation and IVUS vessel wall segmentation demonstrated
very promising results. The proposed method is applicable to $2$-D and $3$-D.

Figure 10: U-net architectures. For convolution layers, the number indicates
dimensions of feature maps, and $3\times 3$ convolution kernels are utilized
except the last convolution layer, where $1\times 1$ kernels are chosen; for
max pooling and transposed convolution layers, $2\times 2$ kernels and strides
2 are utilized. A Relu layer is utilized following each convolution layer
except the final $1\times 1$ convolution layer. For all skipping connections,
the features are merged by adding.

## References

  * [1] Balakrishna, C., Dadashzadeh, S., Soltaninejad, S.: Automatic detection of lumen and media in the ivus images using u-net with vgg16 encoder. arXiv preprint arXiv:1806.07554 (2018)
  * [2] Balocco, S., Gatta, C., Ciompi, F., Wahle, A., Radeva, P., Carlier, S., Unal, G., Sanidas, E., Mauri, J., Carillo, X., et al.: Standardized evaluation methodology and reference database for evaluating ivus image segmentation. Comput Med Imaging Graph 38(2), 70–90 (2014)
  * [3] Cardinal, M.H., Meunier, J., Soulez, G., Maurice, R.L., Therasse, É., Cloutier, G.: Intravascular ultrasound image segmentation: a three-dimensional fast-marching method based on gray level distributions. IEEE Trans. Med. Imag 25(5), 590–601 (2006)
  * [4] Cardinal, M.H.R., Soulez, G., Tardif, J.C., Meunier, J., Cloutier, G.: Fast-marching segmentation of three-dimensional intravascular ultrasound images: A pre-and post-intervention study. Med. Phys. 37(7Part1), 3633–3647 (2010)
  * [5] Chen, S., de Bruijne, M.: An end-to-end approach to semantic segmentation with 3d cnn and posterior-crf in medical images. In: Medical Imaging meets NIPS Workshop at NIPS 2018 (2019)
  * [6] Farsiu, S., Chiu, S.J., O’Connell, R.V., Folgar, F.A., Yuan, E., Izatt, J.A., Toth, C.A., Group, A.R.E.D.S..A.S.D.O.C.T.S., et al.: Quantitative classification of eyes with and without intermediate age-related macular degeneration using optical coherence tomography. Ophthalmology 121(1), 162–172 (2014)
  * [7] Garvin, M.K., Abramoff, M.D., Wu, X., Russell, S.R., Burns, T.L., Sonka, M.: Automated 3-d intraretinal layer segmentation of macular spectral-domain optical coherence tomography images. IEEE Trans. Med. Imag 28(9), 1436–1447 (2009)
  * [8] Guo, H.: A simple algorithm for fitting a gaussian function [dsp tips and tricks]. IEEE Signal Process Mag 28(5), 134–137 (2011)
  * [9] Hinton, G., Vinyals, O., Dean, J.: Distilling the knowledge in a neural network. In: Deep Learning and Representation Learning Workshop at NIPS 2014 (2015)
  * [10] Li, K., Wu, X., Chen, D.Z., Sonka, M.: Optimal surface segmentation in volumetric images-a graph-theoretic approach. IEEE Trans Pattern Anal Mach Intell 28(1), 119–134 (2006)
  * [11] Long, J., Shelhamer, E., Darrell, T.: Fully convolutional networks for semantic segmentation. In: CVPR 2015. pp. 3431–3440 (2015)
  * [12] Milletari, F., Rothberg, A., Jia, J., Sofka, M.: Integrating statistical prior knowledge into convolutional neural networks. In: MICCAI 2017. pp. 161–168 (2017)
  * [13] Papernot, N., McDaniel, P., Wu, X., Jha, S., Swami, A.: Distillation as a defense to adversarial perturbations against deep neural networks. In: 2016 IEEE Symposium on Security and Privacy (SP). pp. 582–597 (2016)
  * [14] Ravishankar, H., Venkataramani, R., Thiruvenkadam, S., Sudhakar, P., Vaidya, V.: Learning and incorporating shape models for semantic segmentation. In: MICCAI 2017. pp. 203–211 (2017)
  * [15] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: MICCAI 2015. pp. 234–241 (2015)
  * [16] Shah, A., Abrámoff, M.D., Wu, X.: Optimal surface segmentation with convex priors in irregularly sampled space. Med. Image Anal. (2019)
  * [17] Shah, A., Zhou, L., Abrámoff, M.D., Wu, X.: Multiple surface segmentation using convolution neural nets: application to retinal layer segmentation in oct images. Biomed. Opt. Express 9(9), 4509–4526 (2018)


# What are the Composition and Power of the Jet in Cyg X-1?

Andrzej A. Zdziarski Nicolaus Copernicus Astronomical Center, Polish Academy
of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland<EMAIL_ADDRESS>Elise
Egron INAF-Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047
Selargius, Italy

###### Abstract

We calculate the electron-positron pair production rate at the base of the jet
of Cyg X-1 by collisions of photons from its hot accretion flow using the
measurement of its average soft gamma-ray spectra by the Compton Gamma Ray
Observatory and INTEGRAL satellites. We have found that this rate
approximately equals the flow rate of the leptons emitting the observed
synchrotron radio-to-IR spectrum of the jet, calculated using an extended jet
model following that of Blandford & Königl. This coincidence shows the jet
composition is likely to be pair-dominated. The same coincidences were found
before in the microquasar MAXI J1820+070 and in the radio galaxy 3C 120, which
shows that the considered mechanism can be universal for at least some classes
of relativistic jets. Furthermore, we recalculate the jet power of Cyg X-1.
The presence of pairs can strongly reduces the power in the bulk motion of
ions, which then limits the parameter space at which the jet can power the
$\sim$5-pc nebular structure present in its vicinity.

##  1 Introduction

We study here the issue of the abundance of electron-positron (e±) pairs in
the jet of the microquasar Cyg X-1. While there are strong hints that
extragalactic jets contain substantial numbers of e± pairs (e.g., Ghisellini
2012; Pjanka et al. 2017; Snios et al. 2018; Sikora et al. 2020; Liodakis et
al. 2022), this is much less clear for jets in microquasars. In either case,
the mechanisms producing the putative pairs remain uncertain.

A viable mechanism is e± pair production within the jet base by collisions of
photons emitted by the accretion flow (e.g., Henri & Pelletier 1991;
Beloborodov 1999; Levinson & Rieger 2011; Aharonian et al. 2017; Sikora et al.
2020). Figures 1 and 3 of Zdziarski et al. (2022a) and Zdziarski et al.
(2022b), respectively, give illustrations of the assumed geometry, and so we
do not show it again here. The rate of the pairs produced at the base can be
compared with the flow rate of non-thermal relativistic electrons emitting
synchrotron emission far downstream in the jet. Such comparisons were done for
the microquasar MAXI J1820+070 (Zdziarski et al., 2022b) and for the radio
galaxy 3C 120 (Zdziarski et al., 2022a), in which cases these two rates were
found to be similar, and compatible with being equal.

Here, we perform such a comparison for Cyg X-1, which is probably the most
studied microquasar, discovered in X-rays already in 1964 (Bowyer et al.,
1965). Its radio-to-mm spectrum in the hard state is flat (Fender et al.,
2000), with $\alpha\approx 0$ (defined by the energy flux of
$F_{\nu}\propto\nu^{\alpha}$), see Figure 1(a). This is usually explained by
the jet emission being partially synchrotron self-absorbed with both the
distribution of non-thermal electrons and the magnetic energy flux maintained
along the jet (Blandford & Königl, 1979). This allows us to determine the rate
of the total electron and positron flow through the jet based only on
$F_{\nu}$ (independent of $\nu$ for $\alpha=0$) and the main jet parameters,
see Appendix A, without the need to specify either the location of the
emission along the jet or the break frequency (above which the entire jet
emission is optically thin).

On the other hand, the rate of pair production within the jet base can be
estimated based on the hard X-ray/soft $\gamma$-ray spectrum of the system. It
has been well measured for Cyg X-1, see Figure 1(b). The spectrum above
$\sim$200 keV is relatively well approximated by a power law, which allows us
to calculate this rate based on a fit by Svensson (1987), see Appendix A.

Figure 1: The spectra of Cyg X-1 in the hard state. (a) The radio to IR
region. The black error bars give the fluxes from Pandey et al. (2007) and
Fender et al. (2000), and the red error bars are the IR fluxes (Persi et al.,
1980; Mirabel et al., 1996), which are strongly dominated by the emission of
the donor supergiant. The black dotted line gives a power-law fit to the
2.25–221 GHz measurements from Fender et al. (2000), with $\alpha=0$ and
$F_{\nu}=13$ mJy. The 235 and 610 MHz fluxes lie below that fit, which appears
to be due to electron energy losses in the jet termination region (Zdziarski
et al., 2014). The cyan symbols in the IR region show the broken-power law fit
of Rahoui et al. (2011), claiming the possible position of the break
frequency. (b) The hard X-ray to soft $\gamma$-ray average spectra measured by
different instruments. The spectra from the OSSE, BATSE and COMPTEL detectors
onboard Compton Gamma Ray Observatory are from McConnell et al. (2002), the
spectra from the INTEGRAL ISGRI and PICsIT detectors are from Zdziarski et al.
(2012), and that from the SPI detector is from Jourdain et al. (2012b). The
black solid line give a power law approximation to its high-energy tail with
the spectral index of $\alpha_{\rm X}\approx-2.53$ and $EF_{E}(511\,{\rm
keV})\approx 1.57\times 10^{3}$ eV cm-2 s-1. The dashed lines show an
approximate uncertainty of the power-law normalization (with
$\Delta\log_{10}F_{E}=\pm 0.1$).

##  2 Comparison of the flow rates

The most recent determination of the binary parameters of Cyg X-1 is by
Miller-Jones et al. (2021), who obtained the black-hole (BH) mass of $M_{\rm
BH}\approx 21\pm 2{\rm M}_{\sun}$, the mass of its supergiant donor of
$M_{1}\approx 41\pm 7{\rm M}_{\sun}$ (in agreement with Ziółkowski 2014) and
the distance of $D\approx 2.2\pm 0.2$ kpc. The binary inclination was given by
them with very small uncertainties as $i_{\rm b}\approx
27.5\arcdeg_{-0.6}^{+0.8}$. This, however, is in strong conflict with the
inclination inferred from the X-ray polarization of Cyg X-1 measured by the
Imaging X-ray Polarimetry Explorer (Krawczynski et al., 2022), which yields
much larger values of $\gtrsim\\!45\arcdeg$ (see also Ursini et al. 2022). It
is likely that the viewing angle measured by the X-ray polarimetry corresponds
to an inner hot accretion flow being inclined with respect to the binary axis
but aligned with the BH spin axis. Furthermore, the polarization angle is
aligned with the position angle of the jet, which suggest the jet is
perpendicular to the inner hot flow, in which case its inclination is also
$i\gtrsim\\!45\arcdeg$. Thus, we consider in our analysis two different values
for $i$.

The jet bulk Lorentz factor, $\Gamma$, and the opening angle, $\Theta$, were
measured by Tetarenko et al. (2019) as $\Gamma\approx 2.6^{+0.8}_{-0.6}$ and
$\Theta\approx 0.4\arcdeg$–$1.8\arcdeg$. These values were obtained assuming
$i\approx 27\arcdeg$, while a lower $\Gamma$ and a larger $\Theta$ would be
obtained for a larger value of $i$. Still, we adopt $\Gamma=2.6$ and
$\Theta=1\arcdeg$ as our default parameters. Based on the spectrum of Fender
et al. (2000), shown in Figure 1(a), we use $F_{\nu}=13$ mJy. The steady-state
power-law index of accelerated electrons, $p$, in the jet of Cyg X-1 is
unknown, but the dependence of the electron flow rate, $\dot{N}_{\rm e}$, on
$p$ is relatively weak, see Equation (A1). We take $p=3$, approximately
consistent with the models of Malyshev et al. (2013) and Zdziarski et al.
(2014). As the default parameters, we take the energy density of the
relativistic electrons equal to that of the magnetic field, $\beta_{\rm eq}=1$
(which approximately minimizes the jet power, see Appendix A), the kinetic
energy of ions as equal to that of the electrons, $k_{\rm i}=1$, and the
number of quasi-thermal electrons below the minimum power-law Lorentz factor,
$\gamma_{\rm min}$, equal to that in the power-law electrons, $k_{\rm e}=1$.
The dependence on the maximum Lorentz factor, $\gamma_{\rm max}$, is
negligible for $\gamma_{\rm max}\gg\gamma_{\rm min}$ and $p>2$; we assume
$\gamma_{\rm max}=10^{4}$. At the above values and at $D=2.2$ kpc,
$i=27\arcdeg$, we find the rate of the flow of leptons through the jet from
Equation (A1),

$\dot{N}_{\rm e}\approx 6.4\times 10^{39}\beta_{\rm
eq}^{0.47}\\!\left(\frac{\Gamma}{2.6}\right)^{2.7}\\!\\!\left(\frac{\Theta}{1\arcdeg}\right)^{0.74}\\!\\!\left(\frac{\gamma_{\rm
min}}{3}\right)^{-1.5}\\!\\!{\rm s}^{-1}.$ (1)

This rate is also positively correlated with $i$, $\dot{N}_{\rm
e}\propto(1-\beta\cos i)^{\frac{14+6p}{13+2p}}(\sin i)^{\frac{2-2p}{13+2p}}$,
where the first term dominates the dependence. For $i=45\arcdeg$, the
numerical coefficient above increases to $1.8\times 10^{40}$.

We then compare those values with the pair production rate. The average hard
X-ray/soft $\gamma$-ray spectra of the hard state of Cyg X-1 measured by
various detectors are plotted in Figure 1(b). We have approximated the part of
the spectra at $\gtrsim$200 keV as a power law with the energy index of
$\alpha_{\rm X}\approx-2.53$ (corresponding to the photon index of 3.53) and
the normalization at 511 keV of $EF_{E}\approx 1.57$ keV cm-2 s-1 with the
uncertainty of $\Delta\log_{10}EF_{E}=\pm 0.1$. We then use the same formula
as in Zdziarski et al. (2022a), given here by Equation (A6). This yields the
rate of production of the sum of e+ and e- of

$2\dot{N}_{+}\approx 3^{+2}_{-1}\times 10^{40}\left(\frac{R_{\rm
hot}}{10R_{\rm g}}\right)^{-3}\left(\frac{R_{\rm jet}}{10R_{\rm
g}}\right)^{2}{\rm s}^{-1},$ (2)

where $R_{\rm g}=GM_{\rm BH}/c^{2}\approx 3\times 10^{6}$ cm. The estimate of
the radius of the base of the jet of $R_{\rm jet}\approx 10R_{\rm g}$ follows
from GRMHD simulations, see Fig. 3 in Zdziarski et al. (2022b), which used a
simulation by Tchekhovskoy (2015) performed for the dimensionless spin
parameter of $a_{*}=0.99$. On the other hand, the radius of the part of the
hot disk with most of the high-energy emission is relatively uncertain. The
used estimate corresponds to the case of a rapidly spinning BH, with
$a_{*}\approx 1$, where we used twice the half-power radius, which can be
approximated as $\approx[5+28(1-a_{*})]R_{\rm g}$ (Fabian et al., 2014). If
the spin is lower, both $R_{\rm hot}$ and $R_{\rm jet}$ will be larger. We see
that, for likely jet parameters, it is larger than the obtained values for the
lepton flow through the jet, see Equation (1). We have also calculated the
steady state optical depth of the pairs produced within the jet base, Equation
(A8), and found it to be $\tau_{\rm T}\lesssim 0.14$, where the upper limits
corresponds to static pairs. Then, even a very modest advection velocity
upstream the jet, with $\beta_{\rm adv}\gtrsim 0.03c$, would result in the
annihilation of the produced pairs being negligible. Thus, pair production by
the accretion photons within the volume of the jet base appears to be entirely
capable of providing enough leptons for the radio-mm synchrotron emission of
the jet far away from the BH.

##  3 The jet power

The jet power in the kinetic energy of particles and the magnetic field at the
parameters as used in Equation (1) is, see Equation (A4),

$P_{B{\rm e}}\approx 2.3\times
10^{35}\\!\left(\frac{\Gamma}{2.6}\right)^{3.7}\\!\\!\left(\frac{\Theta}{1\arcdeg}\right)^{0.74}\\!\\!\left(\frac{\gamma_{\rm
min}}{3}\right)^{-0.46}\\!\frac{\rm erg}{\rm s}.$ (3)

At $i=45\arcdeg$, the numerical coefficient above increases to $6.4\times
10^{35}$. For $\beta_{\rm eq}\ll 1$ and $\gg 1$, $P_{B{\rm e}}$ is larger by
the approximate factors $\beta_{\rm eq}^{-0.53}$ and $\beta_{\rm eq}^{0.47}$,
respectively, which exponents follow from Equation (A4) at the assumed $p=3$.

The power in the bulk motion of cold ions at the parameters as used in
Equation (1) and for the hydrogen abundance of $X=0.5$ (Miller-Jones et al.,
2021) is, see Equation (A5),

$\displaystyle P_{\rm i}\approx 2.0\times 10^{37}\beta_{\rm
eq}^{0.47}\\!\left(\frac{\Gamma}{2.6}\right)^{2.7}\frac{\Gamma-1}{1.6}\\!\\!\left(\frac{\Theta}{1\arcdeg}\right)^{0.74}\times$
(4) $\displaystyle\left(\frac{\gamma_{\rm
min}}{3}\right)^{-1.5}\left(1-\frac{2\dot{N}_{+}}{\dot{N}_{\rm
e}}\right)\frac{\rm erg}{\rm s}.$ (5)

At $i=45\arcdeg$, the numerical coefficient above increases to $5.7\times
10^{37}$.

If the $\sim$5 pc-diameter structure discovered close to Cyg X-1 (Gallo et
al., 2005) is indeed powered by one of its jets (which appears not fully
certain, Sell et al. 2015), the required total jet power is
$\approx\\!(1$–$3)\times 10^{37}$ erg s-1 (Russell et al., 2007). This can put
significant constraints on the jet parameter space, see also Malzac et al.
(2009).

We consider separately the cases of $\beta_{\rm eq}\ll 1$, $\approx$1 and
$\gg$1\. For $\beta_{\rm eq}\ll 1$ (a strongly magnetized plasma),
$\dot{N}_{\rm e}$ is reduced by $\approx\beta_{\rm eq}^{0.5}$, which makes it
even easier for pairs produced at the jet base to account for the produced
leptons. $P_{\rm i}$ is then lowered by the same factors, and, if pairs
dominate the flow, is lowered even more. Then, $P_{B{\rm e}}$ is larger by
$\approx\beta_{\rm eq}^{-0.5}$, but, given its normalization, an extremely low
$\beta_{\rm eq}$ would be required for $P_{B{\rm e}}$ to power the nebula.
Furthermore, the 235 and 610 MHz data (Figure 1a) indicate that the
relativistic electrons lose energy in the jet termination region, which would
make the powering of the nebula by $P_{B{\rm e}}$ even more difficult. At
$\beta_{\rm eq}\approx 1$, the jet can also be composed mostly of pairs. Then
$P_{\rm i}\gg P_{B{\rm e}}$ unless there are about 100 pairs per ion. Since
$P_{\rm i}$ decreases with the increasing pair abundance, it can power the
nebula for modest pair abundances only. At $\beta_{\rm eq}\gg 1$, $P_{\rm
i}\gg P_{B{\rm e}}$ in general, and the parameter space at which the base pair
production can account for the synchrotron-emitting leptons is reduced, but it
is still possible for some combinations of the parameters. The nebula can then
be powered by $P_{\rm i}$.

##  4 Discussion

We consider here some effects relevant to our results. Detailed applications
of the model of Blandford & Königl (1979) using the break energy of 0.15 eV
obtained by Rahoui et al. (2011), shown in Figure 1(a), lead to the
determination of the distance of the 15 GHz emission from the BH center,
$z\sim 3\times 10^{6}R_{\rm g}$ (e.g., Zdziarski et al. 2012, 2014), or $z\sim
3A$, where $A\approx 4\times 10^{12}$ cm is the separation between the stellar
components. Such a distance is also in a qualitative agreement with the
presence of a strong orbital modulation of the radio emission (Pooley et al.,
1999), which appears to be due to an orbital-phase dependent free-free
absorption by the stellar wind of the donor. The large depth of the modulation
at 15 GHz of $\approx 30\%$ requires then $z$ to be comparable to $A$ (Szostek
& Zdziarski, 2007; Zdziarski, 2012).

On the other hand, Stirling et al. (2001) and Rushton et al. (2010) found that
only $\sim$1/3–1/2 of the 8.4 GHz emission is within the unresolved beam of
VLBA, with the size of $\approx 10^{14}/\sin i$ cm, and the remainder forms an
extended jet up to the distance $\sim$5 times larger. At 15 GHz, the VLBA beam
is about twice smaller than that at 8.4 GHz, while the resolved fraction is
similar (Rushton, 2009; Zdziarski, 2012). An application of the partially
self-absorbed model of Blandford & Königl (1979) yields then the location of
the bulk of the 8.4 GHz emission at $\approx 10^{14}/\sin i$ cm (Heinz, 2006),
i.e., the distance about 30 times larger than that discussed above. This is in
agreement with the result of Tetarenko et al. (2019), who found that 11 GHz
emission lags111While the 11 GHz/X-rays lag appears relatively certain, we
note that the relationship between the location and the emission frequency of
$z\propto\nu^{-0.4}$ derived in Tetarenko et al. (2019) using lags between
different ratio frequencies cannot be reconciled with the model of Blandford &
Königl (1979) and Königl (1981) given the well-established spectral index of
$\alpha\approx 0$. The two quantities together imply the electron density
increasing fast with the distance, e.g., $a=-11.5$, $b=8.5$ (defined in
Appendix A) for $p=3$, obtained by combining equations (3) and (15) in
Zdziarski et al. (2022b). Thus, either the obtained radio–radio lags are
spurious or the model of partially self-absorbed synchrotron emission does not
apply to Cyg X-1. behind X-rays by $33.5^{+1.9}_{-1.7}$ min, which corresponds
to its source at $z\approx 2.5\times 10^{14}$ cm (assuming $i\approx
30\arcdeg$ and the jet velocity of $\beta\approx 0.9$). However, any orbital
modulation due to wind absorption at this distance would be very tiny. A
resolution of this discrepancy appears to require the presence of two
dissipation regions in the jet, one at distances of the order of the
separation, and one further away, as proposed in Zdziarski (2012). The former
can occur due to the jet-stellar wind interaction (Perucho & Bosch-Ramon,
2012; Yoon & Heinz, 2015; Yoon et al., 2016). The two components have similar
fluxes, as follows from the comparison of the resolved and unresolved VLBA
components. The hypothesis of the two dissipation regions can be verified by
future studies of radio/X-ray time lags. In the case of the study of Tetarenko
et al. (2019), there was a $\approx$2000-s gap in the X-ray coverage starting
just at the onset of the main flare seen in the 11 and 9-GHz light curves.
Then, any time lag present in the data shorter than several hundred s would be
missing in that analysis.

This would affect our results by lowering the lepton flow rate, $\dot{N}_{\rm
e}$, by a factor of $\sim$2\. This is because the electrons produced in the
inner dissipation region could, after losing their energy, be re-accelerated
in the outer dissipation zone. Thus, it would increase the parameter space
within which the pairs produced by photon-photon collisions at the jet base
can account for the synchrotron emission far upstream in the jet. Still, given
the uncertainties of this process we do not take this into account in our
formalism.

Another issue concerns the measured very high polarization degree (reaching
unity at high energies) found in the INTEGRAL SPI data by Jourdain et al.
(2012a). If this is due to the synchrotron emission of the jet, it would be
emitted close to the onset of the electron acceleration, which is inferred to
be at $z\gg R_{\rm g}$ in Cyg X-1 (e.g., Zdziarski et al. 2014). Then, the
emitted soft $\gamma$-rays would not have the density sufficient for efficient
pair production. However, we consider that origin uncertain, and in this work
we assume the soft $\gamma$-rays are emitted within the accretion flow. This
is agreement with the modelling by McConnell et al. (2002), and strongly
suggested by the shape of the spectrum, with the tail smoothly joining the
emission at lower energies.

##  5 Conclusions

The observed hard X-ray/soft $\gamma$-ray emission of Cyg X-1 in its hard
state extends to several MeV, forming a steep high-energy tail originating
close to the peak (in $EF_{E}$) of the X-ray emission at lower energies, see
Figure 1(b). The origin of this X-ray/$\gamma$-ray spectrum is likely to be
Comptonization of some soft seed photons by hybrid (thermal with a high-energy
tail) electrons within the hot accretion flow (McConnell et al., 2002). We
have found that e± pair production by collisions of the photons from the
photon tail can readily supply enough electrons and positrons to account for
the radio-to-mm emission of the jet, under the assumption that the latter
originates from the partially-self-absorbed synchrotron mechanism, which is
the standard explanation for jet spectra with $\alpha\sim 0$ (Blandford &
Königl, 1979). Thus, the composition of the jet in Cyg X-1 can be dominated by
e± pairs. Such conclusions were achieved before for the microquasar MAXI
J1820+070 (Zdziarski et al., 2022b) and the radio galaxy 3C 120 (Zdziarski et
al., 2022a), which shows that the considered mechanism can be quite universal.

The presence of pairs reduces the density of the ions in the jet, and,
consequently, the power in the ion bulk motion (dominating the total power in
the absence of pairs) can be strongly reduced. This makes it less likely that
the jet of Cyg X-1 powers the nebula discovered by Gallo et al. (2005).
However, no reliable value of the jet power vs. the pair abundance can be
obtained due to the present uncertainty regarding the parameters of the jet
(e.g., its equipartition parameter).

## Acknowledgements

We thank the referee for valuable comments. We acknowledge support from the
Polish National Science Center under the grant 2019/35/B/ST9/03944. This
research benefitted from discussion at Team meetings at the International
Space Science Institute (Bern).

##  Appendix A Formulae

We give here the explicit formulae for the electron flow rate and the jet
power based on the formulation of the Blandford & Königl (1979) and Königl
(1981) model presented in Zdziarski et al. (2019), see also Zdziarski et al.
(2022b). We consider power-law dependencies of the electron number density and
the magnetic field strength, $n(z,\gamma)=n_{0}(z/z_{0})^{-a}\gamma^{-p}$
(where $\gamma$ is the electron Lorentz factor), $B(z)=B_{0}(z/z_{0})^{-b}$,
respectively. Here $z_{0}$ is the distance of the onset of the acceleration.
We assume the canonical values of $a=2$, $b=1$, for which $\alpha=0$. In this
case, the rate of the flow of leptons (electrons and positrons) and the jet
power are constant through the jet, and they are thus independent of the value
of $z_{0}$. We give a formulation also suitable for extragalactic sources,
including the dependence on the cosmological redshift, $z_{\rm r}$, and with
the distance to the source, $D$, being the luminosity distance. The flow rate
can be derived as

$\displaystyle\dot{N}_{\rm e}=(1+k_{\rm
e})\frac{3^{\frac{18+2p}{13+2p}}c^{\frac{23+2p}{13+2p}}(m_{\rm
e}/4)^{\frac{5}{13+2p}}f_{N}\beta\Gamma}{e^{2}\delta^{\frac{14+6p}{13+2p}}}\times$
$\displaystyle\left[\frac{\beta_{\rm eq}}{(1+k_{\rm
i})(f_{E}-f_{N})}\right]^{\frac{3+2p}{13+2p}}\left[\frac{\pi C_{2}(p)}{\sin
i}\right]^{\frac{2p-2}{13+2p}}\times$ (A1)
$\displaystyle\left[\frac{5F_{\nu}D^{2}\tan\Theta}{C_{1}(p)(1+z_{\rm
r})\Gamma_{\rm E}\left(\frac{p-1}{p+4}\right)}\right]^{\frac{8+2p}{13+2p}},$

where $k_{\rm e}$ is the relative contribution of quasi-thermal electrons
below $\gamma_{\rm min}$, $e$ is the electron charge, $C_{1,2}$ are functions
of $p$ defined, e.g., in equations (8) and (9) of Zdziarski et al. (2022b),
respectively, $\Gamma_{\rm E}$ is the Euler Gamma function, $\delta$ is the
Doppler factor,

$\displaystyle\delta=\frac{1}{\Gamma(1-\beta\cos i)},\quad\beta_{\rm
eq}={n_{0}m_{\rm e}c^{2}(1+k_{\rm i})(f_{E}-f_{N})\over B_{0}^{2}/8\pi},$ (A2)
$\displaystyle f_{E}\equiv\begin{cases}{\gamma_{\rm min}^{2-p}-\gamma_{\rm
max}^{2-p}\over p-2},&p\neq 2;\cr\ln{\gamma_{\rm max}\over\gamma_{\rm
min}},&p=2,\cr\end{cases}\quad f_{N}\equiv\frac{\gamma_{\rm
min}^{1-p}-\gamma_{\rm max}^{1-p}}{p-1},$ (A3)

and $k_{\rm i}$ is the fractional contribution to the total kinetic energy
density in particles other than the power-law electrons, i.e., ions and
electrons below $\gamma_{\rm min}$.

Then the jet power in leptons and magnetic field of both the jet and the
counterjet is also constant along the distance from the BH,

$\displaystyle P_{B{\rm
e}}=\frac{2^{\frac{3+2p}{13+2p}}3^{\frac{5}{13+2p}}c^{\frac{49+6p}{13+2p}}m_{\rm
e}^{\frac{18+2p}{13+2p}}\beta\Gamma^{2}(3+2\beta_{\rm
eq})}{e^{2}\delta^{\frac{14+6p}{13+2p}}}\times$
$\displaystyle\left[\frac{(1+k_{\rm i})(f_{E}-f_{N})}{\beta_{\rm
eq}}\right]^{\frac{10}{13+2p}}\left[\frac{\pi C_{2}(p)}{\sin
i}\right]^{\frac{2p-2}{13+2p}}\times$ (A4)
$\displaystyle\left[\frac{5F_{\nu}D^{2}\tan\Theta}{C_{1}(p)(1+z_{\rm
r})\Gamma_{\rm E}\left(\frac{p-1}{p+4}\right)}\right]^{\frac{8+2p}{13+2p}}.$

This power is minimized for $\beta_{\rm eq}=15/(3+2p)$. The usable power in
the bulk motion of cold ions is

$P_{\rm i}=\mu_{\rm e}m_{\rm p}c^{2}(\Gamma-1)\left(\dot{N}_{\rm
e}-2\dot{N}_{+}\right),$ (A5)

where $\mu_{\rm e}=2/(1+X)$ is the mean electron molecular weight, $X$ is the
hydrogen mass content and $m_{\rm p}$ is the proton mass. The term
$(\Gamma-1)$ takes into account the fact that rest-energy flow through the
jets is provided by the accretion flow. Then, the minimum of the total jet
power depends on the unknown pair abundance, and we thus do not give it here.

For completeness, we also provide the used formula for pair production rate by
the accretion photons within the jet base, see equations (A10–A11) of
Zdziarski et al. (2022a) and Svensson (1987),

$\displaystyle\dot{N}_{+}\approx 2\pi R_{\rm jet}^{2}R_{\rm
hot}n_{1}^{2}\sigma_{\rm T}c\frac{\ln(E_{\rm c}/511\,{\rm
keV})}{(1-\alpha_{\rm X})^{5/3}(2-\alpha_{\rm X})},$ (A6) $\displaystyle
n_{1}\approx F_{E}(511\,{\rm keV})\frac{4\pi D^{2}}{2\pi R_{\rm hot}^{2}c},$
(A7)

where $n_{1}$ is the differential photon density at 511 keV (in the source
frame), above the hot disk, $F_{E}$ and $E$ are assumed to have the same
energy unit (e.g., eV), $R_{\rm hot}$ and $R_{\rm jet}$ are the characteristic
radii of the hot accretion flow and the jet, respectively, $\sigma_{\rm T}$ is
the Thomson cross section, and $E_{\rm c}$ is the upper cutoff of the photon
power law, which we assume $\approx$3 MeV. The upper limit on the optical
depth of the pairs within the jet base can be roughly estimated as (Zdziarski
et al., 2021)

$\tau_{\rm T,max}\approx\frac{14\ln(E_{\rm c}/511\,{\rm keV})\sigma_{\rm
T}F_{E}(511\,{\rm keV})D^{2}}{(1-\alpha_{\rm X})^{5/3}(2-\alpha_{\rm X})},$
(A8)

which is obtained by equating the local pair production and pair annihilation
rates. The rate of pair advection will be faster than that of the annihilation
for $\beta_{\rm adv}>(3/16)\tau_{\rm T}$ (Zdziarski et al., 2021). Thus, if
$\tau_{\rm T}$ is low, even a very low $\beta_{\rm adv}$ will prevent
substantial annihilation within the jet base. Then, most of the pairs will be
advected upstream, where further annihilation is completely negligible.

## References

  * Aharonian et al. (2017) Aharonian, F. A., Barkov, M. V., & Khangulyan, D. 2017, ApJ, 841, 61, doi: 10.3847/1538-4357/aa7049
  * Beloborodov (1999) Beloborodov, A. M. 1999, MNRAS, 305, 181, doi: 10.1046/j.1365-8711.1999.02384.x
  * Blandford & Königl (1979) Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34, doi: 10.1086/157262
  * Bowyer et al. (1965) Bowyer, S., Byram, E. T., Chubb, T. A., & Friedman, H. 1965, Science, 147, 394, doi: 10.1126/science.147.3656.394
  * Fabian et al. (2014) Fabian, A. C., Parker, M. L., Wilkins, D. R., et al. 2014, MNRAS, 439, 2307, doi: 10.1093/mnras/stu045
  * Fender et al. (2000) Fender, R. P., Pooley, G. G., Durouchoux, P., Tilanus, R. P. J., & Brocksopp, C. 2000, MNRAS, 312, 853, doi: 10.1046/j.1365-8711.2000.03219.x
  * Gallo et al. (2005) Gallo, E., Fender, R., Kaiser, C., et al. 2005, Nature, 436, 819, doi: 10.1038/nature03879
  * Ghisellini (2012) Ghisellini, G. 2012, MNRAS, 424, L26, doi: 10.1111/j.1745-3933.2012.01280.x
  * Heinz (2006) Heinz, S. 2006, ApJ, 636, 316, doi: 10.1086/497954
  * Henri & Pelletier (1991) Henri, G., & Pelletier, G. 1991, ApJ, 383, L7, doi: 10.1086/186228
  * Jourdain et al. (2012a) Jourdain, E., Roques, J. P., Chauvin, M., & Clark, D. J. 2012a, ApJ, 761, 27, doi: 10.1088/0004-637X/761/1/27
  * Jourdain et al. (2012b) Jourdain, E., Roques, J. P., & Malzac, J. 2012b, ApJ, 744, 64, doi: 10.1088/0004-637X/744/1/64
  * Königl (1981) Königl, A. 1981, ApJ, 243, 700, doi: 10.1086/158638
  * Krawczynski et al. (2022) Krawczynski, H., Muleri, F., Dovčiak, M., et al. 2022, arXiv e-prints, arXiv:2206.09972.  https://arxiv.org/abs/2206.09972
  * Levinson & Rieger (2011) Levinson, A., & Rieger, F. 2011, ApJ, 730, 123, doi: 10.1088/0004-637X/730/2/123
  * Liodakis et al. (2022) Liodakis, I., Blinov, D., Potter, S. B., & Rieger, F. M. 2022, MNRAS, 509, L21, doi: 10.1093/mnrasl/slab118
  * Malyshev et al. (2013) Malyshev, D., Zdziarski, A. A., & Chernyakova, M. 2013, MNRAS, 434, 2380, doi: 10.1093/mnras/stt1184
  * Malzac et al. (2009) Malzac, J., Belmont, R., & Fabian, A. C. 2009, MNRAS, 400, 1512, doi: 10.1111/j.1365-2966.2009.15553.x
  * McConnell et al. (2002) McConnell, M. L., Zdziarski, A. A., Bennett, K., et al. 2002, ApJ, 572, 984, doi: 10.1086/340436
  * Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046, doi: 10.1126/science.abb3363
  * Mirabel et al. (1996) Mirabel, I. F., Claret, A., Cesarsky, C. J., Boulade, O., & Cesarsky, D. A. 1996, A&A, 315, L113
  * Pandey et al. (2007) Pandey, M., Rao, A. P., Ishwara-Chandra, C. H., Durouchoux, P., & Manchanda, R. K. 2007, A&A, 463, 567, doi: 10.1051/0004-6361:20065856
  * Persi et al. (1980) Persi, P., Ferrari-Toniolo, M., Grasdalen, G. L., & Spada, G. 1980, A&A, 92, 238
  * Perucho & Bosch-Ramon (2012) Perucho, M., & Bosch-Ramon, V. 2012, A&A, 539, A57, doi: 10.1051/0004-6361/201118262
  * Pjanka et al. (2017) Pjanka, P., Zdziarski, A. A., & Sikora, M. 2017, MNRAS, 465, 3506, doi: 10.1093/mnras/stw2960
  * Pooley et al. (1999) Pooley, G. G., Fender, R. P., & Brocksopp, C. 1999, MNRAS, 302, L1, doi: 10.1046/j.1365-8711.1999.02225.x
  * Rahoui et al. (2011) Rahoui, F., Lee, J. C., Heinz, S., et al. 2011, ApJ, 736, 63, doi: 10.1088/0004-637X/736/1/63
  * Rushton et al. (2010) Rushton, A., Miller-Jones, J. C. A., Paragi, Z., et al. 2010, in Proceedings of Science, Vol. 10, 61
  * Rushton (2009) Rushton, A. P. 2009, PhD thesis, University of Manchester, UK
  * Russell et al. (2007) Russell, D. M., Fender, R. P., Gallo, E., & Kaiser, C. R. 2007, MNRAS, 376, 1341, doi: 10.1111/j.1365-2966.2007.11539.x
  * Sell et al. (2015) Sell, P. H., Heinz, S., Richards, E., et al. 2015, MNRAS, 446, 3579, doi: 10.1093/mnras/stu2320
  * Sikora et al. (2020) Sikora, M., Nalewajko, K., & Madejski, G. M. 2020, MNRAS, 499, 3749, doi: 10.1093/mnras/staa3128
  * Snios et al. (2018) Snios, B., Nulsen, P. E. J., Wise, M. W., et al. 2018, ApJ, 855, 71, doi: 10.3847/1538-4357/aaaf1a
  * Stirling et al. (2001) Stirling, A. M., Spencer, R. E., de la Force, C. J., et al. 2001, MNRAS, 327, 1273, doi: 10.1046/j.1365-8711.2001.04821.x
  * Svensson (1987) Svensson, R. 1987, MNRAS, 227, 403
  * Szostek & Zdziarski (2007) Szostek, A., & Zdziarski, A. A. 2007, MNRAS, 375, 793, doi: 10.1111/j.1365-2966.2006.11346.x
  * Tchekhovskoy (2015) Tchekhovskoy, A. 2015, Launching of Active Galactic Nuclei Jets, ASSL, Vol. 414 (Springer), 45, doi: 10.1007/978-3-319-10356-3_3
  * Tetarenko et al. (2019) Tetarenko, A. J., Casella, P., Miller-Jones, J. C. A., et al. 2019, MNRAS, 484, 2987, doi: 10.1093/mnras/stz165
  * Ursini et al. (2022) Ursini, F., Matt, G., Bianchi, S., et al. 2022, MNRAS, 510, 3674, doi: 10.1093/mnras/stab3745
  * Yoon & Heinz (2015) Yoon, D., & Heinz, S. 2015, ApJ, 801, 55, doi: 10.1088/0004-637X/801/1/55
  * Yoon et al. (2016) Yoon, D., Zdziarski, A. A., & Heinz, S. 2016, MNRAS, 456, 3638, doi: 10.1093/mnras/stv2954
  * Zdziarski (2012) Zdziarski, A. A. 2012, MNRAS, 422, 1750, doi: 10.1111/j.1365-2966.2012.20754.x
  * Zdziarski et al. (2012) Zdziarski, A. A., Lubiński, P., & Sikora, M. 2012, MNRAS, 423, 663, doi: 10.1111/j.1365-2966.2012.20903.x
  * Zdziarski et al. (2022a) Zdziarski, A. A., Phuravhathu, D. G., Sikora, M., Böttcher, M., & Chibueze, J. O. 2022a, ApJ, 928, L9, doi: 10.3847/2041-8213/ac5b70
  * Zdziarski et al. (2014) Zdziarski, A. A., Pjanka, P., Sikora, M., & Stawarz, Ł. 2014, MNRAS, 442, 3243, doi: 10.1093/mnras/stu1009
  * Zdziarski et al. (2019) Zdziarski, A. A., Stawarz, Ł., & Sikora, M. 2019, MNRAS, 485, 1210, doi: 10.1093/mnras/stz475
  * Zdziarski et al. (2022b) Zdziarski, A. A., Tetarenko, A. J., & Sikora, M. 2022b, ApJ, 925, 189, doi: 10.3847/1538-4357/ac38a9
  * Zdziarski et al. (2021) Zdziarski, A. A., Jourdain, E., Lubiński, P., et al. 2021, ApJ, 914, L5, doi: 10.3847/2041-8213/ac0147
  * Ziółkowski (2014) Ziółkowski, J. 2014, MNRAS, 440, L61, doi: 10.1093/mnrasl/slu002