# Non-convex regularization based on shrinkage penalty function

Manu Ghulyani & Muthuvel Arigovindan

###### Abstract

Total Variation regularization (TV) is a seminal approach for image recovery.
TV involves the norm of the image’s gradient, aggregated over all pixel
locations. Therefore, TV leads to piece-wise constant solutions, resulting in
what is known as the ”staircase effect.” To mitigate this effect, the Hessian
Schatten norm regularization (HSN) employs second-order derivatives,
represented by the pth norm of eigenvalues in the image hessian vector, summed
across all pixels. HSN demonstrates superior structure-preserving properties
compared to TV. However, HSN solutions tend to be overly smoothed. To address
this, we introduce a non-convex shrinkage penalty applied to the Hessian’s
eigenvalues, deviating from the convex lp norm. It is important to note that
the shrinkage penalty is not defined directly in closed form, but specified
indirectly through its proximal operation. This makes constructing a provably
convergent algorithm difficult as the singular values are also defined through
a non-linear operation. However, we were able to derive a provably convergent
algorithm using proximal operations. We prove the convergence by establishing
that the proposed regularization adheres to restricted proximal regularity.
The images recovered by this regularization were sharper than the convex
counterparts.

Department of Electrical Engg., Indian Institute of Science, Bengaluru-12,
Karnataka, India.
<EMAIL_ADDRESS>&<EMAIL_ADDRESS>
##  1 Introduction

Total Variation (TV) [24] is widely applicable because of its ability to
preserve edges. But, image reconstruction (restoration) via TV regularization
leads to piece-wise constant estimates. This effect is known as the staircase
effect. A well-known workaround for the above problem is to use higher-order
derivatives [1, 2, 3, 4] of the image rather than only the first-order
derivative. The use of higher-order derivatives leads to smooth intensity
variations on the edges rather than sharp jumps in the intensity, thereby
eliminating staircase artifacts. The above workaround led to many solutions,
important ones being TV-2, Total Generalized Variation (TGV), Hessian-Schatten
norm, and others. Hessian-Schatten (HS) norm regularization is an important
work because of it theoretical properties and good performance for a wide
variety of inverse problems [5, 6, 7]. Although HS norm regularization leads
to good reconstruction quality, it leads to smoothing of solution-images,
which is a common drawback of all convex regularizations. Also, it is well
known that non-convex regularizations [8] lead to sharper images. But,
convergence of non-convex and non-smooth optimization algorithms is difficult
to establish. Due to these difficulties in optimization of non-convex
functionals, there are very few works that explore higher-order derivative
based non-convex regularization functionals, and also have convergence
guarantees. An important work by [9] explores the properties non-convex
potential functionals, and also give an algorithm to solve the image
reconstruction problem. It is also important to mention the work by [10] that
analyses properties of edges of the recovered images via non-convex
regularization functionals. There are many works that explore non-convex
first-order total variation, for e.g. [11, 12]. The work by [12] is important
as it provides convergence as well as recovery guarantees of the proposed
reconstruction algorithm. To the best of our knowledge, there is no non-convex
regularization that exploits the structural information encoded in the
singular values of the image hessian. This is because computing these singular
values involves a non-linear operation without a known closed-form solution.
In this work, we derive a non-convex regularization inspired from the Hessian-
Schatten norm [4] and non-convex shrinkage penalty [12]. In this work, we use
the shrinkage penalty on the singular values of the hessian. Although non-
convex regularizations are designed to better approximate the $l_{0}$ norm,
non-convexity has many drawbacks such as convergence issues and no recovery
guarantees. In addition to being non-convex, the optimization problem for the
non-convex formulation of HS is also non-smooth. This further adds to the
complexity of the problem in terms of optimization. Now, the above problems
are solved by the following contributions in this work:

  1. 1.

Design of a non-convex regularization retaining the theoretical and structural
properties of the original HS-norm,

  2. 2.

algorithm to solve the image restoration problem with the proposed non-convex
functional,

  3. 3.

convergence results for the algorithm, and

  4. 4.

establshing various theoretical properties of the restoration cost with the
proposed regularization.

###  1.1 Organization of the paper

In section 2, we give describe the proposed non-convex regularization. In
section 3, we give details of the image restoration problem and show the
numerical results for the image restoration problem in section 4. Finally, all
the theoretical results and proofs are given in section 5.

##  2 Formulation

###  2.1 Forward model

The degradation model for a linear imaging inverse problem is expressed as
follows:

$\mathbf{m}=\mathcal{T}(\mathbf{u})+\mathbf{\eta},$ (1)

In this work, our approach involves considering images in a
(lexicographically) scanned form, departing from the conventional 2-D array
perspective. Therefore, the measurement image takes on a vector
representation: $\mathbf{m}\in\mathbb{C}^{N}$, while
$\mathbf{u}\in\mathbb{R}^{N}$ signifies the original image specimen.
Additionally, the operator $\mathcal{T}$ is a linear operator representing the
forward model.

This paper focuses on MRI image reconstruction. In MRI reconstruction, the
forward model $\mathcal{T}$ can be understood as the composition of two
operators: $\mathcal{T}=\mathcal{M}\circ\mathcal{F}$. The operator
$\mathcal{M}$ corresponds to the sampling trajectory and can be represented
through multiplication by a diagonal matrix, which embodies a 2D mask
consisting of 1s (where sampling occurs) and zeros (where sampling is absent).
On the other hand, $\mathcal{F}$ symbolizes the 2D Discrete Fourier Transform.
Additionally, $\mathbf{\eta}\in\mathbb{C}^{N}$ denotes the Gaussian
measurement noise.

It is important to note that when referring to a pixel in the square image
$\mathbf{u}$, containing $N$ pixels, at the coordinates $[r_{1},r_{2}]$, the
notation $[\mathbf{u}]_{\mathbf{r}}$ is used, rather than
$u_{r_{1}\sqrt{N}+r_{2}}$. This choice of notation is made to enhance clarity
and conciseness when indicating access to the pixel positioned at coordinate
$\mathbf{r}$. Also, $\sum_{\mathbf{r}}$ denotes the summation over all pixel
locations.

###  2.2 Hessian-Schatten norm regularization

The $q$-Hessian-Schatten-norm [4] ($\mathcal{HS}_{q}(\cdot)$) at any pixel
location of an image is defined as the $l_{q}$ norm of the singular values of
the image Hessian. The corresponding $q$-Hessian-Schatten-norm (HS)
regularization functional for an image is obtained as the sum of these norm
values across all pixel locations. Let $\mathbf{D}_{xx},\ \mathbf{D}_{xy},\
\mathbf{D}_{yx}\,\text{and }\mathbf{D}_{yy.}$ denote the discrete second
derivative operators (i.e, discrete analogue of second-order partial
derivative $\partial(\cdot)/\partial x\partial y$ etc.), then the discrete
Hessian , $\mathcal{H}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{2\times 2\times
N}$ can be defined as:

$[\mathcal{H}(\mathbf{u})]_{\mathbf{r}}=\begin{pmatrix}[\mathbf{D}_{xx}(\mathbf{u})]_{\mathbf{r}}&[\mathbf{D}_{xy}(\mathbf{u})]_{\mathbf{r}}\\\
[\mathbf{D}_{yx}(\mathbf{u})]_{\mathbf{r}}&[\mathbf{D}_{yy}(\mathbf{u})]_{\mathbf{r}}\end{pmatrix}\in\mathbb{R}^{2\times
2},$

for all $N$ pixel locations indexed by $\mathbf{r}$. Now, let
$\sigma_{1}([\mathcal{H}(\mathbf{u})]_{\mathbf{r}})$ and
$\sigma_{2}([\mathcal{H}(\mathbf{u})]_{\mathbf{r}})$ denote the singular
values of $[\mathcal{H}(\mathbf{u})]_{\mathbf{r}}$. With this, the Hessian-
Schatten-norm can be defined as:

$\mathcal{HS}_{q}(\mathbf{u})=\sum_{\mathbf{r}}\big{[}|\sigma_{1}([\mathcal{H}(\mathbf{u})]_{\mathbf{r}})|^{q}+|\sigma_{2}([\mathcal{H}(\mathbf{u})]_{\mathbf{r}})|^{q}\big{]}^{1/q},$

where $q$ is considered to lie in $[1,\infty].$ This is because the above
choice of $q$ makes $\mathcal{HS}_{q}(\cdot)$ convex, and therefore efficient
convex optimization algorithms (e.g. ADMM[5], primal-dual splitting [4] etc.)
can be used to obtain the reconstruction. The original work [4] also proposed
solutions using proximal operators for $q\in\\{1,2,\infty\\}$. Although the
convexity of the HS-norm described above is an advantage with respect to
optimization and convergence, it has been verified theoretically as well as by
numerical experiments that non-convex regularization functionals lead to a
better quality of the recovered image. This motivates the extension of the HS
penalty to non-convex formulation, which can lead to better reconstruction. In
this following section, we describe the non-convex HS (based) regularization
obtained by applying the shrinkage penalty [12] on the singular values of the
Hessian.

###  2.3 Shrinkage Penalty

The shrinkage penalty, as discussed in [12], possesses theoretical
properties—such as the exact recovery of sparse vectors and the convergence of
proximal algorithms—that resemble those of the conventional $l_{1}$ penalty,
despite its non-convex nature.

The shrinkage penalty ($g_{q}(\cdot)$) is not explicitly defined; rather, it
is characterized by its proximal operation. The proximal operation ($s_{q}$)
is the solution to the following cost for $t$:

$\displaystyle\gamma(x,t)=\rho g_{q}(t)+\frac{1}{2}(t-x)^{2}.$ (2)

The function $s_{q}$ is given by:

$s_{q}(x)=\operatorname*{arg\,min}_{t}\gamma(x,t)=\max\left\\{|x|-\rho^{2-q}|x|^{1-q},0\right\\}\cdot\text{sign}(x).$

This expression of $s_{q}(\cdot)$ is known as the $q$-shrinkage operation.
Notably, when $q=1$, $s_{q}(\cdot)$ corresponds to the familiar soft-
thresholding operation (which is the solution to eq. 2 with $g_{q}(t)$
replaced by $|t|$). For a given proximal mapping ($s$), the corresponding cost
($g$) will exist given the conditions outlined in 2.1 are satisfied.

######  Literature Theorem 2.1.

[12] Consider a continuous function $s:[0,\infty)\rightarrow\mathbb{R}$ and
satisfies

$s(x)=\begin{cases}0\ \ x\leq\lambda\\\ \text{strictly increasing}\ \ \
x\geq\lambda,\par\end{cases}$

also $s(x)\leq x.$ With this $s$, define $S:\mathbb{R}\rightarrow\mathbb{R}$
such that $S:x\mapsto s(x)sign(x)$, then $S(\cdot)$ is a proximal mapping of
an even, continuous and strictly increasing function $g.$ Moreover, $g(\cdot)$
is differentiable in $(0,\infty)$, and $g(\cdot)$ is non-differrentiable at
$0$ if and only if $\lambda>0$ with $\partial g(0)=[-1,1]$.

It can be observed that these conditions are satisfied by the $s_{q}$ defined
above. The function ($g_{q}(\cdot)$) derived from the shrinkage function has
some interesting properties:

  * •

$g_{q}(\cdot)$ is coercive for $q\in(0,1)$

  * •

$g_{q}^{\prime\prime}(x)<0$ for all $x\in(0,\infty)$, this means that
$g_{q}^{\prime}:(0,\infty)\rightarrow(0,\infty)$ is invertible and
$(g_{q}^{\prime})^{-1}:(0,\infty)\rightarrow(0,\infty)$ is well defined.

These properties will be used in showing the existence of the solution for the
regularized image reconstruction, and the restricted proximal regularity of
the cost.

###  2.4 QSHS: q-Shrinkage Hessian-Schatten penalty

With this $g_{q}$, we can define the shrinkage-Schatten penalty ($f(\cdot)$)
on the singular values of the image Hessian $\mathcal{H}(\mathbf{u})$ at pixel
location $\mathbf{r}$ as:

$\displaystyle
f([\mathcal{H}\mathbf{u}]_{\mathbf{r}})=g_{q}(\sigma_{1}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))+g_{q}(\sigma_{2}([\mathcal{H}\mathbf{u}]_{\mathbf{r}})).$
(3)

Without the closed form solution for $g_{q}$, we can still use $g_{q}$ for
image restoration as we can solve the following optimization problem (which is
a step in the ADMM algorithm described in section 3) in terms of the shrinkage
operation $s_{q}(\cdot).$

$H^{*}=\operatorname*{arg\,min}_{H}\frac{1}{2}\|M-H\|_{2}^{2}+\rho
f(H),M\in\mathbb{R}^{2\times 2}.$

To solve the above problem, we define $M=USV^{T}$ and $H=U_{1}S_{1}V_{1}^{T}$
to be the singular value decompositions of $M$ and $H$ respectively. Following
the approach by [4], we apply Von Neumann’s trace inequality to obtain:
$\|M-H\|_{2}\geq\|S-S_{1}\|_{2}$. Now, using the result obtained above, we
get: $\rho f(H)+\frac{1}{2}\|M-H\|_{2}^{2}\geq
0.5(\sigma_{2}(H)-\sigma_{2}(M))^{2}+0.5(\sigma_{1}(H)-\sigma_{1}(M))^{2}+\rho
g_{q}(\sigma_{1}(H))+\rho g_{q}(\sigma_{2}(H))$. Since, the problem is
separable we can obtain

$\sigma_{i}^{*}=\operatorname*{arg\,min}_{\sigma}\frac{1}{2}(\sigma-\sigma_{i}(M))^{2}+\rho
g_{q}(\sigma)=s_{q}(\sigma_{i}(M)),$

for $i=1,2$. Therefore,

$\displaystyle{H^{*}}=U\begin{pmatrix}s_{q}(\sigma_{1}(M))&0\\\
0&s_{q}(\sigma_{2}(M))\end{pmatrix}V^{T}.$ (4)

Here, the proposed ${H^{*}}$ will have a sparser set of singular values, when
compared to the original ${H}$. By this formulation it is clear that the
proposed non-convex functional will retain the properties of the original HS
formulation and lead to sharper results.

The above defined $f(\cdot)$ has many propoerties similar to HS norm:

  * •

We prove that the QSHS penalty satisfies a technical condition of the so-
called restricted proximal regularity (please refer to definition 5.1 for
details), which generalizes the concept of convexity. This condition helps us
to show the convergence of the proposed algorithm. We prove this result in
section 5.

  * •

The (continuous analogue of) QSHS penalty is translational and rotational
invariant. We state the result rigorously in the form of the following
proposition. propositionrotationinvariant Let
$u:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be a twice continuously differentiable
function and let $\mathcal{H}$ denote the hessian operator,
$f(u)\stackrel{{\scriptstyle
def}}{{=}}\int_{\mathbf{r}}g_{q}(\sigma_{1}(\mathcal{H}u(\mathbf{r})))+g_{q}(\sigma_{2}(\mathcal{H}u(\mathbf{r})))d\mathbf{r}$,
then $f(u)=f(u\circ R_{\theta})$ for any rotation matrix $R_{\theta}.$ The
proof of the above proposition is similar to the one presented in [4]. The
result given in [4] can be directly extended to our formulation as the
proposed penalty is based on the singular values of the hessian. Therefore, we
skip the proof.

##  3 Image restoration problem

The recovered image ($\mathbf{u}^{*}$) can be obtained by solving the
following optimization problem:

$\displaystyle\mathbf{u}^{*}=\operatorname*{arg\,min}_{u\in
S}\frac{1}{2}\|\mathcal{T}\mathbf{u}-\mathbf{m}\|_{2}^{2}+\rho\sum_{\mathbf{r}}f([\mathcal{H}\mathbf{u}]_{\mathbf{r}}).$
(5)

Here, $f$ is the shrinkage penalty defined in eq. 3 and $S$ is the set where
desired solution lies. For example, one widely used choice for $S$ is the
positive orthant, i.e $\\{\mathbf{u}|[\mathbf{u}]_{\mathbf{r}}\geq 0\ \
\forall\mathbf{r}\\}.$ We prove the following lemma that guarantees the
existence of the solution of the optimization problem given in eq. 5.
lemmacoercive If
$\mathcal{N}(\mathcal{T})\cap\mathcal{N}(\mathcal{H})=\\{\mathbf{0}\\}$, the
image restoration cost
$f(\mathbf{u})=\frac{1}{2}\|\mathcal{T}\mathbf{u}-\mathbf{m}\|^{2}+\rho\sum_{\mathbf{r}}g_{q}(\sigma_{1}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))+g_{q}(\sigma_{2}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))$
is coercive. Also, $f(\cdot)$ is continuous, therefore, existence of the
minimum point is guaranteed by the Weierstrass theorem. For the complete
proof, please refer to section 5.2.

We solve the optimization by ADMM approach. Although ADMM is (conventionally)
guaranteed to converge for convex functions, but there are recent works (e.g.
[13]) that demonstrate the effectiveness of ADMM for non-convex problems. In
order to derive the ADMM algorithm, we first write a constrained formulation
of eq. 5. It can be verified that eq. 5 is equivalent to the following
constrained problem ($I_{S}$ denotes the indicator function on set $S$):

$\displaystyle\mathbf{u}^{*}$ $\displaystyle=\operatorname*{arg\,min}_{u\in
S}\frac{1}{2}\|\mathcal{T}\mathbf{u}-\mathbf{m}\|_{2}^{2}+\rho\sum_{\mathbf{r}}f([\mathbf{H}]_{\mathbf{r}})+I_{S}(\mathbf{v}),$
$\displaystyle\text{ subject to
}[\mathcal{H}\mathbf{u}]_{\mathbf{r}}=[\mathbf{H}]_{\mathbf{r}}\
\forall\mathbf{r}\text{ and }\mathbf{u}=\mathbf{v}.$ (6)

Note that the constrained formulation decouples the two terms in eq. 5. The
ADMM involves minimization of the augmented Lagrangian ($\mathcal{L}()$) which
is given as:

$\displaystyle\mathcal{L}($
$\displaystyle\mathbf{u},\mathbf{H},\mathbf{v},\hat{\mathbf{u}},\hat{\mathbf{H}})=$
(7) $\displaystyle 0.5\
\|\mathcal{T}\mathbf{u}-\mathbf{m}\|_{2}^{2}+\rho\sum_{\mathbf{r}}f([\mathbf{H}]_{\mathbf{r}})+I_{S}(\mathbf{v})$
$\displaystyle+\frac{\beta}{2}\sum_{\mathbf{r}}\|[\mathcal{H}\mathbf{u}]_{\mathbf{r}}-[\mathbf{H}]_{\mathbf{r}}\|_{F}^{2}+\langle[\hat{\mathbf{H}}]_{\mathbf{r}},\|[\mathcal{H}\mathbf{u}]_{\mathbf{r}}-[\mathbf{H}]_{\mathbf{r}}\rangle+\frac{\beta}{2}\|\mathbf{u}-\mathbf{v}\|_{2}^{2}+\langle\hat{\mathbf{u}},\mathbf{u}-\mathbf{v}\rangle.$

The ADMM algorithm is composed of minimization of $\mathcal{L}()$ w.r.t
$\mathbf{u},\mathbf{H}$ and $\mathbf{v}$ cyclically, and then updating
Lagrange multipliers $\hat{\mathbf{H}}$ and $\hat{\mathbf{u}}.$

For any iteration $k\in\mathbb{N}$, the algorithm advances through the
following four steps:

Step 1, minimization w.r.t $\mathbf{v}$: $\mathbf{v}$ is updated as:
$\mathbf{v}^{(k+1)}=\operatorname*{arg\,min}_{\mathbf{v}}\mathcal{L}(\mathbf{u}^{(k)},\mathbf{H}^{(k)},\mathbf{v},\hat{\mathbf{u}}^{(k)},\hat{\mathbf{H}}^{(k)}).$
This reduces to the following minimization on completing the squares:

$\displaystyle\mathbf{v}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{v}}I_{S}(\mathbf{v})+\frac{\beta}{2}\
\|\mathbf{u}^{(k)}-\mathbf{v}+\frac{\hat{\mathbf{u}}^{(k)}}{\beta}\|_{2}^{2}$
$\displaystyle=P_{S}(\mathbf{u}^{(k)}+\frac{\hat{\mathbf{u}}^{(k)}}{\beta})$
(8)

Here, $P_{S}(\cdot)$ is the projection on set $S$.
Step 2, minimization w.r.t $\mathbf{H}$: In step 2, $\mathbf{H}$ is updated
as:

$\mathbf{H}^{(k+1)}=\operatorname*{arg\,min}_{\mathbf{H}\in\mathbb{R}^{2\times
2\times
N}}\mathcal{L}(\mathbf{u}^{(k)},\mathbf{H},\mathbf{v}^{(k+1)},\hat{\mathbf{u}}^{(k)},\hat{\mathbf{H}}^{(k)}).$

Since, the above minimization is separable for each pixel location
$\mathbf{r},$ we solve the minimization for a fixed $\mathbf{r}.$ This
minimization reduces to

$\operatorname*{arg\,min}_{[\mathbf{H}]_{\mathbf{r}}\in\mathbb{R}^{2\times
2}}\frac{\beta}{2}\|[\mathcal{H}\mathbf{u}^{(k)}]_{\mathbf{r}}-[\mathbf{H}]_{\mathbf{r}}+\frac{[\hat{\mathbf{H}}^{(k)}]_{\mathbf{r}}}{\beta}\|_{F}^{2}+\rho
f([\mathbf{H}]_{\mathbf{r}}).$

The solution to the above problem has already been done in the previous
section (eq. 4), where $M$ plays the role of
$[\mathcal{H}\mathbf{u}^{(k)}]_{\mathbf{r}}+[\hat{\mathbf{H}}^{(k)}]_{\mathbf{r}}.$
Step 3, minimizing w.r.t $\mathbf{u}$: Updating $\mathbf{u}$ is essentially
minimizing $\frac{1}{2}\
\|\mathcal{T}\mathbf{u}-\mathbf{m}\|_{2}^{2}+\frac{\beta}{2}\sum_{\mathbf{r}}\|[\mathcal{H}\mathbf{u}]_{\mathbf{r}}-[\mathbf{H}^{(k+1)}]_{\mathbf{r}}+\frac{[\hat{\mathbf{H}}^{(k)}]_{\mathbf{r}}}{\beta}\|_{F}^{2}+\frac{\beta}{2}\|\mathbf{u}-\mathbf{v}^{(k+1)}+\frac{\hat{\mathbf{u}}^{(k)}}{\beta}\|_{2}^{2}$
w.r.t $\mathbf{u}.$ The minimizer of the above cost can be written as:

$\displaystyle\big{[}\mathcal{T}^{*}\mathcal{T}+\beta$
$\displaystyle(\mathbf{D}_{xx}^{T}\mathbf{D}_{xx}+\mathbf{D}_{xy}^{T}\mathbf{D}_{xy}+\mathbf{D}_{yx}^{T}\mathbf{D}_{yx}+\mathbf{D}_{yy}^{T}\mathbf{D}_{yy}+\mathcal{I})\big{]}\mathbf{u}^{(k+1)}=$
(9)
$\displaystyle\mathcal{T}^{*}\mathbf{m}+\beta\mathbf{v}^{(k+1)}-\mathbf{u}^{(k)}+\beta(\mathbf{D}_{xx}^{T}\bar{\mathbf{H}}_{11}+\mathbf{D}_{xy}^{T}\bar{\mathbf{H}}_{12}+\mathbf{D}_{yx}^{T}\bar{\mathbf{H}}_{21}+\mathbf{D}_{yy}^{T}\bar{\mathbf{H}}_{22}).$
(10)

Here, $\left[\begin{array}[]{cc}\bar{\mathbf{H}}_{11}&\bar{\mathbf{H}}_{12}\\\
\bar{\mathbf{H}}_{21}&\bar{\mathbf{H}}_{22}\end{array}\right]=\mathbf{H}^{(k+1)}-\hat{\mathbf{H}}^{(k)}/\beta.$The
above problem is linear, and as all the operators involved are block
circulant, the equation can be solved efficiently by using 2D-FFTs.

Step 4, updating multipliers: After the above steps multipliers can be updated
as:

  * •

$\hat{\mathbf{u}}^{(k+1)}=\hat{\mathbf{u}}^{(k)}+\beta(\mathbf{u}^{(k+1)}-\mathbf{v}^{(k+1)})$,
and

  * •

$\forall\mathbf{r}:[\hat{\mathbf{H}}^{(k+1)}]_{\mathbf{r}}=[\hat{\mathbf{H}}^{(k)}]_{\mathbf{r}}+\beta[\mathcal{H}\mathbf{u}^{(k+1)}]_{\mathbf{r}}-\beta[\mathbf{H}^{(k+1)}]_{\mathbf{r}}.$

###  3.1 Convergence Guarantees

ADMM was proposed in [14, 15]. ADMM typically converges for convex problems
[16], but can fail to converge for multi-block (3 or more) splitting. The
behaviour of ADMM for non-convex and non-smooth problems was largely unknown
and many questions are still unanswered. But, owing to successful results of
the algorithm in many applications (especially in signal processing
literature, see for e.g. [17, 18, 19]) there has been a lot of interest in
understanding the convergence of the ADMM for non-convex and non-smooth
problems. There are many frameworks that establish the convergence of non-
convex ADMM [20, 21, 13, 22, 23]. The works by [22] and [23] need restrictive
assumptions on the iterates, which are difficult to verify. The [20] work
requires that the hessian of the smooth part (data-fitting) of the cost be
lower-bounded, this is not true as $\mathcal{T}$ has a non-trivial null space
for most of the imaging inverse problems. [21] prove the convergence for only
a special class of optimization problems, and the framework is not general.
[23] provide the most general framework and allow us to prove that the
algorithm is (subsequentially) convergent. We prove the following theorem that
guarantees that any sub-sequential limit of the sequence generated by the
above algorithm is a stationary point of the image restoration cost.
theoremconvergenceIf
$\mathcal{N}(\mathcal{T})\cap\mathcal{N}(\mathcal{H})=\\{\mathbf{0}\\}$ and
$\beta$ is sufficiently large the iterates generated algorithm defined in
section 3 by steps 1-4 are bounded. Moreover, each limit point of the sequence
generated by the iterate is a stationary point of the image restoration cost
$f(\mathbf{u})=\frac{1}{2}\|\mathcal{T}\mathbf{u}-\mathbf{m}\|^{2}+\rho\sum_{\mathbf{r}}g_{q}(\sigma_{1}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))+g_{q}(\sigma_{2}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))$
(defined in eq. 5). The theorem provides the following assurance: when the
sequence produced by the algorithm converges, its convergence will occur at
the point where the sub-gradient of the restoration cost reaches zero. Given
that the sequence is bounded, the existence of a converging subsequence is
guaranteed. Consequently, the limit of this subsequence will correspond to the
point where the sub-gradient of the cost becomes zero. For the proof of the
above theorem, please see section 5.2.

##  4 Simulation Results

To demonstrate the effectiveness of the proposed method, we compare the
reconstruction results with q-Hessian Schatten norm [4] (for $q$=1 and 2) and
TV-1 [24]. Hessian Schatten norm for $q=2$ is popularly known as TV-2. We use
two sampling masks ($\mathcal{M}$) with sampling densities 18 and 9 percent.
We add noise with $\sigma=2.5.$ In numerical simulations, we use a data-set
with 5 typical MRI-images (see fig. 1) of size $256\times 256.$ For the
proposed shrinkage penalty we use $q=.5$ (in eq. 2). The optimal
regularization parameter was tuned (by golden-section method) to obtain
minimum Mean Squared Error (MSE). The SSIM scores of the reconstructions are
given in the table 1. The table clearly shows that the proposed method
performs better than all other methods by a significant margin. To demonstrate
the visual difference between the images, we show result of image 2 for mask-2
(shown in fig. 3 and zoomed view in fig. 4). Clearly, the proposed method
recovers sharper images as dot like structures are much sharper in the
proposed method. Also, from the algorithm it is clear that there is no
significant computational cost associated with the $q-$ shrinkage step.

Im | Mask | TV-1 | TV-2 | HS-1 | Proposed
---|---|---|---|---|---
1 | 1 | 0.915 | 0.920 | 0.925 | 0.958
2 | 0.790 | 0.804 | 0.816 | 0.884
2 | 1 | 0.960 | 0.964 | 0.966 | 0.980
2 | 0.893 | 0.903 | 0.905 | 0.953
3 | 1 | 0.937 | 0.937 | 0.940 | 0.957
2 | 0.815 | 0.815 | 0.820 | 0.849
4 | 1 | 0.936 | 0.938 | 0.941 | 0.966
2 | 0.866 | 0.871 | 0.875 | 0.920
5 | 1 | 0.924 | 0.930 | 0.932 | 0.949
2 | 0.816 | 0.817 | 0.825 | 0.859
Table 1: Table showing SSIM values of reconstructions Figure 1: Test Images
Figure 2: Masks for sampling trajectories Figure 3: Result for Im2 and Mask 2
Figure 4: Zoomed view of the result

##  5 Theoretical Results and Proofs

###  5.1 Properties of QSHS penalty

Now, we formally define the concept of restricted proximal regularity.

######  Definition 5.1.

(Restricted proximal regularity) A lower semi-continuous function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{\infty\\}$ is restricted proximal
regular if for any $M>0$ and any bounded set $\Omega$ there exists
$\gamma\equiv\gamma(M,\Omega)$ such that the following holds for all
$y\in\Omega,x\in\\{\,x\in\Omega\mid\|p\|\leq M\ \ \ \forall p\ \ \ \in\partial
f(x)\,\\}$, and for all $d\in\partial f(x)$:

$f(y)-f(x)-\langle d,(y-x)\rangle\geq-\frac{\gamma}{2}\|y-x\|^{2}.$

The following proposition establishes the restricted proximal regularity of
the proposed QSHS penalty. The proof proceeds by following the methodology
outlined in the proof of restricted proximal regularity for $l_{q}$ norms
($q\in(0,1)$) as presented in [13]. However, a notable challenge in this
context is the lack of a closed-form expression for the penalty. Consequently,
we use the abstract properties of $g_{q}$ to prove the result. propositionrpr
Consider any $\mathbf{H}\in\mathbb{R}^{2\times 2},$ then
$r(\mathbf{H})\stackrel{{\scriptstyle
def}}{{=}}g_{q}(\sigma_{1}(\mathbf{H}))+g_{q}(\sigma_{2}(\mathbf{H}))$ is
restricted proximal regular.

###### Proof.

We use the following result [25] for the sub-gradient: Let
$H\in\mathbb{R}^{2\times 2}$ and $H=[U\ U_{1}]\begin{bmatrix}S&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\\\ \end{bmatrix}[V\ V_{1}]^{T}$ be the singular value
decomposition (SVD) of $H$, $r(H)=\sum_{i=1}^{2}g_{q}(S_{ii}(\mathbf{H})),$
then $UDV^{T}+U_{1}\Theta V_{1}^{T}\in\partial r(H),$ where $D$ is diagonal
matrix with entries $(D)_{ii}=g_{q}^{\prime}(S_{ii}),\text{ and }\Theta$ is
any arbitrary matrix.

Without loss of generality we choose, $\Omega=\\{\,X\in\mathbb{R}^{2\times
2}\mid\|X\|\leq P\,\\}$. Now, for any $M,P>0$, we intend to show that

$r(B)-r(A)-\langle T,B-A\rangle\geq-\frac{\gamma}{2}\|B-A\|^{2}$

for all $B\in\Omega,A\in\Omega_{M}\stackrel{{\scriptstyle
def}}{{=}}\\{\,X\in\Omega\mid\|T\|\leq M,\ \ \ \forall T\ \ \ \in\partial
r(X)\,\\}$, and for all $T\in\partial r(A)$ . We do this in the following
cases:

Case 1:

$\|B-A\|\geq\epsilon_{0}=\frac{1}{3}(g_{q}^{\prime})^{-1}(M)$ Note that the
above condition is equivalent to

$\displaystyle\frac{-\|B-A\|}{\epsilon_{0}}\leq-1.$ (11)

First, it can be observed that,

$\displaystyle r(B)-r(A)-\langle T,B-A$
$\displaystyle\rangle\stackrel{{\scriptstyle a}}{{\geq}}-r(A)+\|T\|\|B-A\|,$
(12) $\displaystyle\stackrel{{\scriptstyle b}}{{\geq}}-R_{max}-M\|B-A\|.$ (13)

Here, $(a)$ is true as $r(\cdot)$ is non-negative and by Cauchy-Scwartz
inequality; while $(b)$ is true as $T$ is bounded and $r(\cdot)$ is a
continuous function on a bounded and closed set, therefore, it attains maximum
$R_{max}.$ Now, using eq. 11 we obtain

$\displaystyle r(B)-r(A)-\langle
T,B-A\rangle{\geq}-\big{(}\frac{R_{max}+M\epsilon_{0}}{\epsilon_{0}^{2}}\big{)}\|B-A\|^{2}.$
(14)

Case 2:

$\|B-A\|<\epsilon_{0}$. To prove this, we first define
$\Omega^{\prime}=\\{\,T\in\mathbb{R}^{n\times n}\mid\|T\|\leq P,\
\min_{i}\sigma_{i}(T)\geq\epsilon_{0}\,\\}$.

Now, we decompose $B=U_{B}\Sigma_{B}V_{B}^{T}$,
$B=U_{B}\Sigma_{B}^{(1)}V_{B}^{T}+U_{B}\Sigma_{B}^{(2)}V_{B}^{T}=B_{1}+B_{2},$
where all singular values of $B_{1}$ are greater than $\epsilon_{0}.$ Clearly,
$B_{1}\in\Omega^{\prime}.$ We show that $A\in\Omega^{\prime}.$ This is proved
by contradiction. Assume the contrary that $A\notin\Omega^{\prime}.$ This
means, $\exists\ i$ such that
$\sigma_{i}(A)<\epsilon_{0}\implies\sigma_{i}(A)<3\epsilon_{0}.$ Now, since
$g_{q}^{\prime}(\cdot)$ is non-increasing, we have
$g_{q}^{\prime}(\sigma_{i}(A))>g_{q}^{\prime}(3\epsilon_{0})=g_{q}^{\prime}((g_{q}^{\prime})^{-1}(M))=M.$
Let $A=U_{A}S_{A}V_{A}^{T}$ be the singular value decomposition of A. Define
$T_{1}\stackrel{{\scriptstyle def}}{{=}}U_{A}S^{\prime}V_{A}^{T},$ where
$\\{S^{\prime}\\}_{kk}=g_{q}^{\prime}(\\{S_{A}\\}_{kk})$ for all $k$. Now, by
lemma $T_{1}\in\partial r(A).$ But, $\|T_{1}\|\geq\|S^{\prime}\|\geq
g_{q}^{\prime}(\sigma_{i}(A))>M.$ This contradicts the fact that
$A\in\Omega_{M}.$ Hence, $A\in\Omega^{\prime}.$ Now, we define a function,
$F:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n\times n}$, which is defined
as

$F:X\mapsto U_{X}D_{X}^{\prime}V_{X}^{T}.$

Here, $X=U_{X}D_{X}V_{X}^{T}$ is the singular value decomposition of $X$ and
$D_{X}^{\prime}$ is a diagonal matrix which is defined as
$(D_{X}^{\prime})_{ii}=g_{q}^{\prime}((D_{X})_{ii})$. Since $F$ is continuous
on compact set on $\Omega^{\prime}$, it is Lipschitz continuous on
$\Omega^{\prime}$ [26], this means $\|F(B)-F(A)\|\leq L\|B-A\|$. Now, by
Taylor’s expansion we have:

$\displaystyle r(B_{1})-r(A)-\langle
B_{1}-A,U_{A}S^{\prime}V_{A}^{T}\rangle\geq\frac{-L}{2}\|B_{1}-A\|^{2}$ (15)

. Now, $\|U_{2}^{T}U_{B}\|\leq\frac{\|A-B_{1}\|}{\epsilon_{0}}$ and
$\|V_{2}^{T}V_{B}\|\leq\frac{\|A-B_{1}\|}{\epsilon_{0}}$ [27]. Now,

$\displaystyle\langle U_{2}^{T}\Theta
V_{2},B_{1}-A\rangle=\langle\Theta,U_{2}U_{B}^{T}\Sigma_{B}^{(1)}V_{B}V_{2}^{T}\rangle\geq-\frac{M^{2}\|B_{1}-A|\|^{2}}{\epsilon_{0}^{2}}.$
(16)

Also, $r(B_{2})-<T,B_{2}>\geq 0$ and by triangle inequality we
get:$\|B_{1}-A\|\leq\|B_{1}-B\|+\|B-A\|\leq 2\|B-A\|.$ Adding eq. 15 and eq.
16 we get,

$\displaystyle r(B)-r(A)-\langle
B-A,T\rangle\geq-(\frac{L}{2}+\frac{4M^{2}}{\epsilon_{0}^{2}})\|B-A\|^{2}.$
(17)

∎

###  5.2 Properties of Restoration cost

The following helps us to establish that the restoration cost is coercive.

######  Claim 0.1.

If $\mathcal{N}(\mathcal{T})\cap\mathcal{N}(\mathcal{H})=\\{\mathbf{0}\\}$ the
function
$\|\mathcal{T}\mathbf{u}\|+\|\mathcal{H}\mathbf{u}\|_{*}\geq\gamma\|\mathbf{u}\|,\text{where
}\gamma>0.$ Here, $\|\mathcal{H}\mathbf{u}\|_{*}$ is
$\mathcal{H}S_{1}(\mathbf{u}),$ the conventional $l-1$ Hessian Schatten norm.

###### Proof.

The above statement is trivial if $\mathbf{u}=\mathbf{0}.$ If
$\mathbf{u}\neq\mathbf{0}$, let
$\hat{\mathbf{u}}=\frac{\mathbf{u}}{\|\mathbf{u}\|}$, then
$\|\mathcal{T}\hat{\mathbf{u}}\|+\|\mathcal{H}\hat{\mathbf{u}}\|_{*}\geq\inf_{\|p\|=1}\|\mathcal{T}{\mathbf{p}}\|+\|\mathcal{H}{\mathbf{p}}\|_{*}$.
Since, $\|\mathbf{p}\|=1$ is a compact set, $\exists\mathbf{p}_{min}$ (with
$\|\mathbf{p}_{min}\|=1$) such that
$\inf_{\|p\|=1}\|\mathcal{T}{\mathbf{p}}\|+\|\mathcal{H}{\mathbf{p}}\|_{*}=\|\mathcal{T}{\mathbf{p}_{min}}\|+\|\mathcal{H}{\mathbf{p}_{min}}\|_{*}$.
Define,
$\gamma=\|\mathcal{T}{\mathbf{p}_{min}}\|+\|\mathcal{H}{\mathbf{p}_{min}}\|_{*}.$
Clearly, $\gamma\neq 0$ as we will get a vector in intersection of the null
spaces, i.e.
$\mathbf{p}_{min}\in\mathcal{N}(\mathcal{T})\cap\mathcal{N}(\mathcal{H})$,
this contradicts the hypothesis. Re-substituting
$\hat{\mathbf{u}}=\frac{\mathbf{u}}{\|\mathbf{u}\|}$ completes the proof. ∎

*

###### Proof.

Without loss of generality, we prove the theorem for $\rho=1.$ Consider the
level set $\mathcal{L}_{\eta}(f)\stackrel{{\scriptstyle
def}}{{=}}\\{\,\mathbf{x}\mid f(\mathbf{x})\leq\eta\,\\}.$ Now, if
$f(\mathbf{x})\leq\eta\implies\frac{1}{2}\|\mathcal{T}\mathbf{x}-\mathbf{m}\|^{2}\leq\eta.$
Now, by triangle inequality

$\displaystyle\|\mathcal{T}\mathbf{x}\|\leq\sqrt{2\eta}+\|\mathbf{m}\|.$ (18)

Since, $\mathbf{u}\in\mathcal{L}_{\eta}(f)\implies
g_{q}(\sigma_{i}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))\leq\eta\ \ \
\forall\mathbf{r}\text{ and }i=1,2.$ As $g_{q}(\cdot)$ is coercive for
$q\in(0,1)$, we have $\sigma_{i}([\mathcal{H}\mathbf{u}]_{\mathbf{r}})\leq M$
for some $M>0.$ This above statement is true because of the fact that any
level set of a coercive function is compact. By Taylor’s series of $g_{q}$
around $0$ we can see that:
$g_{q}(\sigma_{i}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))=g_{q}^{\prime}(\gamma_{i,\mathbf{r}})\sigma_{i}([\mathcal{H}\mathbf{u}]_{\mathbf{r}})$
for $\gamma_{i,\mathbf{r}}\in[0,2M)$. Since, $g_{q}^{\prime}(\cdot)$ is
decreasing, therefore, $g_{q}^{\prime}(\gamma_{i,\mathbf{r}})\geq
g_{q}^{\prime}(2M)\stackrel{{\scriptstyle def}}{{=}}C_{2M}.$ Hence,
$\mathbf{x}\in\mathcal{L}_{\eta}(f)\implies
C_{2M}\|\mathcal{H}\mathbf{u}\|_{*}\leq\sum_{\mathbf{r}}g_{q}(\sigma_{1}([\mathcal{H}\mathbf{u}]_{\mathbf{r}}))+g_{q}(\sigma_{2}([\mathcal{H}\mathbf{u}]_{\mathbf{r}})\leq\eta.$
Now, we use the following 0.1 to show that the level set
$\mathcal{L}_{\eta}(f)$ is bounded. Using, 0.1 we get

$\displaystyle\gamma\|\mathbf{u}\|\leq\|\mathcal{T}\mathbf{u}\|+\|\mathcal{H}\mathbf{u}\|_{*}\leq\frac{\eta}{C_{2M}}+\sqrt{2\eta}+\|\mathbf{m}\|.$
(19)

This means the level set $\mathcal{L}_{\eta}(f)$ is bounded. Combining with
the fact the level set is closed as $f(\cdot)$ is continuous implies that the
level set is compact. As this is true for any level set, therefore $f$ is
coercive.∎

###  5.3 Proof of convergence

$N_{iter}\leftarrow 1000$;

$\mathbf{u}^{(0)}=\mathbf{0}$;

$\mathbf{v}^{(0)}=\mathbf{0}$;

$\mathbf{w}^{(0)}=\mathbf{0}$;

if  _$i\leq N_{iter}$_ then

$\mathbf{u}^{(i+1)}\leftarrow\operatorname*{arg\,min}_{\mathbf{u}}\mathcal{L}_{\beta}(\mathbf{u},\mathbf{v}^{(i)},\mathbf{w}^{(i)})$;

$\mathbf{v}^{(i+1)}\leftarrow\operatorname*{arg\,min}_{\mathbf{v}}\mathcal{L}_{\beta}(\mathbf{u}^{(i+1)},\mathbf{v},\mathbf{w}^{(i)})$;

$\mathbf{w}^{(i+1)}\leftarrow\mathbf{w}^{(i)}+\beta(\mathbf{A}\mathbf{u}^{(i+1)}+\mathbf{B}\mathbf{v}^{(i+1)})$;

$i\leftarrow i+1$;

Algorithm 1 Non-convex ADMM

We will use the following theorem by [13] to show the convergencve of our
algorithm.

######  Literature Theorem 5.1.

([13], theorem 2.2) Consider the minimization of the function
$\phi(\mathbf{u},\mathbf{v})=h(\mathbf{v})+g(\mathbf{v})$ subject to
$\mathbf{A}\mathbf{u}+\mathbf{B}\mathbf{v}=\mathbf{0}$ by non-convex ADMM
algorithm (Algorithm 1). Define the augmented Lagrangian,
$\mathcal{L}_{\beta}(\mathbf{u},\mathbf{v},\mathbf{w})\stackrel{{\scriptstyle
def}}{{=}}\phi(\mathbf{u},\mathbf{v})+\mathbf{w}^{T}(\mathbf{A}\mathbf{u}+\mathbf{B}\mathbf{v})+\frac{\beta}{2}\|\mathbf{A}\mathbf{u}+\mathbf{B}\mathbf{v}\|^{2}.$
If

C1

: $\phi(\mathbf{u},\mathbf{v})$ is coercive on the set
$\\{\,(\mathbf{u},\mathbf{v})\mid\mathbf{A}\mathbf{u}+\mathbf{B}\mathbf{v}=\mathbf{0}\,\\};$

C2

: $Im(A)\subset Im(B)$, where $Im$ denotes the image of the linear operator;

C3

:$\mathbf{A}$ and $\mathbf{B}$ are full column rank;

C4

: $g$ is restricted proximal regular, and

C5

:$h$ is Lipschitz smooth,

then algorithm generates a bounded sequence that has atleast one limit point,
and each limit point is a stationary point of $\mathcal{L}_{\beta}(\cdot)$

Now, we show the convergence of algorithm using the above theorem.

*

###### Proof.

We establish the validity of the aforementioned theorem by fulfilling the
conditions ((C1)-(C5)) of 5.1. In comparison to the splitting presented in
section 3, we set
$h\equiv\frac{1}{2}\|\mathcal{T}(\cdot)-\mathbf{m}\|_{2}^{2}$ and
$g\equiv\rho\sum_{\mathbf{r}}f(\cdot)+I_{S}(\cdot)$. The operator $\mathbf{A}$
is analogous to a linear operator, satisfying
$\mathcal{A}(\mathbf{u})-\mathbf{v}=\mathbf{0}$. For each pixel location
$\mathbf{r}$,
$\mathcal{A}:[\mathbf{u}]_{\mathbf{r}}\rightarrow([\mathcal{H}\mathbf{u}]_{\mathbf{r}},[\mathbf{u}]_{\mathbf{r}})$.
Regarding the constraints, it is evident that $\mathbf{B}$ corresponds to the
negative identity matrix.

(C1) follows from section 3, and (C4) follows from definition 5.1. (C2) is
trivial since $\mathbf{B}$ is the negative identity. (C5) holds true due to
the quadratic nature of the data-fitting cost. To ascertain (C3), we must
demonstrate that $\mathcal{N}(\mathbf{B})=\mathbf{0}$ and
$\mathcal{N}(\mathbf{A})=\mathbf{0}$. Since $\mathbf{B}$ is the negative
identity matrix, $\mathcal{N}(\mathbf{B})=\mathbf{0}$. For $\mathbf{A}$, let
$\mathcal{A}(\mathbf{z})=\mathbf{0}$, implying
$[\mathcal{H}\mathbf{z}]_{\mathbf{r}}=\mathbf{0}$ and
$[\mathbf{z}]_{\mathbf{r}}=\mathbf{0}$ for all $\mathbf{r}$, which concludes
$\mathbf{z}=\mathbf{0}$. Hence, all the conditions are met. The next step is
to establish that the stationary point of the augmented Lagrangian coincides
with that of the restoration cost. Suppose
$(\mathbf{u}^{*},\mathbf{v}^{*},\mathbf{w}^{*})$ is the stationary point of
the augmented Lagrangian; this implies:

  1. 1.

$\mathcal{A}\mathbf{u}^{*}-\mathbf{v}^{*}=\mathbf{0}$,

  2. 2.

$\nabla h(\mathbf{u}^{*})+\mathcal{A}^{T}\mathbf{w}^{*}=0,$ and

  3. 3.

$\partial g(\mathbf{v}^{*})-\mathbf{w}^{*}\ni\mathbf{0}.$

Rearranging 2 and 3 we obtain $\nabla
h(\mathbf{u}^{*})\in-\mathcal{A}^{T}\partial g(\mathbf{v}^{*})$. Now, we use 1
to get

$\displaystyle-\nabla
h(\mathbf{u}^{*})-\mathcal{A}^{T}\partial(\mathcal{A}\mathbf{u}^{*})\ni\mathbf{0}$
(20) $\displaystyle\implies-\partial f(\mathbf{u}^{*})\ni\mathbf{0}$ (21)
$\displaystyle\implies\partial f(\mathbf{u}^{*})\ni\mathbf{0}.$ (22)

∎

##  6 References

## References

  * [1] Scherzer O 1998 Computing 60 1–27 ISSN 0010-485X
  * [2] Bredies K, Kunisch K and Pock T 2010 SIAM Journal on Imaging Sciences 3 492–526
  * [3] Lysaker M, Lundervold A and Tai X C 2003 IEEE Transactions on image processing 12 1579–1590
  * [4] Lefkimmiatis S, Ward J P and Unser M 2013 IEEE transactions on image processing 22 1873–1888
  * [5] Lefkimmiatis S and Unser M 2013 IEEE transactions on image processing 22 4314–4327
  * [6] Liu L, Li X, Xiang K, Wang J and Tan S 2017 IEEE transactions on medical imaging 36 2588–2599
  * [7] Wang Q, Qu G and Ji D 2021 Mathematical Methods in the Applied Sciences 44 1674–1687
  * [8] Arigovindan M, Fung J C, Elnatan D, Mennella V, Chan Y H M, Pollard M, Branlund E, Sedat J W and Agard D A 2013 Proceedings of the National Academy of Sciences 110 17344–17349
  * [9] Nikolova M, Ng M K and Tam C P 2010 IEEE Transactions on Image Processing 19 3073–3088
  * [10] Nikolova M 2005 Multiscale Modeling & Simulation 4 960–991
  * [11] Sidky E Y, Chartrand R and Pan X 2007 Image reconstruction from few views by non-convex optimization 2007 IEEE Nuclear Science Symposium Conference Record vol 5 pp 3526–3530
  * [12] Woodworth J and Chartrand R 2016 Inverse Problems 32 075004
  * [13] Wang Y, Yin W and Zeng J 2019 Journal of Scientific Computing 78 29–63
  * [14] Glowinski R and Marroco A 1975 Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique 9 41–76
  * [15] Gabay D and Mercier B 1976 Computers & mathematics with applications 2 17–40
  * [16] Davis D and Yin W 2016 Convergence rate analysis of several splitting schemes Splitting methods in communication, imaging, science, and engineering (Springer) pp 115–163
  * [17] Chartrand R and Wohlberg B 2013 A nonconvex admm algorithm for group sparsity with sparse groups 2013 IEEE international conference on acoustics, speech and signal processing (IEEE) pp 6009–6013
  * [18] Liavas A P and Sidiropoulos N D 2015 IEEE Transactions on Signal Processing 63 5450–5463
  * [19] Wen Z, Yang C, Liu X and Marchesini S 2012 Inverse Problems 28 115010
  * [20] Li G and Pong T K 2015 SIAM Journal on Optimization 25 2434–2460
  * [21] Hong M, Luo Z Q and Razaviyayn M 2016 SIAM Journal on Optimization 26 337–364
  * [22] Magnússon S, Weeraddana P C, Rabbat M G and Fischione C 2015 IEEE Transactions on Control of Network Systems 3 296–309
  * [23] Wang F, Cao W and Xu Z 2018 Science China Information Sciences 61 1–12
  * [24] Rudin L I, Osher S and Fatemi E 1992 Physica D: nonlinear phenomena 60 259–268
  * [25] Watson G A 1992 Linear Algebra Appl 170 33–45
  * [26] Ding C, Sun D, Sun J and Toh K C 2018 Mathematical Programming 168 509–531
  * [27] Li R C and Stewart G 2000 Linear Algebra and its Applications 313 41–51


# Quasi-periodic dipping in the ultraluminous X-ray source,
NGC 247 ULX-1

W. N. Alston1, C. Pinto2, D. Barret3, A. D’Aì2, M. Del Santo2, H. Earnshaw4
E-mail<EMAIL_ADDRESS>A. C. Fabian5, F. Fuerst1, E. Kara6, P. Kosec6,
M. J. Middleton7, M. L. Parker5 F. Pintore8, A. Robba2,9, T. P. Roberts10, R.
Sathyaprakash11, D. Walton5, E. Ambrosi2
1European Space Agency (ESA), European Space Astronomy Centre (ESAC),
Villanueva de la Cañada, Madrid, E-28691, Spain
2INAF/IASF Palermo, via Ugo La Malfa 153, I-90146 Palermo, Italy
3Institut de Recherche en Astrophysique et Planétologie (IRAP), 9 Avenue du
Colonel Roche, BP 44346, 31028, Toulouse Cedex 4, France
4Cahill Center for Astronomy and Astrophysics, California Institute of
Technology, Pasadena, CA 91125, USA
5Institute of Astronomy, Madingley Rd, Cambridge, CB3 0HA.
6MIT Kavli Institute for Astrophysics and Space Research, Cambridge, MA 02139,
USA
7Department of Physics and Astronomy, University of Southampton, Highfield,
Southampton SO17 1BJ, UK
8INAF-IASF Milano, via A. Corti 12, I-20133 Milano, Italy
9Università degli Studi di Palermo, Dipartimento di Fisica e Chimica, via
Archirafi 36, I-90123 Palermo, Italy
10Centre for Extragalactic Astronomy, Department of Physics, Durham
University, South Road, Durham, DH1 3LE
11Institut de Ciències de l‘Espai, Carrer de Can Magrans, 08193, Cerdanyola
del Vallès, Barcelona

(Accepted 2021 May 18. Received 2021 May 18; in original form 2021 April 20.)

###### Abstract

Most ultraluminous X-ray sources (ULXs) are believed to be stellar mass black
holes or neutron stars accreting beyond the Eddington limit. Determining the
nature of the compact object and the accretion mode from broadband
spectroscopy is currently a challenge, but the observed timing properties
provide insight into the compact object and details of the geometry and
accretion processes. Here we report a timing analysis for an $800$ ks XMM-
Newton campaign on the supersoft ultraluminous X-ray source, NGC 247 ULX-1.
Deep and frequent dips occur in the X-ray light curve, with the amplitude
increasing with increasing energy band. Power spectra and coherence analysis
reveals the dipping preferentially occurs on $\sim 5$ ks and $\sim 10$ ks
timescales. The dips can be caused by either the occultation of the central
X-ray source by an optically thick structure, such as warping of the accretion
disc, or from obscuration by a wind launched from the accretion disc, or both.
This behaviour supports the idea that supersoft ULXs are viewed close to edge-
on to the accretion disc.

###### keywords:

Accretion discs – X-rays: binaries – X-rays: individual: NGC 247 ULX-1

††volume: in press††pubyear: 2021††pagerange: Quasi-periodic dipping in the
ultraluminous X-ray source, NGC 247 ULX-1–Quasi-periodic dipping in the
ultraluminous X-ray source, NGC 247 ULX-1

##  1 Introduction

Ultraluminous X-ray sources (ULXs) are X-ray bright, off-nuclear sources with
luminosity $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}10^{39}$ ${\rm erg}~{}{\rm s}^{-1}$(see
Kaaret et al. 2017 for a review). Proposed explanations for these luminous
sources include i. sub-Eddington accreting stellar mass compact objects whose
radiation is beamed into our line of sight (King et al. 2001); ii. stellar-
mass BHs in a super-Eddington regime (Poutanen et al. 2007) or, iii.
intermediate-mass black holes accreting at sub-Eddington rates from low-mass
companion stars (Colbert & Mushotzky 1999, Farrell et al. 2009).

Figure 1: Background subtracted light curves for 8 new XMM-Newton observations
with $dt=200$ s. The soft ($0.3-1.0$ keV) and hard ($1.0-4.0$ keV) bands are
shown in blue and green respectively. The grey curve shows the estimated
background contribution. Observations are separated by a vertical dotted line
for plotting purposes, and are denoted by their XMM-Newton revolution number
(e.g. r$3660$) and OBSID (where $*=084486$).

The recent detection of coherent pulsations in several ULXs indicates they are
powered by neutron stars (NS) accreting at very high Eddington rates, with
luminosities up to $\sim 500\hbox{$L_{\rm Edd}$}$ (Bachetti et al. 2014,
Israel et al. 2017a, b, Carpano et al. 2018; Sathyaprakash et al. 2019,
Rodriguez Castillo et al. 2019, Chandra et al. 2020). Pulsations are not
always detectable as high count rates and pulsed fraction are required. In the
absence of pulsations, spectral analysis alone makes it difficult to
distinguish between a BH or NS candidate (e.g. NGC 1313 X-1, Walton et al.
2020), which means that NSs are likely numerous among the compact objects of
ULXs (King & Lasota 2019; Pintore et al. 2017; Walton et al. 2018; Middleton &
King 2017; Wiktorowicz et al. 2019).

Fourier timing techniques are powerful tools for investigating in a model-
independent way the geometry of accretion discs, the size scales involved, and
the nature of the compact objects (e.g. Middleton et al. 2007; Heil et al.
2009; Pinto et al. 2017; Kara et al. 2020). For instance, the presence of time
lags between certain energy bands identifies different regions of X-ray
emission in the accretion disc. The soft time lags found in three ULXs
indicates that substantial Compton down-scattering is occurring in the outer
disk (e.g., De Marco et al. 2013; Pinto et al. 2017; Kara et al. 2020). Power
spectral density (PSD) analysis allows us to distinguish different accretion
processes occurring at various accretion rate and identify different zones
within the accretion disc.

Dipping behaviour in the X-ray light curve have also been observed in some
ULXs, e.g. NGC 55 (Stobbart et al. 2004), M 94 (Lin et al. 2013), NGC 628 (Liu
et al. 2005), NGC 5408 X-1 (Pasham & Strohmayer 2013; Grisé et al. 2013 and M
51 ULX-7 (Hu et al. 2021; Vasilopoulos et al. 2021). The number of dips
observed in these sources has so far been limited, however they suggest we are
viewing these sources at high inclination or could be evidence for a
precessing disc.

NGC 247 ULX-1 is a supersoft ultraluminous X-ray source. It is unique in that
it switches between the soft ultraluminous X-ray state (SUL) and the supersoft
ultraluminous state (ULS or SSUL, e.g. Feng et al. 2016). Two previous short
($\sim 30$ ks), XMM-Newton observations revealed strong dips in source flux in
the light curve. These dips lasted several ks with the flux decreasing by an
order of magnitude (Feng et al. 2016). In this paper we report a timing study
of the dipping in the light curves of NGC 247 ULX-1 using a recently awarded
$800$ ks XMM-Newton campaign (PI: Pinto). PSD analysis reveals multiple
preferred periodic-like timescales for the dips. A detailed spectral analysis
is provided in two companion papers: Pinto et al. 2021 on the high-resolution
spectroscopy and D’Aì et al. in prep on the broadband spectroscopy.

##  2 Observations and data analysis

NGC 247 ULX-1 was observed 8 times by XMM-Newton over 1 month, starting
2019-12-03 (see Pinto et al. 2021 for observation details). Seven observations
have duration $\sim 115$ ks and one lasts for $\sim 55$ ks. Two short ($\sim
35$ ks) observations from 2016 are not included in this analysis as they are
too short to probe Fourier frequencies below $\sim 1\times 10^{-4}$ Hz.

The EPIC-pn Observation Data Files (ODFs) were processed following standard
procedures using the XMM-Newton  Science Analysis System (SAS; v19.0.0), using
the conditions PATTERN 0–4 and FLAG = 0. Source counts were extracted from a
circular region with radius 20 arc seconds. The background is extracted from a
large rectangular region on the same chip, avoiding the Cu ring on the outer
parts of the chip. Following Alston et al. (2013, 2019) for flares with
duration $\leq 200$ s, the source light curve was cut out and interpolated
between the gap by adding Poisson noise using the mean of neighbouring points.
The interpolation fraction was typically $<0.3$ %. High background segments,
typically towards the end of the revolutions (e.g. revolutions (revs)
$3633,3675$), were excluded from the analysis.

Figure 1 shows the background subtracted light curves for the $0.3-1.0$ keV
(hereafter soft) and $1.0-4.0$ keV (hereafter hard) bands. These bands are
motivated by broadband spectral analysis in Pinto et al. 2021 and D’Ai et al.
in prep. Above $4.0$ keV the light curve is dominated by the background so we
exclude it from the PSD analysis in Sec. 3.1, however these bands are used
where appropriate in Sec. 3.6.

Clear dipping behaviour is apparent in four of the observations (revs
$3663,3675,3676,3677$), with the dip duration ranging from $~{}\sim 1-20$ ks.
Four observations do not show any obvious dipping behaviour: revs
$3660,3674,3678,3680$. We therefore divide the observations into two data sets
- the ‘dip’ and ‘non-dip’ segments. With the exception of rev $3663$, the dip
segments have higher source flux and are spectrally harder (see also D’Ai et
al. in prep, Pinto et al. 2021). In the following analysis, we use a segment
length of $110$ ks for the dip epochs and $55$ ks for the no-dip epochs, owing
to the shorter duration of rev $3680$.

##  3 Timing analysis

###  3.1 The power spectrum

Figure 2: The Poison-noise subtracted PSDs in $\nu P(\nu)$ units for the
‘dips’ (blue circles) and ‘no dips’ (orange open squares) segments. The upper
panel shows the hard band ($1.0-4.0$ keV) and the middle panel shows the soft
band ($0.3-1.0$ keV). The lower panel shows the coherence between the hard and
soft bands. A simple cubic spline has been overlaid to highlight the frequency
dependence.

The PSD was estimated using the standard method of calculating the periodogram
(e.g. Priestley 1981), with an ${\rm[rms/mean]}^{2}$ normalisation, then
averaging over $M$ segments at each Fourier frequency. Non-overlapping
neighbouring frequency bins are then averaged (see e.g. Vaughan et al. 2003).
In order to pick out features over a broad range of frequencies we used a
geometrically averaged frequency bins, each bin spanning a factor
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}1.3$ in frequency larger than the previous
bin. Figure 2 shows the soft and hard band noise-subtracted PSDs in units of
$\nu P(\nu)$ for both the dipping and non-dipping segments (the noise is
estimated from standard formulae, e.g. Vaughan et al. 2003).

The dipping observations show clear peaks in power at $\sim 1\times 10^{-4}$
Hz and $\sim 2\times 10^{-4}$ Hz. A further broader peak can be seen at $\sim
8\times 10^{-4}$ Hz. At frequencies below $\sim 2\times 10^{-5}$ Hz the PSD is
rolling over towards zero power. With the exception of the appearance of
narrow peaks in power, the shape of the dip and no-dipping segments in the
hard band looks similar. The non-dip PSD has a smooth power-law like frequency
dependence, with a steeper slope in the hard band.

###  3.2 Coherent variations

We study the correlated variations between the two bands using the Fourier
coherence. This provides a measure of the linear correlation between two time
series $x(t),y(t)$ as a function of Fourier frequency, in the range [0,1]
(Vaughan & Nowak 1997). We estimated cross-spectrum, $C_{xy}(f)$, by first
averaging the complex values over $M$ non-overlapping segments, then averaging
over neighbouring frequency bins, in a similar manner to the PSD. This
amplitude part of the cross-spectrum gives the coherence, $\gamma^{2}$, which
is then corrected for Poisson noise, see e.g. Uttley et al. 2014 for more
details.

The coherence is shown in Figure 2 for both the dipping and non-dipping light
curves. A simple cubic spline has been fit to the data to illustrate the
frequency dependence of the data. The coherence is high ($\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.9$) at frequencies
where the soft and hard band PSDs peak, and drops at frequencies between these
peaks in power. In contrast, a high coherence at low frequencies is observed
in the non-dipping light curves, which falls off sharply to zero at $\sim
2\times 10^{-4}$ Hz, most likely caused by the onset of Poisson noise at high
frequencies.

The phase part of the cross-spectrum provides an estimate of the frequency
dependent delay between $x(t),y(t)$, which can be converted to a time delay in
seconds (Vaughan & Nowak 1997). We estimated this for both the dipping and
non-dipping segments and found the delay is consistent with zero lag at all
Fourier frequencies.

###  3.3 Modelling the power spectra

To characterise the observed peaks in variability power in the dip segments we
modelled the PSD using multiple Lorentzian functions;

$P(\nu)=\frac{N_{\rm l}(\sigma_{\rm l}/2\pi)}{(\nu-{\nu}_{\rm
c})^{2}+(\sigma_{\rm l}/2)^{2}}$ (1)

where ${\nu}_{\rm c}$ and $\sigma_{\rm l}$ are the line centroid and width
parameters, respectively, and $N_{\rm l}$ is a normalisation. The quantity
$Q=\nu_{c}/\sigma_{l}$ gives a measure of how narrow, or ‘coherent’ the
Lorentzian is. Here, instead of subtracting the Poisson (‘white’) noise in
each PSD band we model this with an additive constant, $C$. The hard and soft
band have similar shape, so we jointly model the data, allowing certain
parameters untied in model variants described below. The models are fit by
minimising the $\Large\chi^{2}$ statistic within the package xspec. The
parameter and $\Large\chi^{2}$ space were explored with Goodman & Weare (2010)
Markov Chain Monte Carlo (MCMC) sampling. Quoted parameter errors are from the
$90\%$ credible intervals. We used a linear binning factor of $5$ to improve
the frequency resolution in the modelling. With a time bin $dt=50$ s this
results in a combined total of $1602$ data bins for the soft and hard PSDs.
The resulting PSDs and model fits along with the model standard residuals are
shown in Figure 3.

Figure 3: PSD modelling of the dipping observations using multiple Lorentzian
components. The best fit model (BB+3Lor($Q_{\rm l_{1}}=5$) with a zero-
centered Lorentzian (dashed) and three additional narrower Lorentzians
(dotted) are shown for the hard and soft bands in the upper and middle panels
respectively. The horizontal dot-dashed line is the Poisson noise parameter.
The lower three panels show the standard residuals for this and two less
complex models.

The Akaike information criterion was used for model selection to compare
models with additional complexity (Akaike 1974). The corrected form of this
statistic takes account for the bias from the small sample size, which can be
expressed as

${\rm AICc}=2k+{\chi}^{2}+\frac{2k^{2}+2k}{n-k-1}$ (2)

where $k$ is the number of free model parameters, $n$ is the number of data
bins in the PSD and $\chi^{2}$ is used to represent the negative of twice the
logarithm of the likelihood of the model. The model with the lowest ${\rm
AICc}$ is the most preferred over all models. The ${\rm\Delta AICc}={\rm
AICc}-\min({\rm AICc})$, with a ${\rm\Delta AICc}<2$ suggesting weak evidence
in favour of the model. A ${\rm\Delta AICc}\sim 3-7$ suggests moderate
evidence for the model and $>10$ implies the model is strongly favoured (e.g.
Burnham & Anderson 2007).

We initially model the PSDs using a single broad Lorentzian component (model:
Broadband), with all parameters free, giving $\Large\chi^{2}$/dof (degrees of
freedom) $=1884.9/1594$. The $\Large\chi^{2}$ and AICc are given in Table 1
and the residuals of this model are shown in Figure 3. The next more complex
model has the addition of a narrower Lorentzian ($l_{2}$) at $\sim 2\times
10^{-4}$ Hz (model: BB+1Lor). With ${\nu}_{\rm c}$ and $\sigma_{\rm l}$ tied
between the two PSDs and $N_{\rm hard}$ and $N_{\rm soft}$ free, this improves
the fit by $\Delta\Large\chi^{2}=22.7$ for 4 additional dof. Compared to the
‘Broadband’ model the ${\rm\Delta AICc}=14.6$, indicating a strong preference
for the more complex model. Allowing any of the parameters to be untied
between the two PSDs does not notably improve the model fit.

The addition of a Lorentzian ($l_{3}$) with ${\nu}_{\rm c}\sim 8\times
10^{-4}$ Hz (model BB+2Lor) reduces the $\Large\chi^{2}=14.9$ for 4 additional
dof. Untying the parameters of the additional Lorentzian between the two PSDs
does not significantly improve the fit. Compared to the previous model the
${\rm\Delta AICc}=10.9$, again suggesting the more complex model is strongly
supported.

The most complex model we test involves the addition of a narrow Lorentzian
($l_{1}$) with ${\nu}_{\rm c}\sim 1\times 10^{-4}$ Hz (model: BB+3Lor). With
all but the normalisation parameters tied across the two energy bands, the
best fit model gives $\Large\chi^{2}$/dof $=1824.1/1584$, with 4 additional
parameters compared to the previous model. The ${\rm\Delta AICc}=14.98$, again
suggesting the more complex model is strongly supported. The low sampling
resolution at low frequencies means the fit is largely insensitive to
${\nu}_{\rm c}$ of the broadband Lorentzian, so we fix this to zero (model:
BB+3Lor($\nu_{\rm c_{\rm BB}}=0$)). This gives no notable change to the fit,
but helps remove some parameter degeneracies between ${\nu}_{\rm c_{\rm BB}}$
and $\sigma_{\rm BB}$. Equally, ${\nu}_{\rm c,l_{1}}$ and $\sigma_{\rm l_{1}}$
(i.e. the Lorentzian at ${\nu}_{\rm c}\sim 1\times 10^{-4}$ Hz also display
degeneracies, so we fix these at approximately the median MCMC values of
$Q=\nu_{c}/\sigma_{l}=5$ (model: BB+3Lor($Q_{\rm l_{1}}=5$)). The parameters
for this model are shown in Table 2. The addition of more Lorentzian
components provides an improvement of ${\rm\Delta AICc}<2$, so we do not
consider those models further.

For the no-dip segments, a single zero-centered Lorentzian provides an
adequate description of the hard band PSD. The additional parameters are
$\sigma_{\rm no-dip,hard}=4.55$${}_{-2.28}^{+1.25}$ $\times 10^{-5}$ Hz and
$N_{\rm no-dip,hard}=0.29$${}_{-0.07}^{+0.06}$. Additional components do not
significantly improve the AICc so we do not explore this further.

Model | $N_{\rm pars}$ |  $\Large\chi^{2}$ / d.o.f | ${\rm AICc}$ | ${\rm\Delta AICc}$
---|---|---|---|---
$(1)$ | $(2)$ | $(3)$ | $(4)$ | $(5)$
Broadband (BB) | 8 | 1884.9 / 1594 | 1900.9 | 40.5
BB+1Lor | 12 | 1862.2 / 1590 | 1886.3 | 25.9
BB+2Lor | 14 | 1847.2 / 1588 | 1875.5 | 15.1
BB+3Lor | 18 | 1824.1 / 1584 | 1860.5 | 0.1
BB+3Lor ($\nu_{\rm c_{\rm BB}}=0$) | 17 | 1826.4 / 1585 | 1860.8 | 0.4
BB+3Lor ($Q_{\rm l_{1}}=5$) | 16 | 1828.0 / 1586 | 1860.4 | *min*
Table 1: Fit statistics for the PSD model to the hard and soft band dip segments. Columns (1) states the model, (2) the number of free parameters, (3) is the $\Large\chi^{2}$ / d.o.f, (4) is the corrected Akaike criterion, ${\rm AICc}$, and (5) is the ${\rm\Delta AICc}$ compared to the minimum model. Comp. | $\nu_{\rm c}$ | $\sigma_{\rm l}$ | $N_{\rm hard}$ | $N_{\rm soft}$ | Q
---|---|---|---|---|---
| $\times 10^{-4}$ | $\times 10^{-4}$ | $\times 10^{-1}$ | $\times 10^{-2}$ | -
$(1)$ | $(2)$ | $(3)$ | $(4)$ | $(5)$ | $(6)$
$l_{\rm BB}$ |  $0$ fixed |  $0.37$${}_{-0.28}^{+0.25}$ |  $1.80$${}_{-0.70}^{+0.90}$ |  $7.95$${}_{-2.89}^{+0.04}$ | -
$l_{1}$ |  $0.89$${}_{-0.14}^{+0.04}$ | $0.18$ |  $0.89$${}_{-0.47}^{+0.37}$ |  $0.49$${}_{-0.28}^{+1.02}$ | 5f
$l_{2}$ |  $2.16$${}_{-0.23}^{+0.08}$ |  $1.20$${}_{-0.54}^{+0.84}$ |  $1.05$${}_{-0.46}^{+0.31}$ |  $2.44$${}_{-1.06}^{+0.01}$ | 2
$l_{3}$ |  $6.95$${}_{-6.01}^{+1.43}$ |  $7.07$${}_{-2.0}^{+4.2}$ |  $0.28$${}_{-.12}^{+.53}$ |  $2.0$${}_{-1.1}^{+1.0}$ | 1
Table 2: Table for the best fitting model parameters for the BB+3Lor ($Q_{\rm
l_{3}}=5$) model. Errors are from the MCMC posterior $90\%$ credible
intervals. This model consists of a broad zero-centered Lorentzian (BB) plus 3
narrower Lorentzians ($l_{1},l_{2},l_{3}$) plus an additive constant $P_{\rm
n}$ for Poisson noise. All components are tied across energy bands except
$N_{l}$ and the noise levels, $P_{\rm n,hard}=23.95$${}_{-0.45}^{+0.34}$, and
$P_{\rm n,soft}=9.53$${}_{-0.16}^{+0.14}$

###  3.4 Periodogram of dipping observations

To check that no one particular observation is driving the peaks in the
average PSD we estimate the unbinned periodogram. Figure 4 shows the
periodogram for the four dipping observations. It is well known that the
periodogram estimate at each frequency follows a ${\Large\chi}^{2}_{\nu=2}$
distribution so large scatter around the true underlying power spectrum is
expected (e.g. Priestley 1981). We therefore also show a Gaussian smoothed
with periodogram over $n=3$ neighbouring frequency bins. Peaks in power at
$\sim 1\times 10^{-4}$ Hz and $\sim 2\times 10^{-4}$ Hz are found in each
periodogram. Another peak or bend can be seen at $\sim 0.4\times 10^{-4}$ Hz
in all but rev $3663$, which shows power at lower frequencies.

We refit these periodograms using the maximum likelihood methods (MLE) for
unbinned periodograms outlined in Vaughan (2010) and Barret & Vaughan (2012).
The power spectral modelling results are consistent with those found in
Section 3.3, so we do not show those here. For rev $3676$ the Lorenztzian
centroids $l_{1}$ and $l_{2}$ appear to shift to higher frequencies, however
modelling reveals these to be consistent within error bars with the other
revs. The only notable difference is the increase in power at low frequencies
in rev $3663$, which requires a zero-centered Lorentzian with width
$\sigma_{l_{\rm BB}}\sim 1\times 10^{-5}$ Hz. Removing rev $3663$ from the
combined PSD modelling in Section 3.3 does not significantly alter the
results, so we choose to leave it in to increase the number of estimates in
each frequency bin.

Figure 4: Periodograms for the hard band dipping observations (without Poisson
noise subtraction). The green curve is the Gaussian smoothed periodogram over
$n=3$ neighbouring bins.

###  3.5 Change point analysis

In order to study what aspect of the dips is causing the peaks in the power
spectra, we modelled the abrupt light curve changes in the dip segments using
change point analysis. The R package
bcp111https://cran.r-project.org/web/packages/bcp/ was used which implements
the Bayesian change point method of Wang & Emerson 2015. This partitions the
light curve into constant mean ‘blocks’ between change points. MCMC sampling
is used to model the posterior probability at each change point, $p_{\rm ch}$,
defined between [0,1]. Given the sharp changes in the light curves when
dipping occurs, this method is particularly well suited to the data.

The Poisson noise starts to dominate at $\sim 2\times 10^{-3}$ Hz so a time
bin of $dt=500$ s was used in order to smooth over any noise variations. A
posterior probability $p_{\rm ch}>0.4$ was used to identify significant change
point locations. From these we determine the dip onset time $t_{\rm 0}$, dip
duration $t_{\rm d}$, and dip midpoint, $t_{m}=(t_{\rm 0}+t_{\rm d})/2$. The
period $P$ is then taken as the time between successive dip midpoints.

Figure 5 shows the distribution of $t_{\rm d}$ and $P$ for all four dipping
observations. A Gaussian kernel density estimate is overlaid on the observed
distributions. The distribution of $t_{\rm d}$ peaks at $\sim 3$ ks and falls
off approximately log-normally with increasing dip duration, i.e. there is a
much larger probability for shorter dips. The distribution of $P$ is multi-
modal with peaks at $\sim 5$ and $10$ ks, with a smaller probability of
periods at $\sim 17$ ks. A similar distribution is observed if we represent
the $P$ as the time between successive $t_{\rm 0}$, suggesting no relationship
between the duration of the dips and the dip period.

Figure 5: Density of dip duration, $t_{\rm d}$, (left) and dip period, $P$,
(right) as estimated using Bayesian change point analysis. The Gaussian kernel
density estimate is shown in orange.

###  3.6 Variability spectra

We study the energy dependence of the variability at a given timescale using
frequency resolved covariance spectra (Wilkinson & Uttley 2009). Using a high
S/N reference band, the correlated variability is picked out in a given
comparison band, thus improving the S/N of the variability spectrum compared
to traditional rms spectra. We compute the covariance spectra in the Fourier
domain (Uttley et al. 2011; Cassatella et al. 2012; Uttley et al. 2014). The
$0.3-4.0$ keV band is used as the broad reference band, but the choice of this
band does not affect the shape of the resultant covariance spectra.

We estimate covariance spectra in three frequency bands: $f_{1}=2.0-4.0\times
10^{-5}$ Hz, $f_{2}=0.9-1.1\times 10^{-4}$ Hz and $f_{3}=1.9-2.4\times
10^{-4}$ Hz. These timescales probe the long timescale spectral variability as
well as the spectral variability associated with the two prominent dipping
timescales. Figure 6 shows the resulting variability spectra. The top panels
shows these data in absolute units (${\rm cts}~{}{\rm s}^{-1}~{}{\rm
keV}^{-1}$) whilst the lower panels show the variability spectra in fractional
($\%$) units. In the dipping segments, the long timescale spectral variability
is approximately flat, with a $10\%$ rms amplitude. The dipping timescales,
$f_{2}$ and $f_{3}$, show an approximately linear dependence with energy, with
larger variability amplitude at higher energies. This is true of all frequency
bands up to higher frequencies, but we do not show those here.

For the dipping segments, as an additional check, we calculate the variability
amplitude of each of the components in the preferred model (BB+3Lor $Q_{\rm
l_{3}}=5$) by integrating the rms in each Lorentzian component. This shows the
same trend in energy space as the frequency dependent covariance spectra. For
the no-dip segments, the fractional variability amplitude is almost below $1$
keV, with a sharp increase in rms for energies above this.

Figure 6: Variability spectra for the dip (left) and no-dip (right) segments.
The frequency ranges correspond to the low frequency, the first and the second
dip peaks respectively. The upper panels show the time averaged and
variability spectra in absolute units. The lower panels show the fractional
covariance spectra for each frequency band. Missing bins are due to high
Poisson noise in a given energy band.

##  4 Discussion

###  4.1 The dipping behaviour

The light curves of NGC 247 ULX-1 display deep and quasi-periodic dips
imprinted on top of stochastic noise variations. Modelling the PSD with
multiple Lorentzians reveals three preferred timescales for the dipping
behavior. The peak in power of a Lorentzian occurs at frequency $\nu_{\rm
max}=\sqrt{\nu_{c}^{2}+\sigma_{l}^{2}}$ (see e.g. Belloni et al. 2002). For
the preferred model, this gives $\nu_{\rm max,BB}\sim 0.4\times 10^{-4}$ Hz,
$\nu_{\rm max,l_{1}}\sim 0.9\times 10^{-4}$ Hz, $\nu_{\rm max,l_{2}}\sim
2\times 10^{-4}$ Hz and $\nu_{\rm max,l_{3}}\sim 9\times 10^{-4}$ Hz. The
coherence between the hard and soft bands is high ($\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.9$) at the
Lorentzian peaks and drops where the Lorentzians overlap, providing further
evidence for multiple dipping timescales.

Time series change point analysis reveals the dip periods, $P$, cluster around
$\sim 5$ and $10$ ks. This distribution of $P$ is consistent with the peak
timescales observed in the power spectra (see Figure 3). The Lorentzian $Q$
values tell us that there is some scatter around a preferred periodicity, as
is observed in the change point period analysis. The dip durations $t_{d}$
preferentially take on values between $1-5$ ks, which will add some additional
broadening to the measured Lorentzian $Q$ values. There are several broad dips
with $t_{d}>15$ ks, which clearly stand out from the other measured dip
durations in Figure 5. This could be suggesting an alternative mechanism or
timescale for these longer duration dips.

The change point analysis does not find any periods with $P<1000$ s suggesting
this analysis is either insensitive to this dips, e.g. they could be of lower
amplitude and therefore harder to detect, in agreement with the lower
amplitude of the Lorentzian $l_{3}$. Alternatively, this high frequency
component is related to an increase in the broadband noise continuum, which
appears when the source enters the higher-flux / dipping epochs.

The time averaged and frequency resolved covariance spectra in Figure 6 show
that the variations are larger at higher energies in both the dip and no-dip
epochs. The exception being the $f_{1}=2.0-4.0\times 10^{-5}$ Hz frequency
band in the dip epochs. Here, the fractional covariance is independent of
energy. The coherence between the soft and hard bands (see Figure 2) also
drops to $\sim 0.6$ at these frequencies. This suggests that the long
timescale variability is different to the dip generating mechanism at shorter
timescales.

###  4.2 Origin of the dips in NGC 247 ULX-1

Assuming a black hole with $M_{\rm BH}=10M_{\odot}$, the peak timescales
$\nu_{\rm max}\sim 0.5-9\times 10^{-4}$ Hz correspond to an orbital radius
(${}^{3}\sqrt{GMt^{2}/4\pi^{2}}$) of $10^{11}$ cm, or $2-10\times 10^{4}$
$R_{G}$ (where $R_{G}=GM/c^{2}$). A similar value is obtained for the case of
a NS compact object. In both cases, this region is very distant from the
compact object, suggesting an origin in the outer accretion flow. Two likely
scenarios for the dips are warping of an accretion disc seen close to edge on,
periodically obscuring the central source each orbit, or, the launching of a
wind following changes in local $\dot{M}$. Alternatively, a combination of
both processes could be causing the dips at different timescales.

The radius associated with Compton scattering through a medium around the
spherisation radius can be expressed as $R_{\rm C}=\sqrt{ct_{\rm
esc}/\sigma_{\rm T}n_{\rm e}}$ (where $c$ is the speed of light, $\sigma_{\rm
T}$ is the Thomson scattering cross-section, $n_{\rm e}$ is the electron
density, and $t_{\rm esc}$ is the timescale for photons to escape) would be of
the order of $10^{9}$ cm, assuming a disc-like $n_{\rm e}$ density ($10^{20}$
cm-3), while such radius would be comparable ($10^{11-12}$ cm) if we adopt a
wind-like density ($10^{14}$ cm-3). On the other hand, the wind radius,
$R_{\rm wind}=L_{\rm ion}/n_{\rm e}\xi$ (where $L_{\rm ion}$ is the ionising
luminosity and $\xi$ the ionisation parameter), that we estimate from the
best-fit parameters of the XMM-Newton RGS-EPIC joint fit is also around
$10^{11-13}$ cm (see Pinto et al. 2021). The lower density assumption is
favoured as this agrees with both the outflow and dip timescale properties.

The mid plane of the accretion disc, i.e. the spherisation radius where the
disc height is comparable to the distance from the accretor, is instead
$R_{\rm sph}=5/3\dot{M}R_{\rm in}\sim 100$ $R_{G}$ (assuming a standard
$R_{in}=6$ $R_{G}$ and $\dot{M}\sim 10$ $\dot{M}_{\rm Edd}$, which would be
sufficient to achieve $L_{X}\sim 5\times 10^{39}\hbox{${\rm erg}~{}{\rm
s}^{-1}$}$ for $M_{\rm BH}=10M_{\odot}$). This means that the timescales of
the variability detected in NGC 247 ULX-1 are likely associated with the outer
disc and, most likely, with the region where the wind becomes optically thin
(e.g. $n_{\rm e}\lesssim 10^{14}$ cm-3) and rather than with the inner
accretion flow or the optically thick wind ($10^{20}$ cm-3) which would be
expected near the spherisation radius. For a comparison, the wind photosphere
in the Galactic super-Eddington source SS 433 is $\sim 10^{12}$ cm (see, e.g.
Fabrika 2004), which further indicates that NGC 247 ULX-1 quasi-periodic dips
might be associated with some structure present in the outer disc where the
photosphere becomes optically thin.

The presence of dips and their long timescales argue in favour of a moderate-
to-high inclination in this source, which agrees with the scenario in which
soft ULXs are primarily seen along the edge of the wind. The dips seen in both
BH and NS X-ray Binaries (XRBs) are believed to be caused by periodic
obscuration of the central X-ray source by structure located in the outer
regions of a disc. These dips are observed in sources at high inclination
($i>60$ deg, e.g. Frank et al. 1987). For NGC 247 ULX-1, Yao & Feng (2019)
derived an $i>67$ deg via modeling the UV/optical spectral energy distribution
(SED) with an irradiation model, providing further evidence for a high
inclination viewing angle in this source. Additionally, winds observed in BH
or NS X-ray binaries (XRBs) have been mostly observed in sources in which the
disc is inclined at a large angle to the line of sight (e.g. Ponti et al.
2012; Díaz Trigo & Boirin 2016).

Given the difficulty to access the inner regions of the accretion flow, mostly
obscured by the scale height of the disc at the spherisation radius and the
wind, we are not surprised that we did not detect the time lags between the
soft and hard band that have been found in other, harder, ULXs. At these
angles, the photons would undergo loads of Compton scattering, losing the
original information.

If the dips are caused by a wind structure, what is not clear is whether the
short-term variations in optical depth are just stochastic variability (at a
given $\dot{M}$) due to the clumpy nature of the wind. Alternatively, each
variation corresponds to a change in the underlying $\dot{M}$ variability in
the local accretion rate where the wind is launched. This variability in the
wind launching would thereby cause these quasi periodic dips, rather than
being pre-existing structure with the appropriate scale height. This idea is
supported by the onset of dipping behaviour in the higher flux observations.
In fact, the bolometric luminosity of SSUL sources slightly increases when the
blackbody temperature decreases and the blackbody radius increase. This
behaviour is common to several other SSUL (e.g. Feng et al. 2016; Urquhart &
Soria 2016a).

The broadband spectral model decomposition in D’Aì et al,  in prep reveals two
broad thermal components. These are interpreted as emission from an extended,
optically thick, cold photosphere, and a hotter emission component,
originating from regions closer to the compact object. In the ‘out of dip’
epochs (i.e. when a dip is not occurring), the short term variations are
predominantly driven by changes in the flux of the hotter component. During
the dips, both broad spectral components show significant decrease in their
flux, particularly the hotter component. One scenario, where the dipping
mechanism is caused by occultation further out than the typical radii where
both components originate, would be consistent with what we infer from the
timing properties in the present paper.

An alternative explanation for the dips is that of intrinsic luminosity
variations due to the propeller effect switching on-and-off accretion in a
magnetised NS (see e.g. Revnivtsev & Mereghetti 2015 for a review). This was
suggested as the origin of large amplitude variations in the X-ray light curve
of the ULX M 82 X-1 by Tsygankov et al. (2016), however these occur on
timescales of $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}$ weeks. Interestingly, the ‘hiccuping’
X-ray light curve of the accreting millisecond X-ray pulsar, IGR J18245-2452,
displays a remarkable resemblance to those in NGC 247 ULX-1 (Ferrigno et al.
2014), However, the dips in this source show complex spectral hardening,
whereas in NGC 247 ULX-1 the dips are spectrally harder (see Pinto et al. 2021
Figure 2). The detection of pulsations in NGC 247 ULX-1 would provide more
evidence for this scenario, which is scope for future work.

###  4.3 Dips in other ULXs

For NGC 55, Stobbart et al. (2004) report several dips in the XMM-Newton light
curve. These are fewer in number, but have a similar dip duration and
separation between consecutive dips. The authors argued the dips were likely
due to material at $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}2\times 10^{11}$ cm from the central
object (assumed to be a BH). At this radius the orbital timescale is $\approx
10^{4}$ s, consistent with our PSD peak timescales. This is in agreement with
the idea that a disc warp or extended wind structure is causing the dips each
orbit in NGC 247 ULX-1.

The two broad dips with $t_{d}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}15$ ks in observations rev $3663$ and rev
$3677$ are similar to the broad deep dips seen in M 51 ULX-1 (Urquhart & Soria
2016b). These authors argued that these dips are caused by eclipsing events
and placed limits on the orbital period. CG X-1, a ULX in the galaxy Circinus,
has also shown multiple eclipsing events (Qiu et al. 2019), with the authors
finding an orbital period of $\sim 7$ hrs. It is possible that we are seeing a
similar eclipses in NGC 247 ULX-1, meaning the separation of $\sim 19$ days
for these deep dips may correspond to some integer multiple of the orbital
period. If this is the case, then the dipping activity observed in NGC 247
ULX-1 may be triggered by the companion at some point in the orbital phase.
This remains scope for future work.

A tentative $\sim 2$ hr ($\sim 1\times 10^{-4}$ Hz) quasi-periodic oscillation
(QPO) was claimed in NGC 628 (Liu et al. 2005). The light curve of this source
displays a similar dipping process to NGC 247 ULX-1. The authors scaled this
frequency to QPOs typically seen in black hole X-ray binaries (XRBs; see e.g.
Remillard & McClintock 2006) and inferred a black hole mass $M_{\rm BH}\sim
10^{3}M_{\odot}$. However, the light curve properties of the ULX dips are
strikingly different to QPOs seen in XRBs and active galaxies (e.g. Ingram &
Motta 2020, Alston et al. 2014; Alston et al. 2015, 2016). If the dips are
caused by a similar process to NGC 247, then this would argue in favour of a
stellar mass compact object in NGC 628.

##  5 Conclusions

A long XMM-Newton observation of NGC 247 ULX-1 has revealed dipping in the
light curves on $\sim 1-10$ ks timescales. This is consistent with the idea
that supersoft ULXs are viewed from high inclination, with the dipping caused
by an obscuration process from the outer accretion region. These findings will
motivate future studies on the source geometry and emission mechanisms at
super-Eddington accretion.

## Acknowledgments

WNA is supported by an ESA research fellowship. AD, MDS and FP acknowledge
financial contribution from the agreement ASI-INAF n.2017-14-H.0 and INAF
main-stream. HPE acknowledges support under NASA contract NNG08FD60C. TPR
gratefully acknowledges support from the Science and Technology Facilities
Council (STFC) as part of the consolidated grant award ST/000244/1. This paper
is based on observations obtained with XMM-Newton, an ESA science mission with
instruments and contributions directly funded by ESA Member States and the USA
(NASA).

## Data Availability

All of the data underlying this article are publicly available from ESA’s XMM-
Newton Science Archive (XSA; https://www.cosmos.esa.int/web/xmm-newton/xsa)
and NASA’s HEASARC archive (https://heasarc.gsfc.nasa.gov/). Light curves are
available from the author upon request.

## References

  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Alston et al. (2013) Alston W. N., Vaughan S., Uttley P., 2013, MNRAS, 435, 1511
  * Alston et al. (2014) Alston W. N., Markevičiūtė J., Kara E., Fabian A. C., Middleton M., 2014, MNRAS, 445, L16
  * Alston et al. (2015) Alston W. N., Parker M. L., Markevičiūtė J., Fabian A. C., Middleton M., Lohfink A., Kara E., Pinto C., 2015, MNRAS, 449, 467
  * Alston et al. (2016) Alston W., Fabian A., Markevičiutė J., Parker M., Middleton M., Kara E., 2016, Astronomische Nachrichten, 337, 417
  * Alston et al. (2019) Alston W. N., et al., 2019, MNRAS, 482, 2088
  * Bachetti et al. (2014) Bachetti M., Harrison F. A., Walton D. J., Grefenstette B. W., Chakrabarty D., et al., 2014, Nature, 514, 202
  * Barret & Vaughan (2012) Barret D., Vaughan S., 2012, ApJ, 746, 131
  * Belloni et al. (2002) Belloni T., Psaltis D., van der Klis M., 2002, ApJ, 572, 392
  * Burnham & Anderson (2007) Burnham K., Anderson D., 2007, Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer New York, https://books.google.es/books?id=IWUKBwAAQBAJ
  * Carpano et al. (2018) Carpano S., Haberl F., Maitra C., Vasilopoulos G., 2018, MNRAS, 476, L45
  * Cassatella et al. (2012) Cassatella P., Uttley P., Maccarone T. J., 2012, MNRAS, 427, 2985
  * Chandra et al. (2020) Chandra A. D., Roy J., Agrawal P. C., Choudhury M., 2020, MNRAS, 495, 2664
  * Colbert & Mushotzky (1999) Colbert E. J. M., Mushotzky R. F., 1999, ApJ, 519, 89
  * De Marco et al. (2013) De Marco B., Ponti G., Miniutti G., Belloni T., Cappi M., Dadina M., Muñoz-Darias T., 2013, MNRAS, 436, 3782
  * Díaz Trigo & Boirin (2016) Díaz Trigo M., Boirin L., 2016, Astronomische Nachrichten, 337, 368
  * Fabrika (2004) Fabrika S., 2004, ASPR, 12, 1
  * Farrell et al. (2009) Farrell S. A., Webb N. A., Barret D., Godet O., Rodrigues J. M., 2009, Nature, 460, 73
  * Feng et al. (2016) Feng H., Tao L., Kaaret P., Grisé F., 2016, ApJ, 831, 117
  * Ferrigno et al. (2014) Ferrigno C., et al., 2014, A&A, 567, A77
  * Frank et al. (1987) Frank J., King A. R., Lasota J. P., 1987, A&A, 178, 137
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Grisé et al. (2013) Grisé F., Kaaret P., Corbel S., Cseh D., Feng H., 2013, MNRAS, 433, 1023
  * Heil et al. (2009) Heil L. M., Vaughan S., Roberts T. P., 2009, MNRAS, 397, 1061
  * Hu et al. (2021) Hu C.-P., Ueda Y., Enoto T., 2021, ApJ, 909, 5
  * Ingram & Motta (2020) Ingram A., Motta S., 2020, arXiv e-prints, p. arXiv:2001.08758
  * Israel et al. (2017a) Israel G. L., Belfiore A., Stella L., Esposito P., Casella P., et al., 2017a, Science, 355, 817
  * Israel et al. (2017b) Israel G. L., Papitto A., Esposito P., Stella L., Zampieri L., et al., 2017b, MNRAS, 466, L48
  * Kaaret et al. (2017) Kaaret P., Feng H., Roberts T. P., 2017, Ann. Rev. A&A, 55, 303
  * Kara et al. (2020) Kara E., Pinto C., Walton D. J., Alston W. N., Bachetti M., et al., 2020, MNRAS, 491, 5172
  * King & Lasota (2019) King A., Lasota J.-P., 2019, MNRAS, 485, 3588
  * King et al. (2001) King A. R., Davies M. B., Ward M. J., Fabbiano G., Elvis M., 2001, ApJ, 552, L109
  * Lin et al. (2013) Lin D., Irwin J. A., Webb N. A., Barret D., Remillard R. A., 2013, ApJ, 779, 149
  * Liu et al. (2005) Liu J.-F., Bregman J. N., Lloyd-Davies E., Irwin J., Espaillat C., Seitzer P., 2005, ApJ, 621, L17
  * Middleton & King (2017) Middleton M. J., King A., 2017, MNRAS, 470, L69
  * Middleton et al. (2007) Middleton M., Done C., Gierliński M., 2007, MNRAS, 381, 1426
  * Pasham & Strohmayer (2013) Pasham D. R., Strohmayer T. E., 2013, ApJ, 764, 93
  * Pinto et al. (2017) Pinto C., Alston W., Soria R., Middleton M. J., Walton D. J., et al., 2017, MNRAS, 468, 2865
  * Pinto et al. (2021) Pinto C., et al., 2021, arXiv e-prints, p. arXiv:2104.11164
  * Pintore et al. (2017) Pintore F., Zampieri L., Stella L., Wolter A., Mereghetti S., Israel G. L., 2017, ApJ, 836, 113
  * Ponti et al. (2012) Ponti G., Fender R. P., Begelman M. C., Dunn R. J. H., Neilsen J., Coriat M., 2012, MNRAS, 422, L11
  * Poutanen et al. (2007) Poutanen J., Lipunova G., Fabrika S., Butkevich A. G., Abolmasov P., 2007, MNRAS, 377, 1187
  * Priestley (1981) Priestley M., 1981, Spectral analysis and time series. Academic Press, London
  * Qiu et al. (2019) Qiu Y., et al., 2019, ApJ, 877, 57
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, Ann. Rev. A&A, 44, 49
  * Revnivtsev & Mereghetti (2015) Revnivtsev M., Mereghetti S., 2015, Space Sci. Rev., 191, 293
  * Rodriguez Castillo et al. (2019) Rodriguez Castillo G. A., Israel G. L., et al., 2019, arXiv e-prints, p. arXiv:1906.04791
  * Sathyaprakash et al. (2019) Sathyaprakash R., Roberts T. P., Walton D. J., Fuerst F., et al., 2019, MNRAS, p. L104
  * Stobbart et al. (2004) Stobbart A. M., Roberts T. P., Warwick R. S., 2004, MNRAS, 351, 1063
  * Tsygankov et al. (2016) Tsygankov S. S., Mushtukov A. A., Suleimanov V. F., Poutanen J., 2016, MNRAS, 457, 1101
  * Urquhart & Soria (2016a) Urquhart R., Soria R., 2016a, MNRAS, 456, 1859
  * Urquhart & Soria (2016b) Urquhart R., Soria R., 2016b, ApJ, 831, 56
  * Uttley et al. (2011) Uttley P., Wilkinson T., Cassatella P., Wilms J., Pottschmidt K., Hanke M., Böck M., 2011, MNRAS, 414, L60
  * Uttley et al. (2014) Uttley P., Cackett E. M., Fabian A. C., Kara E., Wilkins D. R., 2014, A&ARv, 22, 72
  * Vasilopoulos et al. (2021) Vasilopoulos G., Koliopanos F., Haberl F., Treiber H., Brightman M., Earnshaw H. P., Gúrpide A., 2021, ApJ, 909, 50
  * Vaughan (2010) Vaughan S., 2010, MNRAS, 402, 307
  * Vaughan & Nowak (1997) Vaughan B. A., Nowak M. A., 1997, ApJ, 474, L43
  * Vaughan et al. (2003) Vaughan S., Edelson R., Warwick R. S., Uttley P., 2003, MNRAS, 345, 1271
  * Walton et al. (2018) Walton D. J., Fürst F., Heida M., Harrison F. A., Barret D., et al., 2018, ApJ, 856, 128
  * Walton et al. (2020) Walton D. J., et al., 2020, MNRAS, 494, 6012
  * Wang & Emerson (2015) Wang X., Emerson J. W., 2015, arXiv e-prints, p. arXiv:1509.00817
  * Wiktorowicz et al. (2019) Wiktorowicz G., Lasota J.-P., Middleton M., Belczynski K., 2019, ApJ, 875, 53
  * Wilkinson & Uttley (2009) Wilkinson T., Uttley P., 2009, MNRAS, 397, 666
  * Yao & Feng (2019) Yao Y., Feng H., 2019, ApJ, 884, L3