A Suzaku, NuSTAR, and XMM-Newton view on variable absorption and relativistic reflection in NGC 4151

TL;DR

This study uses Suzaku, NuSTAR, and XMM-Newton data to disentangle complex absorption and relativistic reflection in NGC 4151, revealing variable absorption features and a vertically extended corona.

Contribution

It introduces a physically motivated lamp-post reflection model with simultaneous multi-instrument data, highlighting the role of a vertically extended corona and variable absorbers in NGC 4151.

Findings

Detection of variable neutral and ionized absorbers on short and long timescales.

Evidence for a vertically extended corona rather than two compact sources.

Spectral variability driven by changes in absorber column densities.

Abstract

We disentangle X-ray disk reflection from complex line-of-sight absorption in the nearby Seyfert NGC 4151, using a suite of Suzaku, NuSTAR, and XMM-Newton observations. Extending upon earlier published work, we pursue a physically motivated model using the latest angle-resolved version of the lamp-post geometry reflection model relxillCp_lp together with a Comptonization continuum. We use the long-look simultaneous Suzaku/NuSTAR observation to develop a baseline model wherein we model reflected emission as a combination of lamp-post components at the heights of 1.2 and 15.0 gravitational radii. We argue for a vertically extended corona as opposed to two compact and distinct primary sources. We find two neutral absorbers (one full-covering and one partial-covering), an ionized absorber ($\log \xi = 2.8$), and a highly-ionized ultra-fast outflow, which have all been reported previously. All analyzed spectra are well described by this baseline model. The bulk of the spectral variability between 1 keV and 6 keV can be accounted for by changes in the column density of both neutral absorbers, which appear to be degenerate and inversely correlated with the variable hard continuum component flux. We track variability in absorption on both short (2 d) and long ($\sim$1 yr) timescales; the observed evolution is either consistent with changes in the absorber structure (clumpy absorber at distances ranging from the broad line region (BLR) to the inner torus or a dusty radiatively driven wind) or a geometrically stable neutral absorber that becomes increasingly ionized at a rising flux level. The soft X-rays below 1 keV are dominated by photoionized emission from extended gas that may act as a warm mirror for the nuclear radiation.

Figures

Tables

Table 1. Table 1: XMM-Newton  and Suzaku  observations in 2011/2012 with screened exposure times. Listed are the satellite, the observation ID of the observation, the date when the observation started, the exposure time after screening, and the number of counts detected. The check-symbol ( √ square-root \surd ) denotes observations that are considered for the data analysis. Observations labeled with a cross (✗) are excluded. The observation Nu Suz only contains data that are fully simultaneous to Suz 3.

abbrv.	obsid	det	time	exp [ks]	cnts [$\times {10}^4$]
(XMM)	657840101	pn	2011-05-11			✗
		MOS 1				✗
		MOS 2				✗
XMM 1	657840201	pn	2011-06-12	2.3	1.4	$\surd$
		MOS 1		0.2	0.07	✗
		MOS 2		0.1	0.04	✗
XMM 2	657840301	pn	2011-11-25	5.7	2.8	$\surd$
		MOS 1		6.9	1.3	$\surd$
		MOS 2		7.1	1.3	$\surd$
XMM 3	657840401	pn	2011-12-09	6.6	3.0	$\surd$
		MOS 1		8.8	1.5	$\surd$
		MOS 2		9.0	1.5	$\surd$
XMM 4	679780101	pn	2012-05-13	6.3	2.6	$\surd$
		MOS 1		8.6	1.3	$\surd$
		MOS 2		8.7	1.3	$\surd$
XMM 5	679780201	pn	2012-06-10	8.7	1.6	$\surd$
		MOS 1		12.5	0.9	$\surd$
		MOS 2		12.5	0.9	$\surd$
XMM 6	679780301	pn	2012-11-14	3.8	2.4	$\surd$
		MOS 1		1.8	0.2	$\surd$
		MOS 2		2.0	0.3	$\surd$
XMM 7	679780401	pn	2012-12-10	6.6	2.7	$\surd$
		MOS 1		9.5	0.5	$\surd$
		MOS 2		9.6	1.4	$\surd$
(XMM)	679780501	pn	2012-12-10			✗
		MOS 1				✗
		MOS 2				✗
Suz 1	906006010	XIS 0	2011-11-17	61.7	18.5	$\surd$
		XIS 1		61.7	17.7	$\surd$
		XIS 3		61.7	19.0	$\surd$
		HXD		54.6	8.1	$\surd$
Suz 2	906006020	XIS 0	2011-12-18	55.1	21.1	$\surd$
		XIS 1		55.1	21.3	$\surd$
		XIS 3		55.1	23.4	$\surd$
		HXD		32.1	7.3	$\surd$
Suz 3	707024010	XIS 0	2012-11-11	150.2	32.5	$\surd$
		XIS 1		150.2	32.8	$\surd$
		XIS 3		150.2	35.2	$\surd$
		HXD		139.9	18.9	$\surd$
NuSuz	60001111002/3/5	FPM A	2012-11-12/14	106.1	71.3	$\surd$
		FPM B		106.2	67.4	$\surd$

Table 2. Table 2: Best-fit parameters of the blend of Gaussian lines for Suz 3. All Gaussian lines are fitted with zero width. Centroid energies marked with ∞ \infty are adopted from Schurch et al. ( 2004 ) , Ogle et al. ( 2000 ) , and Vainshtein & Safronova ( 1978 ) , line blends marked with † † \dagger from the Chandra /ACIS study of Wang et al. ( 2011a ) . Frozen parameters are denoted with an asterisk ( ∗ ∗ \ast ). When marked with the symbol ♓ ♓ \pisces , the line identification is uncertain and can be confused with an instrumental edge. The lines are consistent with the study on XMM-Newton /RGS data in this work. Note that the O vii line falls outside the sensitive energy range of Suzaku /XIS. For the possible line-blends Ne ix  RRC / Mg xi  r and Si xiii  (r/f), we additionally fit for the centroid energies, because of unclear line identifications. Also, we require two narrow Gaussians with zero width at 1.58 keV and 1.72 keV, which may be due to intrinsic Si i  K α 𝛼 \alpha and Mg xi   1 ​ s ​ 3 ​ p − 1 ​ s 2 1 𝑠 3 𝑝 1 superscript 𝑠 2 1s3p-1s^{2} or due to calibration effects at the instrumental Al K and Si K edges.

Line	E∗ [keV] / ${\lambda }^{\ast }$ [Å]	flux [Ph s-1 cm-2]
O vii (r,i,f) †	0.57 / 21.80	$\left(1.1\pm 0.3\right)\times {10}^{-3}$
O viii Ly $\alpha$ ∞	0.65 / 18.97	$\left(2.67\pm 0.19\right)\times {10}^{-4}$
O vii RRC ∞	0.74 / 16.77	$\left(9.4\pm 1.0\right)\times {10}^{-5}$
O viii RRC ∞	0.87 / 14	$\left({7.6}_{-1.0}^{+0.7}\right)\times {10}^{-5}$
Ne ix (r,i,f) †	0.91 / 13.55	$\left(1.02\pm 0.09\right)\times {10}^{-4}$
Ne x Ly $\alpha$ †	1.03 / 12.04	$\left(1.8\pm 0.4\right)\times {10}^{-5}$
Ne ix RRC / Ne x Ly $\beta$ †	1.2 / 10.33	$\left(3.0\pm 0.4\right)\times {10}^{-5}$
Ne x RRC / Mg xi (r) ∞	${1.351}_{-0.005}^{+0.004}$ / $9.18\pm 0.03$	$\left(4.4\pm 0.3\right)\times {10}^{-5}$
Mg vii Ly $\alpha$ ∞	1.47 / 8.43	$\left(2.2\pm 0.3\right)\times {10}^{-5}$
Mg xi $1s3p-1s^2$ ${}_{}{}^{\infty }♓$	1.58 / 7.85	$\left(2.8\pm 0.4\right)\times {10}^{-5}$
Si i K$\alpha$ ${}_{}{}^{\infty }♓$	${1.7200}_{-0.0000}^{+0.0017}$ / ${7.2084}_{-0.0071}^{+0.0000}$	$\left(3.6\pm 0.4\right)\times {10}^{-5}$
Si xiii (r,f) ∞	${1.8473}_{-0.0013}^{+0.0008}$ / ${6.712}_{-0.003}^{+0.005}$	$\left(4.4\pm 0.4\right)\times {10}^{-5}$
Si xiv Ly $\alpha$ ∞	2.0 / 6.20	$\left(2.3\pm 0.4\right)\times {10}^{-5}$
Continuum	$\Gamma$	norm [Ph keV-1 s-1 cm-2]
	1.72∗	$(1.33\pm 0.05)\times {10}^{-3}{}^{\ast }$

Table 3. Table 3: Best-fit parameters of the baseline model for the simultaneous 150 ks data of Suz 3 and Nu Suz . The black hole spin of both lamp-post components is set to its maximum value of a = 0.998 𝑎 0.998 a=0.998 , while the radii of the inner and outer disk are kept at the default values of 1 ​ r EH 1 subscript 𝑟 EH 1\,r_{\mathrm{EH}} and 499 ​ r g 499 subscript 𝑟 g 499\,r_{\mathrm{g}} , respectively. Parameters marked with the symbol † † \dagger are tied amongst one another while those marked with an asterisk ( ∗ ∗ \ast ) are frozen. The reflection fraction is no free parameter here, as the incident continuum ( nthcomp ) and both reflection continua (LP 1 , LP 2 ) are fitted independently ( refl_frac=-1 ). The normalization of nthcomp is defined at unity for a norm of 1 at 1 keV. The normalization of xillver and relxill is defined in the Appendix of Dauser et al. ( 2016 ) .

Model component	Parameter	Value
${\chi }^2$ (dof) = 1094 (762)
Detconst	XIS 0	$0.998\pm 0.005$
	XIS 1	$0.954\pm 0.005$
	XIS 3	$1.012\pm 0.005$
	HXD	$1.215\pm 0.009$
	FPMA$\ast$	1
	FPMB	$1.030\pm 0.004$
CAGal	$N_{\text{H,Gal}}^{\ast }$ [${10}^{22}{\mathrm{cm}}^{-2}$]	0.023
XSTAR 1	$N_{\text{H}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	$1.2\pm 0.3$
	$\log \xi [\mathrm{erg}\mathrm{cm}{s}^{-1}$]	${2.82}_{-0.11}^{+0.10}$
CA1	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${21}_{-2}^{+4}$
cov. factor	f${}_{\text{cvr}}$	${0.46}_{-0.05}^{+0.06}$
CA2	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${8.8}_{-0.6}^{+0.5}$
Abs. line	$E\left[\mathrm{keV}\right]$	${8.17}_{-0.09}^{+0.11}$
	$\sigma [\mathrm{keV}$]	${0.34}_{-0.07}^{+0.14}$
	Depth [$2\pi \sigma {\tau }_{\mathrm{line}}$]	${0.0339}_{-0.0020}^{+0.0124}$
HXCOMP	norm	$0.047\pm 0.002$
	${\Gamma }^{\ast }$	1.72
	k $T_e^{\ast }\left[\mathrm{keV}\right]$	399
	k $T_{\mathrm{bb}}^{\ast }\left[\mathrm{keV}\right]$	0.05
Fe K$\alpha$	norm [$\mathrm{Ph}{s}^{-1}{\mathrm{cm}}^{-2}$]	$\left({2.27}_{-0.14}^{+0.13}\right)\times {10}^{-4}$
	$E\left[\mathrm{keV}\right]$	${6.394}_{-0.006}^{+0.005}$
Fe K$\beta$	norm [$\mathrm{Ph}{s}^{-1}{\mathrm{cm}}^{-2}]{}^{\ast }$	$0.12\times {\mathrm{norm}}_{\mathrm{Fe}K\alpha }$
	$E^{\ast }\left[\mathrm{keV}\right]$	7.1
LP1	norm	${6.9}_{-4.9}^{+1.2}$
	height [$r_{\mathrm{EH}}$]	${1.1000}_{-0.0000}^{+0.0013}$
	$i^{†}$	$3_{-0}^{+6}$
	${\Gamma }^{\ast }$	1.72
	$\log \xi$ [$\mathrm{erg}\mathrm{cm}{s}^{-1}$]	${2.835}_{-0.048}^{+0.016}$
	$Z_{\mathrm{Fe}}^{†}$	${2.5}_{-0.5}^{+0.6}$
	k$T_e^{\ast }\left[\mathrm{keV}\right]$	399
LP2	norm	$\left(8\pm 2\right)\times {10}^{-4}$
	height [$r_{\mathrm{EH}}$]	${14.1}_{-1.9}^{+3.8}$
	$i^{†}$	$3_{-0}^{+6}$
	${\Gamma }^{\ast }$	1.72
	$\log {\xi }^{†}$ [$\mathrm{erg}\mathrm{cm}{s}^{-1}]$	2.84
	$Z_{\mathrm{Fe}}^{†}$	${2.5}_{-0.5}^{+0.6}$
	k$T_e^{\ast }\left[\mathrm{keV}\right]$	399
SXCOMP	norm	$\left(1.33\pm 0.05\right)\times {10}^{-3}$
	${\Gamma }^{\ast }$	1.72
	k $T_e^{\ast }\left[\mathrm{keV}\right]$	399
	k $T_{\mathrm{bb}}^{\ast }\left[\mathrm{keV}\right]$	0.05

Table 4. Table 4: Values obtained for reflection fraction R f subscript 𝑅 f R_{\mathrm{f}} , the reflection strength R s subscript 𝑅 s R_{\mathrm{s}} and normalization of the two lamp-post components LP 1 and LP 2 . Frozen parameters are denoted by an asterisk ( ∗ ∗ \ast ).

	free LP1	free LP2
norm LP1	$4.59\pm 0.13$	$0.306\pm 0.011$
norm LP2	$\left(6.6\pm 0.4\right)\times {10}^{-4}$	$\left(1.822\pm 0.019\right)\times {10}^{-3}$
$R_f^{{\mathrm{LP}}_1}$	$1.50\pm 0.06$	${22.5}^{\ast }$
$R_f^{{\mathrm{LP}}_2}$	${1.2}^{\ast }$	$0.44\pm 0.02$
$R_s^{{\mathrm{LP}}_1}$	$0.69$	${10.5}^{\ast }$
$R_s^{{\mathrm{LP}}_2}$	${0.64}^{\ast }$	$0.24$
$R_s^{{\mathrm{LP}}_2+{\mathrm{LP}}_2}$	$0.68$

Table 5. Table 5: Triplet line parameters of the ions O vii and Ne ix . Both lines are described with Voigt profiles. The thermal velocities are unconstrained and not listed in the table. The continuum parameters are frozen to those of the SXCOMP.

Line	E [keV]	FWHM [keV]	flux [Ph s-1 cm-2]
O vii w	${0.57200}_{-0.00006}^{+0.00013}$	$\left({2.6}_{-1.1}^{+0.5}\right)\times {10}^{-3}$	$\left({1.30}_{-0.10}^{+0.13}\right)\times {10}^{-4}$
O vii x+y	${0.56702}_{-0.00007}^{+0.00022}$	$\left({2.7}_{-0.8}^{+1.7}\right)\times {10}^{-3}$	$\left({7.9}_{-0.7}^{+1.8}\right)\times {10}^{-5}$
O vii z	${0.55946}_{-0.00006}^{+0.00005}$	$\left({2.5}_{-0.5}^{+0.6}\right)\times {10}^{-3}$	$\left({4.51}_{-0.13}^{+0.15}\right)\times {10}^{-4}$
Ne ix w	${0.9191}_{-0.0006}^{+0.0007}$	$\le 9\times {10}^{-3}$	$\left({2.4}_{-0.9}^{+0.5}\right)\times {10}^{-5}$
Ne ix x+y	${0.9132}_{-0.0015}^{+0.0013}$	$\le 0.041$	$\left({2.2}_{-0.7}^{+0.8}\right)\times {10}^{-5}$
Ne ix z	$0.9024\pm 0.0004$	$\le 6\times {10}^{-3}$	$\left(5.8\pm 0.6\right)\times {10}^{-5}$

Table 6. Table 6: Summary of all XMM-Newton -related model parameter that are left free to vary. The normalization of the Comptonized continua HXCOMP and SXCOMP are defined as unity if the norm equals one at 1 keV. Parameters that are marked with an asterisk ( ∗ ∗ \ast ) are frozen to the parameters found for the long-look observation Suz 3.

		XMM 1	XMM 2	XMM 3	XMM 4	XMM 5	XMM 6	XMM 7
Detconst	pn	1∗	1∗	1∗	1∗	1∗	1∗	1∗
	MOS 1	–	$1.29\pm 0.02$	$1.24\pm 0.02$	$1.27\pm 0.03$	$1.27\pm 0.03$	$1.27\pm 0.06$	${1.43}_{-0.04}^{+0.05}$
	MOS 2	–	$1.28\pm 0.02$	$1.28\pm 0.02$	$1.29\pm 0.03$	$1.33\pm 0.04$	$1.32\pm 0.06$	${1.32}_{-0.02}^{+0.03}$
HXCOMP	norm	${0.037}_{-0.007}^{+0.006}$	$0.0645\pm 0.0007$	${0.061}_{-0.005}^{+0.004}$	${0.077}_{-0.007}^{+0.005}$	${0.035}_{-0.006}^{+0.003}$	${0.055}_{-0.014}^{+0.007}$	${0.059}_{-0.005}^{+0.004}$
SXCOMP	norm	$\left({1.34}_{-0.24}^{+0.16}\right)\times {10}^{-3}$	$\left({1.42}_{-0.32}^{+0.06}\right)\times {10}^{-3}$	$\left(9.3\pm 3.0\right)\times {10}^{-4}$	$\left({1.36}_{-0.27}^{+0.10}\right)\times {10}^{-3}$	$\left({1.10}_{-0.04}^{+0.07}\right)\times {10}^{-3}$	$\left({1.3}_{-0.4}^{+0.2}\right)\times {10}^{-3}$	$\left({1.15}_{-0.05}^{+0.00}\right)\times {10}^{-3}$
LP 1	norm	$\le 10$	$\le 0.4$	$\le 2$	$\le 3$	$\le 3$	$\le 6$	$\le 2$
	height [$r_{\mathrm{EH}}$]	1.1∗	"	"	"	"	"	"
LP 2	norm	$\left({1.6}_{-0.9}^{+1.1}\right)\times {10}^{-3}$	$\left(1.06\pm 0.13\right)\times {10}^{-3}$	$\left(1.1\pm 0.5\right)\times {10}^{-3}$	$\left(2.2\pm 0.6\right)\times {10}^{-3}$	$\left({6.7}_{-1.0}^{+3.6}\right)\times {10}^{-4}$	$\left(9\pm 7\right)\times {10}^{-4}$	$\left({1.85}_{-0.13}^{+0.42}\right)\times {10}^{-3}$
	height [$r_{\mathrm{EH}}$]	14.1∗	"	"	"	"	"	"
CA1	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${25}_{-8}^{+12}$	${11.8}_{-0.7}^{+0.8}$	$14\pm 2$	$17\pm 4$	${25.3}_{-5.1}^{+1.9}$	${26}_{-7}^{+10}$	${12.3}_{-2.4}^{+1.4}$
	f${}_{\text{cvr}}$	0.46∗	"	"	"	"		"
CA2	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${12.4}_{-0.8}^{+1.6}$	${6.80}_{-0.26}^{+0.11}$	$5.9\pm 0.3$	${12.8}_{-0.6}^{+0.5}$	${15.6}_{-0.7}^{+0.4}$	${9.2}_{-0.7}^{+0.6}$	${9.05}_{-0.35}^{+0.17}$
Fe K$\alpha$	norm [$\mathrm{Ph}{s}^{-1}{\mathrm{cm}}^{-2}$]	$\left({2.4}_{-0.8}^{+0.7}\right)\times {10}^{-4}$	$\left({2.0}_{-0.4}^{+0.6}\right)\times {10}^{-4}$	$\left(2.5\pm 0.6\right)\times {10}^{-4}$	$\left({2.3}_{-0.7}^{+0.8}\right)\times {10}^{-4}$	$\left({2.8}_{-0.4}^{+0.2}\right)\times {10}^{-4}$	$\left({2.7}_{-0.9}^{+0.4}\right)\times {10}^{-4}$	$\left({1.650}_{-0.005}^{+0.621}\right)\times {10}^{-4}$
	E [keV]	6.4∗	"	"	"	"	"	"
XSTAR 1	$N_{\text{H}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	$\left(6_{-6}^{+34}\right)\times {10}^{21}$	$\left(1.2\pm 0.7\right)\times {10}^{22}$	$\left({1.2}_{-1.2}^{+1.4}\right)\times {10}^{22}$	$\left({3.2}_{-1.8}^{+2.6}\right)\times {10}^{22}$	$\le 8.5\times {10}^{21}$	$\left({1.9}_{-1.8}^{+2.8}\right)\times {10}^{22}$	$\le 9.5\times {10}^{21}$

Table 7. Table 7: Summary of all Suzaku -related model parameter that are left free to vary. The normalization of the Comptonized continua HXCOMP and SXCOMP are defined as unity if the norm equals one at 1 keV. Parameters that are marked with an asterisk ( ∗ ∗ \ast ) are frozen to the parameters found for the long-look observation Suz 3.

		Suz 1	Suz 2	Suz 2-1	Suz 2-2	Suz 3	Suz 3-1	Suz 3-2	Suz 3-3	Suz 3-4
Detconst	XIS 0	1∗	1∗	1∗	1∗	$0.998\pm 0.005$	${1.076}_{-0.018}^{+0.019}$	$1.010\pm 0.010$	${0.961}_{-0.007}^{+0.008}$	$0.957\pm 0.008$
Detconst	XIS 1	$1.017\pm 0.006$	$1.032\pm 0.006$	$0.949\pm 0.009$	$1.047\pm 0.006$	$0.954\pm 0.005$	${1.043}_{-0.017}^{+0.019}$	$1.000\pm 0.011$	$0.910\pm 0.007$	$0.923\pm 0.008$
	XIS 3	$0.995\pm 0.006$	$0.975\pm 0.006$	$0.977\pm 0.009$	${0.986}_{-0.005}^{+0.006}$	$1.012\pm 0.005$	${1.113}_{-0.018}^{+0.020}$	$1.053\pm 0.011$	$0.997\pm 0.008$	${0.993}_{-0.008}^{+0.009}$
	HXD	$1.30\pm 0.02$	$1.29\pm 0.02$	$1.32\pm 0.03$	$1.09\pm 0.04$	$1.215\pm 0.009$	$1.35\pm 0.03$	${1.239}_{-0.018}^{+0.019}$	$1.186\pm 0.013$	$1.197\pm 0.014$
	FPMA	–	–	–	–	1∗	1∗	1∗	1∗	1∗
	FPMB	–	–	–	–	$1.030\pm 0.004$	$1.05\pm 0.02$	$0.998\pm 0.007$	$1.039\pm 0.005$	$1.039\pm 0.007$
HXCOMP	norm	$0.0698\pm 0.0013$	$0.0726\pm 0.0012$	$0.072\pm 0.002$	${0.055}_{-0.004}^{+0.007}$	$0.047\pm 0.002$	${0.0406}_{-0.0010}^{+0.0012}$	$0.0481\pm 0.0005$	$0.0445\pm 0.0004$	$0.0509\pm 0.0005$
SXCOMP	norm	$1.33\times {10}^{-3}$ ∗	"	"	"	$\left(1.33\pm 0.05\right)\times {10}^{-3}$	"	"	"	"
LP1	norm	$\left({5.3}_{-1.2}^{+1.4}\right)\times {10}^{-4}$	$\left({7.1}_{-1.9}^{+2.8}\right)\times {10}^{-4}$	$\left(6_{-3}^{+4}\right)\times {10}^{-4}$	${0.4}_{-0.3}^{+0.4}$	${6.9}_{-4.9}^{+1.2}$	"	"	"	"
	height [$r_{\mathrm{EH}}$]	${24}_{-7}^{+9}$	${14}_{-4}^{+6}$	${14}_{-6}^{+10}$	${1.53}_{-0.15}^{+6.31}$	${1.1000}_{-0.0000}^{+0.0013}$	"	"	"	"
LP2	norm	–	–	–	–	$\left(8\pm 2\right)\times {10}^{-4}$	"	"	"	"
	height [$r_{\mathrm{EH}}$]	–	–	–	–	${14.1}_{-1.9}^{+3.8}$	"	"	"	"
CA1	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${13}_{-2}^{+3}$	${9.2}_{-1.1}^{+1.5}$	${12}_{-4}^{+6}$	${5.6}_{-0.9}^{+1.6}$	${21}_{-3}^{+4}$	${13.7}_{-1.8}^{+4.8}$	${12.2}_{-1.4}^{+2.3}$	${24.2}_{-1.7}^{+1.9}$	${15.7}_{-1.4}^{+2.1}$
	f${}_{\text{cvr}}$	${0.46}_{-0.10}^{+0.11}$	$0.45\pm 0.06$	${0.33}_{-0.08}^{+0.13}$	${0.49}_{-0.10}^{+0.09}$	${0.46}_{-0.05}^{+0.06}$	${0.57}_{-0.18}^{+0.13}$	$0.57\pm 0.14$	$0.51\pm 0.04$	$0.57\pm 0.09$
CA2	$N_{\text{H,int}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	${6.9}_{-0.8}^{+0.7}$	$4.4\pm 0.3$	${5.3}_{-0.6}^{+0.4}$	$3.9\pm 0.3$	${8.8}_{-0.6}^{+0.5}$	${7.7}_{-1.3}^{+1.6}$	${7.3}_{-1.4}^{+1.2}$	${9.1}_{-0.6}^{+0.5}$	${7.7}_{-1.0}^{+0.9}$
Fe K$\alpha$	norm [$\mathrm{Ph}{s}^{-1}{\mathrm{cm}}^{-2}$]	$\left(2.4\pm 0.2\right)\times {10}^{-4}$	$\left(2.5\pm 0.2\right)\times {10}^{-4}$	$\left(2.8\pm 0.3\right)\times {10}^{-4}$	$\left({2.70}_{-0.14}^{+0.16}\right)\times {10}^{-4}$	$\left({2.27}_{-0.14}^{+0.13}\right)\times {10}^{-4}$	$\left({2.00}_{-0.18}^{+0.22}\right)\times {10}^{-4}$	$\left({1.75}_{-0.11}^{+0.16}\right)\times {10}^{-4}$	$\left({2.35}_{-0.09}^{+0.16}\right)\times {10}^{-4}$	$\left(2.31\pm 0.14\right)\times {10}^{-4}$
	E [keV]	6.394∗	"	"	"	${6.394}_{-0.006}^{+0.005}$	"	"	"	"
XSTAR 1	$N_{\text{H}}$ [${10}^{22}{\mathrm{cm}}^{-2}$]	$1.22\times {10}^{22}$	"	"	"	$\left(1.2\pm 0.3\right)\times {10}^{22}$	"	"	"	"

Table 8. Table 8: Line parameters of Voigt profiles fitted to the combined RGS data of all XMM-Newton  observations as part of the monitoring. The parameters of the O vii and Ne ix lines are adopted from Table  5 . The thermal velocities are unconstrained and not listed in the table. The continuum parameters are frozen to those of the SXCOMP.

Line	E [keV]	FWHM [keV]	$k_{B}T$ [keV]	flux [Ph s-1 cm-2]
C vi $\alpha$	$0.36637\pm 0.00008$	$\left({4.3}_{-1.5}^{+1.8}\right)\times {10}^{-3}$	–	$\left(2.7\pm 0.3\right)\times {10}^{-4}$
C v RRC	${0.3901}_{-0.0005}^{+0.0004}$	–	$\left({4.4}_{-1.0}^{+1.4}\right)\times {10}^{-3}$	$\left(1.9\pm 0.3\right)\times {10}^{-4}$
N vi f	$0.41855\pm 0.00009$	$\le 1.5\times {10}^{-3}$	–	$\left({1.35}_{-0.15}^{+0.16}\right)\times {10}^{-4}$
N vi i	$0.4250\pm 0.0004$	$\le 6\times {10}^{-3}$	–	$\left(2\pm 3\right)\times {10}^{-5}$
N vi r	$0.42909\pm 0.00015$	$\left(8_{-8}^{+30}\right)\times {10}^{-4}$	–	$\left({9.6}_{-1.5}^{+1.6}\right)\times {10}^{-5}$
C vi $\beta$	$0.4342\pm 0.0002$	$\left(7_{-7}^{+50}\right)\times {10}^{-4}$	–	$\left({3.6}_{-3.6}^{+1.4}\right)\times {10}^{-5}$
C vi $\gamma$	$0.4580\pm 0.0002$	$\le 4.8\times {10}^{-2}$	–	$\left({3.5}_{-3.5}^{+0.5}\right)\times {10}^{-6}$
C vi RRC	${0.4875}_{-0.0006}^{+0.0005}$	–	${0.010}_{-0.005}^{+0.000}$	$\left({1.5}_{-0.4}^{+0.2}\right)\times {10}^{-4}$
N vii Ly $\alpha$	${0.49874}_{-0.00008}^{+0.00015}$	$\le 2.1\times {10}^{-2}$	–	$\left({1.6}_{-1.6}^{+0.5}\right)\times {10}^{-3}$
O vii f	0.55946	$2.5\times {10}^{-3}$	–	$4.5\times {10}^{-4}$
O vii i	0.56702	$2.7\times {10}^{-3}$	–	$7.9\times {10}^{-5}$
O vii r	0.57200	$2.6\times {10}^{-3}$	–	$1.3\times {10}^{-4}$
O viii Ly $\alpha$	${0.65158}_{-0.00011}^{+0.00013}$	$\left({3.4}_{-1.8}^{+2.1}\right)\times {10}^{-3}$	–	$\left({1.58}_{-0.12}^{+0.11}\right)\times {10}^{-4}$
N vii RRC	${0.6629}_{-0.0007}^{+0.0004}$	–	$\left({10}_{-0}^{+1.2}\right)\times {10}^{-3}$	$\left({2.5}_{-0.5}^{+0.6}\right)\times {10}^{-5}$
O vii RRC	$0.7366\pm 0.0005$	–	$\left({3.0}_{-0.7}^{+0.8}\right)\times {10}^{-3}$	$\left(5.5\pm 0.7\right)\times {10}^{-5}$
O viii Ly $\beta$	$0.7723\pm 0.0006$	$\le 4.5\times {10}^{-3}$	–	$\left({1.6}_{-0.4}^{+0.5}\right)\times {10}^{-5}$
O viii Ly $\gamma$	$0.8228\pm 0.0008$	$\le 6.3\times {10}^{-2}$	–	${0.094}_{-0.094}^{+0.007}$
O viii RRC	${0.8677}_{-0.0010}^{+0.0008}$	–	$\left({4.0}_{-1.4}^{+1.9}\right)\times {10}^{-3}$	$\left({4.2}_{-0.6}^{+0.7}\right)\times {10}^{-5}$
Ne ix f	0.9024	$\le 6\times {10}^{-3}$	–	$5.8\times {10}^{-5}$
Ne ix i	0.9132	$0.041$	–	$2.2\times {10}^{-5}$
Ne ix r	0.9191	$\le 9\times {10}^{-3}$	–	$2.4\times {10}^{-5}$
Ne x Ly $\alpha$	$1.0192\pm 0.0006$	$\le 9\times {10}^{-3}$	–	$\left({2.7}_{-1.1}^{+0.5}\right)\times {10}^{-5}$
Continuum	$\Gamma$	norm [Ph keV-1 s-1 cm-2]
	1.72∗	$\left(1.33\pm 0.05\right)\times {10}^{-3\ast }$

Equations

$$\ell =4\pi \frac{m_{\mathrm{p}}{L}_{\mathrm{int}}}{d{m}_{\mathrm{e}}{L}_{\mathrm{Edd}}},$$

$$R=3.6\times 10^{17}\frac{M_{\mathrm{BH}}}{10^7{M}_{\text{astrosun}}}{\left(\frac{n_{\mathrm{H}}}{10^9{\mathrm{cm}}^{-3}}\right)}^2{\left(\frac{\Delta t}{2\mathrm{d}}\right)}^2{\left(\frac{N_{\mathrm{H}}}{10^{22}{\mathrm{cm}}^{-2}}\right)}^{-2}\mathrm{cm},$$

$$R_{\mathrm{d}}=0.13\mathrm{pc}{\left(\frac{L_{\mathrm{bol}}}{10^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.5}{\left(\frac{T_{\mathrm{d}}}{1500\mathrm{K}}\right)}^{-2.6},$$

Full text

1 Introduction

Active galactic nuclei (AGN) efficiently return energy to their environment both via accretion onto supermassive black holes (SMBHs) resulting in broad-band radiation and via the ejection of matter through collimated jets and outflows. Jets can form by extracting energy from the inner edge of an accretion disk around a rotating Kerr black hole (Blandford & Znajek, 1977; Tchekhovskoy et al., 2011) while a maximally spinning black hole allows the innermost stable circular orbit (ISCO) of the disk to lie very close to the black hole. A hot electron plasma of yet unclear geometry, the so-called “corona”, is supposed to upscatter soft thermal seed photons from the accretion disk (Haardt, 1993; Dove et al., 1997a, b; Belmont et al., 2008) approximately resulting in an X-ray power-law spectrum that can be observed both directly and reflected off the disk. For radio-quiet AGN in the absence of strong jets, the X-ray spectrum is entirely dominated by emission processes from the compact X-ray-emitting regions. The prominent fluorescent Fe K$\alpha$ line at a rest-frame energy of 6.4 keV as well as a broad Compton hump peaking around 20–30 keV are well explained by distant reflection (George & Fabian, 1991; García et al., 2013) off a standard optically thick and geometrically thin $\alpha$-disk (Shakura & Sunyaev, 1973). Many “bare” and unabsorbed Seyfert 1 galaxies that are observed under shallow inclination angles are supposed to allow a direct line-of-sight to the inner disk and reflected radiation from this region, where strong relativistic Doppler shifts, gravitational redshifts, and light bending result in a significantly broadened iron line feature (e.g., Walton et al., 2013, and references therein).

The diagnostic power of broad iron lines for studying the accretion physics in the presence of strong gravity has been used in a number of models that are designed for fitting data in XSPEC and other fitting packages. A model that allows to fit for the black hole spin as a free parameter and to convolve a broad-band reflection continuum with a relativistic kernel is provided with relconv by Dauser et al. (2010), who also include an overview of previous models in this field. A number of studies using any of these models consistently find larger disk emissivities closer to the black hole, while the outer disk can be well described with a power-law emissivity law of an $\alpha$-disk (e.g., Wilms et al., 2001; Fabian et al., 2002; Brenneman & Reynolds, 2006; Larsson et al., 2007, for MCG$-$6-30-15 as well as a number of other AGN studied, e.g., by Ponti et al. 2010, Brenneman et al. 2011, Wilkins & Fabian 2011, Dauser et al. 2012 or Risaliti et al. 2013). These steep inner emissivities strongly motivate a “lamp-post” geometry for the corona as opposed to a corona that sandwiches the accretion disk (Svoboda et al., 2012, and references therein). In the lamp-post geometry (Martocchia & Matt, 1996), the source of hard photons that irradiate the accretion disk is located above the black hole, on the symmetry axis of the system. A possible physical realization would be that the corona is the base of a jet (Markoff et al., 2005; Wilkins & Gallo, 2015; King et al., 2017). The first model allowing to directly fit the broad Fe K $\alpha$ line in a self-consistent lamp-post geometry has been published by Dauser et al. (2013). As a further improvement, García et al. (2014) link the reflected spectrum xillver (García et al., 2013) from each point of an ionized disk with the correct relativistic transfer function. Their model relxill_lp is the only model to date that allows self-consistent fits of reflection features under consideration of a primary source at a certain height above the disk in an angle-dependent way. A coverage of high signal-to-noise (S/N) data above 10 keV, as provided by NuSTAR, is important in such models where the model parameters are very sensitive to changes in the spectral shape (Walton et al., 2016; Dauser et al., 2016). A dedicated study of the sample of bare AGN (Walton et al., 2013) with relxill_lp will follow by Fink et al. (2017, in prep.).

Many continua of bare AGN, where most of the relativistically blurred features have been detected, are free of strong and neutral absorption but show signs of ionized warm absorption or outflows. Still, Compton-thick absorption has been claimed to explain the broad iron-line features of AGN (e.g., Miller et al., 2008; Turner et al., 2009). Walton et al. (2010) provide strong arguments against this scenario and in favor of inner-disk reflection.

The only intermediate-class Seyfert galaxy where a variable cold and clumpy absorber (Risaliti et al., 2005, 2007, 2009b; Maiolino et al., 2010) has been observed in conjunction with clear evidence for relativistic reflection is NGC 1365 (Risaliti et al., 2009a, 2013; Brenneman et al., 2013; Walton et al., 2013, 2014). For this source, Risaliti et al. (2009a) are able to disentangle a partial-covering absorber with low covering fraction from similarly broad spectral features of blurred reflection by making use of the variability of the absorber. Mrk 766 may be an additional example, however, the detection of a relativistically broadened iron line is not yet clear (Miller et al., 2007; Patrick et al., 2012).

In Seyfert galaxies, variability in line-of-sight absorption across a wide range of timescales has been observed. One common interpretation of such variability is the passage of discrete clouds across the line-of-sight (Nenkova et al., 2008a, b). Clouds in Cen A (Rivers et al., 2011) and NGC 3227 (Lamer et al., 2003) are observed via centrally peaked column-density absorption profiles on timescales of weeks to months and are inferred to reside in the outer broad line region (BLR) or inner dusty torus. Meanwhile, much shorter absorption events ($\lesssim 1-3$ d) have been detected, for example, in Mrk 766 (Risaliti et al., 2011), Fairall 9 (Lohfink et al., 2012), NGC 1365 (Risaliti et al., 2007, 2009b), Swift J2127.4+5654 (Sanfrutos et al., 2013), or NGC 3227 (Beuchert et al., 2015). These patterns are consistent with transiting clumps at the distance of the BLR.

NGC 4151 is a close-by ($z=0.003319$; de Vaucouleurs et al. 1991) and well-studied Seyfert 1.5 galaxy. Complex and variable line-of-sight absorption has been reported from X-ray spectra for over three decades (e.g., Holt et al., 1980; Yaqoob et al., 1989; Fiore et al., 1990). The absorption has been modeled in a variety of ways across various X-ray missions; successful models have typically incorporated combinations of approximately two absorbers, for example, sometimes cold and/or warm components. Having at least one partial-covering (covering fraction typically $\sim 30$–70%) component is common. Variability in absorption structure has been observed across a wide range of timescales, but positively identifying the responsible component(s) remains difficult (de Rosa et al., 2007; Puccetti et al., 2007; Wang et al., 2010).

Keck et al. (2015) have recently provided solid evidence for relativistic reflection off the inner disk in NGC 4151 with additional signs for absorption variability. The goal of this paper is to revisit the spectral modeling of NGC 4151 and to pursue a physically-motivated model of its broad-band X-ray spectrum and spectral variability. We consider multiple datasets with high count statistics. We use the most updated relativistic reflection model code, relxillCp_lp, and remain mindful of potential degeneracies between relativistic reflection and complex absorption (e.g., Keck et al., 2015, for NGC 4151). The choice of this model is strongly motivated by the detection of non-relativistic radio jets at a velocity of \sim$$0.05\,c (Wilson & Ulvestad, 1982; Mundell et al., 2003; Ulvestad et al., 2005) and by related Comptonization models in the jet-base of microquasars (Markoff & Nowak, 2004; Markoff et al., 2005). In Sect. 2 we provide an overview of the observations we consider for the data analysis presented in Sect. 3. For the analysis, we carefully motivate a baseline model, which we then apply to all observations. We also investigate the inherent spectral variability as well as the sensitive parameters of the relativistic reflection components. The soft X-rays are separately investigated. The results and implications are discussed in Sect. 4 and concluded in Sect. 5.

2 Observations and data reduction

The observing log in Table 1 lists all examined XMM-Newton, Suzaku, and NuSTAR observations. They divide into the long-look Suzaku observation Suz 3 with the simultaneous long-look NuSTAR observation NuSuz. These observations are complemented with two additional Suzaku observations of significant exposure, Suz 1 and Suz 2, from roughly one and two years before the joint Suzaku/NuSTAR campaign. We also add a number of shorter XMM-Newton  observations in order to probe variability on timescales from months to years.

2.1 Suzaku

In the analysis of the Suzaku observations we make use of data from the X-ray Imaging Spectrometer (XIS; Koyama et al., 2007) and the Hard X-ray Detector (HXD; Takahashi et al., 2007). In particular, we use the front (XIS 0,3) and back-illuminated (XIS1) chips in the $3\times 3$ and $5\times 5$ editing modes and reprocess the unfiltered event lists by applying standard procedures for event file screening and attitude correction using aeattcor2. We use the calibration releases 2015-10-05 for Suz 1 and Suz 2 and the later version from 2016-02-04 for Suz 3. A previous investigation of Suz 3 by Keck et al. (2015) has shown discrepancies between the unfolded spectra among the different XIS detectors below 2.5 keV. With the calibration from 2016-02-04, this effect seems to be reduced. We therefore consider data in the full range between 0.6–10 keV. We find that both releases result in identical event files. For Suzaku/HXD, we use the latest calibration release from 2011-09-15. Spectra of both modes and all detectors are extracted from circular regions of $\sim$90″ radius, centered on the point source. We exclude pixels above a threshold of 4% pile-up as estimated with the tool pileest. The resulting spectra of the $3\times 3$ and $5\times 5$ modes are merged using phaadd for each XIS. In order to guarantee sufficient statistics, all spectra are binned to a minimal S/N of 10 but at least 11 channels per bin in presence of spectral lines and at least 20 channels per bin for the continuum. The size of each bin is larger than the resolution of the response grid at $\sim$6 keV111https://heasarc.gsfc.nasa.gov/docs/astroe/prop_tools/suzaku_td/node10.html. Calibration uncertainties are prominent around the Si and Au K edges at the energies of 1.8 keV and $\sim$2.2 keV, leading us to exclude counts from the energy intervals 1.72–1.88 keV and 2.19–2.37 keV.

Figure 1 shows the hardness ratio evolution for the 55 ks and 150 ks observations Suz 2 and Suz 3, respectively, revealing significant variability.

We use the Bayesian block analysis (Scargle et al., 1998, 2013) to divide both observations Suz 2 and Suz 3 into single sub-spectra (Suz 2A, Suz 2B and Suz 3A–Suz 3D) that are highlighted in color in the figure. This method uses Bayesian statistics to identify time intervals that are, given the Poissonian likelihood, compatible with constant hardness ratio. We extract observations from these blocks with the same procedures as described above. The exposures divide into 25.0, 25.6, 55.9, and 43.7 ks for the observations Suz 3A–Suz 3D resulting in $27.3\times 10^{4}$, $19.0\times 10^{4}$, $33.3\times 10^{4}$ , and $30.9\times 10^{4}$ counts for the combination (XIS 0+XIS 1+XIS 3). Similarly, we find 23.0 and 32.1 ks for the observations Suz 2A and Suz 2B with $23.1\times 10^{4}$ and $37.2\times 10^{4}$ counts.

The non-imaging HXD-PIN data are extracted for the whole $34\arcmin\times 34\arcmin$ field of view. The HXD data are binned to a minimum S/N of 40, and a S/N of 20 for resolved sub-spectra. Where simultaneous NuSTAR data are present, cross-calibration constants are fitted relative to the NuSTAR focal plane module A (FPM A) and compared to those found by Madsen et al. (2015). Otherwise spectra are normalized relative to XIS0222http://heasarc.gsfc.nasa.gov/docs/suzaku/analysis/abc/.

2.2 NuSTAR

NuSTAR (Harrison et al., 2013) is the first instrument to focus hard X-rays. X-rays are focused on the two focal-plane modules FPM A and FPM B. We extract data from within 3–78 keV with the standard NuSTAR Data Analysis Software NuSTARDAS-v.1.5.1, which is part of HEASOFT-v.6.18. Due to the high flux of NGC 4151, we extract source counts from within a relatively large region of 90″ radius on both chips FPM A and FPM B, and background counts from a region of the same size located $\sim$340″ off-source but close enough not to introduce much bias due to the spatial dependence of the background (Wik et al., 2014). We explicitly extract NuSTAR data from a time interval that is fully simultaneous to Suz 3 (NuSuz) and from multiple intervals that are subsets of Suz 3A–Suz 3D (NuSuz,A–NuSuz,D) with exposures of 106 ks, 2.5 ks, 25.8 ks, 50.2 ks, and 27.6 ks, respectively. This translates to a range of 3–$62\times 10^{4}$ counts for the observations NuSuz,A to NuSuz,D. Spectra are binned to a minimum S/N of 100 for the integrated spectrum and to 20 for the four individual spectra, which only leaves data below 50 keV. Due to irregularities in the cross-calibration between XMM-Newton, Suzaku, and NuSTAR, only data taken above 5 keV are considered for NuSTAR.

2.3 XMM-Newton

There exist nine observations of NGC 4151 with the EPIC camera (Strüder et al., 2001; Turner et al., 2001) between 2011-05-11 and 2012-12-10. In addition, the XMM-Newton observation with the ID 679780501 has to be excluded due to poor data quality. All observations are in small-window mode. We follow the standard procedure to extract data of all detectors (EPIC-pn, EPIC-MOS1, EPIC-MOS2, and RGS) using the SAS v.14 and the most recent calibration files. After creating calibrated event lists with filtered hot and bad pixels, events between 10 and 12 keV are screened for particle flaring with a threshold of $0.4\,\mathrm{cnts}\,\mathrm{s}^{-1}$. We extract all counts within a maximum possible radius of 40–43″ for EPIC-pn, and within $\sim$120″ for EPIC-MOS. While the region size is physically limited by the chip border in the case of EPIC-pn, we are able to extract counts from nearly 100% of the encircled energy fraction of the on-axis PSF for EPIC-MOS. Background counts are extracted from an off-source spot on the chip within 45″ and 89″ for EPIC-pn and EPIC-MOS, respectively. We detect significant pileup in all three cameras. In the case of the MOS, we exclude the central pixels within 20″ of the source position in all cases. For the EPIC-pn, we exclude data from the inner 15″ for all observations, except for 0679780201, where 17.5″ needed to be excluded due to the higher count-rate of 8.7 cnts s*-1*. The EPIC-pn data are binned by a factor of 2 between 0.5–1.0 keV, 4 between 1–5 keV, 6 between 5–8 keV, and 10 above that. The EPIC-MOS data are binned to a minimum S/N of 10, with additional geometrical binning of 3, 5, and 12 channels/bin in the 0.5–1.0, 1.0–3.0, and 3.0–10 keV bands, respectively. This choice guarantees at least 20–25 cnts bin*-1* and provides an optimal trade-off between a decently binned continuum and sufficient data bins around line features in the spectrum. The first XMM-Newton observation, XMM 1, will not be further used due to strong particle flaring. For XMM 2, we exclude the EPIC-MOS data due to the small amount of net-exposure after filtering the event-files.

From the RGS (den Herder et al., 2000) data, we extract both orders with the task rgsproc and combine the spectra of the detectors RGS 1 and RGS 2 for each order with the task rgscombine. When fitting the RGS data, we simultaneously include both orders and all considered observations. We choose a geometrical binning of a factor of 3 for individual spectra to limit the oversampling of the theoretical RGS energy resolution as suggested by Kaastra & Bleeker (2016). Given the lack of sufficient counts per bin, we choose Cash-statistics for the further data analysis. Due to the low effective area at short wavelengths, we only consider data below 1.3 keV (9.5Å).

3 X-Ray spectral analysis

In the following, we examine all XMM-Newton, Suzaku and NuSTAR  observations. These observations were taken over a period of more than one year. We provide a detailed investigation of the spectral components as well as their variability. We apply the Galactic column of $N_{\mathrm{H,Gal}}=2.3\times 10^{20}\,\mathrm{cm}^{-2}$ (Kalberla et al., 2005) in all cases. In Fig. 2, we show the spectra of all observations in the top panel. The data imply a lack of obvious variability of the soft X-rays below $\sim$1 keV. The source is moderately variable above 6 keV within the range of 1–$1.5\times 10^{-3}\,\mathrm{Photons}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{keV}^{-1}$ at 10 keV with XMM 5 catching the source in an exceptionally low flux state. Strong spectral variability, in contrast, is apparent for the range between $\sim$1 and 6 keV, both in spectral shape and normalization. This argues strongly against a physical partial coverer to explain the continuum including the soft emission below $\sim$1 keV (see earlier work by Holt et al., 1980; Perola et al., 1986; Fiore et al., 1990) but in favor of diffuse emission from optically thin and ionized gas on larger spatial scales (e.g., Wang et al., 2011a, and references therein).

3.1 Motivation for a blurred-reflection component

As a first approach to model the X-ray data, we simultaneously fit all XMM-Newton and Suzaku spectra with an incident and reflected power law (xillver; García et al., 2013, reflection off an optically thick, geometrically thin, ionized accretion disk) that is absorbed by near-neutral intrinsic material (zxipcf) and Galactic foreground gas (tbnew). We fit the diffuse emission below 1 keV phenomenologically with a partial-covering version of zxipcf with a covering fraction of $\sim$98%. We emphasize that the partial coverer has no physical meaning in this case and reflects the combination of a fully absorbed with an attenuated power law, both of identical slope. A physical interpretation of the latter would be nuclear emission scattered off distant and large-scale gas. Figure 2 shows the spectra with corresponding fits. The modeling confirms a strongly changing column density ranging between 2.6 and $14.0\times 10^{22}\,\mathrm{cm}^{-2}$ as well as a variable power-law flux. The ionization states of the reflecting material are found to lie between $\log\xi\sim-3$ and $-1.1$.

For simplicity, we freeze the cutoff energy of the incident continuum to $E_{\mathrm{cut}}=300\,\mathrm{keV}$ and the inclination angle to $i=30^{\circ}$. The Galactic column is considered as foreground absorption in all spectral fits. The residuals indicate decent fits to the continua of the individual observations and the ionized reflection component accounts well for a narrow Fe K$\alpha$ and Fe K$\beta$ line contribution but leaves line-like residuals below 2 keV, strengthening the notion of line-emitting plasma. Broad line-like residuals between $\sim$5 and 6.4 keV, and a hard excess above 8 keV are reminiscent of an additional component of blurred reflection. See also Fig. 3 for a zoom into these features.

3.2 Building on a recent investigation of the Suzaku data:

application of an improved model for relativistic reflection

Keck et al. (2015) have presented a detailed study of the joint Suzaku/NuSTAR campaign (Suz 3/$\mathrm{Nu}_{\mathrm{Suz}})$ for data above 2.5 keV. They discuss two possible models to explain the data: (1) their best-fit model dominated by relativistic inner-disk reflection and (2) an entirely absorption-dominated model. The first is given by a convolution of the unblurred and initially angle-resolved reflection continuum xillver with the relativistic code relconv (Dauser et al., 2010) describing reflection off the inner parts of an accretion disk. The authors also test for the more self-consistent model relxill (García et al., 2013), which links the relativistic transfer-function and the angle-resolved disk-reflection spectrum xillver at each point of the disk. Although García et al. (2013) predict deviations of up to 20% for parameters of relconv and relxill, Keck et al. (2015) find a negligible statistical difference. They attempt to fit for a lamp-post geometry with the convolution code relconv_lp (Dauser et al., 2013) and require two distinct lamp-post components at different heights ($\sim 1.3\,r_{\mathrm{EH}}$ and $\sim 14\,r_{\mathrm{EH}}$, where $r_{\mathrm{EH}}$ is the event horizon). The authors, however, find significant S-shaped residuals between 3–5 keV and therefore reject this solution on statistical grounds. The immediate aim of our study is to investigate this model description using a lamp-post geometry that is physically motivated not only by steep inner emissivities (Svoboda et al., 2012) but also by independent reverberation studies (e.g., Kara et al., 2013, and references therein). As a first improvement, we replace the convolution model with the self-consistent and fully angle-resolved version relxill_lp and adopt the parameters found by Keck et al. (2015). In Fig. 4, we show this model evaluated for the data of Suz 3 and NuSuz in the same energy range considered by Keck et al. (2015). We find similar statistics and residuals as compared to the convolution model. The S-shaped residuals around the turnover of the absorber between 3 and 5 keV imply a yet unmodeled partial coverer, which the authors investigated as part of their independent model (2).

This model combines two partial coverers, one with a column of $\sim 6\times 10^{23}\,\mathrm{cm}^{-2}$ and low covering fraction ($\sim$40%) and a second one with $\sim 1.3\times 10^{23}\,\mathrm{cm}^{-2}$ and a near-maximum covering fraction ($\sim$94%). The model provides a decent description of the continuum above 2.5 keV, but leaves residuals reminiscent of a broad iron line and is statistically less preferred. This demonstrates that modeling the complex spectrum of NGC 4151 is not straightforward but combined with degeneracies between these two solutions.

In the following we build a baseline model using the long-look 150 ks Suz 3 and simultaneous NuSuz observations. We extend the work by Keck et al. (2015) by considering the entire energy range covered by the instruments. In order to flatten the S-shaped residuals in our broad-band X-ray continuum, we combine their models (1) and (2) and apply the latest self-consistent relativistic reflection code relxillCp_lp333The model relxillCp_lp extends on relxill_lp (Dauser et al., 2013) taking into account the Comptonization continuum nthcomp as primary continuum., which is, similar to the latest version of xillver, calculated using the Comptonization continuum nthcomp (Zdziarski et al., 1996). In the following, we fix the seed photon temperature to 50 eV, which is used to calculate the xillver tables.

3.3 A re-investigation of the 150 ks simultaneous Suzaku  and NuSTAR observations with relxillCp_lp

3.3.1 Derivation of a baseline model

The residuals at each step are shown in Fig. 5. First, we fit the data with the Comptonization continuum nthcomp, hereafter referred to as HXCOMP. It provides a more physically motivated primary continuum with an intrinsic cutoff as opposed to a power law with external cutoff (García et al., 2015). This continuum is absorbed by fully-covering neutral gas (CA, tbnew_simple_z). The residuals are shown in Fig. 5a and in particular indicate an unmodeled soft continuum.

In previous studies using Chandra data, the soft continuum was accounted for with a bremsstrahlung component (Ogle et al., 2000; Wang et al., 2011a). The latter authors, however, emphasize their lack of a physical motivation for this component. We therefore adopt a simple scenario, in which the nuclear Comptonized continuum is scattered off distant and large-scale gas and model this component with an unabsorbed soft Comptonization continuum (SXCOMP, Fig. 5b). Here, we adopt the same parameters as for the HXCOMP but leave the normalization free to vary. This way of modeling introduces a minimal set of additional degrees of freedom. Also, at CCD energy resolution and given the blend of emission lines present (see below), we cannot constrain the exact form of the soft continuum. Fits using Comptonization or bremsstrahlung continua yield statistically identical fits.

On top of the soft continuum, a number of lines appear in the residuals, reminiscent of emission from the ionized large-scale gas component that has been extensively studied by Wang et al. (2011a), for example, using Chandra and physical emission codes. We focus mainly on the nuclear properties of the X-ray spectrum and instead fit this emission with a phenomenological blend of Gaussians (Fig. 5c). The centroid energies are adopted from those lines that are significantly detected in gratings data of XMM-Newton/RGS (Schurch et al., 2004) and Chandra/LETG (Ogle et al., 2000, see also Vainshtein & Safronova 1978), or as line-blends by Chandra/ACIS (Wang et al., 2011a). Table 2 lists the parameters of all fitted lines with their ion identification, line flux and centroid energies. For deriving the final photon line fluxes and their uncertainties, we freeze all continuum parameters including the SXCOMP normalization to the best-fit value found in the final step. This also reduces degeneracies between the line-blend and the continuum. The remaining residuals show a yet unmodeled swing in the continuum around 2 keV, strong and broad emission forming an iron-line complex at approximately 6 keV and extra curvature above 10 keV.

Narrow emission components at the centroid energies of the Fe K$\alpha\,(\beta)$ lines at $\sim$6.4 (7.1) keV can be modeled with either two narrow Gaussians (Warwick et al., 1989; Zdziarski et al., 2002; Schurch et al., 2003; Wang et al., 2010) with a frozen flux-ratio of 12% (Fig. 5d) or a component of distant reflection (xillver) in Fig. 5e. Both options leave a broad emission feature between 5 keV and 6 keV, as previously seen by Wang et al. (1999), as well as equal continuum residuals. We are therefore not able to confirm a strong statistical need for a distant reflection component. We continue to use two Gaussian components as a phenomenological description of the iron line emission that likely originates from Compton-thin gas. This way, we only add a minimal set of additional degrees of freedom to the model.

The continuum results in S-shaped residuals below 5 keV, similar to what Keck et al. (2015) found. We replace the fully covering absorber (CA) with the combination of a partially covering absorber (CA1) and a fully covering absorber (CA2), both neutral (see Fig. 5f). CA1 requires a covering fraction of $\sim$40–50%. A similar dual neutral absorber has frequently been applied before (e.g., Wang et al., 2010; Keck et al., 2015, and references therein). Further tests reveal that the introduced spectral curvature can not be reproduced with an ionized warm absorber.

The broad pattern just below 6 keV remains, even after flattening the continuum below 5 keV. It likely features the red wing of an extremely blurred Fe K$\alpha$ line (e.g., Dauser et al., 2010). We attempt to account for this feature with the previously introduced model relxillCp_lp. We call this component LP1. Its height hits the lower limit at $\sim 1.1\,r_{\mathrm{EH}}$, which corresponds to $1.2\,r_{\mathrm{g}}$ for the spin fixed at its maximum value. Note that we fit the incident and reflected continua independently at this point. Due to the lack of data above 50 keV, we are unable to constrain the cutoff energy (expressed via the electron temperature in our model) and freeze it at $kT_{\mathrm{e}}=399\,\mathrm{keV}$ or $E_{\mathrm{cut}}\approx 1000\,\mathrm{keV}$. We note, however, that Malizia et al. (2014) find $E_{\mathrm{cut}}=196_{-32}^{+47}\,\mathrm{keV}$ with INTEGRAL, Swift, and XMM-Newton data between 2 and 100 keV. We expect significant bias for their measurement of the cutoff due to: the lack of soft X-rays (García et al., 2015); the use of a lower photon index of $\Gamma\sim$1.63; and the background dominance of INTEGRAL as opposed to the NuSTAR data used in this work. We emphasize that for primary sources close to the black hole, the cutoff energy has to be corrected for the gravitational redshift as outlined by Niedźwiecki et al. (2016). The resulting residuals (Fig. 5g) illustrate that the underlying broad feature as well as parts of the Compton hump above 10 keV are successfully fitted by LP1.

Dips around 7 keV suggest an additional column of highly ionized absorption. We can greatly improve the fit using a XSTAR absorption component for dense, coronal gas with $\log\xi\sim 2.8$. Warm absorber components with similar ionization have previously been found by Weaver et al. (1994b), Schurch & Warwick (2002), and Keck et al. (2015). The XSTAR model removes line-like residuals close to the centroid energies of Fe xxv He $\alpha$, Fe xxvi Ly $\alpha,$ and Fe xxv He $\beta$, at zero velocity offset relative to systemic. There remain narrow absorption-like residuals around 8 keV; they can be modeled with a broad Gaussian absorption component that has an energy centroid of $8.17^{+0.11}_{-0.09}\,\mathrm{keV}$ (rest frame). We identify this feature as a blueshifted Fe xxv He $\alpha$ or Fe xxvi Ly $\alpha$ line, although the possibility of unmodeled contributions from Fe xxvi Ly $\beta$ (8.25 keV rest frame) associated with the systemic warm absorber cannot be ruled out. The implied velocities are $0.22^{+0.02}_{-0.01}\,\mathrm{c}$ for He-like Fe or $0.17^{+0.02}_{-0.01}\,\mathrm{c}$ for H-like Fe, suggesting an ultra-fast outflow (UFO), as previously constrained with a highly ionized XSTAR component (see Tombesi et al., 2011, 2013, for a physical modeling of the UFO with XSTAR). The flattened residuals are shown in Fig. 5h.

An excess around $\sim$6 keV implies a still broad iron line with a much less smeared red wing. We therefore add a second lamp-post component (LP2) with a primary source at the larger height of $\sim 17\,r_{\mathrm{EH}}$, which also takes care of remaining excess residuals above 20 keV. While both components LP1 and LP2 already describe the entire Compton hump above 10 keV, further tests can exclude the need for a third, distant reflection component to model also the narrow iron line. We find a best-fitting baseline model with $\chi^{2}/\mathrm{dof}=1094/762$ and overall flat residuals in Fig. 5i.

3.3.2 The best-fit baseline model

We found a solid model for the broad continuum from the soft to the hard X-rays, which, below, we apply to the remaining observations. This baseline model combines a set of four complex absorbers with a physical description of blurred reflection as part of the self-consistent lamp-post geometry. We fit both LP components with a tied slab-ionization and find a common value of $\log\xi\sim 2.8$. The fit is hitting the lower-limit for the inclination of 3.1∘ with 90% uncertainties allowing values as high as 10∘ (see the contours between disk inclination and spin in Fig. 6 where we unfroze the spin parameter). The wide range originates in the $\cos i$-dependence of the model. This inclination is consistent with the constraint of $\theta<30^{\circ}$ that Cackett et al. (2014) found with reflection-component reverberation mapping. We find no obvious correlation with the spin parameter and fix the spin at its maximum value. The relatively flat turnover between 3 and 6 keV is well described by the dual neutral absorber rather than a single full-covering absorber. The column densities of the two neutral absorbers are found to be $N_{\mathrm{H,1}}\sim 21\times 10^{22}\,\mathrm{cm}^{-2}$ with $f_{\mathrm{cov}}=0.46$ and $N_{\mathrm{H,2}}\sim 9\times 10^{22}\,\mathrm{cm}^{-2}$ for the full-covering absorber. The HXCOMP photon index is well constrained to $\Gamma=1.72\pm 0.01$ for the long-look observation and will therefore be kept fixed for the remaining observations of lower count statistics in order to reduce degeneracies within the complex model composite. We list the best-fit parameters and their uncertainties in Table 3 and show all model components in Fig. 7.

We test for degeneracies inherent to the complex continuum for the case of Suz 3 and study the resulting contours in $\Delta\chi^{2}$-space (Fig. 8). The contours and derived uncertainties suggest that we are able to constrain and clearly separate between both lamp-post continua with primary sources at different heights and a two-fold absorber, composed of a partial-covering and full-covering neutral column. We only observe a tentative correlation between the LP2-normalization and the column density of CA1. This correlation is likely not physical but caused by both models describing a similar spectral shape around 6 keV. We can show that the additional absorber of highly ionized gas is free of degeneracies with the continuum and can be kept separate throughout the analysis.

3.3.3 The reflection fraction and reflection strength

In our best-fit baseline model for Suz 3, we fit the primary continuum independently from the reflection continua LP1 and LP2. The normalizations of the reflection continua have no geometrical interpretation in this case. In other words, this does not allow us to infer the reflection fraction $R_{\mathrm{f}}$, which is intrinsic to the lamp-post geometry and corresponds to the ratio of photons escaping the system towards the observer over the coronal photons incident on the disk (Dauser et al., 2014, 2016). The reflection fraction therefore mainly depends on the effect of light-bending, which is stronger for primary lamp-post sources close to the disk, that is, predicting larger values of reflection fraction $R_{\mathrm{f}}$ (Miniutti et al., 2004). In the following, we use a workaround in order to find access to the intrinsic reflection fraction for both LP1 and LP2. We re-define the baseline model, remove the independent HXCOMP continuum and thaw $R_{\mathrm{f}}$ with refl_frac$>0$. This model describes the same spectral components but now provides two primary continua, instead of one, that are intrinsically linked to the reflected continua via the lamp-post geometry, that is, it allows to directly fit for $R_{\mathrm{f}}$ as a more meaningful parameter.

If we fit for both primary continua and both reflection fractions at a time, we find strong degeneracies and large uncertainties with $N_{\mathrm{LP_{1}}}=0.5^{+6.6}_{-0.4}$ and $N_{\mathrm{LP_{2}}}=\left(1.68^{+0.11}_{-0.84}\right)\times 10^{-3}$ as well as $R_{\mathrm{f}}^{\mathrm{LP_{1}}}=14^{+36}_{-13}$ and $R_{\mathrm{f}}^{\mathrm{LP_{2}}}=0.43^{+49.58}_{-0.05}$. For that reason, we mutually freeze the reflection fraction of one component to the value predicted for the emission of a point source at the previously determined height above the disk (fixReflFrac=2) and fit for the reflection fraction of the other component as well as both normalizations of the primary continua. The results are shown in Table 4. In the left column, LP2 is fixed to a predicted reflection fraction of 1.2 for the given height of $h_{\mathrm{LP_{2}}}=14.1\,r_{\mathrm{EH}}\sim 15.0\,r_{\mathrm{g}}$. A fit of the reflection fraction of LP1 is now well constrained as $R_{\mathrm{f}}^{\mathrm{LP_{1}}}=1.50\pm 0.06$, which is rather low with respect to its low height of $h_{\mathrm{LP_{1}}}=1.1\,r_{\mathrm{EH}}$ and the therefore predicted fraction of 22.5. This is absorbed by a much larger normalization of LP1 as opposed to LP2, owing to strong degeneracies between these parameters. In the second case, we set the reflection fraction of LP1 to $R_{\mathrm{f}}^{\mathrm{LP_{1}}}=22.5$ as predicted for its height $h_{\mathrm{LP_{1}}}=1.1\,r_{\mathrm{EH}}=1.17\,r_{\mathrm{g}}$ and freely fit $R_{\mathrm{f}}^{\mathrm{LP_{2}}}$. We find a reasonably low reflection fraction of $R_{\mathrm{f}}^{\mathrm{LP_{2}}}=0.437\pm 0.023$ as expected for a primary lamp-post source at larger height. Also, the LP-normalizations do not diverge as strongly as in the case before. We can still demonstrate that degeneracies make it challenging to interpret the reflection fraction as a probe of the lamp-post geometry for the case of two interacting lamp-post sources.

While the reflection fraction can not simply be inferred from the observed spectra without knowledge of the geometry (Dauser et al., 2014, 2016), the reflection strength $R_{\mathrm{s}}$ is defined as the strength of the Compton hump of the reflection model with respect to the primary continuum, that is, the flux-ratio of the reflected to the incident continuum in the 20–40 keV energy band. The derived numbers in Table 4 imply a behavior very similar to that observed for the reflection fraction. The reflected LP1-spectrum seems to be too weak with $R_{\mathrm{s}}^{\mathrm{LP_{1}}}=0.69$ compared to the value of 10.5 as predicted by the geometry of a point source very close to the black hole, which is again accounted for by degeneracies between the reflection fraction and the LP-normalizations. In contrast, we find a reasonable value of $R_{\mathrm{s}}^{\mathrm{LP_{2}}}=0.24$ at a larger height, featuring a rather weak Compton hump close to the predicted value of 0.64, predicted for a point source at the given height. The reflection strength of the combination of both lamp-post components is 0.68.

3.4 Test for relativistic reflection in other Suzaku and XMM-Newton observations

We have now derived a robust baseline model based on the average 150 ks observation by Suzaku (Suz 3) and NuSTAR  (NuSuz) that can also be applied to the remaining observations of Suzaku and XMM-Newton. We again use the original model description with one primary continuum (HXCOMP) that is fitted independently of the reflection continua (see Table 3). We re-fit the observations Suz 1 and Suz 2 and freeze all parameters except of the two absorbers CA1 and CA2 and the normalizations of the incident continuum and the narrow Fe K$\alpha$ and K$\beta$ lines. This approach yields overall good fits, except for broad excess residuals in the iron band, indicating extra variability of the reflection components. A re-fit of the normalizations of the lamp-post continua results in strong degeneracies. When fitting with both lamp-post components switched off, we find the gray residuals shown in Fig. 9, arguing for the presence of blurred reflection features via broad features in the iron band. We can, however, demonstrate that the data of Suz 1 and Suz 2 do not allow us to disentangle two emitting sources at different heights.

Given the comparatively large amount of counts for Suz 1 and Suz 2 we attempt to constrain a single lamp-post component with variable height. We find overall good fits (Fig. 9, bottom panels). The statistics improve substantially compared to the model lacking relativistic reflection with $\Delta\chi^{2}=55\,(1)$ for Suz 1 and with $\Delta\chi^{2}=47\,(1)$ for Suz 2. The lamp-post heights are fitted with $24^{+7}_{-9}\,r_{\mathrm{EH}}$ for Suz 1 and $14^{+4}_{-6}\,r_{\mathrm{EH}}$ for Suz 2. Due to parameter degeneracies we can neither state variability of this single LP-component between Suz 1 and Suz 2, nor between the time-resolved spectra of Suz 2. A direct comparison with the complex and well constrained double-lamp-post source in Suz 3 are not possible either. We suggest, however, that the single component fitted to Suz 1 and Suz 2 may be a blend of these components LP1 and LP2.

Similar to the Suzaku data, the XMM-Newton data reveal visible broad features in the iron-band residuals of a fit using the unblurred xillver model (Fig. 3). Due to the strong degeneracies arising for a free lamp-post height in the models of Suz 1 and Suz 2, we freeze the heights of both lamp-post components to those derived for Suz 3 and only fit the normalizations. Table 6 shows that LP1 is undetected by XMM-Newton in contrast to LP2, where we can constrain its normalization to within $\sim$12–60% with the exception of a very good constraint of $\sim$4% for XMM 5. Together with the XMM-Newton data, we report significant variability in normalization for LP2 over time with a minimum timescale of 20–30 d.

These results demonstrate the need for a good number of counts to properly constrain one or even multiple components of relativistically blurred reflection. Otherwise, it is challenging to simultaneously probe the stability of both reflection components, LP1 and LP2, over time.

3.5 Spectral variability probed with Suzaku and XMM-Newton

As we have shown in Sect. 3, NGC 4151 shows significant spectral variability between 1 and 6 keV, on which we will focus in the following. We consider all XMM-Newton, Suzaku, and NuSTAR observations between 2011 and 2012 including the time-resolved observations Suz 2A,B as well as Suz 3A-D and NuSuz,A-D.

To address the variability of the spectral components, we apply the baseline model and allow only a few parameters to vary. Besides the cross-normalizing detector constants and the flux normalizations of the model components, these are $N_{\text{H,1}}$, $N_{\text{H,2}}$, and the covering fraction $f_{\text{cvr}}$ of CA 1. We detect strong degeneracies between $N_{\text{H,1}}$ and $f_{\text{cvr}}$ for XMM-Newton data and therefore fix the covering fraction to the weighted mean with respect to the observations Suz 1, Suz 2, and Suz 3, which all lie very closely to $f_{\mathrm{cvr}}=0.46$. Note that the SXCOMP normalization is kept frozen to the value derived for Suz 3 for all Suzaku observations but allowed to vary for XMM-Newton. The derived values are, however, consistent with the frozen value.

The baseline model fits well to all observations with only the few above mentioned free parameters. All parameters and uncertainties are listed in the Tables 6 and 7. Figure 10 shows all spectra with overlaid fits in the top panel and residuals in the bottom panels. Most previous studies have found the X-ray emission below $\sim$1–2 keV to be non-variable (e.g., Yang et al., 2001; de Rosa et al., 2007); Landt et al. (2015) report on weakly-variable coronal O vii emission. In contrast, Wang et al. (2010) find evidence for significant variability of the soft continuum. Our multiple observations yield 0.6–1.0 keV fluxes that remain within $\sim$6% of each444The amount of variability is calculated as the standard deviation over the average. but underlie large uncertainties. This range contrasts with $\sim$20% variability for the 7–10 keV flux of the incident HXCOMP continuum. We find no direct correlation between the SXCOMP and HXCOMP variability. We can also exclude correlated variability between the Ne ix emission line and the incident 7–10 keV flux. Within the statistical uncertainties, the Ne ix line flux remains constant over the probed timescale. Due to the CCD resolution of our spectra and the blend of emission lines, it is not clear if the weak flux variability within 0.6 and 1.0 keV is due to soft continuum, the Ne ix flux, or both. The narrow Fe K$\alpha$/$\beta$ line, in contrast, is stable over the monitored time interval.

We find significant absorption variability with column density changes of a maximum of 50% for CA1 and 20% for CA2. Spectral variability is not only found between single observations but also within the 150 ks-long observation Suz 3 as we show in Fig. 11. The evolution of the parameters of both neutral absorbers and the incident HXCOMP is shown in Fig. 12. We find overall larger columns $N_{\text{H,1}}$ but smaller uncertainties for the partial coverer than for the full-covering column $N_{\text{H,2}}$. Both columns $N_{\text{H,1}}$ ($N_{\text{H,2}}$) seem to show correlated variability over time with an initial decline around MJD 55900 followed by an increase between MJD 56080 and MJD 56200 and a subsequent decline after MJD 56240 back down to a baseline-level of $N_{\mathrm{H,1}}\,(N_{\mathrm{H,2}})\sim 12\,(9)\times 10^{22}\,\mathrm{cm}^{-2}$. The normalizations of the primary continuum (HXCOMP) can be shown to be variable down to a timescale of $\sim$20 d over a normalization range of $\sim 0.03$–$0.08\,\mathrm{Photons}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$ at 1 keV. The time-resolved measurements for Suz 2 and Suz 3 clearly strengthen the presence of variability of these parameters on timescales as short as days. This remains true even though the uncertainties of the parameters found for Suz 2A,B are larger compared to those inferred from the spectra Suz 3A-D. The covering fraction can be assumed as constant within the uncertainties. The same applies to the time-resolved results of Suz 3A-D, where we find covering fractions scattering around 0.55. This value is slightly larger than that derived for the total observation Suz 3, which can be attributed to model degeneracies. The columns $N_{\text{H,1}}$ and $N_{\text{H,2}}$ peak during the observation Suz 3D and cover a dynamic range of $N_{\mathrm{H,1}}\,(N_{\mathrm{H,2}})\sim 13$–$25\,(6$–$10)\times 10^{22}\,\mathrm{cm}^{-2}$. The primary continuum, in contrast, shows a more complex variability pattern with normalizations between $0.04\,\mathrm{Ph}\,\mathrm{keV}^{-1}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$ and $0.05\,\mathrm{Ph}\,\mathrm{keV}^{-1}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$.

The correlated variability between CA1 and CA2 on all probed timescales can likely be attributed to degeneracies arising between both absorbers (Fig. 13) that reveal similar columns. We probably also observe systematics related to the complexity of the dual absorber and to the assumptions made for the baseline model. According to these results, we can therefore not claim the observed variability to be inherent in one or the other absorber for most of the analyzed observations. In contrast, both contributing absorbers are likely well separated by the single 150 ks-long observation Suz 3 (see Fig. 8). We further observe an anticorrelation between the 7–10 keV flux of the pre-absorbed HXCOMP continuum and both columns $N_{\text{H,1}}$ and $N_{\text{H,2}}$, which is highlighted in Fig. 14. This relation is stronger for CA1 with the Pearson correlation coefficient $r_{\mathrm{P,CA_{1}}}=77$% and a very low p-value of $P_{\mathrm{P,CA_{1}}}=0.8$% as compared to CA2 with $r_{\mathrm{P,CA_{2}}}=60$% and $P_{\mathrm{P,CA_{2}}}=7.2$%.

3.6 The soft X-rays and the NLR

NGC 4151 and its ionized environment have been subject to a number of studies that combine high spectral and spatial resolution (e.g., Wang et al., 2011a, and references therein). In Sect. 3.3, we motivate the description of the soft X-rays below 2 keV with the SXCOMP continuum, which is complemented with a blend of Gaussian emission lines that have been measured to persist over decades (see also Ogle et al., 2000; Yang et al., 2001; Schurch et al., 2004; Wang et al., 2011a). While Perola et al. (1986), Weaver et al. (1994a), Weaver et al. (1994b), Warwick et al. (1995) and Wang et al. (2010) detect no signs for major variability in the soft X-rays, we find tentative signs for variability at a low dynamic range in Sect. 3.5. As the S/N of individual RGS spectra are too low for measuring line fluxes over time, we combine the XMM-Newton/RGS data of all observations ranging over more than 1.5 years. The resulting spectrum reveals a series of highly significant emission lines, dominated by the H-like and He-like ions O viii, O vii , and Ne ix (Fig. 15) on top of the weak SXCOMP continuum. The emission lines are modeled with Voigt profiles. Whenever the S/N of a line-feature is too low to constrain its centroid energy, we use the values from Vainshtein & Safronova (1978) to set the starting parameters. Additional radiative recombination continua (RRC) are described with the redge model. The model is simultaneously fitted to all observations with the resulting parameters listed in Table 8. The centroid energies of all fitted lines are consistent with those derived earlier, for example, by Ogle et al. (2000) or Schurch et al. (2004). The poor constraints on the width and the turbulent velocities arise due to the choice of the Voigt profile as the physically motivated description of the emission lines, coming along with more free parameters than a Gaussian approximation. Despite these increased uncertainties, the Voigt profile shape is statistically required by the most prominent lines. In contrast, we can roughly constrain the plasma temperatures of the prominent RRC features.

A description of the soft X-rays using physical emission codes for optically thin and ionized gas is beyond the scope of this paper, whose main focus lies on a solid description of the nuclear continuum. Also, the quality of the XMM-Newton/RGS data does not allow for improvement of results from previous studies with Chandra. We therefore remain with a phenomenological modeling of the soft X-rays and provide a straight-forward gas diagnostic based on the two prominent line triplets of O vii and Ne ix at $\sim$0.56 keV (22 Å) and $\sim$0.91 keV (13.5 Å), respectively (Fig. 16). Both lines are statistically best described by Voigt profiles. We list the parameter constraints in Table 5. In addition to the statistical uncertainties, we also assume non-vanishing systematical uncertainties for the fit of the hardly detected resonance and intercombination lines of Ne ix.

Both triplets consist of a resonance ($w$), an intercombination ($x+y$) and a forbidden line ($z$), while the intercombination line is a blend of the lines $x$ and $y$. Relating the strengths of these lines can be used as a powerful diagnostic for the density and temperature of the emitting gas (Gabriel & Jordan, 1969; Porquet & Dubau, 2000; Bautista & Kallman, 2000). The ratio $R=z/(x+y)$ is sensitive to the gas density in that the rate of collisions increases with $n^{2}$, which suppresses the forbidden line emission. The ratio $G=(z+x+y)/w$, in turn, is sensitive to the temperature, which positions the gas between being dominated by collisions or recombination. A hot plasma gives rise to collisions and therefore a strong resonance line, that is, small values of $G$. For larger values of $G\gtrsim 4$, the plasma is dominated by recombination and the triplet levels with the intercombination and forbidden line have large statistical weight, therefore featuring a plasma dominated by photoionization. For both He-like ions, the ratios are well consistent with those derived using XMM-Newton/RGS by Armentrout et al. (2007). For O vii, the ratios ($G=4.05\pm 0.37$ and $R=5.8\pm 1.0$) indicate gas dominated by photoionization. These ratios imply a low gas temperature of $T_{\mathrm{e}}\lesssim 10^{5}\,\mathrm{K}$, which is consistent with that obtained from the narrow RRC features. In contrast to Schurch et al. (2004), who find a moderate $R$-ratio of 3.9, our value is more consistent with that measured by Landt et al. (2015), arguing for a gas of very low density with $n_{\mathrm{e}}\sim 10^{3}\,\mathrm{cm}^{-3}$. We detect a relatively strong resonance line for Ne ix, resulting in a low $G$-ratio of $3.333\pm 1.029$ and $R$-ratio of $2.6\pm 0.9$. This result implies a hybrid plasma (e.g., Porquet & Dubau, 2000; Bautista & Kallman, 2000) that is likely in pressure equilibrium with a collisional ionization gas phase (Wang et al., 2011a). The temperature ($\sim 10^{6}\,K$) and density (a few $\times 10^{11}\,\mathrm{cm}^{-3}$) are estimated to be significantly higher compared to the O vii gas. Again, the results are overall consistent with earlier work (Schurch et al., 2004; Armentrout et al., 2007; Ogle et al., 2000; Wang et al., 2011a). The suppression of the forbidden line of Ne ix and therefore the $R$-ratio can alternatively arise from a strong UV field, where photoionization is still dominant (Mewe & Schrijver, 1978).

4 Discussion

Extending earlier work by Keck et al. (2015) by including soft X-ray data, in this paper we re-investigated data from the joint Suzaku/NuSTAR campaign of NGC 4151. The resulting broad energy coverage and the use of the lamp-post model relxillCP_lp allowed us to disentangle complex and variable absorption from blurred disk reflection originating in the regime of strong gravity. We described the blurred reflection with two lamp-post components at different heights, while the bulk of the 1–6 keV continuum absorption was attributed to an absorber consisting of partial- and full-absorbing components. We also detected narrow absorption features that can be modeled with a highly ionized ($\log\xi\sim 2.8$) warm absorber and an additional absorption line indicating an UFO. The soft emission was modeled with a SXCOMP scattering continuum and a number of emission lines that are identified with XMM-Newton/RGS.

4.1 Size of the corona and pair production

Our baseline model encompasses two point-like primary lamp-post sources at the heights of $1.2\,r_{\mathrm{g}}$ and $15.0\,r_{\mathrm{g}}$, which result in strongly and moderately blurred spectral reflection components (see also Keck et al., 2015). Individual studies, for example, by Dauser et al. (2012), Parker et al. (2014) and Fink et al. 2017, in prep., that apply the model relxill_lp, confirm the need for relatively low source heights between $\sim$2 and 4$r_{\mathrm{g}}$ and compact coronae in order to explain the observed spectra. This has been independently confirmed by simulations of Svoboda et al. (2012) as well as reverberation studies (Fabian et al., 2009; Zoghbi et al., 2010, for 1H 0707$-$495). In the following we discuss the implications and limitations of this solution.

For the moment, we assume two distinct emission regions as they result from our modeling. Dovčiak & Done (2016) provide estimates on the spatial extent of spherical coronae depending on their height above the BH. For an observed luminosity of $L_{\mathrm{X,obs}}^{0.3-10\,\mathrm{keV}}=0.001\,L_{\mathrm{Edd}}$ and a photon index of $\Gamma=2$, the extent of a corona at a height of 10–$20\,r_{\mathrm{g}}$ grows to values larger than $1\,r_{\mathrm{g}}$. The size is decreasing to $\sim 0.4\,r_{\mathrm{g}}$ at lower heights and again increasing to $\sim 1\,r_{\mathrm{g}}$ for a height close to $1\,r_{\mathrm{g}}$. This evolution of coronal sizes is caused by light-bending close to the BH and gravitational redshifts at larger heights. The models fitted to the XMM-Newton, Suzaku, and NuSTAR data examined in this paper imply the incident power law to have a mean observed luminosity of $L_{\mathrm{X,obs}}^{0.3-10\,\mathrm{keV}}/L_{\mathrm{Edd}}=0.001$–0.002, which is comparable to the luminosity assumed by Dovčiak & Done (2016). We therefore adopt a coronal radius of $1\,r_{\mathrm{g}}$ for our components. If the corona of LP1 was in fact situated at $1.2\,r_{\mathrm{g}}$ above the BH, its radius would have to be even smaller in order not to interfere with the event horizon.

The compactness of a corona with radius $d$ can be expressed as

[TABLE]

where $L_{\mathrm{int}}$ is the intrinsically emitted luminosity, $L_{\mathrm{Edd}}\sim 1.5\times 10^{38}\,(M/M_{\mathrm{\astrosun}})\,\mathrm{erg}\,\mathrm{s}^{-1}=4.5\times 10^{45}\,\mathrm{erg}\,\mathrm{s}^{-1}$ is the Eddington luminosity555The Eddington luminosity is calculated with the black hole mass $M_{\mathrm{BH}}=3\times 10^{7}\,M_{\mathrm{\astrosun}}$ (Hicks & Malkan, 2008). and $d$ is the radius of a simplified spherical corona in units of $r_{\mathrm{g}}$. The compactness can be interpreted as the optical depth of the corona with respect to pair production and has been extensively discussed by Lightman & Zdziarski (1987), Svensson (1987), Dove et al. (1997a) and Dove et al. (1997b). A useful parameter in that regard is $\Theta=k\,T_{\mathrm{e}}/m_{\mathrm{e}}\,c^{2}\sim E_{\mathrm{cut}}/2\,m_{\mathrm{e}}\,c^{2}$ (García et al., 2015), where $E_{\mathrm{cut}}$ is the high-energy cutoff of the Comptonized continuum. The cross-section for electron/positron pair production per photon peaks between $E=m_{\mathrm{e}}\,c^{2}=511\,\mathrm{keV}$ and $\sim$1 MeV. We therefore evaluate the intrinsic luminosity for the model of each observation both between 0.1 and 200 keV (the energy range used by Fabian et al. 2015) and 0.1 and 1000 keV, which includes the peak of the cross-section.

If we replace the model description of our best-fit baseline model (Table 3) with two separate incident continua and fit for both normalizations ($N_{\mathrm{LP_{1}}}$, $N_{\mathrm{LP_{2}}}$) and reflection fractions ($R_{\mathrm{f}}^{\mathrm{LP_{1}}}$, $R_{\mathrm{f}}^{\mathrm{LP_{2}}}$), these four parameters are strongly degenerate (see Sect. 3.3.3 for the numbers). The dynamic range of the normalizations comprises also the solutions of Table 4 and therefore provides the full range of allowed values in a conservative way. The observed source luminosities therefore range between $3.2\times 10^{41}\,\mathrm{erg}\,\mathrm{s}^{-1}\lesssim L_{\mathrm{obs,0.1-200\,keV}}^{\mathrm{LP_{1}}}\lesssim 2.54\times 10^{43}\,\mathrm{erg}\,\mathrm{s}^{-1}$ and $1.17\times 10^{43}\,\mathrm{erg}\,\mathrm{s}^{-1}\lesssim L_{\mathrm{obs,0.1-200\,keV}}^{\mathrm{LP_{2}}}\lesssim 2.50\times 10^{43}\,\mathrm{erg}\,\mathrm{s}^{-1}$. Especially in the case of LP1, where the corona would be situated deep inside the gravitational potential, we need to transform the observed luminosities into the intrinsic frame of the corona. The conversion factor $(1+z_{\mathrm{g}})^{\Gamma}$ depends on the gravitational redshift $z_{\mathrm{g}}=\left(1/\sqrt{1-2\,h/(h^{2}+a^{2})}\right)-1$ between the corona and the observer as well as the photon index $\Gamma$. Fixing the spin at $a=0.998$, we find that the conversion factor increases quickly from $\sim$1.14 ($z_{\mathrm{g}}\sim 0.08$) at $h\sim 15.0\,r_{\mathrm{g}}$ to $\sim$37 ($z_{\mathrm{g}}\sim 7.2$) at $h\sim 1.2\,r_{\mathrm{g}}$. In addition, in the case of LP1 we also need to apply a correction factor to account for photon trapping. Due to strong light-bending, at the height of $1.2\,r_{\mathrm{g}}$ only $\sim$1% of the photons of an isotropically emitting primary source reach the observer, while a significant fraction of the photons ($\sim$13%) reach the accretion disk, leading to a large reflection fraction of $R_{\mathrm{f}}^{\mathrm{LP_{1}}}=22$ for LP1 (see also Table 4). At this reflection fraction, the corresponding luminosity becomes $L_{\mathrm{int,0.1-200\,keV}}^{\mathrm{LP_{1}}}=0.4L_{\mathrm{Edd}}$. Leaving $R_{\mathrm{f}}$ free, the incident continuum can reach observed luminosities as large as $2.54\times 10^{43}\,\mathrm{erg}\,\mathrm{s}^{-1}$. For this value, the corona would intrinsically exceed $L_{\mathrm{Edd}}$ by a factor of ten.

The full range of allowed luminosities translates into a compactness-range of $2500\lesssim\ell_{\mathrm{LP_{1}}}\lesssim 200000$ and $5\lesssim\ell_{\mathrm{LP_{2}}}\lesssim 11$. With these constraints, we can show that the compactness parameter $\ell$ would, within the uncertainties, well cover the parameter space below the limits for pair production (Stern et al., 1995; Dove et al., 1997a, b; Fabian et al., 2015) for LP2. It would, however, exceed this limit in some extreme cases of LP1. Here, we also have to correct the parameter $\Theta$ and the cutoff energy for the gravitational energy shift arising between the compact source and the observer. In the case of LP1, this would shift the primary source above the pair limit for all compactness values independent of the value for $E_{\mathrm{cut}}$. As a result, the corona would be entirely optically thick for photons with energies $E>511\,\mathrm{keV}$ with respect to pair production.

Such extremely compact coronae have not been found in the sample of AGN and XRBs studied by Fabian et al. (2015). Also, the limitations that we outline above and that are also mentioned by Niedźwiecki et al. (2016) challenge the interpretation of LP1 as a distinct component. The component LP1 can therefore only be interpreted in combination with LP2. Although the superposition of the two distinct lamp-post components LP1 and LP2 serves well to describe our data, we hesitate to claim the two point sources to represent a realistic description of the corona. Also, the strong degeneracies between the largely uncertain normalizations of both lamp-post components and the low reflection fraction found for LP1 argue against two extremely compact and distinct coronae. In contrast, a continuous structure enclosing both components seems more likely. The detected spectral signatures of strongly blurred reflection yet indicate in a model-independent way that at least part of the corona must lie very close to the black hole. The low reflection fraction measured for LP1 also points towards a vertically extended structure (Keck et al., 2015). It may be reflected by outflows (see also King et al., 2017) that have been detected for NGC 4151 (Kraemer et al., 2005; Tombesi et al., 2010, 2011), where relativistic aberration (Beloborodov, 1999; Malzac et al., 2001) can effectively reduce the observed fraction of reflected photons. These outflows together with its non-relativistic jets (e.g., Pedlar et al., 1993; Ulvestad et al., 2005) could well fit in this picture and may provide a natural environment for the source of primary X-ray photons (Markoff & Nowak, 2004; Markoff et al., 2005; Wilkins & Gallo, 2015; King et al., 2017).

Reverberation studies have independently suggested horizontally (Wilkins & Fabian, 2012) and, in particular, vertically extended primary sources above the black hole for Ark 564 (Zoghbi et al., 2010), Mrk 335 (Kara et al., 2013), IRAS 13224$-$3809 (Wilkins & Gallo, 2015) and also for NGC 4151 (Zoghbi et al., 2012; Cackett et al., 2014), which indicates a primary source within $\sim$5–$10\,r_{\mathrm{g}}$ from the black hole. This may, in conjunction with the high-quality spectral information of the long-look Suzaku and NuSTAR data, be well in agreement with a jet-based geometry, which can still be relatively radially compact. Magnetic fields (e.g., Baczko et al., 2016, for the extreme case of the radio galaxy NGC 1052) make this region an efficient emitter of synchrotron photons (e.g., Merloni et al., 2000; Markoff et al., 2005), which can act as additional seed photons for Comptonization processes in the corona and therefore at least reduce its transverse extent (Dovčiak & Done, 2016) close to or below the value of $1\,r_{\mathrm{g}}$, which we have previously assumed for our estimates of the compactness.

4.2 Complex absorption variability

We model four layers of absorption: two neutral layers, one partially covering the nucleus with $f_{\mathrm{cvr}}\sim 0.5$ (CA1/$N_{\text{H,1}}$) and one full-covering absorber (CA2/$N_{\text{H,2}}$), a warm absorber with $\log\xi\sim 2.8$ as well as an additional broad absorption line around 8 keV that can either be formed by a Fe xxvi Ly $\beta$ line that is unmodeled by the warm absorber or, more likely, a blueshifted Fe xxv He $\alpha$ or Fe xxvi Ly $\alpha$ line, indicative of an UFO at a speed of $\sim$0.16–0.24 c.

We showed that one or both of the absorbers CA1 or CA2 account for the bulk of the spectral variability between 1–6 keV on timescales from days to years. This agrees with Puccetti et al. (2007) and de Rosa et al. (2007) for a very similar set of absorbers. The partial-covering column $N_{\text{H,1}}$ ranges between $10$ and $25\times 10^{-22}\,\mathrm{cm}^{-2}$. Its covering fraction averages 46% across the Suzaku observations. The full-covering column $N_{\text{H,2}}$, on the other hand, ranges between $5$ and $15\times 10^{-22}\,\mathrm{cm}^{-2}$. The columns $N_{\text{H,1}}$ and $N_{\text{H,2}}$ seem to track one another well; given the data quality and the comparable values of both columns, we can explain this model-dependent result as being due to degeneracies between $N_{\text{H,1}}$ and $N_{\text{H,2}}$. In agreement with earlier analyses using BeppoSAX (Puccetti et al., 2007; de Rosa et al., 2007) and Chandra (Wang et al., 2010), we can therefore not tell if one or the other column dominates the variability. A separation of both neutral absorbers from the inner-disk reflection component is equally challenging, as both describe very similar spectral features at the flat turnover. We could still show that disentangling both is possible with the high count statistics and broad spectral coverage of the long-look Suzaku and NuSTAR observations as well as the variable nature of the absorbers. The latter has also been demonstrated by Risaliti et al. (2009a) and Risaliti et al. (2013) for NGC 1365, which is very similar to NGC 4151 in that regard.

The shortest variability timescale of 2 d has been measured with the time-resolved analysis of Suz 3 (see Fig. 12). The circumnuclear gas in NGC 4151 may be continuous yet non-homogeneous and/or containing (or consisting of) a discrete number of localized clouds. Such scenarios have been explored for this source in the past by Holt et al. (1980), Yaqoob & Warwick (1991), Zdziarski et al. (2002), and Puccetti et al. (2007).

In the case where the variable absorber consists solely of clouds, variability either by $N_{\mathrm{H}}$ or the covering fraction may be associated with single clouds entering or exiting the line-of-sight. We can not exclude variability of the covering fraction due to inherent degeneracies that lead us to fix this parameter. The model by Nenkova et al. (2008b) attempts to explain line-of-sight variability for all inclinations with a Poissonian distribution of clouds that decrease in number density further out. This model has been successfully fitted to infrared (IR) SEDs by Alonso-Herrero et al. (2011) and is able to explain absorption events observed in the X-rays (Markowitz et al., 2014, and, e.g., Beuchert et al. 2015). On average, these studies predict only a few clouds on the line-of-sight for inclinations similar to the the one we measure for NGC 4151. This is consistent with independent estimates for NGC 4151 by Holt et al. (1980), Yaqoob & Warwick (1991) and Zdziarski et al. (2002).

Regardless of whether the variable absorber is CA1 or CA2, we can use the observed variability timescales of $N_{\text{H}}$ to estimate the location of potential clouds that are moving on Keplerian orbits (Risaliti et al., 2002; Puccetti et al., 2007) with

[TABLE]

where $M_{\mathrm{BH}}$ is the black hole mass and where $n_{\mathrm{H}}$ and $N_{\mathrm{H}}$ are the cloud number density and column density, respectively. For $\Delta t$ we use the shortest variability timescale of 2 d. The column density is set to the maximum modeled value of $N_{\mathrm{H,1}}=\sim 25\times 10^{22}\,\mathrm{cm}^{-2}$. The number density is unknown. If we assume number densities for BLR clouds of $n_{\mathrm{H}}=10^{9-10}\,\mathrm{cm}^{-3}$ (Netzer, 1990; Kaspi & Netzer, 1999; Netzer, 2008), we obtain a distance range of $R\sim 5.6\times 10^{-4}$–$5.6\times 10^{-2}\,\mathrm{pc}$ or $3.9\times 10^{2}$–$3.9\times 10^{4}\,r_{\mathrm{g}}$. These distances are consistent with those inferred by Puccetti et al. (2007) and with the distance of the BLR ($\sim 8\times 10^{-3}\,\mathrm{pc}$, Maoz et al., 1991). Note that for the inclination of NGC 4151, the average number of clouds on our line-of-sight stays approximately the same even if we extrapolate the clumpy torus model down to the BLR. In contradiction to this theoretical consideration, Arav et al. (1998) find no signs for distinct BLR clouds in high-resolution optical Keck spectra. This result favors an interpretation of our data with irregular and dynamic absorbing structures. In theory, dust must be entirely or at least partially sublimated at the inferred distances with the dust sublimation radius

[TABLE]

(Nenkova et al., 2008b). We find $R_{\mathrm{d}}\sim 0.13\,\mathrm{pc}=8.9\times 10^{4}\,r_{\mathrm{g}}$ with the assumed bolometric luminosity $L_{\mathrm{bol}}=10^{44}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (Vasudevan & Fabian, 2009) and the dust evaporation temperature of $T_{\mathrm{d}}=1500\,\mathrm{K}$ (Barvainis, 1987). When the inner range of the torus is estimated independently, slightly smaller distances are found; using the 5100Å line luminosity (Kaspi et al., 2005), the outer BLR of NGC 4151 could be constrained to a distance of $\sim 6\times 10^{3}$–$6\times 10^{4}\,r_{\mathrm{g}}$, and the inner rim of the torus to $\sim$0.04 pc using thermal dust reverberation studies by Minezaki et al. (2004) and Burtscher et al. (2009). Schnülle et al. (2015), however, find no signs for dust sublimation in their data and explain this with large graphite dust grains that sublimate at much higher temperatures. If we assumed clouds at a distance of $R=R_{\mathrm{d}}$, the cloud density required by the observed $N_{\mathrm{H}}$ variability pattern would be on the order of $\sim 6\times 10^{10}\,\mathrm{cm}^{-3}$, which can be excluded for a dusty torus (Elitzur, 2007).

On the other hand, we also detect variability over longer timescales of approximately one year that may indicate clouds at larger distances. If we assume, for example, number densities of $10^{7-8}\,\mathrm{cm}^{-3}$, which are typical for the dusty torus (Miniutti et al., 2014; Markowitz et al., 2014), we find $R\sim 2\times 10^{-3}$–$2\times 10^{-1}\,\mathrm{pc}$ or $1.3\times 10^{3}$–$1.3\times 10^{5}\,r_{\mathrm{g}}$. These values put the clouds into the outer BLR or at the inner side of the dusty torus, which is also a common result of Markowitz et al. (2014) for similar timescales. Much larger distances of clouds from well inside the torus would require unrealistically large number densities and are therefore unlikely.

Recently, Couto et al. (2016) published a study that investigates X-ray absorbers with archival long-look Chandra observations. They find a highly ionized column similar to our xstar component as well as an outflowing near-neutral absorber with a speed of $\sim 500\,\mathrm{km}\,\mathrm{s}^{-1}$. This neutral absorber goes back to absorption line features in HST/GHRS/STIS data first mentioned by Weymann et al. (1997) and referred to as the kinematic component “D+Ea” in Table 1 of Kraemer et al. (2001); see also Kraemer et al. (2006) for further usage. Couto et al. fit these absorbers to a number of seven archival Chandra observations with two additional observations by XMM-Newton in 2000 and Suzaku in 2006. They conclude that the bulk of the spectral variability over 14 years is caused by a change in the ionization state of both absorbers as a response to changes in the irradiating luminosity rather than the observed variations in $N_{\text{H}}$. Similar to the absorbers used by Couto et al. (2016), we find the columns of both neutral absorbers, CA1 and CA2, to be anticorrelated with the incident photon flux of the HXCOMP component. A portion of either of these is likely consistent with the component D+Ea. Its outward motion forms a consistent picture with the observed anticorrelation of the incident flux with the column density, in which a still dusty, radiatively driven wind (e.g., Czerny & Hryniewicz, 2011; Dorodnitsyn & Kallman, 2012) causes a decrease of the line-of-sight absorption for stronger radiative driving. We, however, emphasize that the interpretation of the data analyzed in this work strongly depends on the method of modeling. Both the interpretation with orbiting clouds or with an outflowing near-neutral absorber represent structural changes in the absorber. We do not favor one over the other. In a yet different scenario, the ionization state of the absorber can change with varying irradiation, resulting in changes of the equivalent column density. We, however, are not sensitive to this effect with the available count statistics.

On larger scales, Ruiz et al. (2003) and Radomski et al. (2003) detect dusty extended gas in NIR/MIR data, which is fully covering the nucleus similar to dust that has been found within $\sim$4 pc using the Gemini NIR integral field spectrograph Riffel et al. (2009). This dusty component may therefore make up a non-variable portion of our CA.

4.3 The soft X-rays

We model the soft X-rays of all observations with a SXCOMP with the photon index tied to that of the HXCOMP as well as a blend of Gaussian lines that are motivated from high-resolution grating observations of XMM-Newton/RGS and Chandra/HETG. We have shown that the soft flux below 1 keV is, contrary to the hard flux, only mildly variable at a low dynamic range of $\sim$6% at most. The origin of the variability (continuum, emission line(s), or both) is unclear. If we calculate the ratio of the SXCOMP and HXCOMP normalization, we can, on average, infer a low optical depth of the soft emitting gas of $\tau\sim 0.023$.

Due to the comparatively large PSFs of XMM-Newton and Suzaku, we are unable to spatially resolve the line-emitting gas. In contrast, a number of authors have been using Chandra for this purpose, which is both powerful in spatial and spectral resolution. Ogle et al. (2000) and Wang et al. (2011a) show that a considerable part of the soft X-rays is due to distant, extended gas that is spatially coincident with a bi-conical gas distribution (see Storchi-Bergmann et al., 2010, for integral field spectroscopy of [O iii]) and the narrow line region (NLR, Bianchi et al. 2006). A spatially resolved modeling of the extended gas with Cloudy allows Wang et al. (2011a) to conclude a two-phase photoionized medium of intermediate ($\log\xi\sim 1.7$) and high ($\log\xi\sim 2.7$) ionization (see also Armentrout et al., 2007, for similar modeling using Cloudy with XMM-Newton/RGS data) next to the collisionally ionized phase. The higher ionized line-emitting gas phase detected using Cloudy may also be consistent with the highly ionized warm absorber of similar ionization ($\log\xi\sim 2.8$) that we model with XSTAR.

The evidence for a thermal, collisionally ionized gas phase may point towards a contribution of a bremsstrahlung continuum to the soft X-rays, which we, for simplicity, model with a SXCOMP continuum, that is, scattered nuclear Comptonized emission. We are unable to favor either of both options with our data, but provide the reader with a short discussion on the implications of our chosen model. In this picture, the observed that highly ionized phase, which we refer to as a “warm mirror” (Guainazzi et al., 2005; Guainazzi & Bianchi, 2007a; Guainazzi & Bianchi, 2007b) can act as scattering medium for the nuclear emission. This scenario would justify our SXCOMP component to be a long-term average with respect to the variable nuclear HXCOMP (see also Pounds et al., 1986; Yang et al., 2001; Wang et al., 2011b). The low degree of flux-variability that we measure between 0.6 and 1.0 keV may be explained with portions of this warm mirror that are located close enough to the nuclear source to respond to its variability. In fact, around 30%666This number can only be approximate due to the recent improvements on the Chandra PSF. of the soft X-rays originate in a region that is unresolved by Chandra  (Ogle et al., 2000). The response of this mirror at various distances from the source could explain the lack of correlated variability between the soft and hard bands. A correlation of the HXCOMP variations with the prominent O vii and Ne ix line fluxes can neither be claimed nor excluded with respect to the large uncertainties at CCD resolution. Also, we can neither report to be in favor of or against correlated variability of these lines with the column density of the absorbers. This would be expected if the clumpy absorber were to temporarily block the nuclear irradiation onto the diffuse gas and promote recombination of the gas.

Other than the warm mirror, an intrinsic soft excess may also explain the soft emission, which has phenomenologically been modeled with a steep soft power-law (e.g., Yang et al., 2001; Wang et al., 2010). This also involves a fit with an extra bremsstrahlung component (e.g., Warwick et al., 1995). In a more self-consistent picture, the soft excess emission could be provided by a combination of blurred and unblurred ionized disk reflection, that is, by the model components relxillCp_lp and xillver, which we use to describe the relativistic features inherent in the continuum. Although this works relatively well for the soft continuum at CCD resolution, xillver is unable to model the highly ionized species of the H-like and He-like ions O vii/O viii and Ne ix, which we observe with XMM-Newton/RGS. In addition, xillver predicts a number of strong lines from Mg, Si, and S rather than O or Ne, which are not observed with XMM-Newton/RGS. We therefore prefer an independent soft continuum due to extended gas on larger scales.

5 Conclusions and outlook

The unique Seyfert galaxy NGC 4151 allows to both probe the circumnuclear absorber and its strongly variable absorption together with the effects of strong gravity close to the BH via imprints on the reflection spectrum. As part of this work, we conducted a follow-up study based on Keck et al. (2015) and apply the improved model relxillCp_lp that describes blurred reflection in a physically motivated and self-consistent lamp-post geometry together with the complex set of neutral absorbers that has been frequently reported in the literature. We apply the resulting baseline model to all Suzaku, XMM-Newton, and NuSTAR spectra that we consider in this work and perform a time-resolved spectral analysis of the neutral absorbers CA1 and CA2 that are variable on timescales from days to years. We find the soft X-rays below 1 keV to be only mildly variable within a maximum of $\sim$6% as opposed to $\sim$20% for the HXCOMP. Strong spectral variability is apparent between 1 and 6 keV. As a result of our dedicated modeling, we come up with the following conclusions:

We identify two separate point-like lamp-post components LP1 and LP2 in simultaneous long-look Suzaku and NuSTAR spectra at heights of $h\sim 1.2\,r_{\mathrm{g}}$ and $h\sim 15.0\,r_{\mathrm{g}}$, respectively. We applied the most recent model relxillCp_lp that combines the reflection spectrum off each point on an ionized disk (xillver) with the appropriate relativistic transfer function. This model uses the Comptonization continuum nthcomp as primary continuum. The normalizations and reflection fractions of both lamp-post components are highly degenerate. In particular we measured a low reflection fraction for LP1 and find that runaway pair production would dominate for a single and compact corona close to the BH. We therefore propose a vertically extended corona as opposed to two distinct and compact coronae. We emphasize that our results possibly reflect Comptonization processes in a jet-base and emphasize the presence of non-relativistic jets in NGC 4151. An outflowing corona would flatten the emissivity profile, which may explain the low observed reflection fraction. We must additionally consider the jet-base and its magnetic field as a source of synchrotron photons. These additional seed photons may cause the corona to be relatively compact in the horizontal direction. The relevance of magnetic fields for coronae close to the BH is strongly implied and needs to be carefully investigated in the future. 2. 2.

Thanks to the high count statistics provided by the long-look Suzaku/NuSTAR campaign, we are able to constrain a complex system of four separate layers of absorption, that is, two neutral absorbers, one of which is partially covering the nucleus with only 40%–50% (CA1) and one fully covering the same (CA2), a third layer of highly ionized absorption ($\log\xi\sim 2.8$) as constrained with absorption features of Fe xxv and Fe xxvi as well as an UFO with an outflow velocity of $\sim$0.16–$0.24\,c$ from a broad absorption feature around $\sim$8 keV. We also showed that we are able to distinguish the relatively flat turnover of the neutral partial coverer between 3 and 6 keV from the broad and blurred reflection components, which both describe a similar spectral shape in this energy range. 3. 3.

We observe both columns CA1 and CA2 to be strongly variable both on short timescales of 2 d (probed with the long-look Suzaku observation) and long timescales of approximately one year (probed with two additional Suzaku observations as well as a XMM-Newton monitoring). Both absorbers are responsible for the bulk of the spectral variability observed between 1 and 6 keV. Their variability patterns are similar at all timescales, which is likely caused by the observed degeneracies between both columns. We are therefore unable to tell from our modeling, if one or even both absorbers are intrinsically variable. The observed column density evolution could be interpreted as a clumpy absorber, where one or more single clouds transit the line-of-sight at a distance as close as the BLR up to the inner side of the torus. While no BLR clouds have yet been detected in the BLR of NGC 4151, a clumpy nature of the inner dusty torus may still be a valid explanation. On the other hand, the anticorrelation of the irradiating photon flux with the column densities $N_{\text{H,1}}$ and $N_{\text{H,2}}$ offers the alternate explanation of a radiatively driven dusty wind or changes in the ionization degree of the near-neutral absorber. 4. 4.

The soft X-rays below 1 keV are only mildly variable. They likely originate in extended gas, which is included by the large PSF of Suzaku and XMM-Newton. We observe a blend of unresolved emission lines in the CCD spectra together with a weak continuum that we model with the SXCOMP. We analyzed the combined XMM-Newton/RGS spectrum of all seven XMM-Newton pointings and the underlying soft emission with a phenomenological blend of Gaussian line profiles. Line-ratio diagnostics on the dominant He-like triplets O vii and Ne ix suggest that the gas must primarily be photoionized with a minor contribution of a collisionally ionized phase. Related to the choice of the SXCOMP as underlying continuum, we discuss the extended and highly ionized gas in terms of a “warm mirror”, scattering nuclear Comptonized continuum emission into our line of sight. We outline observational evidence for this gas phase, that is, the $\log\xi\sim 2.8$ warm absorber and the soft line emission. The mild degree of variability in the soft X-rays may originate in gas that is located close enough to the nucleus to be able to respond to changes in the hard X-ray continuum within the probed timescale. We can exclude intrinsic soft excess emission due to blurred, ionized reflection.

Acknowledgements.

We thank the anonymous referee for providing us with a number of detailed comments that greatly improved the clarity of this manuscript. We made use of ISIS functions provided by ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/) as well as the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We thank J. E. Davis for the development of the slxfig module that has been used to prepare the figures in this work. This work used data obtained with the Suzaku satellite, a collaborative mission between the space agencies of Japan (JAXA) and the USA (NASA) as well as XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This work was supported under NASA contract No. NNG08FD60C, and made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). A.G.M. acknowledges support from NASA grant NNX15AE64G. A.A.Z. has been supported in part by the Polish National Science Centre grants 2013/10/M/ST9/00729 and 2015/18/A/ST9/00746.

Appendix A Continuum and line parameters