#  Coupling of radiation and magnetospheric accretion flow in ULX pulsars:
radiation pressure and photon escape time

Caitlyn Flexer1 and Alexander A. Mushtukov2
1 Castilleja School, 1310 Bryant St, Palo Alto, CA 94301 Palo Alto, USA
2 Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK
E-mail<EMAIL_ADDRESS>(CF)E-mail:
<EMAIL_ADDRESS>(AAM)

###### Abstract

The accretion flow within the magnetospheric radius of bright X-ray pulsars
can form an optically thick envelope, concealing the central neutron star from
the distant observer. Most photons are emitted at the surface of a neutron
star and leave the system after multiple reflections by the accretion material
covering the magnetosphere. Reflections cause momentum to be transferred
between photons and the accretion flow, which contributes to the radiative
force and should thus influence the dynamics of accretion. We employ Monte
Carlo simulations and estimate the acceleration along magnetic field lines due
to the radiative force as well as the radiation pressure across magnetic field
lines. We demonstrate that the radiative acceleration can exceed gravitational
acceleration along the field lines, and similarly, radiation pressure can
exceed magnetic field pressure. Multiple reflections of X-ray photons back
into the envelope tend to amplify both radiative force along the field lines
and radiative pressure. We analyze the average photon escape time from the
magnetosphere of a star and show that its absolute value is weakly dependent
on the magnetic field strength of a star and roughly linearly dependent on the
mass accretion rate being $\sim 0.1\,{\rm s}$ at $\dot{M}\sim 10^{20}\,{\rm
g\,s^{-1}}$. At high mass accretion rates, the escape time can be longer than
free-fall time from the inner disc radius.

###### keywords:

accretion – accretion discs – X-rays: binaries – stars: neutron – stars:
oscillations

††pubyear: 2023††pagerange:  Coupling of radiation and magnetospheric
accretion flow in ULX pulsars: radiation pressure and photon escape time –A

##  1 Introduction

X-ray pulsars (XRPs) are strongly magnetized, accreting neutron stars (NSs) in
close binary systems (see Mushtukov & Tsygankov 2022 for review). The typical
magnetic field strength at the NS surface in XRPs is estimated to be $\sim
10^{12}\,{\rm G}$, which was determined from the detection of cyclotron lines
in the X-ray energy spectra of some objects (see Staubert et al. 2019 for
review). The apparent luminosity of XRPs covers about nine orders of magnitude
from $\sim 10^{32}\,{\rm erg\ \rm s^{-1}}$ to $\sim 10^{41}\,{\rm erg\ \rm
s^{-1}}$. The lower limit of accretion luminosity is due to the onset of the
centrifugal barrier and the “propeller" effect of accretion (Illarionov &
Sunyaev, 1975; Tsygankov et al., 2016). The brightest XRPs have luminosities
that exceed the Eddington limit for a NS: $L_{\rm Edd}\simeq{4\pi GMm_{\rm
p}c}/[\sigma_{\rm T}(1+{\rm x}_{H})]\simeq 2\times 10^{38}(M/M_{\odot})\,{\rm
erg\ \rm s^{-1}}$, where $m_{\rm p}$ is a mass of proton, $\sigma_{\rm T}$ is
the Thomson scattering cross-section, and ${\rm x}_{H}$ is the hydrogen mass
fraction in accreting material. These XRPs belong to the recently discovered
class of pulsating ultra-luminous X-ray sources (ULXs, see e.g., Fabrika et
al. 2021; King et al. 2023). Currently, six pulsating ULXs have been detected:
M82 X-2 (Bachetti et al., 2014), NGC 7793 P13 (Fürst et al., 2016; Israel et
al., 2017b), NGC 5907 X-1 (Israel et al., 2017a), NGC 300 X-1 (Carpano et al.,
2018), M51 X-7 (Rodríguez Castillo et al., 2020), NGC 1313 X-2 (Sathyaprakash
et al., 2019). Luminosity values comparable to the luminosity of ULX pulsars
can be detected in some Be X-ray transients at the peaks of the outbursts (see
Reig 2011 for review).

The strong magnetic field shapes the geometry of accretion flow in XRPs. The
accretion flow in form of stellar wind or accretion disc is interrupted by the
strong magnetic field at a distance called the magnetospheric radius, where
the accretion material starts to move along magnetic field lines towards NS
surface. By comparing the ram pressure of accreting material $P_{\rm
ram}\propto\rho v^{2}$ (where $\rho$ and $v$ are material mass density and
velocity) and the pressure of the NS magnetic field $P_{B}\propto B^{2}\propto
B^{2}_{0}\,r^{-6}$ (where $B$ is magnetic field strength at distance $r$ from
a NS and $B_{0}$ is field strength at the NS surface at the NS surface), we
can estimate the magnetospheric radius (see, e.g., Chapter 6.3 in Frank et al.
2002). In the case of a magnetic field dominated by the dipole component, the
magnetospheric radius can be estimated as

$\displaystyle R_{\rm m}\simeq 5\times 10^{7}\,\Lambda
B_{12}^{4/7}\dot{M}_{19}^{-2/7}m^{-1/7}R_{6}^{12/7}\,\,{\rm cm},$ (1)

where $\Lambda\leq 1$ is dimensionless coefficient, $B_{12}$ is the field
strength at the NS magnetic pole in units of $10^{12}\,{\rm G}$,
$\dot{M}_{19}$ is the mass accretion rate in units of $10^{19}\,{\rm
g\,s^{-1}}$, $m$ is the mass of a NS in units of Solar masses, and $R$ is NS
radius in units of $10^{6}\,{\rm cm}$. The coefficient $\Lambda$ is typically
taken to be $\Lambda=0.5$ for the case of accretion from the disc and
$\Lambda=1$ for the case of accretion from the wind (Ghosh & Lamb, 1978,
1979). However, in the case of accretion from the disc, it is dependent on the
mass accretion rate and physical conditions at the inner disc radius (see
Psaltis & Chakrabarty 1999; Chashkina et al. 2019 for details and discussion).
At the magnetospheric radius, the accretion flow penetrates into the NS
magnetosphere due to the development of instabilities, including reconnection,
magnetic Rayleigh-Taylor (Arons & Lea, 1976; Kulkarni & Romanova, 2008) and
Kelvin-Helmholtz instabilities (Burnard et al., 1983). Within the
magnetosphere of a NS, material moves toward the NS along magnetic field lines
and finally reaches the stellar surface in a small regions (area $\sim
10^{9}\,{\rm cm^{2}}$) near the magnetic poles of a NS. The kinetic energy of
the accretion flow turns into heat and is emitted in the form of X-rays by hot
spots or accretion columns, depending on the mass accretion rate and surface
magnetic field strength and structure (Basko & Sunyaev, 1976; Mushtukov et
al., 2015).

The NS magnetic field in XRPs is typically assumed to be dominated by the
dipole component. However, evidence of strong non-dipole magnetic field was
reported in a few objects including Her X-1 (Shakura et al., 1991; Mönkkönen
et al., 2022), bright Be X-ray transients SMC X-3 (Tsygankov et al., 2017) and
Swift J0243.6+6124 (Kong et al., 2022), and ULX pulsar NGC 5907 X-1 (Israel et
al., 2017a).

Typically, it is assumed that the accretion flow between the inner disc radius
and NS surface does not modify the X-ray flux and spectra in XRPs because the
magnetospheric flow is optically thin. Sunyaev 1976, however, noticed that
even at low mass accretion rates, the flow within the magnetosphere of a NS
can absorb and reprocess photons of low energies, which can result in
relatively low pulsed fraction observed in soft X-ray energy band. At high
mass accretion rates $\dot{M}\gtrsim 10^{18}\,{\rm g\,s^{-1}}$, the accretion
flow inside the magnetospheric radius is expected to be optically thick
(Mushtukov et al., 2017) due to the Compton scattering. Under this condition,
accretion should affect the processes of spectra formation, pulse profiles and
features of the aperiodic variability (Mushtukov et al., 2019; Koliopanos et
al., 2019; Brice et al., 2023). In particular, optically thick envelopes in
ULX pulsars can shield central NS from an observer, reprocess a large fraction
of X-ray flux and re-emit it in softer multi-color blackbody radiation
(Koliopanos et al., 2019; Brice et al., 2023). Brice et al. 2023 demonstrated
that appearance of optically thick envelopes in XRPs results in one-peak pulse
profile instead of pulse profile with two peaks of lower mass accretion rates.
The numerous absorptions and re-emissions of X-ray photons by the envelope
contributes to the typical escape time of photons from the system. Therefore,
one would expect suppression of X-ray aperiodic variability on time scales
smaller than the typical escape time (Mushtukov et al., 2019).

The dynamics of the flow covering NS magnetosphere is governed by the geometry
of magnetic field lines, gravitational force, centrifugal force and, in the
case of high mass accretion rates and luminosities, the radiative force. The
flow dynamics within the magnetospheric radius was already simulated
numerically by Romanova et al. 2003, 2004, but these simulations were limited
to the case of low mass accretion rates (where one can safely neglect the
influence of radiative force) and relatively small magnetospheres of a NS.
Self-consistent simulations of accretion flow dynamics in bright XRPs requires
calculations of radiative force, which was not performed up to date.

In this paper, we investigate the radiative force applied to the material
covering the magnetosphere of a NS, the role of multiple reflections of X-ray
photons, and the typical time scale of photon escape from the magnetosphere
covered by an optically thick accretion flow. To estimate the radiative force,
we use pre-calculated maps of material distribution over the magnetospheric
surface, and perform Monte Carlo simulations under the assumption that the
magnetic field geometry is dominated by the dipole component, while the
opacity is due to the scattering only.

##  2 Model set up

Figure 1:  The schematic illustration of dipole magnetic field lines geometry.
In general, the gravitational and the radiative forces applied to a given
point at the magnetospheric surface are not strictly in opposite directions
due to the possibility of multiple photon reprocessing by the accretion flow.

We consider a NS with mass $M=1.4M_{\odot}$ and radius $R=10^{6}\,{\rm cm}$ in
a binary system. We assume that the accretion flow forms a disc outside of the
magnetospheric radius (1). The magnetic and rotational axis of a NS are
assumed to be orthogonal to the disc plane. At a distance $R_{\rm m}$ from the
central object, the accretion flow is interrupted by the magnetic field of the
NS. Starting from the inner disc radius, the accretion flow follows magnetic
field lines and moves towards NS surface, forming an envelope that covers the
magnetosphere of the star. In the case of strongly magnetised NSs, $R_{\rm
m}\gg R$ and the majority of energy is released at the NS surface. The
luminosity related to the mass accretion rate is

$\displaystyle L=\frac{GM\dot{M}}{R}\simeq 1.3\times
10^{39}\dot{M}_{19}mR_{6}^{-1}\,{\rm erg\ \rm s^{-1}}.$ (2)

Due to the cylindrical symmetry, the coverage of the magnetosphere with matter
is expected to be independent of the azimuthal angle. The motion of material
along magnetic field lines starts in the equatorial plane with the velocity
close to the thermal velocity in the accretion disc. Estimating the inner disc
radius, we assume that the parameter $\Lambda$ in (1) is not dependent on the
mass accretion rate and equals $0.5$. However, we account for the dependence
of the accretion disc geometrical thickness at the inner radius on the mass
accretion rate. Estimating the geometrical thickness of a disc we follow
Suleimanov et al. 2007 and get

$\displaystyle H_{\rm d}\approx 2.7\times
10^{6}\,\alpha^{-1/10}\Lambda^{9/8}\,\dot{M}_{19}^{-6/35}B_{12}^{9/14}\,m^{-15/28}\,R_{6}^{27/14}\,{\rm
cm}$ (3)

for the case of relatively low mass accretion rates

$\dot{M}_{19}\lesssim 0.04\,\Lambda^{21/20}m^{-1/2}R_{6}^{9/5},$

when the disc is gas pressure dominated and the opacity is given by the
Kramers’ law. At higher mass accretion rates, the disc is assumed to be
radiation pressure dominated with the opacity determined by the electron
scattering. In this case, the geometrical thickness of the disc is expected to
be independent of the radial coordinate (Shakura & Sunyaev, 1973) and can be
estimated as

$\displaystyle H_{\rm d}\approx 1.3\times 10^{7}\,\dot{M}_{19}\,\mbox{cm}.$
(4)

We consider the magnetic field of a NS dominated by a dipole component (see
Fig. 1). The geometry of magnetic field lines in spherical coordinates is
described by

$\displaystyle r(\theta,\varphi)=R_{\rm m}\sin^{2}\theta=R_{\rm
m}\cos^{2}\lambda.$ (5)

The surface area of a segment on the dipole surface is given by

$\displaystyle{\rm d}S=R_{\rm
m}^{2}\cos^{4}\lambda(1+3\sin^{2}\lambda)^{1/2}{\rm d}\lambda{\rm d}\varphi,$
(6)

while the length along the field lines $x$ is related to the latitude:

$\displaystyle{\rm d}x=R_{\rm m}\cos\lambda(1+3\sin^{2}\lambda)^{1/2}{\rm
d}\lambda.$ (7)

The angle between the position vector of the point at the dipole surface and
the tangent line to the dipole magnetic field line (see Fig. 1) is

$\displaystyle\chi={\rm atan}\left(\frac{0.5}{{\rm tan}\lambda}\right).$ (8)

The magnetic field strength in the case dominant dipole component drops with
the distance as $B\propto r^{-3}$ and the field strength at the magnetosheric
radius can be estimated as

$\displaystyle B(r=R_{\rm m})\approx 3\times
10^{6}\,\Lambda^{-3}B_{12}^{-5/7}L_{39}^{6/7}m^{-3/7}R_{6}^{-9/7}\,\,{\rm G}.$
(9)

Note that the stronger the field strength at the NS surface, the smaller it is
at the magnetospheric radius.

The accretion flow moves across the magnetospheric surface strictly along
magnetic field lines under the influence of gravitational, centrifugal and the
radiative forces. The magnitude of acceleration due to the gravitational force
along the magnetic field lines is

$\displaystyle\left|a_{{\rm grav},||}\right|=\cos\chi\,\frac{GM}{r^{2}}\simeq
1.33\times 10^{10}\cos\chi\,\frac{m}{r_{8}^{2}}\,\,{\rm cm\,s^{-2}}.$ (10)

In the case of dipole magnetic field axis orthogonal to the disc plane, the
acceleration due to the centrifugal force is

$\displaystyle a_{{\rm cen},||}(\alpha=0)=\omega^{2}R_{\rm
m}\cos^{3}\lambda\cos(\chi-\lambda),$ (11)

where $\omega=2\pi/P_{\rm spin}$ is the angular velocity of a NS, and $\alpha$
is the angle between the magnetic and rotation axis of a NS.

The radiative force is dependent on the luminosity of the central object, the
beam pattern, the surface density of the material covering the magnetosphere,
the optical thickness, and details of the radiative transfer (i.e., the
dominant processes of radiation and matter interaction). We assume that the
opacity of matter covering the magnetospheric surface is due to the scattering
of photons by electrons (Compton scattering), and the opacity is

$\displaystyle\kappa_{e}=0.34\,{\rm cm^{2}\,g^{-1}}.$ (12)

Free-free absorption can be approximatelly described by the Kramers formula:

$\displaystyle\kappa_{\rm ff}\approx 0.0136\,\rho T_{\rm keV}^{-3.5}\,\,\,{\rm
cm^{2}\,g^{-1}},$ (13)

because the temperature of the magnetospheric accretion flow is expected to be
$T\sim 1\,{\rm keV}$ and the mass density in accretion channel $\rho\ll
1\,{\rm g\,cm^{-3}}$, $\kappa_{\rm ff}\ll\kappa_{e}$.

If one considers only scattering and neglects the absorption of photons in the
magnetospheric flow, each photon that reaches the magnetosphere either

  1. 1.

passes through the layer without interacting with it,

  2. 2.

passes through the layer after several scatterings, or

  3. 3.

is reflected back inside the magnetosphere, having scattered once or multiple
times.

The fraction of photons that pass through the layer without interacting with
the accretion flow is

$f=e^{-\tau_{e}/\cos\theta_{\rm in}},$

where $\theta_{\rm in}$ is the angle between the momentum of incident photon
and local normal to the magnetospheric surface. These photons do not transfer
any momentum to the flow and thus do not contribute to the radiative force.
The photons passing through the layer after one or more scatterings contribute
to the radiative force because the direction of their momentum changes after
each reflection. Moreover, the photons reflected back into the NS
magnetosphere contribute to the radiative force after each reflection. In the
case of high mass accretion rates and optically thick accretion columns,
photons can experience many reflections on average, and the reflection process
can heavily influence the radiative force and radiative pressure.

Figure 2:  Optical thickness of accretion flow covering magnetosphere of a NS
as a function of $\theta$-coordinate. Different curves are given for different
accretion luminosity: $10^{38}$ (lower), $10^{39}$ (middle), $10^{40}\,{\rm
erg\ \rm s^{-1}}$ (upper). Red solid curves are calculated for the case of
non-rotating magnetosphere, while the black dashed ones account for NS
rotation with spin period $P=1\,{\rm s}$. Magnetic field at the NS surface is
fixed at $B=3\times 10^{12}\,{\rm G}$. Calculation are performed according
Mushtukov et al. 2017, 2019.

In the case of azimuthal symmetry of magnetospheric coverage by the accreting
material, the local optical thickness of the accretion flow due to the
scattering can be roughly estimated as

$\displaystyle\tau_{e}\approx\frac{70L_{39}^{6/7}B_{12}^{2/7}}{\beta}\left(\frac{\cos\lambda_{0}}{\cos\lambda}\right)^{3},$
(14)

where $\beta=v/c$ is local dimensionless velocity (see Section 2.1 in
Mushtukov et al. 2017 and Fig. 2 illustrating more accurate calculations), and
$\lambda_{0}\approx{\pi}/{2}-(R/R_{\rm m})^{{1}/{2}}$ is the angle at the
accretion channel base at the NS surface. If the magnetospheric flow is
optically thick, the photons that originated from the central source
experience several reflections back into the NS magnetosphere before they
finally leave the system. The fraction of photons that penetrate through the
layer is dependent on the optical thickness $\tau_{e}$ of a layer and initial
angle between photon momentum and local normal (see Fig. 3). For the case of
an optically thick layer (i.e., $\tau_{e}\gtrsim 1$), where the radiative
transfer is dominated by the isotropic scattering, and photons incident
perpendicular to the surface (i.e., $\theta_{\rm in}=0$), the fraction of
particles penetrating through the layer can be approximated by (see Appendix C
in Mushtukov et al. 2019 and Nagirner et al. 2022 for detailed analytical
discussion) 111Note, that there were a typo in this approximate expression in
Mushtukov et al. 2019.

$\displaystyle f_{\rm out}\simeq\frac{\sqrt{3}}{\tau_{e}+1.42}.$ (15)

For larger $\theta_{\rm in}$, the fraction of photons penetrating through the
layer is smaller.

The construction of a self-consistent model of magnetospheric accretion falls
outside the scope of this paper. To assess the impact of multiple
reprocessings of X-ray photons in the magnetosphere of a NS, we utilize
precalculated maps of matter distribution over the magnetospheric surface.
These maps are based on the models proposed by Mushtukov et al. 2017, 2019
(examples of calculations are illustrated in Fig.,2).

Figure 3:  The fraction of photons that penetrate through the layer of optical
thickness $\tau_{e}$. Red lines show the fractions calculated for different
initial angles between the photon momentum and normal to the layer:
$\theta_{\rm in}=0\degr$, $60\degr$ and $86\degr$ (top to bottom, see Appendix
A). Black dots show approximation calculated according to (15).

##  3 Numerical model

To calculate the radiative force and momentum exchange between the accretion
flow and photons, we perform Monte Carlo simulations of momentum transfer to
the accretion envelope. Each simulation is composed of a few steps:

  1. 1.

We start a photon from the central source. The initial direction of a photon
is described by two angles - $\theta_{i}$ and $\varphi_{i}$. The angles are
determined by random numbers $X_{1},X_{2}\in[0;1]$. In the case of central
source producing isotropic X-ray emission, the initial direction of photon
momentum is given by

$\displaystyle\theta_{i}^{\rm(iso)}$ $\displaystyle=$ $\displaystyle{\rm
acos}(1-2X_{1})\in\left[0;\pi\right],$ (16)
$\displaystyle\varphi_{i}^{\rm(iso)}$ $\displaystyle=$ $\displaystyle 2\pi
X_{2}\in[0;2\pi].$ (17)

Each photon in this case is assigned the same statistical weight $w_{i}=1$. By
changing the statistical weight $w_{i}$ of a photon depending on the initial
direction of its motion we can take into account the possible non-isotropy of
the radiation from the central NS. In the particular case of

$\displaystyle\frac{{\rm d}L}{{\rm d}\Omega_{i}}=g(\theta_{i},\varphi_{i})L,$
(18)

where $\Omega_{i}$ is the solid angle, $g$ is a function describing initial
angular distribution of X-ray photons taken so that $\int_{0}^{\pi}{\rm
d}\theta_{i}\int_{0}^{2\pi}{\rm
d}\varphi_{i}g(\theta_{i},\varphi_{i})\sin\theta_{i}=1$, the photons’
statistical weight is taken as

$\displaystyle w_{i}=g(\theta_{i},\varphi_{i}).$ (19)

Each photon participating in the simulation carries effective momentum
$p_{i}\propto w_{i}$, or more precisely

$\displaystyle p_{i}=\frac{w_{i}}{N_{\rm tot}}\frac{L}{c},$ (20)

where $N_{\rm tot}$ is the total number of photons used in a simulation.

  2. 2.

The photon started from a given coordinate in a certain direction is traced as
it travels to the magnetospheric surface, under the assumption that it is
propagating along a straight-line trajectory. We determine the coordinates of
the point at the magnetospheric surface at which the photon eventually
reaches.

  3. 3.

Using (14), we estimate the local optical thickness $\tau_{\rm e}$ of the
accretion flow at the point where the photon reaches the magnetosphere. We use
the local optical thickness and angle between the normal to the magnetospheric
surface and photon momentum to calculate the probability $f_{\rm out}$ for the
photon to penetrate through the envelope (see approximation 15 and Fig. 3). We
generate a random number $X_{3}$ between 0 and 1 and compare it to $f_{\rm
out}$ in order to determine whether the photon will penetrate through the
envelope or reflect back into the dipole magnetosphere. If the photon
penetrates through the flow, we stop tracing its trajectory and proceed to
step 4. If the photon is reflected back into the NS magnetosphere, we
determine the new direction of its motion, account for the local momentum
transfer to the material along and across magnetic field lines, and continue
tracing its trajectory. Each reflection is associated with momentum transfer
between the accretion flow and X-ray photons. The momentum effectively
transferred to the accretion flow due to a reflection is

$\displaystyle\Delta\textbf{{p}}=\textbf{{p}}_{\rm ini}-\textbf{{p}}_{\rm
fin},$ (21)

where $\textbf{{p}}_{\rm ini}$ and $\textbf{{p}}_{\rm fin}$ are vectors of the
photon’s momentum before and after the reflection respectively, and $w_{i}$ is
a photon effective weight (19). The new direction of photon motion after
reflection is calculated under the assumption of multiple conservative,
isotropic and coherent scattering in the layer of a given optical thickness.
We return to step 2.

  4. 4.

We account for the photon time travel before its escape in our estimations of
the average photon escape time.

The momentum transferred from photons to the accretion flow determines the
radiative force applied to the flow and the corresponding acceleration along
magnetic field lines:

$\displaystyle a_{\rm rad,||}\geq a^{*}_{\rm
rad,||}=\left(\frac{\sum_{j}\Delta p_{j,||}}{N_{\rm
tot}}\right)\frac{L}{c}\frac{1}{\Sigma\Delta S},$ (22)

where $\Sigma$ is the local surface density of the accretion flow, and $\Delta
S$ is the surface area of the segment. The projection of the transferred
momentum onto the direction of the local magnetic field lines is given by the
scalar production

$\displaystyle\Delta p_{j,||}=\Delta\textbf{{p}}_{j}\cdot\textbf{{n}}_{b}$
(23)

where

$\displaystyle\textbf{{n}}_{b}=\left(\begin{array}[]{c}-\cos(\chi-\lambda)\cos\varphi\\\
-\cos(\chi-\lambda)\sin\varphi\\\ \sin(\chi-\lambda)\end{array}\right)$ (27)

is the unit vector tangential to the magnetic field lines.
Note that the estimation (22) assumes that the momentum transfer is
distributed homogeneously over the geometrical thickness of the accretion
flow, which is not necessarily the case.

The radiative pressure on the magnetosphere is related to the photon momentum
transfer across magnetic field lines:

$\displaystyle P_{\rm rad,\perp}=\left(\frac{\sum_{j}\Delta
p_{j,\perp}}{N_{\rm tot}}\right)\frac{L}{c}\frac{1}{\Delta S},$ (28)

where

$\displaystyle\Delta p_{j,\perp}=\Delta\textbf{{p}}_{j}\cdot\textbf{{n}}$ (29)

is the projection of the transferred momentum on the local normal to the
magnetospheric surface

$\displaystyle\textbf{{n}}=\left(\begin{array}[]{c}{\rm
sgn}\lambda\,\sin(\lambda-\chi)\cos\varphi\\\ {\rm
sgn}\lambda\,\sin(\lambda-\chi)\sin\varphi\\\ -{\rm
sgn}\lambda\,\cos(\lambda-\chi)\end{array}\right).$ (33)

##  4 Numerical results

In our simulations, we explore the luminosity range of $10^{38}\,{\rm erg\ \rm
s^{-1}}\lesssim L\lesssim 10^{40}\,{\rm erg\ \rm s^{-1}}$. The lower limit
corresponds to the luminosity level where one would expect an optically thick
magnetospheric accretion flow (Mushtukov et al., 2017). Through Monte Carlo
simulations, we investigate how the mass accretion rate and the multiple
reprocessings of X-ray photons influence acceleration along magnetic field
lines (Section 4.2) and pressure across field lines (Section 4.3) due to the
radiative force. By simulating multiple reprocessings of photons inside the
magnetosphere, we can estimate the average photon escape time from the
envelope (Section 4.4).

###  4.1 Radiative transfer through a layer

Each photon interaction with the magnetospheric surface is associated with
either a penetration through the layer or a reflection back into the NS
magnetosphere. The fraction of photons that penetrate through the layer,
$f_{\rm out}$, is dependent on the optical thickness $\tau_{\rm e}$ and the
incident angle of a photons (see Fig. 3). The larger the thickness and the
angle between the local normal and initial photon momentum, the smaller the
fraction of photons that penetrate through the layer. For the case of photons
that initially move along the local normal, the fraction can be well
approximated by (15) if $\tau_{\rm e}\gg 1$ (compare black dotted line and the
upper red solid line in Fig. 3). It is expected that the most of momentum
along the field lines is transferred to the accretion flow within the layer of
optical thickness $\sim 1$. After the first scattering, the photons “forget"
the initial direction of their motion and contribute to the photon energy flux
across magnetic field lines.

###  4.2 Radiative force and acceleration along the field lines

The radiative force acting along magnetic field lines is proportional to the
momentum transferred along the field lines (see 22 and 23). According to our
simulations, the radiative acceleration along magnetic field lines is always
in the opposite direction of the projection of the gravitational acceleration
along magnetic field lines. A larger mass accretion rate and luminosity
corresponds with a larger radiative force and acceleration along the field
lines (see Fig. 4a). This occurs because there is a larger energy release at
the NS surface, and a greater optical thickness of the accretion flow results
in a larger number of photon reflections (and hence, more momentum is
transferred to the accreting material) before photons escape the system. To
estimate the contribution of the radiative force to the acceleration, we
assume that the momentum of reprocessed photons is distributed homogeneously
over the local optical thickness of the accretion flow (i.e. there is no
gradient of radiative acceleration along $\tau_{\rm e}$). Multiple reflections
of X-ray photons affects the acceleration, especially in regions located close
to the disc plane, and could increase the ratio $a^{*}_{\rm rad}/a_{\rm grav}$
by nearly an order of magnitude (see Fig. 4b). It is expected that in the case
of $a^{*}_{\rm rad}/a_{\rm grav}>1$, the radiative force will begin to
decelerate the accretion flow. Remarkably, the distribution of the
acceleration due to the radiative force becomes less sensitive to the specific
beam pattern at higher mass accretion rates and luminosity. It happens because
at high mass accretion rates and larger optical thickness of accretion flow
the radiative force is largely due to the photons that have already
experienced one or a few reprocessing.

Figure 4:  The ratio of the radiative to gravitational accelerations along the
field lines as a function of coordinate $\theta_{B}$. (a) Different lines
correspond to different luminosities: $L=10^{38},\,10^{39},\,10^{40}\,{\rm
erg\ \rm s^{-1}}$ (blue dashed-dotted, solid red, and black dashed lines
respectively). (b) The accretion luminosity is fixed at $10^{40}\,{\rm erg\
\rm s^{-1}}$. Black (red) lines correspond to the surface field strength of
$B=10^{14}\,{\rm G}$ ($B=3\times 10^{12}\,{\rm G}$). The solid lines account
for multiple reflection of X-ray photons while the dashed ones account for the
first reflection only. The dips at $\theta_{B}\sim 90^{\circ}$ appear due to
the accretion disc that makes the effective surface density extremely large
(see expression 22). One can see that under the condition of sufficiently high
luminosity, the acceleration along the field lines due to the radiative force
can exceed the acceleration due to the gravitational force. Note that the
ratio is considered with a minus sign, and positive values indicate scenarios
where the acceleration resulting from the radiative force opposes the
acceleration due to gravity along the field lines.

###  4.3 Radiative pressure across the field lines

Reflections of X-ray photons lead to momentum transfer across magnetic field
lines, resulting in the anticipated radiative pressure. The radiative pressure
across the field lines increases with luminosity and may surpass the local
pressure of the magnetic field, denoted as $P_{B}=B^{2}/(8\pi)$, particularly
in regions near the accretion disc plane, where the magnetic field is weaker
(refer to Fig. 5a). When $P_{\rm rad}\gtrsim P_{\rm B}$, the geometry of the
magnetic field lines can be influenced by the radiation pressure. Note that
the local sharp maxima around $\theta_{\rm d}\approx 90\degr$ in the simulated
curves, corresponding to $L=10^{38}\,{\rm erg\ \rm s^{-1}}$ and $10^{39}\,{\rm
erg\ \rm s^{-1}}$ in Fig. 5a, arise from the presence of the accretion disc,
which reprocesses all photons reaching it.

The radiative pressure is significantly impacted by multiple reflections of
X-ray photons, causing the pressure to increase by more than an order of
magnitude (see Fig. 5b). Similarly to the radiate acceleration, the pressure
becomes less sensitive to the specific beam pattern of higher mass accretion
rates.

Figure 5:  The ratio of the radiation pressure $P_{\rm rad}$ and the magnetic
field pressure $P_{B}$ at the magnetospheric surface as a function of
$\theta_{B}$ coordinate. (a) Different lines correspond to different accretion
luminosities: $L=10^{38},\,10^{39},\,10^{40}\,{\rm erg\ \rm s^{-1}}$ (from
bottom to top). One can see that under condition of sufficiently high
accretion luminosity, the radiative pressure can be larger than the local
pressure of magnetic field (the field strength is fixed and is taken to be
$B=3\times 10^{12}\,{\rm G}$ at the NS surface). (b) Solid lines show the
result accounting for multiple reflections, while the dashed lines illustrate
the case when only the first reflection is taken into account. Red and black
lines correspond to the case of surface magnetic field strength $3\times
10^{12}\,{\rm G}$ and $B=10^{14}\,{\rm G}$ respectively, luminosity is fixed
at $L=10^{40}\,{\rm erg\ \rm s^{-1}}$.

###  4.4 Photon escape time

Photon escape time is affected by the mass accretion rate, as it determines
the optical thickness of the flow and the geometrical size of the NS
magnetosphere. The larger the mass accretion rate and luminosity, the longer
the average time photons will travel within the magnetosphere (see Fig. 6).
The dependence of the average escape time on the mass accretion rate is
approximately linear. Note, however, that according to Chashkina et al. 2019,
the inner disc radius stabilizes and is weakly dependent on the mass accretion
rate when it is advective or dominated by the radiation pressure (see Fig. 26
in Mushtukov & Tsygankov 2022). If it is indeed the case, the dependence of
the photon travel time on the mass accretion rate should be stronger. The
absolute value of the escape time tends to be longer in the case of a stronger
magnetic field, due to a larger magnetosphere of a NS (see Eq. 1). The
dependence of the escape time on the NS field strength is not strong because
(i) the magnetospheric radius depend weakly on the field trength: $R_{\rm
m}\propto B^{4/7}$, and (ii) the larger the magnetosphere, the smaller the
optical thickness of the magnetospheric flow. At an accretion luminosity of
$\sim 10^{40}\,{\rm erg\ \rm s^{-1}}$, the photon escape time is $\sim
0.1\,{\rm s}$ (see Fig. 6a).

The characteristic time of photon escape from the magnetosphere determines the
time scale on which the radiation force and pressure are established. If the
typical timescale of the accretion flow variations are smaller than (or about)
the characteristic time of photon escape, the radiation field inside the
magnetosphere is determined not only by the current configuration of the
accretion flow, but also by the history of its variations. At sufficiently
high mass accretion rates, the escape time of X-ray photons can be longer than
the typical time of free-fall from the inner radius of accretion disc (see
Fig. 6b):

$\displaystyle t_{\rm ff}\simeq 0.2\,m^{-1/2}\left(\frac{R_{\rm
m}}{10^{8}\,{\rm cm}}\right)^{3/2}\,{\rm s},$ (34)

where $R_{\rm m}$ is given by (1).

It is expected that the escape time of photons from the magnetosphere affect
features of aperiodic variability in XRPs. The variability on timescales
$t<t_{\rm esc}$ is expected to be suppressed (see, e.g., Mushtukov et al.
2019).

Figure 6:  The typical escape time of X-ray photons from the magnetosphere of
a NS covered by optically thick envelope in units of seconds (a), and in units
of free-fall time from the inner radius of accretion disc (b). Black and red
lines illustrate the results for a surface magnetic field strength of
$B=3\times 10^{12}\,{\rm G}$ and $10^{14}\,{\rm G}$, respectively. The breaks
in both curves are due to accretion disc transition from the geometrically
thin to the geometrically thick state at the magnetospheric radius.

##  5 Summary and discussion

We have analyzed interactions of photons with the magnetospheric accretion
flow in bright XRPs ($L\gtrsim 10^{38}\,{\rm erg\ \rm s^{-1}}$) including ULXs
powered by accretion onto a NS. In these objects, the accretion flow covering
NS magnetosphere is expected to be optically thick due to the photons
scattering by electrons and influence spectra and pulse profile formation
(Mushtukov et al., 2017; Brice et al., 2023). Using precalculated maps of
material distribution over the magnetospheric surface (Mushtukov et al., 2019)
and applying Monte Carlo simulations, we have investigated the distribution of
radiative force along magnetic field lines and radiation pressure across the
magnetic field lines over the magnetospheric surface. We have estimated the
typical time scale of photon escape from the magnetosphere of a NS at
different mass accretion rates. Our conclusions are:

  * •

Multiple reflections of X-ray photons increase (sometimes by more than an
order of magnitude) both the radiative force along magnetic field lines (see
Fig. 4b) and the radiation pressure across the field lines (see Fig. 5b).

  * •

At sufficiently high luminosities, the magnitude of acceleration along the
magnetic field lines due to the radiative force can exceed the acceleration
along the field lines due to the gravitational force, particularly in regions
near the equatorial plane of the magnetic dipole (see Fig. 4a). In these
regions one would expect significant deceleration of accretion flow by the
radiative force.

  * •

At an accretion luminosity of $L\gtrsim 5\times 10^{39}\,{\rm erg\ \rm
s^{-1}}$, the radiative pressure across magnetic field lines can exceed the
local magnetic field pressure (see Fig. 5a). The ratio $P_{\rm rad}/P_{B}$
tends to be larger in the regions located close to the equator of the magnetic
dipole. In the case of $P_{\rm rad}/P_{B}>1$ one would expect local disruption
of the accretion flow and possibly radiation driven mass losses.

  * •

We have demonstrated that increasing the mass accretion rate and luminosity
results in a larger average escape time for photons inside the magnetosphere
(see Fig. 6). The average time of photon escape depends on the mass accretion
rate and size of NS magnetosphere (i.e., on the NS magnetic field strength and
structure). The dependence of the photon escape time on the mass accretion
rate is approximately linear, while the dependence on the stellar field
strength is weak. At accretion luminosity $L\approx 10^{40}\,{\rm erg\ \rm
s^{-1}}$, the avarage photon escape time $\sim 0.1\,{\rm s}$. The typical time
scale of photon escape determines the time scale of radiative force and
pressure establishing, which has to be taken into account in calculations of
magnetosphere accretion flow dynamics.

The significant radiative force along magnetic field lines, not previously
considered in earlier publications (Mushtukov et al., 2017, 2019; Brice et
al., 2023), can profoundly impact the dynamics of the accretion flow within
the magnetospheric radius. In this scenario, the analysis of material
distribution over the magnetospheric surface, its influence on spectra
formation, and even the stability of the accretion flow must account for the
effects of the radiative force. As the radiative acceleration along magnetic
field lines can surpass the acceleration due to gravitational force, the
accretion flow towards the neutron star surface may be impeded by radiation.
This, in turn, results in a decrease in accretion luminosity, a reduction in
radiative force, and the subsequent resumption of the accretion process.
Consequently, the appearance of oscillations in mass accretion rate within the
magnetosphere is anticipated (Mushtukov et al., 2024).

In line with previous findings (Mushtukov et al., 2017), the radiation
pressure on the accretion flow can surpass the local magnetic field pressure.
The radiation pressure has the potential to distort magnetic field lines,
thereby influencing the geometry and dynamics of the magnetospheric accretion
flow.

The extended timescale for photon escape from the envelope at high mass
accretion rates and luminosities (see Fig. 6) is anticipated to influence the
timing properties of X-ray aperiodic variability. Strong suppression is
expected in the power density spectra on timescales shorter than the average
time of photon escape (see Mushtukov et al. 2019 for detailed discussions).
Given that the typical free-fall time from the inner disc radius in X-ray
pulsars is approximately $\sim 0.1\,{\rm s}$ (comparable to the timescale of
photon escape during high luminosity states) and the accretion flow can be
rendered unstable due to the influence of radiative force and pressure
(Mushtukov et al., 2024), the structure of the magnetospheric flow and
radiative field becomes coupled. Under these conditions, their temporal
evolution cannot be calculated separately in numerical models.

## Acknowledgements

AAM thanks UKRI Stephen Hawking fellowship. We are grateful to Simon Portegies
Zwart for discussions and our anonymous referee for useful comments and
suggestions which helped us improve the paper.

## Data availability

The calculations presented in this paper were performed using a private code
developed and owned by the corresponding author. All the data appearing in the
figures are available upon request.

## References

  * Arons & Lea (1976) Arons J., Lea S. M., 1976, ApJ, 207, 914
  * Bachetti et al. (2014) Bachetti M., et al., 2014, Nature, 514, 202
  * Basko & Sunyaev (1976) Basko M. M., Sunyaev R. A., 1976, MNRAS, 175, 395
  * Brice et al. (2023) Brice N., Zane S., Taverna R., Turolla R., Wu K., 2023, MNRAS, 525, 4176
  * Burnard et al. (1983) Burnard D. J., Arons J., Lea S. M., 1983, ApJ, 266, 175
  * Carpano et al. (2018) Carpano S., Haberl F., Maitra C., Vasilopoulos G., 2018, MNRAS, 476, L45
  * Chashkina et al. (2019) Chashkina A., Lipunova G., Abolmasov P., Poutanen J., 2019, A&A, 626, A18
  * Fabrika et al. (2021) Fabrika S. N., Atapin K. E., Vinokurov A. S., Sholukhova O. N., 2021, Astrophysical Bulletin, 76, 6
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: Third Edition
  * Fürst et al. (2016) Fürst F., et al., 2016, ApJ, 831, L14
  * Ghosh & Lamb (1978) Ghosh P., Lamb F. K., 1978, ApJ, 223, L83
  * Ghosh & Lamb (1979) Ghosh P., Lamb F. K., 1979, ApJ, 234, 296
  * Illarionov & Sunyaev (1975) Illarionov A. F., Sunyaev R. A., 1975, A&A, 39, 185
  * Israel et al. (2017a) Israel G. L., et al., 2017a, Science, 355, 817
  * Israel et al. (2017b) Israel G. L., et al., 2017b, MNRAS, 466, L48
  * King et al. (2023) King A., Lasota J.-P., Middleton M., 2023, New Astron. Rev., 96, 101672
  * Koliopanos et al. (2019) Koliopanos F., Vasilopoulos G., Buchner J., Maitra C., Haberl F., 2019, A&A, 621, A118
  * Kong et al. (2022) Kong L.-D., et al., 2022, ApJ, 933, L3
  * Kulkarni & Romanova (2008) Kulkarni A. K., Romanova M. M., 2008, MNRAS, 386, 673
  * Mönkkönen et al. (2022) Mönkkönen J., Tsygankov S. S., Mushtukov A. A., Doroshenko V., Suleimanov V. F., Poutanen J., 2022, MNRAS, 515, 571
  * Mushtukov & Tsygankov (2022) Mushtukov A., Tsygankov S., 2022, arXiv e-prints, p. arXiv:2204.14185
  * Mushtukov et al. (2015) Mushtukov A. A., Suleimanov V. F., Tsygankov S. S., Poutanen J., 2015, MNRAS, 447, 1847
  * Mushtukov et al. (2017) Mushtukov A. A., Suleimanov V. F., Tsygankov S. S., Ingram A., 2017, MNRAS, 467, 1202
  * Mushtukov et al. (2019) Mushtukov A. A., Ingram A., Middleton M., Nagirner D. I., van der Klis M., 2019, MNRAS, 484, 687
  * Mushtukov et al. (2024) Mushtukov A. A., Ingram A., Suleimanov V. F., DiLullo N., Middleton M., Tsygankov S. S., van der Klis M., Portegies Zwart S., 2024, arXiv e-prints, p. arXiv:2402.12965
  * Nagirner et al. (2022) Nagirner D. I., Milanova Y. V., Dementyev A. V., Volkov E. V., 2022, Astrophysics, 65, 560
  * Psaltis & Chakrabarty (1999) Psaltis D., Chakrabarty D., 1999, ApJ, 521, 332
  * Reig (2011) Reig P., 2011, Ap&SS, 332, 1
  * Rodríguez Castillo et al. (2020) Rodríguez Castillo G. A., et al., 2020, ApJ, 895, 60
  * Romanova et al. (2003) Romanova M. M., Ustyugova G. V., Koldoba A. V., Wick J. V., Lovelace R. V. E., 2003, ApJ, 595, 1009
  * Romanova et al. (2004) Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E., 2004, ApJ, 610, 920
  * Sathyaprakash et al. (2019) Sathyaprakash R., et al., 2019, MNRAS, 488, L35
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shakura et al. (1991) Shakura N. I., Postnov K. A., Prokhorov M. E., 1991, Pisma v Astronomicheskii Zhurnal, 17, 803
  * Staubert et al. (2019) Staubert R., et al., 2019, A&A, 622, A61
  * Suleimanov et al. (2007) Suleimanov V. F., Lipunova G. V., Shakura N. I., 2007, Astronomy Reports, 51, 549
  * Sunyaev (1976) Sunyaev R. A., 1976, Soviet Astronomy Letters, 2, 111
  * Tsygankov et al. (2016) Tsygankov S. S., Lutovinov A. A., Doroshenko V., Mushtukov A. A., Suleimanov V., Poutanen J., 2016, A&A, 593, A16
  * Tsygankov et al. (2017) Tsygankov S. S., Doroshenko V., Lutovinov A. A., Mushtukov A. A., Poutanen J., 2017, A&A, 605, A39

##  Appendix A ON THE FRACTION OF PHOTONS PENETRATING THROUGH THE LAYER

In the case of a layer illuminated by photons, the photons can either
penetrate through, be absorbed, or be reflected by the layer. The
probabilities of these different scenarios depend on the optical thickness of
the layer, as well as the nature of opacity resulting from the scattering and
absorption processes.

Consider a plane-parallel layer with a given optical thickness, denoted as
$\tau_{e}$. Assume that isotropic scattering, without photon absorption (i.e.,
conservative scattering), is the sole process that governs photon interaction
within the layer. Under this condition, photons are either reflected or
penetrate through the layer, and the total number of photons is conserved. The
fraction of photons that penetrate through the layer depends on both
$\tau_{e}$ and the initial direction of the photon momentum relative to the
normal to the layer, represented by the angle $\theta_{\rm in}$. For the
specific case where $\theta_{\rm in}=0$, there exists an approximate
expression (15) for the fraction of photons capable of penetrating through the
layer. For arbitrary angles $\theta_{\rm in}$, the fraction of photons that
penetrate through the layer can be determined numerically. To calculate the
fraction of photons that penetrate through the layer, we conduct Monte Carlo
simulations using the following algorithm, which can be divided into several
steps:

  1. 1.

Initiate a photon with a given initial direction of motion, determined by the
angle between the momentum and the normal, denoted as $\theta$, from a
boundary of a layer (i.e., at optical depth $\tau_{\rm ini}=0$).

  2. 2.

Generate a random number $X_{1}$ to determine the free path of a photon inside
a layer in units of mean free path:

$\displaystyle\Delta\tau=-\ln X_{1}.$ (35)

  3. 3.

Use the obtained free path of the photon to determine the optical depth of the
next scattering event:

$\displaystyle\tau_{\rm fin}=\tau_{\rm ini}+\Delta\tau\,\cos\theta.$ (36)

If $\tau_{\rm fin}<0$ ($>\tau_{e}$), the photon is reflected (penetrates
through) the layer. In this case, we account for the photon and return to step
1. If $0<\tau_{\rm fin}<\tau_{e}$, the photon undergoes the next scattering
inside the layer, and proceed to step 4.

  4. 4.

Set the optical depth of reemitted photons as $\tau_{\rm ini}=\tau_{\rm fin}$.
Re-emit the photon isotropically, and calculate a new direction of photon
motion:

$\displaystyle\theta={\rm arccos}(1-2X_{2}),$ (37)

where $X_{2}$ is a random number. Return to step 2.

The fraction of photons that penetrate through a layer is

$\displaystyle f_{\rm out}\approx\frac{N_{\rm out}}{N_{\rm tot}},$ (38)

where $N_{\rm tot}$ is a total number of photons in a simulation, and $N_{\rm
out}$ is a number of photons penetrated through a layer. Results of numerical
simulations are shown in Fig. 3. The larger the optical thickness and angle
between the normal and photon momentum, the smaller the fraction of photons
penetrating through the layer. At $\tau_{e}\gg 1$, $f_{\rm
out}\propto\tau_{e}^{-1}$. For $\tau_{e}>50$ and $\theta<1.5\,{\rm rad}$, one
can use approximation

$\displaystyle f_{\rm out}(\tau_{e},\theta)\approx f_{\rm
out}(\tau_{e},\theta=0)-\frac{\theta^{2}}{2.2\tau_{e}}.$ (39)


# Deep Nonparametric Convexified Filtering for Computational Photography,
Image Synthesis and Adversarial Defense

Jianqiao Wangni
Zhipu AI
<EMAIL_ADDRESS>
###### Abstract

We aim to provide a general framework of for computational photography that
recovers the real scene from imperfect images, via the Deep Nonparametric
Convexified Filtering (DNCF). It is consists of a nonparametric deep network
to resemble the physical equations behind the image formation, such as
denoising, super-resolution, inpainting, and flash. DNCF has no
parameterization dependent on training data, therefore has a strong
generalization and robustness to adversarial image manipulation. During
inference, we also encourage the network parameters to be nonnegative and
create a bi-convex function on the input and parameters, and this adapts to
second-order optimization algorithms with insufficient running time, having
10X acceleration over DIP. With these tools, we empirically verify its
capability to defend image classification deep networks against adversary
attack algorithms in real-time.

##  1 Introduction

Computational photography aims to recover real scenes from imperfect images
captured by cameras. Understanding the physics behind the image formation,
like gain control, aperture, exposure time and depth of focus, etc, are
fundamental works in the area. Besides the physics principles that are clear
and universal, lots of factors like object depth, are treated as a random
variable. For computational photography, we usually solve an inverse problem
with a statistical prior of these factors. The characterization of the prior
is thus an important task and contributes significantly to the selection of a
suitable algorithm. Take image denoising, for example, the Wiener filter works
better for Gaussian white noise, and the median filter fits the salt-and-
pepper noise [6]. Besides, any estimation or pre-calibration of the noise
level function [20] [15] helps with denoising. Meanwhile, some other imaging
factors are more difficult to estimate from limited images. Take the single
image blurring task as an example; there is no proper prior knowledge of the
essential variables like object depth, which directly affects the point spread
function (PSF) of out-of-focus blur kernel; object and camera movement are
also unknown. There were some works in the line of deep learning that learn a
statistical prior from massive training data for photography tasks, e.g,
super-resolution [9] , image dehazing [3], deblurring [26] and denoising
[32][30]. However, rigorously, such approaches have to rely on the assumption
that testing images resemble training images to ensure generalization.

Figure 1:  Image denoising. Left to right: clean, noisy, TV [5], DIP (92.7
sec), DeepRED[19] +TV (68.1 sec), DNCF (11.9 sec).

We view the aforementioned computational photography algorithms as
nonparametric, that their parameters are specific to each image, and do not
depend on training data. On a very different track of applications, where the
images are used in semantic tasks like classification with deep neural
networks, the researchers found the networks are extremely easy to be fooled,
and the images seem to be deceptive in the eyes of deep learning [2] [25][14]
[17][4]. Without precaution, an adversary party can manually create images to
make the deep network having incorrect predictions with high-confidence, even
those images were only manipulated in details that human eyes are unable to
notice. We understand this from the intrinsic nature of the deep networks,
which are mostly over-parameterized, and most math operations are
differentiable. It is easier to get better performance for deep learning, as
these parameters are practical to optimize. This can be a double-edged sword,
since if one party can easily train the parameters of the network, then the
adversary party can easily manipulate the images to the wrong side. This
inspires us to think that any deep network that has no specific parameters
will leave no chances of attack to the adversary, by being exactly the
opposite of conventional deep learning approaches.

A recent method named Deep Image Prior (DIP)[27][7][19][13] addresses the
problem by describing an image prior implicitly by the network structure, and
it no longer needs any training. DIP assumes that each image has an underlying
parameterization, which is specific to the individual image itself. The target
image $I_{t}$ is synthesized by a neural network $f$ with $\theta$ as its
parameters, and it is generated from a random vector variable $z$.

$\displaystyle I_{t}=f_{\theta^{\star}}(z).$ (1)

The network may consist of convolution, downsampling, deconvolution, and
upsampling layers specifically for different tasks. Simply defining the
structure of $f$, the target images $I_{t}$ is reconstructed through gradient-
based optimization w.r.t. $\theta$, based on the observed source image
$I_{s}$,

$\displaystyle\theta^{\star}=\arg\min_{\theta}\mathcal{L}\left(I_{s},f_{\theta}(z)\right).$
(2)

where the loss function $\mathcal{L}$ can be square $\ell_{2}$ distance or
total variation. The interpretation behind the formulation resembles that the
target image should be an optimal point of a regularized energy function
$E(I;I_{s})+\mathcal{S}(I)$, the regularization $\mathcal{S}(I)=0$ for
$\exists z,I=f_{\theta}(z)$ and $\mathcal{S}(I)=+\infty$ otherwise, which
states that $I$ is constrained by the network structure, and the prior is
defined by its expressiveness.

We study a better network structure of $f$ that better fits the physical
process and a faster inference pipeline. Start with considering the simplicity
of the loss function driven by physics, e.g.; we observe a deblurring model
that adopts a denoising feature as

$\displaystyle\widehat{I}_{t},\widehat{\mathcal{K}}=\min_{I_{t},\mathcal{K}}\|I_{t}\ast\mathcal{K}-I_{s}\|^{2},\quad$
(3)

where $\mathcal{K}$ is the unknown point spread function (PSF). This blind
deconvolution model is though nonconvex w.r.t the unknown variables
$\\{I_{t},\mathcal{K}\\}$, once the PSF is known, the loss function is convex
w.r.t the target image. The neural networks are typically composed of many
layers, being inherently highly nonlinear and nonconvex w.r.t to the network
parameter and inputs, and may consist of up to millions of parameters, so they
lose the simplicity of physical equations, many of which are linear or convex.
In this paper, we try to explore a framework that are as expressive as deep
networks and also being explainable and nonparametric as blind deconvolution.

Figure 2:  Visualization of Super-resolution with limited running time. Left
to right: high resolution, low resolution, DIP (353 sec), DeepRED + TV (112.1
sec), DNCF (12.0 sec).

##  2 Approach

This section presents our nonparametric deep network approach to solve
computational photography. Our first consideration is to have a decomposed
neural network of two meaningful parts: a physical interference model $f$
sequentially connects to the generative model $g$, thus
$I_{s}=f_{\theta}(g_{\theta}(z))$. Here $f$ and $g$ have different set of
network parameters, but for convenience we use $\theta$ to represent a
concatenation of both network, which also indicates the two networks could be
jointly trained. Ideally, the intermediate feature $g_{\theta}(z)$ simulates
the real scene $I_{t}$ and $f_{\theta}$ simulate the physical interference
such as blur and noise to approximate the observed image $I_{s}$, though
$I_{t}$ may not directly exist in any layer of $f$. The first initiative is
the following objective with a clear definition of the convolution kernel
$\mathcal{K}$:

$\displaystyle
I_{t},\theta^{\star}=\arg\min_{y,\theta}\|y\ast\mathcal{K}-I_{s}\|+\beta\|y-g_{\theta}(z)\|.$
(4)

We wish to reserve the highly strong expressivity of $g$ to simulate a complex
enough scene and simultaneously make $f$ towards a simple function to mimic a
physical corruption process. To further make the physical interference more
complex, we substitute the convolution by $f$ and obtain:

$\displaystyle
I_{t},\theta^{\star}=\arg\min_{y,\theta}\|f_{\theta}(y)-I_{s}\|+\beta\|y-g_{\theta}(z)\|.$
(5)

Figure 3:  Visualization of Super-resolution with limited running time. Left
to right: high resolution, low resolution, DIP (353 sec), DeepRED + TV (112.1
sec), DNCF (12.0 sec).

Here we pursue a smart initialization of $g_{\theta}(z)$, such as
$y=G_{\rho}(I_{s})+\epsilon$, where $G$ could be Gaussian smoothing filters of
standard deviation $\rho$ or bicubic interpolation, for specific tasks. This
is to suggest that in the degrade case of $f$ being an identity mapping, the
gradient of the objective function are nonzero. By initializing the
intermediate feature at $y_{0}=G_{\rho}(I_{s})$, to approximate a feasible
solution of image, but not too close to $I_{s}$ as to having a degenerated
solution like $I_{t}=I_{s}$ and $f$ is an identity mapping. We therefore add
additional regularization $\mathcal{R}$ and further refine the results by
gradient descent, by jointly optimizing the following:

$\displaystyle\mathcal{L}(\theta,y)=\|f_{\theta}(y)-I_{s}\|^{2}+\beta\|y-G_{\rho}(I_{s})\|^{2}+\mathcal{R}(\theta,y).$

We refer to our methods as the Deep Nonparametric Nonnegative Nonlinear
Filters (DNCF). As it degenerates to a simple filter $G$ as a trivial
solution, and if $\beta$ is small enough, it could lead to a fixed point
solution as $y=f_{\theta}(y)=I_{s}$. We generate the real scene by choosing a
better result (OPT) according to the photometric measurement, via heuristic
selection. This prevents the hard cases of DIP being degenerated, e.g. the
current result $f_{\theta}(y)$ within restricted running time still resembles
random signal (as $y$ is randomly initialized), and in this case we take
passive prediction $y$, or try a aggressive prediction
$f_{\theta}(f_{\theta}(y^{\star}))$ for luck,

$\displaystyle
I_{t}=OPT(y^{\star},f_{\theta}(y^{\star}),f_{\theta}(f_{\theta}(y^{\star})),\quad
y^{\star}=\arg\min_{y}\min_{\theta}\mathcal{L}(y,\theta).$

Figure 4:  Visualization of image inpainting. From left to right: the original
image, masked image, recovered images with DIP, and the results with D3NF.
(inpaint1 on the top and inpaint2 on the bottom).

Inspired by the previous work on regularization for blind deconvolution, we
also seek a proper formulation of $\mathcal{R}$ ideally with convexity, and
easier inference. We note that the second quadratic term $\|y-g(I_{s})\|^{2}$
contributes to the strong convexity to the objective function, and maintains
consistency between two branches of prediction. We also include a idea that a
non-negative combination of convex functions is still convex. By mathematical
reduction, we know that a sequential network $f$ is convex w.r.t. the random
vector $z$, given input data $I_{s}$, if all weights $\theta$ have non-
negative elements, and the activation function like ReLU is convex and non-
decreasing. The pioneering work from the input-convex network [1] uses this
rule and puts a strong nonnegative constraint on the parameters, which also
limits the expressiveness of the network and therefore the final performance.
In this paper, we use an interpolation between the convex network and an
arbitrary one, through a soft regularization on the parameter weights. We use
$\mathcal{R}(\theta)$ denotes the additional regularization to the negative
weights

$\displaystyle\mathcal{R}(\theta)=\gamma\|\max(-\theta,0)\|.$ (6)

The norm may change to $\ell_{1}$ or $\ell_{2}$. In an extreme case that
$\gamma$ goes to infinity, then this becomes a constrained problem that
requires all weights are nonnegative as [1]. Another middle way that goes
between the rigorous input convexity and nonconvexity is to limit partial
channels, or partial layers to have nonnegative weights:

$\displaystyle\mathcal{R}^{\prime}(\theta)=\sum\gamma\|\max(-\theta[indx],0)\|,$
(7)

where $indx$ is the indices of a subset of parameters. We typically only
select those layers above a predefined depth to regularize, to keep lower
layer being expressive.

We use alternating optimization between the network parameters $\theta$ and
the intermediate variables $y$, similiar to the Expectation-Maximization (EM).
Although EM is nonconvex in general and is hard to optimize, a part of it
might have an analytical solution. E.g., gradient descent can achieve a linear
convergence rate on a strongly convex and smooth function; that is, to obtain
a suboptimality gap decreases at an exponential speed[22]. However, without
the strong convexity, the decreasing rate may slow down inversely linear (for
smooth and convex functions)[21]. Unfortunately, the (strong) convex property
does not hold for the popular approaches based on neural networks.

##  3 Computational photography

Denoising. Blind denoising is to recover the underlying clean image from a
noisy observation without knowledge of the noise formation process. Assuming
the additive noise model as $I_{s}=I_{t}+\epsilon$, where $I_{s}$ is source
image for this task, i.e the noisy image and $I_{t}$ is the target image, i.e.
clean image, and $\epsilon$ follows Gaussian distribution. DNCF optimizes the
following objective function ($G$ as Gaussian smoothing):

$\displaystyle\mathcal{L}(\theta,y)=\|f_{\theta}(y)-I_{s}\|^{2}+\beta\|y-G(I_{s})\|^{2}+\mathcal{R}(\theta,y).$

Texture synthesis and inpainting. This technique [24][16][11][10] is
complementary to computational photography, if there is a defected region in
pictures, and the synthesis and inpainting technique is hereafter to repair by
filling in using textures from visually similar images. Some further
applications include artistic manipulation of removing an object or individual
from a picture while keeping the other pixels as photo-realistic as possible.
There were several lines of classical approaches. For the inpainting problem
with a picture $I_{s}$ that is removed pixels according to a mask
$M\in\\{0,1\\}^{H\times W}$, DIP optimize the parameterization for the
following objective function,

$\displaystyle\mathcal{L}_{inpaint}(\theta,y)=\|M\odot f_{\theta}(y)-M\odot
I_{s}\|^{2}+\mathcal{R}(\theta,y)$

Single image super-resolution (SR). This is to generate higher resolution
images from a lower resolution image restricted by camera sensors. If the
resolution degradation problem is known to be from physical factors like blur
effects, i.e. the observed image is a simple smoothed interpolation of the
high-resolution target image, then SR algorithms could resemble image
deblurring by alternating the estimation the blur kernel and deconvolution,
such as blind image deconvolution as Eq.(3). By denoting $I_{s}$ as the source
image of lower resolution, $D$ as the downsampling filter, $G$ as a one step
SR filter, e.g. bicubic interpolation, DNCF-SR optimizes the following
objective:

$\displaystyle\mathcal{L}_{SR}{(\theta,y)}=\|D(f_{\theta}(y))-I_{s}\|^{2}+\beta\|y-G(I_{s})\|^{2}+\mathcal{R}(\theta,y).$

Flash/no-flash photography. Under some lousy lighting conditions, there will
be a flash lighting control trade-off challenging to find. On the one hand,
flash images contain unrealistic colors as photons from ambient illumination
is comparably under-estimated, and this reduces the brightness of the natural
color of objects, on the other hand, no-flash images have to rely on strong
gains of cameras as photons are insufficient, and makes a stronger noise being
inevitable as a byproduct. A representation algorithm [23] is to use a joint
bilateral filter on the no-flash image, where the flash image serves to
provide edges information due to its advantage of less noise. To apply DIP on
this task, we denote the flash image as $I_{f}$ and the no-flash image as
$I_{nf}$, and let $\epsilon$ follows multivariate Gaussian distribution and
$\gamma$ to be the noise regularization, $G$ as the Gaussian smoothing filter,
and DNCF-Flash optimizes the following:

$\displaystyle\mathcal{L}_{flash}(\theta,y)=\|f_{\theta}(y)-I_{nf}\|^{2}+\beta\|y-G(I_{f}+\epsilon)\|^{2}+\mathcal{R}(\theta,y).$

Figure 5:  Visualization of the optimization procedure. In left-right, top-
down order, flash image, initialized $I_{t}$, $5$ middle results by every 200
iterations, and then the original non-flash image.

##  4 Adverarial defense for machine learning

There are several ways to categorize adversarial attack methods. One
perspective to classify is based on the final goal of the adversary. For
example, poisoning attacks refer to add fake images and labels to make the
overall trained models be poor performance [2], and this of course, requires
access to the training sets and are plausible for web-crawled data. The
evasion attacks give up on misleading the models but try to generate samples
that are hard to recognize, and this resembles an extreme version of data
augmentation, e.g. to put a mask on a traffic sign [12]. Some adversary
attacks aim for fix classes of samples $I$, say with label $t$, and they try
to generate visually similar samples $I^{\prime}$ but with false label
$t^{\prime}$, and these are referred to as targeted attack, and the
generalized attacks without specific label designation are non-targeted
attacks. From another perspective, we could categorize the attacks by whether
the adversary has access to the machine learning models configurations and
parameters, whose leakage will boost the capability of attacks. White-box
attacks refer to the methods with knowledge of the machine learning models.
such as L-BFGS [25], FGSM [14], BIM [17] and CW [4].

Fast gradient sign method (FGSM) This was proposed in [14]. We denote the
classifier model as $C$, and the misleading images $I^{adv}$ are generated
from real samples $I$ and the true label $t$ as

$\displaystyle
I^{adv}=I+\epsilon\operatorname*{sign}(\nabla_{I}\mathcal{L}(C(I),t)).$ (8)

Basic iterative method (BIM) This [17] could be viewed as a multi-step variant
of FGSM of smaller step size $\alpha$ with clipping to $\epsilon$-ball each
iteration, which is a projecting operation to restrict the adversary examples
being within $\epsilon$ distance from the original image $I_{0}$

$\displaystyle
I^{adv}_{n}=\textit{Proj}_{\epsilon}(I_{n-1}+\alpha\operatorname*{sign}(\nabla_{I}\mathcal{L}(C(I_{n-1}),t))).$
(9)

Carlini $\&$ Wagner The L-BFGS attacks was proposed in[25] which formulated
the attack as an optimization algorithm for finding an optimal adversarial
example $I$ that is predicted as an designated label $t$, which can be relaxed
as an regularized optimization problem

$\displaystyle\min_{I}\|I_{t}-I\|+\mathcal{L}(I,t),\quad s.t.\quad
I\in[0,1]^{m}.$ (10)

and solving this with L-BFGS. On top of the objective function, Carlini and
Wagner proposed the following margin loss in [4], to maximize the relative
confidence of the label $t$ comparing to all other labels,

$\displaystyle\min_{I}\|I_{t}-I\|+\max(\max_{t^{\prime}}\mathcal{L}(I,t^{\prime})-\mathcal{L}(I,t),0),$
(11)

DNCF as adversary defense The difference between the physical interference and
adversary attacks is that the later case could be as bad as possible for
classification. However, DNCF is insensitive to those pixel-wise changes, and
will generate results randomly each time. DNCF optimizes the following
function with $G$ as Gaussian smoothing and $\epsilon$ as random perturbation:

$\displaystyle\mathcal{L}(\theta,y)=\|I^{adv}-f_{\theta}(y)\|^{2}+\beta\|G(I^{adv})+\epsilon-y\|^{2},$

and perform the classification model on the denoising result
$f_{\theta^{\star}}(y^{\star})$. Comparing to computational photography, DNCF
might be easier to remove adversary noises without the need to fixate on the
pixel-wise accuracy for visual effects.

Figure 6:  ImageNet data attacked by CW, with predictions of attacked images
on the top and defended predictions on the bottom.

##  5 Optimization perspectives

This partial channel-wise convexification in Eq.(7) aims to provide a weaker
definition of strong convexity (SC), e.g., the constant nullspace strong
convexity (CNSC) [31] [33], and try to pursue a faster optimization algorithm
under these conditions. CNSC rises naturally from the notation of the general
convexity: as a function may not be strongly convex on the entire linear space
$\mathcal{Y}$, it is still possible that the function is strongly convex on a
subspace $\Phi$ of $\mathcal{Y}$; however, it may not be convex on the
orthogonal space $\Phi^{\top}$. For notation, we use $\textit{Proj}_{\Phi}$ to
represent the projection to a subspace $\Phi$.

######  Definition 1.

A function $f(y)$ that is twice-differentiable w.r.t. $y\in\mathbb{R}^{K}$, is
defined as having constant nullspace strong convexity of constant $m\geq 0$ if
there exists a subspace $\Phi\in\mathbb{R}^{K^{\prime}}$ on which the function
$f(y)$ only depends, as $z=\textit{Proj}_{\Phi}(y)$ and $f^{\prime}(z)=f(y)$,
and the Hessian matrix $H(z)$ of $f^{\prime}(z)=f(y)$ has the following
property:

$\displaystyle v^{\top}H(z)v\geq m\|v\|^{2},\quad v^{\top}H(z)=0,\quad\forall
v\in\Phi.$ (12)

We notice that if the subspace $\Phi$ is equivalent to $\mathbb{R}^{K}$, then
CNSC equals to SC (Strong Convexity). Generally speaking, for the problem that
takes much fewer data samples, technically, only one sample during inference,
than the number of elements in that sample, CNSC is a more reasonable
condition to satisfy since the inference procedure, since the gradient vector
w.r.t. input lies in a much lower-dimensional subspace. For one last layer in
the decomposition,

$\displaystyle\mathcal{L}(y)=\|\sigma(\theta_{l}^{\top}y+\theta^{\top}_{l}f_{h}(z)+b)-I_{s}\|^{2}_{2}+\lambda\|y\|^{2},$
(13)

then the gradient vectors $\partial\mathcal{L}(y)$ lies in the subspace
spanned by $\theta_{l}$, whose rank might very likely be much lower than $y$,
which is typically high-dimensional by design.

######  Lemma 2.

If the objective function lies in subspace $\Phi$, that there exists an
alternative expression of $L$, say $L^{\prime}$ only depends on the subspace
$\Phi$, such that
$\mathcal{L}(y)=\mathcal{L}^{\prime}(\textit{Proj}_{\Phi}(y))$, then the
gradient vector $p(y)$ and Hessian matrix $H$ also lies in the subspace
$\Phi$.

We propose to use a quasi-Newton optimization algorithm to update the
intermediate features $y$. The basic idea of this second-order algorithms is
to minimize the second-order expansion around each iteration. Based on the the
last iteration $y_{t}$, the gradient vector
$p_{t}=\partial\mathcal{L}/\partial{y^{t}}$, and an approximation
$B_{t}\in\mathbb{R}^{D\times D}$ to the Hessian matrix, we give the ascent
direction $d_{t}$ of this iteration by

$\displaystyle d_{t}=\arg\min_{d}p_{t}^{\top}d+\frac{1}{2}d^{\top}B_{t}d.$
(14)

The step size for this iteration $\alpha_{t}$ is obtained by line search,
which is iterative tried over $\\{b^{c}\\},c\in[0,1,\cdots],b\in(0,1)$, from
the largest to the smallest, to meet the Armijo rule of

$\displaystyle\mathcal{L}(y^{t}+\alpha_{t}d_{t})\leq\mathcal{L}(y^{t})+\alpha_{t}\xi
p_{t}^{\top}d_{t}.$ (15)

where $\xi\in(0,1)$ is a constant. Applying BFGS algorithm [8] [29], the
approximated Hessian $B_{t+1}$ for the next iteration is:

$\displaystyle
B_{t+1}=B_{t}-\frac{B_{t}s_{t}s_{t}^{\top}B_{t}}{s_{t}^{\top}B_{t}s_{t}}+\frac{o_{t}o_{t}^{\top}}{o_{t}^{\top}s_{t}},$
(16)

where $s_{t}=y_{t+1}-y_{t},\quad o_{t}=p_{t+1}-p_{t}$. However, only updating
$y$ leads to insufficient diversity in the result, as $y$ is constrained in
the subspace $\Phi$, so we heuristically update $\theta$ and $y$
alternatively.

Table 1: Performance of DNCF for adversary defense: the running time (seconds) of attacks. the accuracy ($\%$) of the original network, the accuracy after being attacked, and the accuracy after defense. The running time of DNCF is 126 seconds for 200 images. Attack | time | Ac(orig) | Ac(attack) | Ac(defense)
---|---|---|---|---
CW | 1646.2 | 75.48 | 0.00 | 64.42
BIM | 32.26 | 67.79 | 0.00 | 21.63
FGSM | 5.60 | 72.60 | 18.27 | 33.65
PGD | 30.21 | 74.04 | 0.00 | 15.87
FFGSM | 5.45 | 75.00 | 16.83 | 30.29

##  6 Experiments

Visualization of convexification. We first visualize the convexification
effect of Eq.(7) on the neural networks for smaller tasks. We give a simple
example of synthetic data as the fully input convex network (FICNN) as
described in [1], with a $2$-dimensional regression problem. This network
takes both the data $z\in\mathbb{R}^{2}$ and a label proposal $y\in\\{0,1\\}$
as the network input and predicts the compatibility score of the proposal. The
layers within the network have a special structure, in the $l$th hidden layer,
we require $z_{l}$ of $f$ has the following formulation,

$\displaystyle z_{l}=\sigma(\theta_{l}^{z}z_{l-1}+\theta_{l}^{y}y+b_{l}),\quad
z_{l}\in\mathbb{R}^{K_{l}},$ (17)

where $\sigma$ is ReLU function. The last layer projects to a scalar variable
and the loss function is square distance. We name the network with our
proposed regularization as CVXR-Net. FICNN uses the same network architecture
with CVXR-Net, they both have two fully connected layers of $200$-dimensional
latent features. Beside Partially ICNN (PICNN) is built on top of ICNN, except
that latent features $z_{l}$ have a branch connecting $y$ that allows positive
weights and therefore has relatively stronger expressive power. FICNN and
PICNN use projected gradient descent for optimization to keep the network
weights strictly non-negative. After training both networks to convergence, we
demonstrate the decision boundary in figure(7) for several synthetic data
sets, where red and blue points are from two different classes. We see that
CVXR-Net produces a better boundary for harder examples and maintains a
similar behavior on easier examples.

Table 2: Performance comparison between DNCF (less than 700 iterations) and DIP (with unconstrained running time). Here $\downarrow$ means better images have smaller values and vice versa.  task | DIP | DNCF | Figures
---|---|---|---
Denoise | 9.9e-4$\downarrow$ | 9.4e-3 $\downarrow$ | fig1
Flash | 6.887e-3 $\downarrow$ | 6.705e-3 $\downarrow$ | fig5
Inpaint1 | 7.08e-4$\downarrow$ | 5e-4 $\downarrow$ | fig4
Inpaint2 | 9.521e-3$\downarrow$ | 8.4e-3$\downarrow$ | fig4
SR | 19.067$\uparrow$ | 19.447 $\uparrow$ | fig2

Figure 7:  Visualization of prediction regions. 3 columns on the top:
Nonnegative regularization; 3 columns on the bottom: ICNN.

Computational photography. We conduct the experiments of DNCF with simpler
version of the original implementation of DIP by smaller network. The visual
results are in figures (1 2 4). We use the skip network for all experiments,
an encoder-decoder architecture with skip connections. The network consists of
several blocks, and each block has an equal number of channels in each of the
convolution layers. Each block has four convolution layers, four batch
normalization layers, and four ReLU layers, where the second convolution layer
has a stride of $\times 2$, therefore, it reduces the resolution and an
upsampling layer follows to expand the resolution by $\times 2$. The last
block also consists of an extra convolution layer to project the last feature
tensor to $3$-channel that matches the images. For general applications, the
output of the network is in the same resolution with the input, except that
for super-resolution, there is a $\times 8$ resolution expansion effect, since
we set to generate images in this resolution. We compare against total
variation (TV) denoising [5], bicubic interpolation, and DeepRED [19] combined
to We restricted the algorithm to run for at most 400 iterations on NVIDIA GTX
1060 with 6GB graphics memory, as the original parameter of DIP runs
considerably slow, as long as 300 seconds for a best result of super-
resolution. We put the details of results for different tasks in Table.(2).We
start with $\gamma$ in the table as the regularization coefficient. Each time
we found that the convexification regularization is larger than the mean
square error (MSE) loss in forward-propagation, we reduce $\gamma$ to $1/4$ of
its last value, to prevent that the regularization dominates the total loss.
This adaptive regularization also saves us from extensive tuning $\gamma$ for
each task individually. Our essence is to make the DIP practical on personal
computers, therefore despite that, we inherit most hyperparameters, e.g, the
number of layers, we reduce numbers of channels for those applications to fit
in $4GB$ GPU memory, and reduce the number of extra iterations without visible
changes. We report the performance in Table.(2). In the table, we show MSE of
both DIP and DNCF; except for the super-resolution (SR) task, the metric is
the peak SNR (PSNR) of the ground-truth against the SR version. Here
$\uparrow$ means the higher the better and $\downarrow$ means reversely. On
the flash task, we run $800$ iterations of optimizing the whole network and
then only optimizing the last block and its $32$-channel intermediate features
for about $1,200$ iterations using LBFGS, then compare with DIP baseline which
is to optimize the whole network using Adam optimization, for the same period
of time. We plot the convergence of MSE, the transition of the image flash-no
flash photography in Figure(5), i.e. from no-flash to flash version, for every
$200$ iterations. Note that for the $5$ middle results of $I_{t}$, and $600$
iterations on a $32$-channel DNCF give a good enough result, in comparison to
original DIP of $96$-channels with the same amount of iterations. The
convexification improves the visual details and decreases MSE. In Table.(2)
DNCF reduces MSE on all applications. We also see a great amount of
information from the figures, e.g. in Figure.(1,2), DNCF generates images with
weaker blur degradation and clearer quality; in Figure.(4) DNCF fills in the
hole with better patterns.

Adversary Defense. We tested the effect of removing adversary attacks from
images and use GoogleNet Inception V3 as the classifier. We randomly select
100 classes from the validation set of ImageNet ILSVRC 2012. We randomly
tested 200 images for each methods of attacking (so the accuracy on clean
images varies slightly), including the aforementioned methods, and projected
gradient descent (PGD) [18] and Fast FGSM (FFGSM) [28]. We set the $\epsilon$
value (the norm of maximum pertubation on the images), $\alpha$ value (the
step size) and $n$ vlaue (number of steps) as follows:
$\epsilon=8/256,\quad\alpha=2/256,\quad n=7$. We use a simple DNCF of only 3
convolution layers. In this case, we could use a batch of image and run them
in parallel, so they actually share a network parameterization as we have
lower standard in photometric measurement. We put the classification results
on original images, attacked images and defended images on Table.(1). We show
the running time comparison between attacking and defense as well. We notice
that DNCF is extremely effective against CW, and have different improvement
against others.

##  7 Conclusion

In this paper we propose a deep nonparametric filter (DNCF) as a general
framework of modelling physical interference and adversary attacks on images,
and lead to a solution of computational photography and adversary defense.

## References

  * [1] Brandon Amos, Lei Xu, and J Zico Kolter.  Input convex neural networks.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 146–155. JMLR. org, 2017.
  * [2] Battista Biggio, B Nelson, and P Laskov.  Poisoning attacks against support vector machines.  In 29th International Conference on Machine Learning, pages 1807–1814. ArXiv e-prints, 2012.
  * [3] Bolun Cai, Xiangmin Xu, Kui Jia, Chunmei Qing, and Dacheng Tao.  Dehazenet: An end-to-end system for single image haze removal.  IEEE Transactions on Image Processing, 25(11):5187–5198, 2016.
  * [4] Nicholas Carlini and David Wagner.  Towards evaluating the robustness of neural networks.  In 2017 ieee symposium on security and privacy (sp), pages 39–57. IEEE, 2017.
  * [5] Antonin Chambolle.  An algorithm for total variation minimization and applications.  Journal of Mathematical imaging and vision, 20(1):89–97, 2004.
  * [6] Tony F Chan and Jianhong Jackie Shen.  Image processing and analysis: variational, PDE, wavelet, and stochastic methods, volume 94.  Siam, 2005.
  * [7] Zezhou Cheng, Matheus Gadelha, Subhransu Maji, and Daniel Sheldon.  A bayesian perspective on the deep image prior.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5443–5451, 2019.
  * [8] John E Dennis, Jr and Jorge J Moré.  Quasi-newton methods, motivation and theory.  SIAM review, 19(1):46–89, 1977.
  * [9] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang.  Image super-resolution using deep convolutional networks.  IEEE transactions on pattern analysis and machine intelligence, 38(2):295–307, 2015.
  * [10] Alexei A Efros and William T Freeman.  Image quilting for texture synthesis and transfer.  In Proceedings of the 28th annual conference on Computer graphics and interactive techniques, pages 341–346, 2001.
  * [11] Alexei A Efros and Thomas K Leung.  Texture synthesis by non-parametric sampling.  In Proceedings of the seventh IEEE international conference on computer vision, volume 2, pages 1033–1038. IEEE, 1999.
  * [12] Kevin Eykholt, Ivan Evtimov, Earlence Fernandes, Bo Li, Amir Rahmati, Chaowei Xiao, Atul Prakash, Tadayoshi Kohno, and Dawn Song.  Robust physical-world attacks on deep learning visual classification.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1625–1634, 2018.
  * [13] Matheus Gadelha, Rui Wang, and Subhransu Maji.  Shape reconstruction using differentiable projections and deep priors.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 22–30, 2019.
  * [14] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572, 2014.
  * [15] Michael D Grossberg and Shree K Nayar.  Modeling the space of camera response functions.  IEEE transactions on pattern analysis and machine intelligence, 26(10):1272–1282, 2004.
  * [16] David J Heeger and James R Bergen.  Pyramid-based texture analysis/synthesis.  In Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, pages 229–238, 1995.
  * [17] Alexey Kurakin, Ian Goodfellow, and Samy Bengio.  Adversarial machine learning at scale.  arXiv preprint arXiv:1611.01236, 2016.
  * [18] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  arXiv preprint arXiv:1706.06083, 2017.
  * [19] Gary Mataev, Peyman Milanfar, and Michael Elad.  Deepred: Deep image prior powered by red.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV) Workshops, Oct 2019.
  * [20] Calvin S McCamy, Harold Marcus, James G Davidson, et al.  A color-rendition chart.  J. App. Photog. Eng, 2(3):95–99, 1976.
  * [21] Yurii Nesterov and Arkadii Nemirovskii.  Interior-point polynomial algorithms in convex programming, volume 13.  Siam, 1994.
  * [22] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [23] Georg Petschnigg, Richard Szeliski, Maneesh Agrawala, Michael Cohen, Hugues Hoppe, and Kentaro Toyama.  Digital photography with flash and no-flash image pairs.  ACM transactions on graphics (TOG), 23(3):664–672, 2004.
  * [24] Javier Portilla and Eero P Simoncelli.  A parametric texture model based on joint statistics of complex wavelet coefficients.  International journal of computer vision, 40(1):49–70, 2000.
  * [25] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  arXiv preprint arXiv:1312.6199, 2013.
  * [26] Xin Tao, Hongyun Gao, Xiaoyong Shen, Jue Wang, and Jiaya Jia.  Scale-recurrent network for deep image deblurring.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8174–8182, 2018.
  * [27] Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Deep image prior.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 9446–9454, 2018.
  * [28] Eric Wong, Leslie Rice, and J Zico Kolter.  Fast is better than free: Revisiting adversarial training.  arXiv preprint arXiv:2001.03994, 2020.
  * [29] Stephen Wright and Jorge Nocedal.  Numerical optimization.  Springer Science, 35(67-68):7, 1999.
  * [30] Junyuan Xie, Linli Xu, and Enhong Chen.  Image denoising and inpainting with deep neural networks.  In Advances in neural information processing systems, pages 341–349, 2012.
  * [31] Ian En-Hsu Yen, Cho-Jui Hsieh, Pradeep K Ravikumar, and Inderjit S Dhillon.  Constant nullspace strong convexity and fast convergence of proximal methods under high-dimensional settings.  In Advances in Neural Information Processing Systems, pages 1008–1016, 2014.
  * [32] Kai Zhang, Wangmeng Zuo, Shuhang Gu, and Lei Zhang.  Learning deep cnn denoiser prior for image restoration.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3929–3938, 2017.
  * [33] Kai Zhong, Ian En-Hsu Yen, Inderjit S Dhillon, and Pradeep K Ravikumar.  Proximal quasi-newton for computationally intensive l1-regularized m-estimators.  In Advances in Neural Information Processing Systems, pages 2375–2383, 2014.