# Single Bell inequality to detect genuine nonlocality in three-qubit
genuinely entangled states

Ignacy Stachura  Owidiusz Makuta  Remigiusz Augusiak  Center for Theoretical
Physics, Polish Academy of Sciences, Aleja Lotników 32/46, 02-668 Warsaw,
Poland

###### Abstract

It remains an open question whether every pure multipartite state that is
genuinely entangled is also genuinely nonlocal. Recently, a new general
construction of Bell inequalities allowing the detection of genuine
multipartite nonlocality (GMNL) in quantum states was proposed in [F. J.
Curchod, M. L. Almeida, and A. Acín, New J. Phys. 21, 023016 (2019)] with the
aim of addressing the above problem. Here we show how, in a simple manner, one
can improve this construction to deliver tighter Bell inequalities for
detection of GMNL. Remarkably, we then prove one of the improved Bell
inequalities to be powerful enough to detect GMNL in every three-qubit
genuinely entangled state. We also generalize some of these inequalities to
detect not only GMNL but also nonlocality depth in multipartite states and we
present a possible way of generalizing them to the case of more outcomes.

##  I Introduction

Using entangled quantum particles one may perform an experiment which gives
rise to correlations that can by no means be explained by a local realistic
theory. The existence of such correlations, referred to as Bell nonlocality,
is one of the most characteristic features of quantum mechanics. Moreover,
nonlocal correlations have been harnessed as a resource for numerous device-
independent applications such as quantum cryptography [1], certification of
randomness [2] or self-testing [3, 4].

One of the fundamental problems in quantum information and foundations of
quantum physics is understanding the relationship between the Bell nonlocality
and quantum entanglement, another key feature of quantum theory. While the
former implies the latter in the sense that a state that exhibits Bell non-
locality must be necessarily entangled, the opposite implication does not
always hold true. Even though a considerable effort has been devoted to
understanding which entangled states do or do not violate Bell inequalities
(see, e.g., Refs. [5, 6, 7, 8]), the relationship between entanglement and
nonlocality for general composite quantum systems remains unclear even in case
of pure states. In fact, while this connection has been elucidated for two-
particle systems, the exploration of the interplay between entanglement and
nonlocality in systems involving more than two particles remains an active
area of investigation. Indeed, it was demonstrated in Refs. [9, 10] that every
pure entangled state violates local realism. Subsequently, in Ref. [11],
Popescu and Rohrlich generalized this observation to the multipartite setting
proving that any pure entangled state displays some form of nonlocality.

While one could conclude that these results definitely resolve the problem,
this is not entirely the case because the truly relevant question in the
multipartite setting is whether genuine entanglement [12] implies genuine
nonlocality [13, 14]. There have been attempts to settle this problem. In
particular, in Ref. [15], using a three-partite form of Hardy-type argument,
it was proven that any three-qubit genuinely entangled state is also genuinely
nonlocal, whereas in Ref. [16] the same conclusion was drawn for arbitrary
$n$-qubit symmetric states. Despite this progress, the question of whether any
genuinely entangled state is genuinely nonlocal, beyond the above simple
cases, remains open.

A possible way to tackle this problem goes through the construction of
suitable Bell inequalities whose violation reveals GMNL. While many such Bell
inequalities have been derived to date (see, e.g., Refs. [13, 14, 17, 18, 19,
20]), none of them has been proven to be powerful enough to demonstrate
genuine nonlocality in every genuinely entangled pure multipartite state.

Another technique to construct Bell inequalities witnessing GMNL which seems
promising in this context was proposed recently by Curchod and collaborators
in Ref. [21]. It takes a bipartite Bell inequality, called a seed inequality,
that obeys certain requirements and reforges it into a multipartite one
detecting GMNL. Importantly, it was analytically proven in Ref. [21] that the
inequalities obtained in this way with the seed inequality being the CHSH Bell
inequality [22] are capable of revealing GMNL of any three-qubit pure state
which is symmetric under permutation of two chosen parties as well as of GHZ-
like states of any number of qubits. These analytical results were also
supplemented by numerical evidence that one of these inequalities is violated
by any four-qubit GME state. It is worth mentioning that more recently,
suitable modifications of these Bell inequalities were used in Ref. [23] to
characterize GMNL in quantum networks.

Our aim in this work is to show that a slightly improved version of one of the
inequalities of Ref. [21] is powerful enough to definitely settle the problem
of equivalence of genuine nonlocality and genuine entanglement for three-qubit
states. In other words, we introduce a single Bell inequality that is capable
of detecting genuine nonlocality in any three-qubit genuinely entangled. We
thus not only recover and confirm the result of Ref. [15], but we also
introduce a universal test of genuine nonlocality for three-qubit states
which, unlike the Hardy-type arguments, is capable of tolerating some amount
of experimental imperfections and noises. On the way of achieving the above
result we show how to modify the construction of Ref. [21] to provide tighter
Bell inequalities detecting GMNL which improves that of Ref. [21]. We also
generalize some of the inequalities of Ref. [21] to detect the nonlocality
depth of multipartite systems. Finally, we explore the possibility of
generalizing the approach of Ref. [21] to Bell scenarios involving more
outcomes, concentrating on a particular case of three-outcome measurements.
Indeed, we show how to suitably modify the CGLMP Bell inequality [24] to
derive an inequality detecting GMNL in tripartite systems which is analogous
to inequalities of Ref. [21]

##  II Preliminaries

Before presenting our main results, let us first introduce some relevant
terminology. We will introduce the notions of genuine multipartite
entanglement (GME) and genuine multipartite nonlocality (GMNL) and then in
order to make the paper self-contained we will recall the constructions of
Bell inequalities detecting GMNL from Ref. [21].

###  II.1 Genuinely multipartite entanglement and nonlocality

Let us consider $n$ observers $A_{1},\ldots,A_{n}$ sharing some $n$-partite
pure state $|\psi\rangle$ from the Hilbert space
$\mathcal{H}=\mathbb{C}^{d_{1}}\otimes\ldots\otimes\mathbb{C}^{d_{n}}$. Let us
then consider a partition of all the parties into two disjoint and nontrivial
sets $S$ and $\bar{S}$ such that
$S\cup\bar{S}=\mathcal{A}:=\\{A_{1},\ldots,A_{n}\\}$. We say that
$|\psi\rangle$ is entangled across this bipartition if
$|\psi\rangle\neq|\phi_{S}\rangle\otimes|\phi_{\bar{S}}\rangle$ for any two
vectors $|\phi_{S}\rangle$ and $|\phi_{\bar{S}}\rangle$ corresponding to the
parties in $S$ and $\bar{S}$, respectively. If the state $|\psi\rangle$ is
entangled with respect to all bipartitions, we call it genuinely multipartite
entangled (GME) [12]. An excellent example of a pure GME state is the
$n$-qudit Greenberger-Horne-Zeilinger (GHZ) state [25]

$|\mathrm{GHZ}_{n,d}\rangle=\frac{1}{\sqrt{d}}\sum_{i=0}^{d-1}|i\rangle^{\otimes
n}.$ (1)

The above definition of entanglement can further be extended to mixed states.
Namely, a state $\rho$ is entangled across a bipartition $S|\bar{S}$ if it
does not admit the following representation

$\rho=\sum_{\lambda}q_{\lambda}\,\rho_{S}^{\lambda}\otimes\rho_{\bar{S}}^{\lambda},$
(2)

where $\rho_{S}^{\lambda}$ and $\rho_{\bar{S}}^{\lambda}$ are quantum states
corresponding to the sets of parties $S$ and $\bar{S}$, whereas $q_{\lambda}$
is some probability distribution. A state that admits the above representation
is called separable across $S|\bar{S}$. We then say that $\rho$ is GME [12] if
it cannot be represented as a convex combination of states that are separable
across any possible bipartition,

$\rho=\sum_{S|\bar{S}}r_{S|\bar{S}}\sum_{\lambda}q_{\lambda}^{S|\bar{S}}\,\rho_{S}^{\lambda}\otimes\rho_{\bar{S}}^{\lambda},$
(3)

where $r_{S|\bar{S}}$ is some probability distribution.

Let us now move on to the notions of genuine multipartite nonlocality and
nonlocality depth. We again consider $n$ parties sharing a state $\rho$ and we
assume that each party performs one of $m$, $d$-outcome measurements on their
share of the state. We label the measurements choices and outcomes of party
$A_{i}$ as $x_{i}=1,\ldots,m$ and $a_{i}=0,\ldots,d-1$, respectively. After
many rounds of repeated measurements, the parties can estimate the underlying
joint probabilities $p(a_{1},\ldots,a_{n}|x_{1},\ldots,x_{n})$ of obtaining
results $a_{1},\ldots,a_{n}$ upon performing the measurements labeled by
$x_{1},\ldots,x_{n}$.

The correlations that are produced in this experiment are described by a
collection of probability distributions that can be ordered into a vector

$\vec{p}=\\{p(\vec{a}|\vec{x}):=p(a_{1},\ldots,a_{n}|x_{1},\ldots,x_{n})\\},$
(4)

where $\vec{a}$ and $\vec{x}$ are shorthands for $a_{1},\ldots,a_{n}$ and
$x_{1},\ldots,x_{n}$, respectively.

To consider correlations as local in a given scenario, we require that the
outcome obtained by one party depends neither on the outcomes nor on the
measurements of the other parties. It is still allowed that the probabilities
can be correlated by some shared information represented by the variable
$\lambda$, and thus we define local correlations as those that admit the
following decomposition:

$p(\vec{a}|\vec{x})=\sum_{\lambda}q_{\lambda}p_{A_{1}}(a_{1}|x_{1},\lambda)\ldots
p_{A_{n}}(a_{n}|x_{n},\lambda),$ (5)

where $q_{\lambda}$ is some probability distribution and
$p_{A_{i}}(a_{i}|x_{i},\lambda)$ are probability distributions representing
the action of each party $A_{i}$. Correlations that cannot be written in the
form (5) are termed Bell nonlocal. Then, the state $\rho$ is called nonlocal
if it gives rise to nonlocal correlations in some Bell scenario.

A more suitable notion in the multipartite case is that of genuine
nonlocality. Consider again a bipartition $S|\bar{S}$ with nonempty sets $S$
and $\bar{S}$. We call a correlation $\vec{p}$ local across the bipartition
$S|\bar{S}$ if

$p(\boldsymbol{a}|\boldsymbol{x})=\sum_{\lambda}q_{\lambda}^{S|\bar{S}}p_{S}(\boldsymbol{a}_{S}|\boldsymbol{x}_{S},\lambda)p_{\bar{S}}(\boldsymbol{a}_{\bar{S}}|\boldsymbol{x}_{\bar{S}},\lambda),$
(6)

where $\boldsymbol{a}_{S}$ and $\boldsymbol{x}_{S}$ are vectors of outcomes
and measurement choices of parties belonging to the set $S$, and
$p_{S}(\boldsymbol{a}_{S}|\boldsymbol{x}_{S},\lambda)$ are probability
distributions that satisfy the non-signaling principle [26]. Then we say that
$\vec{p}$ is bilocal if it can be represented as

$\displaystyle p(\boldsymbol{a}|\boldsymbol{x})$
$\displaystyle=\sum_{S|\bar{S}}r_{S|\bar{S}}p_{S|\bar{S}}(\boldsymbol{a}|\boldsymbol{x}),$
(7)

where $p_{S|\bar{S}}(\boldsymbol{a}|\boldsymbol{x})$ are probability
distributions local across $S|\bar{S}$ (6). Finally, we refer to any $\vec{p}$
that does not admit the above decomposition as genuinely multipartite nonlocal
(GMNL) [13, 14].

Interestingly, one can introduce a more fine-grained classification of
nonlocality in the multipartite case to more accurately describe different
forms of nonlocality in correlations that are neither local nor GMNL. Notions
that suit this purpose are those of $k$-producibility of nonlocality and
nonlocality depth.

Let us consider a partition of the set $\mathcal{A}=\\{A_{1},\ldots,A_{n}\\}$
into $L$ disjoint nonempty subsets $S_{i}$ such that $S_{1}\cup\ldots\cup
S_{L}=\mathcal{A}$ and the size of each of them does not exceed $k$, i.e.,
$|S_{i}|\leqslant k$ for some positive integer $k=1,\ldots,n$. We refer to
this partition as an $L_{k}$ partition of $\mathcal{A}$. Consider now
correlations that admit the following decomposition:

$p(\boldsymbol{a}|\boldsymbol{x})=\sum_{\lambda}q_{\lambda}p_{S_{1}}(\boldsymbol{a}_{S_{1}}|\boldsymbol{x}_{S_{1}},\lambda)\cdot\ldots\cdot
p_{S_{L}}(\boldsymbol{a}_{S_{L}}|\boldsymbol{x}_{S_{L}},\lambda),$ (8)

where as before we assume that
$p_{S_{i}}(\boldsymbol{a}_{S_{i}}|\boldsymbol{x}_{S_{i}})$ satisfy non-
signaling conditions between parties in $S_{i}$. We call correlations
$\vec{p}$ that can be written in this way $k$-producible with respect to the
given $L_{k}$ partition. We then call $\vec{p}$ $k$-producible if it can be
decomposed into the following convex combination

$p(\boldsymbol{a}|\boldsymbol{x})=\sum_{L_{k}}q_{L_{k}}p_{L_{k}}(\boldsymbol{a}|\boldsymbol{x})$
(9)

of correlations $p_{L_{k}}(\boldsymbol{a}|\boldsymbol{x})$ that are
$k$-producible with respect to various $L_{k}$ partitions. Moreover, the
minimal $k$ for which $p(\boldsymbol{a}|\boldsymbol{x})$ can be expressed in
the form (9) is referred to as the nonlocality depth. In other words, we say
that $p(\boldsymbol{a}|\boldsymbol{x})$ has nonlocality depth $k$ if $k$ is
the minimal integer $k=1,\ldots,n$ for which the decomposition (9) exists.
Hence, for $k=N$, $p(\boldsymbol{a}|\boldsymbol{x})$ is GMNL.

###  II.2 Class of Bell inequalities for detection of GMNL

Let us now briefly recall some of the Bell inequalities constructed in Ref.
[21]. To this end, we first need to introduce the well-known CHSH Bell
inequality [22], which serves as a seed inequality for the construction. We
state it here in its quite unusual form

$I:=p(00|00)-p(01|01)-p(10|10)-p(00|11)\leqslant 0,$ (10)

which can be obtained from the standard form by using the nonsignaling
conditions.

Then, let us denote by $I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}$ the CHSH Bell
expression between the parties $A_{i}$ and $A_{j}$ lifted to $n$ parties in
such a way that the measurement choices and outcomes of the remaining $n-2$
parties are $x_{i}=0$ and $a_{i}=0$, that is,

$\displaystyle I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}=$ $\displaystyle
p(\boldsymbol{0}|\boldsymbol{0})-p(0\ldots 01_{i}0\ldots 0|0\ldots
01_{i}0\ldots 0)$ (11) $\displaystyle-p(0\ldots 01_{j}0\ldots 0|0\ldots
01_{j}0\ldots 0)$ $\displaystyle-p(\boldsymbol{0}|0\ldots 01_{i}0\ldots
01_{j}0\ldots 0),$

where $\boldsymbol{0}:=0\ldots 0$ is a shortcut notation for $n$ consecutive
zeros, and the subscripts $i,j$ refer to parties $A_{i}$ and $A_{j}$.
Importantly, for correlations $\vec{p}$ that are local across any bipartition
$S|\bar{S}$ (5) such that $A_{i}\in S$ and $A_{j}\in\bar{S}$, the expression
(11) satisfies $I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant 0$.

In this notation, the Bell inequalities detecting GMNL from Ref. [21] can be
stated as

$\sum_{i=1}^{n}\sum_{i<j=1}^{n}I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant\binom{n-1}{2}p(\boldsymbol{0}|\boldsymbol{0})$
(12)

and

$\sum_{j=2}^{n}I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant(n-2)p(\boldsymbol{0}|\boldsymbol{0}).$
(13)

Notice that the first inequality is invariant under permutation of any pair of
parties, whereas second one involves the lifted CHSH Bell expressions between
a fixed party $A_{1}$, and the remaining parties $A_{i}$ with $i=2,\ldots,n$.

It is worth mentioning that, as indicated in Ref. [21], any other, in general
$m$-partite Bell inequality with $m<n$ in the simplest Bell scenario with two
dichotomic measurements per party, which can be written as

$p(\boldsymbol{0}|\boldsymbol{0})-\sum_{\mathbf{a},\mathbf{x}\neq\mathbf{0}}\beta_{\mathbf{x}}^{\mathbf{a}}p(\mathbf{a}|\mathbf{x})\leqslant
0,$ (14)

where $\beta_{\mathbf{x}}^{\mathbf{a}}\geqslant 0$, can be used as a seed
inequality for the construction. In Ref. [21] a few examples of Bell
inequalities in the bipartite and multipartite scenarios were presented which
admit the above decomposition such as for instance the one found in Ref. [14]
which reads

$\displaystyle I^{1,2,3}_{\mathrm{tri}}:=$ $\displaystyle
p(000|000)-p(010|111)-p(000|011)$ (15)
$\displaystyle-p(001|001)-p(100|110)-p(010|010)$
$\displaystyle-p(100|100)\leqslant 0$

and detects GMNL in three-partite quantum systems. Taking this inequality as a
seed inequality for the construction of [21], the following inequality was
produced

$\sum_{i>2}I_{\mathrm{tri}}^{1,2,i}\leqslant(n-3)p(\mathbf{0}|\mathbf{0}).$
(16)

##  III Results

###  III.1 General construction

In this section, we show a general technique of formulating inequalities that
detect GMNL correlations, and in particular, we show how this technique can be
used to improve the inequalities of Ref. [21], in particular (13). To this
end, let us for a moment take a broader perspective and consider a general
expression composed of various $m$-partite Bell expressions

$\sum_{g\in G}\mathcal{I}^{g},$ (17)

where $G$ is an arbitrarily chosen set of subsets $g\subset\\{1,\ldots,n\\}$
such that $|g|=m$, and $\mathcal{I}^{g}$ are in general $m$-partite Bell
expressions detecting GMNL between parties from the set $\\{A_{i}\\}_{i\in g}$
lifted to the $n$-partite scenario for $n>m$ such as for instance the
expressions in (11).

Our technique requires that the expressions $\mathcal{I}^{g}$ fulfill a couple
of conditions.

  * (i)

First, for all correlations $\vec{p}$ that are local across the bipartition
$S|\bar{S}$ (6), we require that $\mathcal{I}^{g}\leqslant 0$ if there exists
a pair $i,j\in g$ such that $A_{i}\in S$ and $A_{j}\in\bar{S}$.

  * (ii)

To formulate the second condition, we first decompose $\mathcal{I}^{g}$ as

$\mathcal{I}^{g}=\mathcal{I}^{g}_{+}-\mathcal{I}^{g}_{-},$ (18)

where in both expressions $\mathcal{I}^{g}_{\pm}$ probabilities enter with
positive weights. In other words, we decompose each Bell expression into two
terms containing probabilities that enter the Bell expression with positive
and negative weights, respectively. Then, the second condition on expressions
$\mathcal{I}^{g}$ requires that the positive parts of every Bell expression
$\mathcal{I}^{g}$ are the same, that is,

$\forall_{g\in G}\qquad\mathcal{I}^{g}_{+}=\mathcal{I}_{+}.$ (19)

Let us also denote by $T$ the common terms appearing in all the negative parts
$\mathcal{I}_{-}^{g}$. All together this implies that for all $g\in G$ we have

$\mathcal{I}^{g}\leqslant\mathcal{I}_{+}-T$ (20)

Our aim now is to determine an upper bound for the expression (17) for bilocal
models (7). For this purpose, we consider a particular bipartition
$S|\bar{S}$. Then, for all $\vec{p}$ local across the bipartition $S|\bar{S}$
(6) we have

$\sum_{g\in G}\mathcal{I}^{g}\leqslant\sum_{g\in G^{\prime}}\mathcal{I}^{g},$
(21)

where $G^{\prime}\subseteq G$ is a set of such $g$ that either
$\\{A_{i}\\}_{i\in g}\subset S$ or $\\{A_{i}\\}_{i\in g}\subset\bar{S}$. Next,
using the decomposition (19), and Ineq. (20), the above can further be upper
bounded by

$\sum_{g\in
G}\mathcal{I}^{g}\leqslant|G^{\prime}|\left(\mathcal{I}_{+}-T\right).$ (22)

This leads us to our first result which we formulate as the following theorem.

######  Theorem 1.

Consider a Bell expression (17) such that the conditions (i) and (ii) are
satisfied by $I^{g}$ for each $g\in G$. Then, violation of the following Bell
inequality

$\sum_{g\in G}\mathcal{I}^{ij}\leqslant\gamma\left(\mathcal{I}_{+}-T\right)$
(23)

certifies GMNL, where $\gamma=\max_{S|\bar{S}}|G^{\prime}|$ and $T$ is a
common term appearing in all the negative parts $\mathcal{I}_{-}^{g}$ under
the decomposition (18).

###### Proof.

To prove the above theorem first notice that Ineq. (23) is derived from Ineq.
(22) by maximizing over all bipartitions $S|\bar{S}$, meaning that Ineq. (23)
holds true for all $\vec{p}$ local across any $S|\bar{S}$. Next, notice that
all $\vec{p}$ admitting Eq. (7) are a convex combination of probability
distributions $p_{S|\bar{S}}(\boldsymbol{a}|\boldsymbol{x})$ local across
bipartitions $S|\bar{S}$. And since both sides of Ineq. (23) are linear
functions of $\vec{p}$, the violation of Ineq. (23) implies that the
correlation $\vec{p}$ does not admit Eq. (7), or in other words, that
$\vec{p}$ is GMNL. ∎

Example 1. Let us now show how the above considerations lead to an improvement
of Ineq. (13) derived in Ref. [21]. In this case $\mathcal{I}^{g}$ are simply
the lifted CHSH Bell expressions $I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}$
(11) and $G=\\{\\{1,i\\}\\}_{i=2}^{n}$ is the set of all pairs of parties
$A_{1}A_{i}$ for $i>1$. Then, it directly follows that
$\mathcal{I}^{+}=p(\boldsymbol{0}|\boldsymbol{0})$ and at the same time, one
finds that in this case

$T=p(1\mathbf{0}_{n-1}|1\mathbf{0}_{n-1}),$ (24)

where $\mathbf{0}_{n-1}$ is a shorthand for a string of $0$ of length $n-1$.
All this implies that the improved inequality is

$\displaystyle\sum_{j=2}^{n}I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant(n-2)[p(\mathbf{0}|\mathbf{0})-p(1\mathbf{0}_{n-1}|1\mathbf{0}_{n-1})].$
(25)

In order to demonstrate that the above inequality is tighter than that in (13)
we consider a mixture of the $n$-qubit GHZ state (1) and white noise,

$\rho(q)=(1-q)|\mathrm{GHZ}_{n,2}\rangle\\!\langle\mathrm{GHZ}_{n,2}|+q\frac{\mathds{1}}{2^{n}}$
(26)

with $q\in[0,1]$. The numerical comparison between noise levels needed for
violation of either inequality has been presented in Fig. 1. The numerical
analysis shows that, in the case of a three-qubit GHZ state, the violation of
Ineq. (13) is achieved for noise levels of $q<0.05$ while for the new
inequality (25) we can observe violation for $q<0.106$.

Figure 1: The numerical comparison of noise robustness to the white noise of
inequalities (13) (black triangles) and (25) (blue disks), for a noisy GHZ
state (26). The variable $q$ corresponds to the maximal amount of white noise
for which a given inequality detects GMNL of the state $\rho(q)$ (26), while
$n$ denotes the number of qubits. We see that even though the noise robustness
of both inequalities drops with $n$, the new inequality is still outperforming
the old one for every $n$.

The maximal violations of the inequalities were found assuming projective
measurements. Each was parameterized by two real angles and then the standard
methods in Mathematica were used to find the maxima.

Example 2. In a similar vein one can improve the inequality (16). Indeed, by
following the prescription above one realizes that the following inequality is
satisfied by all non-GMNL correlations

$\displaystyle\sum_{i>2}I_{\mathrm{tri}}^{1,2,i}$
$\displaystyle\leqslant(n-3)[p(\boldsymbol{0}|\boldsymbol{0})-p(01\boldsymbol{0}_{n-2}|01\boldsymbol{0}_{n-2})$
(27)
$\displaystyle-p(10\boldsymbol{0}_{n-2}|10\boldsymbol{0}_{n-2})-p(10\boldsymbol{0}_{n-2}|11\boldsymbol{0}_{n-2})],$

where $\mathcal{I}_{+}=p(\boldsymbol{0}|\boldsymbol{0})$ and

$\displaystyle T=$ $\displaystyle
p(01\boldsymbol{0}_{n-2}|01\boldsymbol{0}_{n-2})-p(10\boldsymbol{0}_{n-2}|10\boldsymbol{0}_{n-2})$
(28) $\displaystyle-p(10\boldsymbol{0}_{n-2}|11\boldsymbol{0}_{n-2}).$

###  III.2 All tripartite, pure GME states are GMNL

Let us now show that the improved inequality (25) is capable of detecting
genuine nonlocality in any three-qubit pure genuinely entangled. Whereas, the
equivalence between genuine entanglement and genuine nonlocality for three-
qubit pure states was already established with a Hardy-type argument [15],
here we reproduce this result with a single Bell inequality.

######  Theorem 2.

Every triqubit, pure, genuinely multipartite state violates Ineq. (25).

While the detailed proof can be found in Appendix A, here we outline the main
idea behind it. To detect GMNL we use the improved inequality (25) which for
$n=3$ can be rewritten as

$\displaystyle p(0$ $\displaystyle 00|000)-p(100|100)-p(010|010)$ (29)
$\displaystyle-p(001|001)-p(000|110)-p(000|101)\leqslant 0.$

To simplify our considerations, we make use of the generalized Schmidt
decomposition for three-qubit states introduced in Ref. [27], which for the
case $n=3$, states that every state $|\psi\rangle\in(\mathbb{C}^{2})^{\otimes
3}$ can be decomposed as

$|\psi\rangle=a\mathrm{e}^{\mathbbm{i}\varphi}|000\rangle+b|011\rangle+c|101\rangle+d|110\rangle+e|111\rangle,$
(30)

where $\mathrm{e}$ is the Euler’s number, $\varphi\in\mathbb{R}$, and
$a,b,c,d,e$ are non-negative coefficients satisfying $a\geqslant x$ for all
$x\in\\{b,c,d,e\\}$, and $a^{2}+b^{2}+c^{2}+d^{2}+e^{2}=1$ . Moreover, notice
that the violation of Ineq. (13) implies the violation of Ineq. (25).
Therefore, since the violation of Ineq. (13) was used to show that all pure
GME states, symmetric with respect to a permutation of the second and third
parties are GMNL [21], we can limit ourselves only to states which are not
invariant under such permutation.

To show a violation of Ineq. (29) we chose measurements parametrized by a
coefficient $\alpha\in[0,\pi/2]$ such that for all $\alpha$

$p(001|001)=p(010|010)=p(000|101)=0.$ (31)

Then, we show that for all non-symmetric states (30) there exists $\alpha$ for
which

$\displaystyle p(0$ $\displaystyle 00|000)-p(100|100)-p(000|110)>0,$ (32)

which implies the violation of Ineq. (29).

Interestingly, since this holds for all tripartite, pure states that are not
symmetric with respect to a permutation of two parties, it implies that all
such states are GMNL, and so also GME.

###  III.3 Detection of nonlocality depth

Interestingly, we can easily generalize the approach of Theorem 1 to construct
inequalities detecting nonlocality depth. In fact, the following theorem holds
true.

######  Theorem 3.

Consider a Bell expression (17) such that the conditions (i) and (ii) stated
in Sec. III.1 are satisfied by $I^{g}$ for each $g\in G$. Then, violation of
the following Bell inequality

$\sum_{g\in
G}\mathcal{I}^{g}\leqslant\gamma_{k}\left(\mathcal{I}_{+}-T\right),$ (33)

by correlations $p(\boldsymbol{a}|\boldsymbol{x})$, where
$\gamma_{k}=\max_{S_{1},S_{2},\ldots,S_{L}}|G^{\prime}|$ for all
$|S_{i}|\leqslant k$, implies those correlations have nonlocality depth $k+1$.
As before, $T$ is a common term appearing in all the negative parts
$\mathcal{I}_{-}^{g}$ under the decomposition (18).

###### Proof.

This theorem can be proven in the same way as Theorem 1 and so we do not
repeat it here. ∎

Let us now illustrate Theorem 3 with a couple of examples.

Example 3. First, let us consider the expression in Eq. (13), that is,

$\sum_{j=2}^{n}I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}$ (34)

as well as an $L_{k}$-partition of the set $\mathcal{A}$. Let us then denote
by $S_{1}$ the set such that $A_{1}\in S_{1}$. Then, we can bound each term
$I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}$ in a similarly way as in the proof
of inequality (25): for all $A_{j}\notin S_{1}$ we have
$I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}\leq 0$, while for $A_{i}\in S_{1}$ we
have $I^{1,i}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant
p(\mathbf{0}|\mathbf{0})-p(1\boldsymbol{0}_{n-1}|1\boldsymbol{0}_{n-1})$.
Since the size of $S_{1}$ is upper-bounded by $k$, the following inequality
holds true that for all $k$-producible $p(\boldsymbol{a}|\boldsymbol{x})$:

$\sum_{j=2}^{n}I^{1,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant(k-1)[p(\mathbf{0}|\mathbf{0})-p(1\mathbf{0}_{n-1}|1\mathbf{0}_{n-1})].$
(35)

Violation of this inequality by some $\vec{p}$ implies that the nonlocality
depth of this correlation is at least $k+1$.

Example 4. Interestingly, a similar idea can be applied to Ineq. (12) to
formulate another inequality capable of detecting nonlocality depth. In this
case, however, the derivation of the maximal value over the $k$-producible
correlations is slightly more difficult.

The task is to find the $L_{k}$ partition that maximizes the number of
$I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}$ terms between parties from the same
subset, which can be upper-bounded by
$I_{\mathbf{0}|\mathbf{0}}^{i,j}\leqslant p(\mathbf{0}|\mathbf{0})$. The
optimal strategy for maximizing this number is to take the partition for which
$|S_{i}|=k$ for all $i\in\\{1,\ldots,L-1\\}$ and $|S_{L}|=n\mod k$, where
$L=\lceil n/k\rceil$. We thus obtain the following inequality

$\sum_{i}\sum_{j>i}I^{i,j}_{\boldsymbol{0}|\boldsymbol{0}}\leqslant\left[\left\lfloor\frac{n}{k}\right\rfloor\binom{k}{2}+(n_{k})\binom{n_{k}}{2}\right]p(\boldsymbol{0}|\boldsymbol{0}),$
(36)

where $n_{k}\equiv n\mod k$, violation of which detects nonlocality depth
$k+1$.

###  III.4 A possible generalization to more outcomes

In this section we provide a possible way of generalizing the above
inequalities detecting GMNL to Bell scenarios involving more outcomes,
however, concentrating on the simplest case of $d=3$. The most natural
candidate for the seed inequality is the CGLMP Bell inequality [24], which in
the general $d$-outcome can be written in a form highly resembling the one in
Eq. (10) as [28, 29],

$\displaystyle P(A_{0}$ $\displaystyle<B_{0})-P(A_{0}<B_{1})$ (37)
$\displaystyle-P(B_{1}<A_{1})-P(A_{1}<B_{0})\leqslant 0,$

where

$P(A_{i}<B_{j})=\sum_{\begin{subarray}{c}a,b=0\\\
a<b\end{subarray}}^{d-1}p(ab|ij)\qquad(i,j=0,1).$ (38)

In the particular case of $d=3$, after relabeling the outcomes of $B_{0}$ such
that $0\leftrightarrow 2$, the above expands as

$\displaystyle J_{3}$ $\displaystyle=$ $\displaystyle
p(01|00)+p(00|00)+p(10|00)$ (39) $\displaystyle-p(01|01)-p(02|01)-p(12|01)$
$\displaystyle-p(10|11)-p(20|11)-p(21|11)$
$\displaystyle-p(01|10)-p(00|10)-p(10|10)\leqslant 0.$

We work in the tripartite scenario and hence the $\mathcal{I}^{g}$ expressions
will be denoted more conveniently by $\mathcal{I}^{AB}$, $\mathcal{I}^{AC}$
and $\mathcal{I}^{BC}$. To satisfy the condition that all $\mathcal{I}^{g}$’s
entering the expression (17) have the same positive part we will add two
expressions together, one lifted by $1|0$ and another by $0|0$. Yet, the
expression (39) has a different number of positive terms with outcomes $1$ and
$0$ on a given position. In order to equalize these terms we will add
$\tilde{J}_{3}$ which is obtained from $J_{3}$ by permuting the outcomes $0$
and $1$, that is,

$\displaystyle\tilde{J}_{3}$ $\displaystyle=$ $\displaystyle
p(10|00)+p(11|00)+p(01|00)$ (40) $\displaystyle-p(10|01)-p(12|01)-p(02|01)$
$\displaystyle-p(01|11)-p(21|11)-p(20|11)$
$\displaystyle-p(10|10)-p(11|10)-p(01|10).$

We then consider the following expression

$\mathcal{I}^{i,j}=J_{3}^{i,j}({1|0})+\tilde{J}_{3}^{i,j}({0|0}),$ (41)

where $J_{3}(1|0)$ is the expression $J_{3}$ lifted to the third party with
measurement $0$ and outcome $1$, whereas $\tilde{J}_{3}(0|0)$ is obtained by
lifting $\tilde{J}_{3}$ to the third party with measurement choice and outcome
set to 0. In this way we achieved that the positive terms in
$\mathcal{I}^{AB}$, $\mathcal{I}^{BC}$ and $\mathcal{I}^{AC}$ are the same,

$\displaystyle\mathcal{I}_{+}$ $\displaystyle=$
$\displaystyle\mathcal{I}_{+}^{AB}=\mathcal{I}_{+}^{AC}=\mathcal{I}_{+}^{BC}$
$\displaystyle=$ $\displaystyle p(001|000)+p(010|000)+p(011|000)$
$\displaystyle+p(100|000)+p(101|000)+p(110|000).$

Additionally $\mathcal{I}^{AB}$ and $\mathcal{I}^{AC}$ share some common
negative terms

$\displaystyle T$ $\displaystyle=$ $\displaystyle
p(001|100)+p(010|100)+p(011|100)$ (43)
$\displaystyle+p(100|100)+p(101|100)+p(110|100).$

Now, we have everything to construct Bell inequalities which analogously to
inequalities (12) and (13) witness GMNL in $(3,2,3)$ scenario,

$\mathcal{I}^{AB}+\mathcal{I}^{AC}+\mathcal{I}^{BC}\leqslant\mathcal{I}^{+}$
(44)

and

$\mathcal{I}^{AB}+\mathcal{I}^{AC}\leqslant\mathcal{I}_{+}-T.$ (45)

Both these inequalities can be violated using a three-qutrit GHZ state
$|\mathrm{GHZ}_{3,3}\rangle$ [cf. Eq. (1)]. Inequality (45) was further tested
on pure states drawn randomly with the Haar measure (notice that biseparable
states form a set of measure zero), which are invariant under a permutation of
two parties. We tested $100$ states and in each case, we managed to find
measurements such that obtained correlations violated (45). This may suggest
that the analytical results from Ref. [21] that every pure tripartite GME
state which is invariant under permutation of two the parties $B$ and $C$ is
GMNL, generalize to the tripartite qutrit states.

It is also worth pointing out that the numerically found violations are much
more noise-sensitive as compared to the qubit scenario. For example, the
3-qutrit GHZ state violated inequality (45) with maximal $q\approx 0.02\%$
(26), whereas in the case of three-qubit GHZ more than $10\%$ of white noise
could be added to still detect GMNL.

##  IV Conclusions

In this work, we showed how one can easily improve one of the classes of Bell
inequalities witnessing genuine nonlocality in the multipartite scenario
introduced recently in Ref. [21]. We have shown that using the improved
inequality one can prove the genuine nonlocality of every tripartite, pure,
genuinely entangled state. We also demonstrated how the inequalities of Ref.
[21] can be turned into ones detecting nonlocality depth. We finally presented
a possible approach to generalizing the inequalities of Ref. [21] to the case
of a higher number of outcomes.

There is a couple of ways in which the above research line can be continued.
First, generalizing our proof for tripartite genuinely entangled states, the
inequalities we introduced can perhaps be used to prove that all pure,
genuinely entangled states are genuinely nonlocal which is a long-standing
conjecture. Next, there is a question of how to maximally increase the common
negative term giving a greater noise-robustness of the inequality. Here, a
couple of possible approaches are possible. One could, for example, rewrite
the ”seed” inequalities using non-signaling assumption, or perhaps lift these
inequalities using different measurement settings and results, akin to what we
did for $3$-outcome inequality (45).

Lastly, there is a possibility of adapting this technique for different
scenarios. One example of such would be a network scenario, either with
independent sources (see e.g. [30]) or in the local operations and shared
randomness model [31]. For both of these, such a technique could prove highly
desirable, as derivation of suitable Bell inequality is a notoriously
difficult problem in the network scenario.

## References

  * Acín et al. [2007] A. Acín, N. Brunner, N. Gisin, S. Massar, S. Pironio, and V. Scarani, Physical Review Letters 98 (2007), 10.1103/physrevlett.98.230501.
  * Pironio et al. [2010] S. Pironio, A. Acín, S. Massar, A. B. de la Giroday, D. N. Matsukevich, P. Maunz, S. Olmschenk, D. Hayes, L. Luo, T. A. Manning, and C. Monroe, Nature 464, 1021–1024 (2010).
  * Mayers and Yao [2004] D. Mayers and A. Yao, Quantum Info. Comput. 4, 273–286 (2004).
  * Šupić and Bowles [2020] I. Šupić and J. Bowles, Quantum 4, 337 (2020).
  * Werner [1989] R. F. Werner, Phys. Rev. A 40, 4277 (1989).
  * Augusiak et al. [2015] R. Augusiak, M. Demianowicz, J. Tura, and A. Acín, Phys. Rev. Lett. 115, 030404 (2015).
  * Augusiak et al. [2014] R. Augusiak, M. Demianowicz, and A. Acín, Journal of Physics A: Mathematical and Theoretical 47, 424002 (2014).
  * Bowles et al. [2016] J. Bowles, J. Francfort, M. Fillettaz, F. Hirsch, and N. Brunner, Physical Review Letters 116 (2016), 10.1103/physrevlett.116.130401.
  * Gisin [1991] N. Gisin, Physics Letters A 154, 201 (1991).
  * Cabello [2002] A. Cabello, Physical Review A 65 (2002), 10.1103/physreva.65.062105.
  * Popescu and Rohrlich [1992] S. Popescu and D. Rohrlich, Physics Letters A 166, 293 (1992).
  * Seevinck and Uffink [2001] M. Seevinck and J. Uffink, Phys. Rev. A 65, 012107 (2001).
  * Svetlichny [1987] G. Svetlichny, Phys. Rev. D 35, 3066 (1987).
  * Bancal et al. [2013] J.-D. Bancal, J. Barrett, N. Gisin, and S. Pironio, Phys. Rev. A 88, 014102 (2013).
  * Yu and Oh [2013] S. Yu and C. H. Oh, “Tripartite entangled pure states are tripartite nonlocal,” (2013), arXiv:1306.5330 [quant-ph] .
  * Chen et al. [2014] Q. Chen, S. Yu, C. Zhang, C. Lai, and C. H. Oh, Physical Review Letters 112 (2014), 10.1103/PhysRevLett.112.140404.
  * Bancal et al. [2009] J.-D. Bancal, C. Branciard, N. Gisin, and S. Pironio, Physical Review Letters 103 (2009), 10.1103/physrevlett.103.090503.
  * Aolita et al. [2012] L. Aolita, R. Gallego, A. Cabello, and A. Acín, Phys. Rev. Lett. 108, 100401 (2012).
  * Augusiak et al. [2019] R. Augusiak, A. Salavrakos, J. Tura, and A. Acín, New Journal of Physics 21, 113001 (2019).
  * Pandit et al. [2022] M. Pandit, A. Barasiński, I. Márton, T. Vértesi, and W. Laskowski, New Journal of Physics 24, 123017 (2022).
  * Curchod et al. [2019] F. J. Curchod, M. L. Almeida, and A. Acín, New Journal of Physics 21, 023016 (2019).
  * Clauser et al. [1969] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Phys. Rev. Lett. 23, 880 (1969).
  * Contreras-Tejada et al. [2021] P. Contreras-Tejada, C. Palazuelos, and J. I. de Vicente, Physical Review Letters 126 (2021), 10.1103/physrevlett.126.040501.
  * Collins et al. [2002] D. Collins, N. Gisin, N. Linden, S. Massar, and S. Popescu, Physical Review Letters 88 (2002), 10.1103/physrevlett.88.040404.
  * Greenberger et al. [1989] D. M. Greenberger, M. A. Horne, and A. Zeilinger, “Going beyond bell’s theorem,” in Bell’s Theorem, Quantum Theory and Conceptions of the Universe, edited by M. Kafatos (Springer Netherlands, Dordrecht, 1989) pp. 69–72.
  * Barrett et al. [2005] J. Barrett, N. Linden, S. Massar, S. Pironio, S. Popescu, and D. Roberts, Physical Review A 71 (2005), 10.1103/physreva.71.022101.
  * Carteret et al. [2000] H. A. Carteret, A. Higuchi, and A. Sudbery, Journal of Mathematical Physics 41, 7932–7939 (2000).
  * Zohren and Gill [2008] S. Zohren and R. D. Gill, Phys. Rev. Lett. 100, 120406 (2008).
  * Zohren et al. [2010] S. Zohren, P. Reska, R. D. Gill, and W. Westra, EPL (Europhysics Letters) 90, 10002 (2010).
  * Branciard et al. [2012] C. Branciard, D. Rosset, N. Gisin, and S. Pironio, Physical Review A 85 (2012), 10.1103/physreva.85.032119.
  * Coiteux-Roy et al. [2021] X. Coiteux-Roy, E. Wolfe, and M.-O. Renou, Phys. Rev. Lett. 127, 200401 (2021).

##  Appendix A Proof of Theorem 2

In this appendix, we present a detailed proof of Thm. 2, which we recall below

######  Theorem 2.

Every tripartite, pure, genuinely multipartite state violates Ineq. (25).

###### Proof.

By the virtue of [27, Theorem 3], every pure, $3$-partite, qubit state can be
written as

$|\psi\rangle_{ABC}=a\mathrm{e}^{\mathbbm{i}\varphi}|000\rangle+b|011\rangle+c|101\rangle+d|110\rangle+e|111\rangle,$
(46)

where $a,b,c,d,e$ are nonnegative and $a\geqslant x$ for $x\in\\{b,c,d,e\\}$,
and by $\mathrm{e}=2.718\dots$ we denote the Euler’s number. Notice, that by
permuting parties, we can rearrange $b,c,d$ without affecting other terms.
Therefore, without a loss of generality, we can take $b\geqslant c\geqslant
d$.

Furthermore, in Ref. [21] it was shown that the inequality (13) is violated by
all pure, genuinely multipartite entangled states that are invariant under the
permutation of a pair of qubits. As the new inequality (25) is tighter than
(13), it directly follows that it also has to be violated by such states.
Therefore, we can restrict ourselves to states that are not invariant under
the permutation of any two qubits, which implies that $b>c>d$.

For the case of $n=3$, Ineq. (25) can be rewritten as

$p(000|000)-p(100|100)-p(010|010)-p(001|001)-p(000|110)-p(000|101)\leqslant
0.$ (47)

Therefore, to complete the proof we have to show that for all parameters
$a,b,c,d,e$ satisfying

$a\geqslant b>c>d\geqslant 0,\qquad a\geqslant e\geqslant 0,\qquad
a^{2}+b^{2}+c^{2}+d^{2}+e^{2}=1,$ (48)

there exist measurements which performed on the state (46) lead to the
violation of Ineq. (47).

To this end, we assume that each party performs projective measurements which
we denote by
$M_{i|j}^{(k)}=\big{|}M_{i|j}^{(k)}\big{\rangle}\\!\big{<}M_{i|j}^{(k)}\big{|}$,
where $i$ is the outcome of the measurement $j$ performed by the party
$A_{k}$, and we take the states $\big{|}M_{i|j}^{(k)}\big{\rangle}$ to be

$\displaystyle\big{|}M_{0|0}^{(1)}\big{\rangle}$
$\displaystyle=\cos{\alpha}|0\rangle+\sin{\alpha}|1\rangle,$
$\displaystyle\hskip 42.67912pt\big{|}M_{0|1}^{(1)}\big{\rangle}$
$\displaystyle=\frac{1}{\eta_{1}}\left(-c^{2}\sin\alpha|0\rangle+a^{2}\cos\alpha|1\rangle\right),$
(49) $\displaystyle\big{|}M_{0|0}^{(2)}\big{\rangle}$
$\displaystyle=|0\rangle,$ $\displaystyle\big{|}M_{0|1}^{(2)}\big{\rangle}$
$\displaystyle=\frac{1}{\eta_{2}}\left(c\sin\alpha|0\rangle+\left(b\cos\alpha+e\sin\alpha\right)|1\rangle\right),$
$\displaystyle\big{|}M_{0|0}^{(3)}\big{\rangle}$ $\displaystyle=|1\rangle,$
$\displaystyle\big{|}M_{0|1}^{(3)}\big{\rangle}$
$\displaystyle=\frac{1}{\eta_{3}}\left(\mathrm{e}^{\mathbbm{i}\phi}a\cos\alpha|0\rangle+c\sin\alpha|1\rangle\right),$

where $\alpha\in[0,\pi/2]$ is a free parameter, and $\eta_{k}$ for $k=1,2,3$
are normalization constants

$\eta_{1}^{2}=c^{4}\sin^{2}\alpha+a^{4}\cos^{2}\alpha,\qquad\eta_{2}^{2}=c^{2}\sin^{2}{\alpha}+\left(b\cos{\alpha}+e\sin{\alpha}\right)^{2},\qquad\eta_{3}^{2}=a^{2}\cos^{2}\alpha+c^{2}\sin^{2}\alpha.$
(50)

Notice, that measurement operators for outcome $1$ can be trivially derived
from the relation $M_{0|j}^{(k)}+M_{1|j}^{(k)}=\mathbb{1}$ for all $j,k$. It
is easy to check, that for this particular choice of measurements, we have

$p(001|001)=p(010|010)=p(000|101)=0,$ (51)

while the remaining probabilities in Ineq. (47) are given by

$\begin{gathered}p(000|000)=c^{2}\sin^{2}\alpha,\qquad
p(100|100)=\frac{1}{\eta_{1}^{2}}c^{6}\sin^{2}\alpha,\\\
p(000|110)=\frac{1}{\eta_{1}^{2}\eta_{2}^{2}}\left[a^{2}c^{2}\sin\alpha\cos\alpha+\left(a^{2}e\cos\alpha-
bc^{2}\sin\alpha\right)\left(b\cos\alpha+e\sin\alpha\right)\right]^{2}.\end{gathered}$
(52)

And so, after substituting the above probabilities to Ineq. (47) we get

$c^{2}\sin^{2}\alpha-\frac{1}{\eta_{1}^{2}}c^{6}\sin^{2}\alpha-\frac{1}{\eta_{1}^{2}\eta_{2}^{2}}\left[a^{2}c^{2}\sin\alpha\cos\alpha+\left(a^{2}e\cos\alpha-
bc^{2}\sin\alpha\right)\left(b\cos\alpha+e\sin\alpha\right)\right]^{2}\leqslant
0.$ (53)

We need to consider two separate cases to prove the violation of this
inequality. First, let us assume that $e=0$. Under this assumption we can take
$\alpha=\pi/4$ which yields

$c^{2}\frac{1}{2}-\frac{1}{2\eta_{1}^{2}}c^{6}\
-\frac{1}{4\eta_{1}^{2}\eta_{2}^{2}}\left(a^{2}c^{2}-b^{2}c^{2}\right)^{2}\leqslant
0.$ (54)

Next, we multiply the above by $4\eta_{1}^{2}\eta_{2}^{2}$, and we use Eq.
(50) to express $\eta_{1}$ and $\eta_{2}$ in terms of the coefficients
$a,b,c,e,\alpha$, which gives us

$\frac{1}{2}(b^{2}+c^{2})(a^{2}-c^{2})(a^{2}+c^{2})-c^{2}\left(a^{2}-b^{2}\right)^{2}\leqslant
0.$ (55)

However, from Ineq. (48) it follows that

$\frac{1}{2}(b^{2}+c^{2})>c^{2}\quad\textrm{and}\quad
a^{2}-c^{2}>a^{2}-b^{2}\quad\textrm{and}\quad a^{2}+a^{2}>a^{2}-b^{2},$ (56)

implying that the left-hand side of Ineq. (55) is always positive, which
implies the violation of Ineq. (47).

For the case $e>0$ let us first focus on the first two terms of Ineq. (54)

$c^{2}\sin^{2}\alpha-\frac{1}{\eta_{1}^{2}}c^{6}\sin^{2}\alpha.$ (57)

It is easy to see that the above is positive for all $\alpha\in]0,\pi/2[$, and
so to show the violation of Ineq. (54) it is sufficient to show that there
exists $\alpha\in]0,\pi/2[$ for which the third term from Ineq. (54) equals
$0$

$a^{2}c^{2}\sin\alpha\cos\alpha+\left(a^{2}e\cos\alpha-
bc^{2}\sin\alpha\right)\left(b\cos\alpha+e\sin\alpha\right)=0.$ (58)

Observe, that the left-hand side is strictly positive for $\alpha=0$ and
strictly negative for $\alpha=\pi/2$, and since the function is continuous for
all $\alpha\in[0,\pi/2]$, there has to exist $\alpha\in]0,\pi/2[$ for which
the above equation holds true. Therefore, Ineq. (47) is violated.

Lastly, let us note that this proof implies that all states of the form (46)
that are not invariant under the permutation of a pair of qubits are GMNL, and
so also GME. And since the GME states that are invariant under such operations
were proven to be GMNL in [21], we have that all pure GME states are GMNL. ∎


11institutetext: Tuorla Observatory, Department of Physics and Astronomy,
20014 University of Turku, Finland
11email<EMAIL_ADDRESS>22institutetext: Anton Pannekoek Institute
for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam,
the Netherlands

# Polarized radiation from an accretion shock in accreting millisecond pulsars
using exact Compton scattering formalism

Anna Bobrikova 11 Vladislav Loktev 11 Tuomo Salmi 2211 Juri Poutanen 11

(Received 8 May 2023 / Accepted 11 August 2023)

Pulse profiles of accreting millisecond pulsars can be used to determine
neutron star (NS) parameters, such as their masses and radii, and therefore
provide constraints on the equation of state of cold dense matter. Information
obtained by the Imaging X-ray Polarimetry Explorer (IXPE) can be used to
decipher pulsar inclination and magnetic obliquity, providing ever tighter
constraints on other parameters. In this paper, we develop a new emission
model for accretion-powered millisecond pulsars based on thermal
Comptonization in an accretion shock above the NS surface. The shock structure
was approximated by an isothermal plane-parallel slab and the Stokes
parameters of the emergent radiation were computed as a function of the zenith
angle and energy for different values of the electron temperature, the Thomson
optical depth of the slab, and the temperature of the seed blackbody photons.
We show that our Compton scattering model leads to a significantly lower
polarization degree of the emitted radiation compared to the previously used
Thomson scattering model. We computed a large grid of shock models, which can
be combined with pulse profile modeling techniques both with and without
polarization included. In this work, we used the relativistic rotating vector
model for the oblate NS in order to produce the observed Stokes parameters as
a function of the pulsar phase. Furthermore, we simulated the data to be
produced by IXPE and obtained constraints on model parameters using nested
sampling. The developed methods can also be used in the analysis of the data
from future satellites, such as the enhanced X-ray Timing and Polarimetry
mission.

###### Key Words.:

methods: numerical – polarization – stars: neutron – techniques: polarimetric
– X-rays: binaries

##  1 Introduction

Rapidly rotating accreting neutron stars (NSs) can be used to study the
properties of the extremely dense matter inside their cores. Some of these are
accretion-powered millisecond pulsars (AMPs), which show X-ray pulsations at
the stellar spin frequency. Their emission profiles can be modeled, accounting
for light bending and other relativistic effects, and from these profiles, the
information about the mass and radius of the NS can be extracted(see e.g.,
Pechenick et al. 1983; Miller & Lamb 1998; Poutanen & Gierliński 2003;
Poutanen & Beloborodov 2006; Morsink et al. 2007; Lo et al. 2013; Watts et al.
2016; Salmi et al. 2018; Bogdanov et al. 2019). The mass and radius
constraints can be translated to the equation of state (EOS) constraints for
matter inside the NS core (e.g., see Lindblom 1992; Lattimer 2012; Hebeler et
al. 2013; Baym et al. 2018).

This technique was recently applied to the case of rotation-powered
millisecond pulsars using data from the Neutron Star Interior Composition
Explorer (NICER, Miller et al. 2019; Riley et al. 2019; Miller et al. 2021;
Riley et al. 2021; Salmi et al. 2022). In these sources, the hot regions on
the NS surface are heated by a magnetospheric return current and the emission
escaping the NS surface can be computed using self-consistent atmosphere
models (see e.g., Zavlin et al. 1996; Ho & Lai 2001; Heinke et al. 2006; Ho &
Heinke 2009; Haakonsen et al. 2012; Salmi et al. 2020). For AMPs, these models
are not adequate because the emitted radiation is Comptonized in an accretion
“shock” formed above the hot region. In reality, a standard gasdynamical shock
might not exist at all, but the kinetic energy of incoming particles is
released in the surface layer through Coulomb collisions and plasma
instabilities (Zel’dovich & Shakura 1969; Suleimanov et al. 2018; Salmi et al.
2020). Comptonization of soft photons generated within the layer as well as by
the underlying NS surface is then responsible for the power-law spectra
extending up to 100 keV as observed from AMPs (Poutanen & Gierliński 2003;
Gierliński & Poutanen 2005; Falanga et al. 2005a, b, 2007, 2011, 2012).
Electron scattering also causes the radiation to be significantly polarized.
Variations of polarization with pulsar phase can be used to determine
geometrical parameters such as inclination angle and magnetic obliquity
(Viironen & Poutanen 2004), which in their turn allow us to improve
constraints on the NS mass and radius.

Previously, polarized emission models for AMPs used the Thomson scattering
approximation in an optically thin NS atmosphere (see Sunyaev & Titarchuk
1985; Viironen & Poutanen 2004; Salmi et al. 2021). However, the polarization
degree due to Compton scattering on hot electrons strongly depends on the
electron temperature (Nagirner & Poutanen 1994; Poutanen 1994); for example,
for 100 keV electrons, polarization is less than 50% of that for Thomson
scattering. On the other hand, the models for nonpolarized radiation —in the
context of pulse profile modeling— typically employ approximate formulas for
the anisotropy of the radiation and empirical models for the Comptonized
spectra (see e.g., Poutanen & Gierliński 2003; Leahy et al. 2008; Steiner et
al. 2009; Salmi et al. 2018). Self-consistent accretion-heated atmosphere
models have also been developed (see e.g., Zampieri et al. 1995; Deufel et al.
2001; Suleimanov et al. 2018), although these are relatively computationally
expensive and have a large number of free parameters.

In this work, we applied the formalism for Compton scattering in a hot slab
(Nagirner & Poutanen 1993; Poutanen & Svensson 1996) to compute the Stokes
parameters of the emergent radiation as a function of energy and emission
angle using only three model parameters: the electron temperature and optical
depth of the hot slab on top of the NS surface, and the temperature of the
seed blackbody photons coming from the NS. By employing pulse profile modeling
for polarized radiation from rapidly rotating oblate NSs (Poutanen 2020;
Loktev et al. 2020), we also simulated data for the Imaging X-ray Polarimeter
Explorer (IXPE, Weisskopf et al. 2022) and updated the NS geometry parameter
constraints previously predicted in Salmi et al. (2021) using the Thomson
scattering model. These methods can also be applied when analyzing the
observations from IXPE111IXPE was launched in December 2021, but has not yet
been able to observe an AMP in outburst. or from future X-ray polarimetric
missions such as the enhanced X-ray Timing and Polarimetry mission (eXTP;
Zhang et al. 2019; Watts et al. 2019). The atmosphere look-up tables produced
here can also be combined with any AMP pulse profile modeling codes.222The
tables are available at https://github.com/AnnaBobrikova/ComptonSlabTables.

The remainder of this paper is structured as follows. In Sect. 2, we present
the theory that describes the formation of polarized radiation in a
Comptonizing slab above the NS surface. We then apply this theory to obtain
the intensity and polarization of the radiation escaping from the slab. We
then describe the method to obtain pulse profiles using our new radiation
model. We compare our resulting spectra and pulse profiles to those obtained
using previous models. We also discuss pulse profiles obtained for different
NS parameters. In Sect. 3, we generate synthetic data and apply our fitting
routine to determine the NS parameters from the data. We discuss the
applications of our model and predictions for the upcoming observations in
Sect. 4. In Sect. 5, we conclude and summarize our findings.

##  2 Emission model

###  2.1 Radiative transfer equation

We consider a simple NS atmosphere model in which the atmosphere is a plane
parallel slab consisting of electrons and lying above an optically thick
source that radiates as a blackbody. The hot slab has the Thomson optical
depth $\tau_{\mathrm{T}}=\sigma_{\mathrm{T}}n_{\mathrm{e}}H$, where $H$ is the
vertical height, $\sigma_{\mathrm{T}}$ is the Thomson cross section, and
$n_{\mathrm{e}}$ is the electron concentration. The electron gas in the
atmosphere is considered to be isotropic and isothermal with the electron
temperature $T_{\mathrm{e}}$. The momentum distribution of electrons is given
by relativistic Maxwellian distribution characterized by the dimensionless
temperature $\Theta_{\mathrm{e}}=kT_{\rm e}/m_{\rm e}c^{2}$:

$f_{\mathrm{M}}(\gamma)=\frac{e^{-\gamma/\Theta_{\mathrm{e}}}}{4\pi\Theta_{\mathrm{e}}\
K_{2}(1/\Theta_{\mathrm{e}})},$ (1)

where $\gamma$ is the electron Lorentz factor and $K_{2}$ is the modified
Bessel function of the second kind. This distribution is normalized to unity:

$\int_{0}^{\infty}f_{\mathrm{M}}(\gamma){\rm
d}^{3}p=4\pi\int_{0}^{\infty}f_{\mathrm{M}}(\gamma)p^{2}{\rm d}p=1,$ (2)

where $p=\sqrt{\gamma^{2}-1}$ is the electron dimensionless momentum. The
radiation field can be described by the Stokes vector $\bm{I}(\tau,x,\mu)$,
where $\mu$ is the cosine of the zenith angle, an angle between the normal to
the slab and the direction of photon propagation, and $x=E/m_{\rm e}c^{2}$ is
the photon energy measured in the units of the electron rest mass.

The Stokes vector usually contains four components $I$, $Q$, $U$, and $V$, but
because of the azimuthal symmetry and the absence of sources of circular
polarization, the latter two are equal to zero. We therefore use just two
Stokes parameters $\bm{I}(\tau,x,\mu)=(I,Q)^{\rm T}$, where superscript T
denotes the transposed vector.

The photon distribution at the bottom of the slab is considered to be
Planckian of the temperature $T_{\mathrm{bb}}$. The incident radiation (for
zeroth scattering order, $n=0$) at the bottom of the slab is given by the
Stokes vector

$\@vec{I}_{n=0}(\tau=0,x,\mu)=\frac{2m_{\mathrm{e}}^{4}c^{6}}{h^{3}}\frac{x^{3}}{e^{x/\Theta_{\mathrm{bb}}}-1}\left(\begin{array}[]{c}1\\\
0\end{array}\right),$ (3)

where $\Theta_{\mathrm{bb}}=kT_{\mathrm{bb}}/m_{\rm e}c^{2}$. In order to
describe the propagation of polarized radiation through the hot electron slab,
we solve the radiative transfer equation (RTE) (Nagirner & Poutanen 1994;
Poutanen & Svensson 1996):

$\mu\frac{{\rm d}\@vec{I}(\tau,x,\mu)}{{\rm
d}\tau}=-\sigma(x)\@vec{I}(\tau,x,\mu)+\@vec{S}(\tau,x,\mu),$ (4)

where ${\rm d}\tau=\sigma_{\rm T}n_{\rm e}dz$ is the Thomson optical depth,
$\@vec{S}$ is the source function (also a Stokes vector), and $\sigma(x)$ is
the dimensionless Compton scattering cross-section (in units of the Thomson
cross-section $\sigma_{\rm T}$). We solve the RTE using the iterative
scattering method of Poutanen & Svensson (1996), where the intensity is
represented as a series expansion in scattering orders. The Stokes vector of
unscattered radiation at all optical depths is

$\@vec{I}_{n=0}(\tau,x,\mu)=\@vec{I}_{n=0}(0,x,\mu)\ {\rm
e}^{-\sigma(x)\tau/\mu}.$ (5)

From the known Stokes vector at the $n$-th scattering order, we find the
source function for the following scattering order as

$\@vec{S}_{n+1}(\tau,x,\mu)=x^{2}\\!\\!\int_{0}^{\infty}\\!\\!\frac{{\rm
d}x_{1}}{x_{1}^{2}}\int_{-1}^{1}\\!\\!\\!{\rm d}\mu_{1}\
\hat{\@vec{R}}(x,\mu;x_{1},\mu_{1})\@vec{I}_{n}(\tau,x_{1},\mu_{1}),$ (6)

where $\hat{\@vec{R}}$ is the $2\times 2$ azimuth averaged redistribution
matrix describing Compton scattering by isotropic electrons (see Appendix A1
in Poutanen & Svensson 1996). If the source function is known, the Stokes
vector is found via the formal solution of the RTE:

$\@vec{I}_{n+1}(\tau,x,\mu)\\!=\\!\left\\{\begin{array}[]{ll}\displaystyle\\!\\!\int_{0}^{\tau}\\!\\!\frac{{\rm
d}\tau^{\prime}}{\mu}\ \@vec{S}_{n+1}(\tau^{\prime},x,\mu)\ {\rm
e}^{-(\tau-\tau^{\prime})\sigma(x)/\mu},&\mu>0,\\\ &\\\
\displaystyle\\!\\!\int_{\tau}^{\tau_{\rm T}}\\!\\!\frac{{\rm
d}\tau^{\prime}}{(-\mu)}\ \@vec{S}_{n+1}(\tau^{\prime},x,\mu)\ {\rm
e}^{-(\tau^{\prime}-\tau)\sigma(x)/(-\mu)},&\mu<0.\end{array}\right.$ (7)

Iterations proceed until the required accuracy of the total Stokes vector
$\@vec{I}=\sum_{n}\@vec{I}_{n}$ is achieved so that the maximal contribution
of the next scattering is less than $1\%$ of the total spectrum in all
energies and angles.

Figure 1: Intensity of the emergent radiation (left) and the PD (right) as
functions of the cosine of the zenith angle. The atmosphere parameters are
$\tau_{\mathrm{T}}$ = 1.0, electron temperature $T_{\rm e}$=50 keV, and seed
photon temperature $T_{\rm bb}$=1 keV. Black, blue, green, orange, and red
solid lines show the model for photon energies of 2, 5, 8, 12, and 18 keV,
respectively. The angular dependence of the intensity is normalized so that
$\int_{0}^{1}\mu I(\mu){\rm d}\mu=1/2$.

Figure 2: Intensity of emergent radiation (left) and the PD (right) as
functions of photon energy. The same atmosphere parameters are used as in Fig.
1. Black, blue, green, orange, red, and magenta solid lines show the model for
$\mu$ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 respectively.

Finally, we can compute the polarization degree (PD) of the emergent radiation
as

$P(x,\mu)=100\%\frac{Q(\tau_{\mathrm{T}},x,\mu)}{I(\tau_{\mathrm{T}},x,\mu)}.$
(8)

The positive values of $P$ imply that the polarization vector lies in the
meridional plane defined by the slab normal and the photon momentum, while for
the negative values, the polarization vector is perpendicular to that plane,
as in the case of optically thick electron scattering atmosphere
(Chandrasekhar 1960; Sobolev 1963).

###  2.2 Spectro-polarimetric properties of accretion shocks

We can now apply the derived formalism to predict the beaming patterns, energy
spectra, and polarization of the accretion shock at the NS surface. We take as
a fiducial set of parameters: $\tau_{\mathrm{T}}=1$, $T_{\mathrm{e}}=50$ keV,
$T_{\mathrm{bb}}=1$ keV. Figures 1 and 2 show the spectrum and the PD as
functions of the cosine of the zenith angle and energy, respectively, for this
set of parameters. The parameters of the source are the same as those used in
the Thomson model of Salmi et al. (2021) (see their Figs. 1 and 2), and so
these figures can be used to illustrate the impact of using an exact Compton
scattering description instead of the Thomson scattering model in the slab of
hot electrons.

When comparing our results with those of Salmi et al. (2021), we note that the
PD differs significantly while the spectrum remains mainly unchanged. This
effect is most obviously seen in the angular dependence of the PD. The
absolute values of the PD are reduced by a factor of two for the photon
energies of 2, 5, and 8 keV, and increased at higher energies. For Compton
scattering, the PD changes the sign already at the photon energy of 8 keV and
is positive at higher values, while for Thomson scattering the sign changes at
photon energies above 18 keV. Last but not least, the maximum in the absolute
value of the PD for the photon energies higher than 8 keV has shifted from
$\mu=0$ to higher values of the cosine of the zenith angle. Keeping in mind
the energy range of the IXPE satellite (2–8 keV), we can say that compared to
the previously developed model for AMPs that used Thomson scattering, the
Compton scattering model predicts a lower PD. This result is in agreement with
our expectations, as for hot electrons the relativistic aberration becomes
important (Poutanen 1994). In the case of cold electrons, the PD of the
Thomson singly scattered radiation is highly dependent on the scattering
angle, reaching 100% at 90°. For relativistic electrons (but still in the
Thomson scattering regime), the same happens in the electron rest frame;
however, because of the relativistic motion of the electron, the photon in the
external frame is preferentially scattered in the direction of the electron
motion, resulting in zero final PD independent of the scattering angle
(Bonometto et al. 1970; Nagirner & Poutanen 1993). In our case of 50 keV
electrons, the electrons are not yet relativistic and the maximum PD is
reduced to 60% from the Thomson case. We can predict here that with the higher
electron temperature, the PD will decrease further. This will be examined in
the following section.

###  2.3 Grid of models

Figure 3: Angular dependency of the PD for different atmosphere parameters.
For each plot, we vary only one parameter (marked on the plot) while other
parameters are from the fiducial set ($\tau_{\mathrm{T}}$ = 1,
$T_{\mathrm{e}}$ = 50 keV, $T_{\mathrm{bb}}$ = 1 keV). Cyan, blue, green,
orange, and red solid lines show the PD on photon energies of 2, 5, 8, 12, and
18 keV, respectively.

The model we consider for the accretion shock has three parameters: the
optical depth of the Comptonizing plasma $\tau_{\mathrm{T}}$, the
dimensionless electron temperature $\Theta_{\rm e}$, and the dimensionless
characteristic photon temperature $\Theta_{\rm bb}$ (associated with $T_{\rm
e}$ and $T_{\rm bb}$, respectively). To determine these parameters from the
upcoming IXPE data, we need a dense grid of models. For the optical depth of
the shock $\tau_{\mathrm{T}}$, we chose the lower and higher value of 0.5 and
3.5, respectively, and varied the parameter in steps of 0.1. For the electron
temperature, we varied $\Theta_{\mathrm{e}}$ from 0.04 to 0.2 (i.e.,
$T_{\mathrm{e}}$ from $\sim$20 to 100 keV) with a step of 0.004 (i.e., $\sim$2
keV). For the seed photon temperature, the dimensionless
$\Theta_{\mathrm{bb}}$ was varied in the interval (1–3)$\times 10^{-3}$ (i.e.,
$T_{\mathrm{bb}}=$0.5–1.5 keV) in steps of $2\times 10^{-4}$ (i.e., $\sim$0.1
keV). Altogether, we computed 13981 models.

The impact of the NS atmosphere parameters on the PD can be seen in Fig. 3.
For the electron temperature, we see the expected decrease in the PD with
$T_{\rm e}$ as electrons become more relativistic. For the two other
parameters, the main feature is that with increasing the parameter value, the
change of sign of the PD happens at higher energies. Within the mentioned
intervals of parameter values, the impact is most significant with increasing
the optical depth of the atmosphere and least significant with changes to the
seed photon energy. For the IXPE observational range of 2–8 keV, this means
that, in general, PD increases with the increase in these parameters.

In order to use the calculated Stokes parameters more effectively, we wrote a
3D interpolation routine over this grid of parameters, so that for any
arbitrary combination of these parameters within the intervals mentioned, we
can obtain the intensity and PD of the emission as functions of energy and
zenith angle. We tested that the results of this interpolation match the
direct calculation results for any combination of the parameters within the
grid to within 8% (within 1% in the IXPE energy range).

###  2.4 Pulse profiles and polarization

We calculated the pulse profiles from the known beaming pattern and PD exactly
as in Salmi et al. (2021). After obtaining the beaming pattern and PD of the
emission on top of the NS atmosphere, we transport them from the hot emitting
regions of the star (called hot spots) to the observer following the formalism
described in Morsink et al. (2007), AlGendy & Morsink (2014), Poutanen (2020),
and Loktev et al. (2020). We consider an oblate, rapidly rotating star in the
Schwarzschild metric (hence, the effects of rotation on space-time are
neglected). Gravitational redshift and the light-bending effect, as well as
the effect of the fast rotation of the star, are accounted for in the
polarization angle (PA) transfer. The only difference in approach is the
reference radius used in time-delay calculation. While Loktev et al. (2020)
use the exact radius of the star at the center of the spot, we follow the path
described in Poutanen & Beloborodov (2006) and applied in Salmi et al. (2021)
and use the equatorial radius as a reference radius.

Figure 4: Pulse profiles of the observed flux, PD, and PA for two antipodal
spots shown for three different energies (2, 5, $8\leavevmode\nobreak\
\mathrm{keV}$). The solid black curves correspond to the total flux (PD, PA),
the dashed blue lines correspond to the contribution from the primary spot,
and the dotted red lines are for the secondary spot. The dash-dotted gray
lines correspond to the Thomson model from Salmi et al. (2021) (combined from
two spots). The NS parameters in both models are those shown in Table 1 of
Salmi et al. (2021).

We again return to the results from Salmi et al. (2021) and compare the pulse
profiles coming from the model with Compton scattering in the atmosphere with
previously obtained ones using Thomson scattering. Figure 4 presents this
comparison for the case of two antipodal spots, both combined and individually
from each spot, for the NS with atmospheric parameters from the fiducial set
(inclination $i$ = 60° and magnetic co-latitude $\theta$ = 20° as in Salmi et
al. 2021). We see that the flux does not change significantly, which is in
agreement with the beaming patterns not changing either. The PD is now
decreasing from 2 to 5 keV and then slowly increasing towards 8 keV. From
inspection of the corresponding lines in Fig. 1, and accounting for the
gravitational redshift, it is clear that the PD is following the same pattern
for the corresponding energies of the emitted photons. The ‘jump’ by 90° in
the PA for the 8 keV photons comes from the PD changing sign.

The main conclusion here is that we expect a significant drop in the observed
PD compared to that predicted by Salmi et al. (2021). The impact of this
reduction in PD on the possibility of observing the AMPs is further discussed
in Sects. 3 and 4.

##  3 Synthetic data and their analysis

###  3.1 Generating the synthetic data

The developed model can be further used to analyze the data. Polarimetric
measurements of the AMPs are yet to be carried out. In the absence of
observations, we used the ixpeobssim framework (Baldini et al. 2022) (version
16) with the pcube algorithm to generate several samples of synthetic data
imitating the data to be obtained from the IXPE observations of an AMP. We
substituted the source with the model described above and obtained the event
list for a specified observing time of 600 ks for a source with a flux of 100
mCrab. In this part of our study, we do not investigate the energy dependency
of PD, and collect all the data in one energy bin. These data sets are in
exactly the same form as publicly available IXPE data products, so we will be
able to apply the developed routines to the real data as soon as they become
available.

As mentioned above, applying the Compton scattering formalism reduces the
expected PD as compared to the previous studies where Thomson scattering in
the NS atmosphere was considered. We synthesized the data for the same
parameter set as used in Fig. 4 (see Table 3, Model 1) and compared the PD to
the so-called minimal detectable polarization (MDP) at 99% significance (see
Weisskopf et al. 2010 for a definition of the MDP and Kislat et al. 2015 for
further details) values for the corresponding phase bins. Figure 5a
illustrates that the obtained values of the PD are below the MDP. According to
Weisskopf et al. (2010), this means that we cannot have a proper measurement
of the PD in this case. From this point, there are three main ways to improve
the quality of the potential measurement: (i) we can increase the number of
counts to reduce the MDP values, (ii) study different energy bands to see if
the PD values are higher there, and/or (iii) choose a more promising source,
which in terms of simulated data means selecting different source parameters.

Figure 5: Simulated PD and normalized Stokes $q$ and $u$ profiles for the NS
parameter set 1 (panels (a) and (b) are for the 2–8 keV and 2–5 keV energy
range, respectively) and set 4 (c) from Table 3. The exposure time is 600 ks
in all three cases. The blue lines show the theoretical model, the purple dots
are the synthesized observed data, and the red dots are the MDP values for the
corresponding phase bins. Table 1: Parameter sets for all the computed models
and the most probable values obtained in the fitting.

Parameter | Injected values |  | Fitting results
---|---|---|---
| 1 | 2 | 3 | 4 |  | 1aaa$a$Columns labeled 1a and 3a correspond to fitting the whole energy band of IXPE, 2–8 keV. | 1bbb$b$Columns labeled 1b and 3b correspond to the data realization generated and fitted in the 2–5 keV energy band.  | 2 | 3aaa$a$Columns labeled 1a and 3a correspond to fitting the whole energy band of IXPE, 2–8 keV. | 3bbb$b$Columns labeled 1b and 3b correspond to the data realization generated and fitted in the 2–5 keV energy band.  | 4
$i$ [deg] | 60 | 10 | 60 | 10 |  | $25_{-12}^{+17}$ | $40_{-11}^{+18}$ | $10_{-3}^{+4}$ | $49_{-12}^{+15}$ | $54_{-10}^{+15}$ | $9_{-3}^{+3}$
$\theta$ [deg] | 20 | 105 | 20 | 105 |  | $21_{-12}^{+16}$ | $15_{-5}^{+7}$ | $101_{-6}^{+6}$ | $20_{-6}^{+7}$ | $17_{-3}^{+4}$ | $104_{-6}^{+4}$
$\chi$[deg] | 0 | 0 | 0 | 0 |  | $-2_{-30}^{+33}$ | $-7_{-9}^{+9}$ | $-6_{-14}^{+14}$ | $2_{-41}^{+8}$ | $-3_{-3}^{+3}$ | $-7_{-11}^{+11}$
$\Delta\phi/2\pi$ | 0.0 | 0.0 | 0.0 | 0.0 |  | $0.01_{-0.24}^{+0.22}$ | $-0.02_{-0.11}^{+0.13}$ | $-0.02_{-0.04}^{+0.04}$ | $0.03_{-0.05}^{+0.05}$ | $-0.02_{-0.03}^{+0.03}$ | $-0.02_{-0.03}^{+0.03}$
$T_{\mathrm{e}}$ [keV] | 50 | 50 | 50 | 50 |  | $58_{-24}^{+27}$ | $62_{-26}^{+26}$ | $59_{-25}^{+28}$ | $63_{-27}^{+26}$ | $65_{-27}^{+24}$ | $64_{-26}^{+24}$
$T_{\mathrm{bb}}$ [keV] | 1.0 | 1.0 | 1.0 | 1.0 |  | $1.0_{-0.3}^{+0.3}$ | $1.0_{-0.3}^{+0.3}$ | $0.9_{-0.3}^{+0.4}$ | $1.0_{-0.3}^{+0.3}$ | $1.1_{-0.3}^{+0.3}$ | $0.9_{-0.3}^{+0.4}$
$\tau_{\mathrm{T}}$ | 1.0 | 1.0 | 1.6 | 1.6 |  | $1.4_{-0.6}^{+1.1}$ | $1.9_{-0.8}^{+1.0}$ | $1.6_{-0.6}^{+0.9}$ | $2.0_{-0.8}^{+0.9}$ | $2.2_{-0.9}^{+0.8}$ | $2.1_{-0.7}^{+0.8}$

333 NS mass, radius, and spin frequency are fixed at
$M=1.4\leavevmode\nobreak\ {M}_{\sun}$, $R$ = 12 km, $\nu$ = 401 Hz. These
parameters, as well as Model 1 parameters, are selected based on the SAX
J1808.4$-$3658 pulsar and expected values for canonical NSs. The errors
represent the 68% credible interval. Additional parameters $\chi$ and
$\Delta\phi$ are the position angle of the pulsar rotation axis and the phase
shift, respectively. As before, $i$ is inclination, $\theta$ is co-latitude of
the spot, $T_{\mathrm{e}}$ and $T_{\mathrm{bb}}$ are the temperatures of the
electron in the slab and seed photons, respectively, and $\tau_{\mathrm{T}}$
is the Thomson optical depth of the slab. Values are averaged over two data
realizations for each scenario.

Increasing the number of counts by prolonging the exposure time is impossible
with AMPs, as they are transient sources observable only during the outbursts,
which typically last for a few weeks (see e.g., Patruno & Watts 2021). Also,
observations lasting longer than 600 ks are unrealistic for scheduling
reasons. The two other options are explored below.

Studying different energy bands is relevant as, according to our predictions
from Sect. 2.1, the PD changes sign within the 2–8 energy band. This means
that averaging over all energies decreases the total PD value. Figure 5b
illustrates the results of analyzing the 2–5 keV energy band instead of the
2–8 keV represented in Fig. 5a. The blue line corresponding to the model
prediction is now significantly higher and closer to the MDP values. However,
the improvement is slight, and we note that the PD values are still below the
MDP values in this case.

Figure 6: Posterior probability distributions for NS parameters for Model 1a
(left panels) and Model 4 (right panels) from Table 3 using one data
realization. The two-dimensional contours correspond to 68% and 95% credible
regions with the corresponding regions colored black and gray. The histograms
show the normalized one-dimensional distribution for a given parameter derived
from the posterior samples. The mean value and the 68% and 95% credible levels
for the parameters are shown by the vertical dashed black, red, and orange
lines, respectively. Blue lines indicate the true values of the parameters
that were used for the data simulations.

In our search for a more optimistic (from the point of view of potential
observations) yet still valid set of both geometrical and atmosphere
parameters, we aimed to reproduce the known spectra of the AMPs with the
parameters that, coming from our previous investigations, are expected to give
a higher PD of the emission. From Fig. 3, we see that the optical depth
affects the PD strongest from all the atmosphere parameters, and so we aimed
to increase it but preserve the overall spectrum behavior presented in
Poutanen & Gierliński (2003). We also attempted to increase the electron
temperature to match the predictions in Poutanen & Gierliński (2003), but the
scenario proposed there gives PD values below the MDP values, and so we chose
not to continue with it. We developed two possible adjustments to the fiducial
scenario (Model 1): we can choose the geometry such that the angle between the
inclination and co-latitude of the spot is almost 90° (as in Model 2) or use
‘more promising’ atmosphere parameters (as in Model 3). We can also combine
these two (Model 4). The parameter values for these four models are given in
the left part of Table 3 in the columns labeled correspondingly 1–4. Figure 5c
shows that for Model 4, the PD now exceeds the MDP values. We also checked
this for Models 2 and 3, and for Model 2 the PD is slightly higher than the
MDP values, while for Model 3 the PD is similar to MDP. In all the models, we
only calculated the emission coming from one spot on the surface of the NS.

In addition, we attempted to improve the quality of the results for Model 1
illustrated in Fig. 5a,b by reducing the number of phase bins, as the PD
varies less in these cases compared to other models. We find that for Model 1,
case a, only reducing the number of bins to one gives us the PD value that
barely exceeds the MDP, meaning that we can detect something from such a
source, but we would lose the phase variability completely. For Model 1, case
b, we can obtain a detection with three phase bins, but only for the bin that
includes the peak (in Fig. 5b the corresponding phase shift is 0.5), and
again, the PD barely exceeds the MDP. We conclude that reducing the number of
bins can support the observations with barely measurable PD, but its impact is
not sufficient for us to rely upon it in cases where the source is otherwise
unobservable with IXPE.

Nevertheless, the most recent development of the IXPE data analysis shows that
for multiple-point observations (i.e., observations split in different phase
bins), information about polarization can be extracted even from the data
points where PD values do not exceed MDP (González-Caniulef et al. 2023;
Suleimanov et al. 2023). As IXPE measures the Stokes parameters $q$ and $u$
rather than PD, we can fit these parameters directly and obtain some
constraints on the geometry of the source. This approach is further explored
in the following section.

###  3.2 Fitting the data

After synthesizing the data, we applied the multinest (Feroz et al. 2009)
multimodal nested sampling algorithm (specifically, the pymultinest package,
see Buchner et al. 2014) to obtain new constraints on both the geometrical and
atmosphere parameters of the source. As in Salmi et al. (2021), we fit the
geometrical parameters, such as inclination, co-latitude, phase shift, and
position angle, but unlike there, we also perform a simultaneous fit over the
parameters of the atmosphere, such as characteristic photon temperature,
electron temperature, and optical depth of the thermalization layer. Table 3
presents the most probable values and 68% credible intervals for these
parameters. Figure 6 shows the results of the fitting presented for the
fiducial set of parameters (Model 1) and for a ‘more promising’ scenario
(Model 4). In both cases, the full energy band of IXPE (2–8 keV) is
considered, and the exposure time is 600 ks.

Figure 6 clearly illustrates the improvement of the fitting results if
comparing the data synthesized from the Model 1-type of source (left) and the
data coming from the Model 4-type of source (right). As we use the uniform
prior for all our fittings, in the former case, the constraints originate from
our direct fitting of the $q$ and $u$ Stokes parameters. This illustrates the
possibility to extract some information even when the MDP values are not
reached. In the latter case, we can clearly see the significant constraints
coming from the PD values exceeding the MDP by approximately a factor of two.

We can also conclude that the constraints on the atmospheric parameters are
notably less restrictive than the ones on the geometrical parameters of the
NS. The main reason for obtaining such broad constraints on atmospheric
parameters lies in the energy range of the IXPE satellite, as we are unable to
detect the whole spectrum and fit it correctly. A solution to this issue could
be to use the observational results from several missions simultaneously; for
instance, fitting the NuSTAR data in its higher energy range to determine the
electron temperature (and consequently the optical depth through the spectral
slope). Then, using the values obtained, we can fit the IXPE data with the
fixed (or constrained) atmospheric parameters. However, we see that even with
large uncertainties on these parameters, we are able to derive constraints on
the geometry of the source from the IXPE data only.

Table 3 also shows that the constraints from the Model 3 data are less
restrictive than the ones from Model 2, and so we can conclude that the
atmospheric parameters are less important in our search for the ‘promising’
source than the geometric parameters of the NS. We also attempted to improve
the constraints by reducing the number of simultaneously fitted parameters
(e.g., we fixed the atmospheric parameters and only fitted the geometrical
parameters). Decreasing the number of fitted parameters reduces the
calculation time, but the credible intervals for the fitted parameters remain
the same.

In the fittings presented here, we fixed the mass and radius of the NS. This
is done for the simplicity of the calculations, as they are quite time-
consuming. Moreover, we only fit the normalized Stokes $q$ and $u$ parameters,
and they are not particularly sensitive to the mass or radius. We could also
fit the Stokes $I$ parameter (i.e., flux) as well and add mass and radius as
free parameters. Adding two more parameters would increase the computational
time significantly, but we do not expect any significant impact of this on the
credible intervals of other parameters.

##  4 Discussion

Our fiducial set of NS parameters is based on the SAX J1808.4$-$3658, and Fig.
4 shows that, for the parameters selected, our predictions for the PD of the
SAX J1808.4$-$3658-like source are low. However, with the uncertainties we
still have on both the geometry and the atmospheric parameters of the source,
there is still a possibility to observe a higher PD. Furthermore, with about
20 AMPs discovered, we can put forward several conclusions about the
perspective of observing measurable PD from this class of sources.

It is known that the product of the electron temperature and optical depth for
the AMPs is almost invariant (Poutanen 2006). Nevertheless, we have learned
that all the atmospheric parameters have their own contribution to the PD
values; see Fig. 3. Optical depth is the most important parameter here; we are
searching for a source with an atmosphere that is sufficiently optically thick
for many Compton scattering events to occur, yet not so thick that the
pulsation pattern is excessively influenced by these scatterings. In this
sense, objects such as IGR J17591$-$2342 and IGR J17511$-$3057 are promising.
The former has an estimated optical depth of 1.59–2.3 and a pulsed fraction of
the emission of 10%–17% (Kuiper et al. 2020; Manca et al. 2023). The latter
has a lower optical depth of 1.34, but higher electron and seed photons
temperatures ($T_{\mathrm{e}}=51$ keV, $T_{\mathrm{bb}}=1.36$ keV) and a
pulsed fraction of 14% (Papitto et al. 2010). Constraints on the geometrical
parameters of the AMPs are among the goals of the polarimetric observations.
So far our main source of information on the possible geometry is the
amplitude of the pulsation, which depends approximately on the product of the
sines of the inclination and the magnetic co-latitude (Poutanen & Beloborodov
2006). Moreover, while for the inclination we usually have some constraints,
little is known about the geometrical configuration of the magnetic field of
the AMP. From the emission patterns shown in Fig. 1, the highest PD is
received when the line of sight and the normal to the emitting area are
orthogonal. One of the most interesting objects to study is Swift
J1749.4$-$2807, for which the inclination is well-constrained,
$74.4\degree<i<77.3\degree$, and the co-latitude is estimated to be
$\theta\approx 50\degree$ (Altamirano et al. 2010). In this case, we can
expect high polarization and we can test whether the inclination obtained from
the polarimetry actually matches the orbital one.

We performed all our fittings with 600 ks of exposure time. It is typical for
AMPs to go into outburst for a couple of weeks, and so in order to collect
enough photons we need either a rather soft spectrum of emission or a very
bright source. We can exclude hard sources such as IGR J18245$-$2452 with a
photon index of $\Gamma\approx 1.4$ (Papitto et al. 2013).

Figure 4 illustrates an additional feature of the PD behavior we expect to
observe: the observed PD is different between cases where the secondary spot
is invisible (e.g., the very beginning of the outburst) and where it is
visible (at the later stages of the outburst). In order to get additional
information about the geometry of the source, it is crucial to observe the
outburst from the earliest to the latest possible stages.

We purposely focused on phase-resolved analysis and did not investigate the
possibilities that energy-resolved analysis could open up. IXPE has an energy
resolution of 0.57 keV (at 2 keV, Ramsey et al. 2021), and so with sufficient
statistics it is possible to study the behavior of PD with energy and obtain
more information about the geometry and atmospheric properties of the source.

The model presented in this paper can be improved further. For instance, the
isothermal (static) slab approximation of the accretion shock can be replaced
with a more physical model where the dynamics of the accreting gas is computed
self-consistently with the electron temperature. Also, the model for seed soft
photons can be modified to account more accurately for the reprocessing of the
hard X-rays at the NS surface, accounting for the polarization properties of
this radiation. We intend to work on these developments in the future.

##  5 Summary

We present a new model to compute the flux and polarization of the emission
coming from the accretion-powered millisecond pulsars. We used the formalism
for Compton scattering in a plane-parallel isothermal hot slab to describe the
propagation of photons in the accretion shock above the NS surface. Figures 1
and 2 show the result of such calculations for the fiducial set of parameters
taken from Salmi et al. (2021). Comparing our results with the previous
studies using Thomson scattering formalism shows that, at least for the same
NS parameters, our model predicts lower PD values.

We calculated the Stokes parameters of the radiation as a function of two
variables, the zenith angle of the emitted photon and its energy, for almost
14,000 combinations of three parameters: electron temperature, the temperature
of the seed blackbody photons, and the Thomson optical depth of the slab. The
tables of the spectral energy distribution of the Stokes parameters for all
the parameters can now be used for further analysis of the accretion-powered
millisecond pulsars, both in studying the pulse profiles and the polarimetric
properties of the emerged radiation. Figure 3 illustrates one of the ways to
use this new data set. We can understand the impact of the NS atmosphere
parameters on the angular dependency of the PD from this particular study.

We then used the relativistic rotating vector model for the oblate NS to
produce observed Stokes parameters as a function of the pulsar phase. Figure 4
shows the observed fluxes coming from the two antipodal spots of an AMP and
the combined total flux as calculated using our model and a total flux
calculated using a previously developed Thomson model. A significant drop in
the PD between these two models is again present here.

We used the developed model of the AMP emission to generate the synthetic data
sets imitating the one that the IXPE satellite could produce after observing
an AMP. We note that the NS with the parameters from our fiducial set would
emit the light with PD values below the MDP values of IXPE (which still does
not exclude the possibility to extract information about the geometry of the
source); however, for a different set of NS parameters, the PD is high enough
to be detectable by IXPE, as shown in Fig. 5. We then applied the multimodal
nested sampling technique to fit the simulated data with our model and obtain
constraints on both the geometrical and atmospheric parameters of the NS. The
results of the fitting are presented in Fig. 6. The main conclusion of this
study is that we can obtain reasonably good constraints on the geometrical
parameters from the IXPE data for some of the AMPs; however, the parameters of
the atmosphere of the NS need to be constrained from the supporting
observations by other X-ray missions such as NuSTAR.

Our findings contribute to the advancement of understanding AMPs and have
significant implications for future research in this field. We are waiting for
the AMPs to be observed by IXPE so that we may apply the developed formalism
to analyze the resulting data. We also hope to see the eXTP mission observing
the AMPs, as it would solve several of the issues discussed in this article:
the larger effective area will lower the MDP, and so the fainter sources will
become observable; the broader energy range will allow us to develop better
constraints on the atmospheric parameters; and finally, the wide field
instrument of eXTP will be able to catch the transients, such as AMPs, and
observe the outbursts from their very beginning.

###### Acknowledgements.

This research was supported by the Academy of Finland grant 333112 and by the
Finnish Cultural Foundation grant 00220175 (AB). TS acknowledges support from
ERC Consolidator Grant (CoG) No. 865768 AEONS (PI: A. Watts). The computer
resources of the Finnish IT Center for Science (CSC) and the FGCI project
(Finland) are acknowledged. We thank Anna Watts, Bas Dorsman, Alessandro Di
Marco, and John Rankin for discussions and testing our models.

## References

  * AlGendy & Morsink (2014) AlGendy, M. & Morsink, S. M. 2014, ApJ, 791, 78
  * Altamirano et al. (2010) Altamirano, D., Cavecchi, Y., Patruno, A., et al. 2010, ApJ, 727, L18
  * Baldini et al. (2022) Baldini, L., Bucciantini, N., Lalla, N. D., et al. 2022, SoftwareX, 19, 101194
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Rep. Prog. Phys., 81, 056902
  * Bogdanov et al. (2019) Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26
  * Bonometto et al. (1970) Bonometto, S., Cazzola, P., & Saggion, A. 1970, A&A, 7, 292
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Radiative transfer (New York: Dover)
  * Deufel et al. (2001) Deufel, B., Dullemond, C. P., & Spruit, H. C. 2001, A&A, 377, 955
  * Falanga et al. (2005a) Falanga, M., Bonnet-Bidaud, J. M., Poutanen, J., et al. 2005a, A&A, 436, 647
  * Falanga et al. (2005b) Falanga, M., Kuiper, L., Poutanen, J., et al. 2005b, A&A, 444, 15
  * Falanga et al. (2011) Falanga, M., Kuiper, L., Poutanen, J., et al. 2011, A&A, 529, A68
  * Falanga et al. (2012) Falanga, M., Kuiper, L., Poutanen, J., et al. 2012, A&A, 545, A26
  * Falanga et al. (2007) Falanga, M., Poutanen, J., Bonning, E. W., et al. 2007, A&A, 464, 1069
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Gierliński & Poutanen (2005) Gierliński, M. & Poutanen, J. 2005, MNRAS, 359, 1261
  * González-Caniulef et al. (2023) González-Caniulef, D., Caiazzo, I., & Heyl, J. 2023, MNRAS, 519, 5902
  * Haakonsen et al. (2012) Haakonsen, C. B., Turner, M. L., Tacik, N. A., & Rutledge, R. E. 2012, ApJ, 749, 52
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11
  * Heinke et al. (2006) Heinke, C. O., Rybicki, G. B., Narayan, R., & Grindlay, J. E. 2006, ApJ, 644, 1090
  * Ho & Heinke (2009) Ho, W. C. G. & Heinke, C. O. 2009, Nature, 462, 71
  * Ho & Lai (2001) Ho, W. C. G. & Lai, D. 2001, MNRAS, 327, 1081
  * Kislat et al. (2015) Kislat, F., Clark, B., Beilicke, M., & Krawczynski, H. 2015, Astroparticle Physics, 68, 45
  * Kuiper et al. (2020) Kuiper, L., Tsygankov, S. S., Falanga, M., et al. 2020, A&A, 641, A37
  * Lattimer (2012) Lattimer, J. M. 2012, Annu. Rev. Nucl. Part. Sci., 62, 485
  * Leahy et al. (2008) Leahy, D. A., Morsink, S. M., & Cadeau, C. 2008, ApJ, 672, 1119
  * Lindblom (1992) Lindblom, L. 1992, ApJ, 398, 569
  * Lo et al. (2013) Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19
  * Loktev et al. (2020) Loktev, V., Salmi, T., Nättilä, J., & Poutanen, J. 2020, A&A, 643, A84
  * Manca et al. (2023) Manca, A., Gambino, A. F., Sanna, A., et al. 2023, MNRAS, 519, 2309
  * Miller & Lamb (1998) Miller, M. C. & Lamb, F. K. 1998, ApJ, 499, L37
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28
  * Morsink et al. (2007) Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244
  * Morsink et al. (2007) Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244
  * Nagirner & Poutanen (1993) Nagirner, D. I. & Poutanen, J. 1993, A&A, 275, 325
  * Nagirner & Poutanen (1994) Nagirner, D. I. & Poutanen, J. 1994, Astrophysics and Space Physics Reviews, 9, 1
  * Papitto et al. (2013) Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517
  * Papitto et al. (2010) Papitto, A., Riggio, A., Di Salvo, T., et al. 2010, MNRAS, 407, 2575
  * Patruno & Watts (2021) Patruno, A. & Watts, A. L. 2021, in Astrophysics and Space Science Library, Vol. 461, Timing Neutron Stars: Pulsations, Oscillations and Explosions, ed. T. M. Belloni, M. Méndez, & C. Zhang, 143–208
  * Pechenick et al. (1983) Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846
  * Poutanen (1994) Poutanen, J. 1994, ApJS, 92, 607
  * Poutanen (2006) Poutanen, J. 2006, Adv. Space Res., 38, 2697
  * Poutanen (2020) Poutanen, J. 2020, A&A, 641, A166
  * Poutanen & Beloborodov (2006) Poutanen, J. & Beloborodov, A. M. 2006, MNRAS, 373, 836
  * Poutanen & Gierliński (2003) Poutanen, J. & Gierliński, M. 2003, MNRAS, 343, 1301
  * Poutanen & Svensson (1996) Poutanen, J. & Svensson, R. 1996, ApJ, 470, 249
  * Ramsey et al. (2021) Ramsey, B. D., Attina, P., Baldini, L., et al. 2021, in Proc. SPIE, Vol. 11821, UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XXII, ed. O. H. Siegmund, 118210M
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27
  * Salmi et al. (2021) Salmi, T., Loktev, V., Korsman, K., et al. 2021, A&A, 646, A23
  * Salmi et al. (2018) Salmi, T., Nättilä, J., & Poutanen, J. 2018, A&A, 618, A161
  * Salmi et al. (2020) Salmi, T., Suleimanov, V. F., Nättilä, J., & Poutanen, J. 2020, A&A, 641, A15
  * Salmi et al. (2022) Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2022, ApJ, 941, 150
  * Sobolev (1963) Sobolev, V. V. 1963, A treatise on radiative transfer (Princeton: Van Nostrand)
  * Steiner et al. (2009) Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279
  * Suleimanov et al. (2023) Suleimanov, V. F., Forsblom, S. V., Tsygankov, S. S., et al. 2023, A&A, in press, arXiv:2305.15309
  * Suleimanov et al. (2018) Suleimanov, V. F., Poutanen, J., & Werner, K. 2018, A&A, 619, A114
  * Sunyaev & Titarchuk (1985) Sunyaev, R. A. & Titarchuk, L. G. 1985, A&A, 143, 374
  * Viironen & Poutanen (2004) Viironen, K. & Poutanen, J. 2004, A&A, 426, 985
  * Watts et al. (2016) Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, Rev. Mod. Phys., 88, 021001
  * Watts et al. (2019) Watts, A. L., Yu, W., Poutanen, J., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 29503
  * Weisskopf et al. (2010) Weisskopf, M. C., Elsner, R. F., & O’Dell, S. L. 2010, in Proc. SPIE, Vol. 7732, Space Telescopes and Instrumentation 2010: Ultraviolet to Gamma Ray, 77320E
  * Weisskopf et al. (2022) Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 026002
  * Zampieri et al. (1995) Zampieri, L., Turolla, R., Zane, S., & Treves, A. 1995, ApJ, 439, 849
  * Zavlin et al. (1996) Zavlin, V. E., Pavlov, G. G., & Shibanov, Y. A. 1996, A&A, 315, 141
  * Zel’dovich & Shakura (1969) Zel’dovich, Y. B. & Shakura, N. I. 1969, Sov. Ast., 13, 175
  * Zhang et al. (2019) Zhang, S., Santangelo, A., Feroci, M., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 29502