Test for Echo: X-Ray Reflection Variability in the Seyfert-2 Active Galactic Nucleus NGC 4388

We report on a study of the narrow Fe Kα line and reflection spectrum in the well-known Seyfert-2 active galactic nucleus (AGN), NGC 4388. X-ray spectra summed from two extensive NICER monitoring campaigns, separated by years, show strong evidence of variation in the direct continuum and reflected emission, but only small variations in the obscuring gas. Fits to the spectra from individual NICER observations find a strong, positive correlation between the power-law photon index, Γ, and direct flux that is commonly observed in unobscured AGN. A search for a reverberation lag between the direct and reflected spectra—dominated by the narrow Fe Kα emission line—measures a timescale of t=16.37−0.38+0.46 days, or a characteristic radius of r=1.374−0.032+0.039×10−2 pc =3.4−0.1+0.1×104GM/c2 . Only one cycle of this tentative lag is observed, but it is driven by a particularly sharp drop in the direct continuum that leads to the subsequent disappearance of the otherwise prominent Fe Kα line. Physically motivated fits to high-resolution Chandra spectra of NGC 4388 measure a line production radius of r=2.9−0.7+1.2×104GM/c2 , formally consistent with the tentative lag. The line profile also prefers a Compton-thick reflector, indicating an origin in the disk and/or thick clumps within a wind. We discuss the strengths and weaknesses of our analysis and methods for testing our results in future observations, and we note the potential for X-ray reverberation lags to constrain black hole masses in obscured Seyferts wherein the optical broad line region is not visible.

. Only one cycle of this tentative lag is observed, but it is driven by a particularly sharp drop in the direct continuum that leads to the subsequent disappearance of the otherwise prominent Fe Kα line.Physically motivated fits to high-resolution Chandra spectra of NGC 4388 measure a line production radius of r GM c 2.9 10 0.7 1.2 4 2

Introduction
NGC 4388 is a heavily obscured Seyfert-2 galaxy.At a distance of d = 18.0 Mpc (Sorce et al. 2014), it is one of the closest examples of this type of active galactic nucleus (AGN) and one of the brightest in X-rays.Nonthermal X-ray emission from the nucleus has been detected out to nearly 500 keV (Beckmann et al. 2004).The obscuration that hides the central engine is not Compton thick, but is only a factor of a few below this key threshold (see, e.g., Miller et al. 2019).Although the mass of the black hole cannot be determined using optical reverberation techniques (see, e.g., Peterson et al. 2004), H 2 O megamaser emission in its disk has enabled a precise mass measurement: M BH = 8.4 ± 0.2 × 10 6 M e (Kuo et al. 2010).
This mass measurement and the high X-ray flux of NGC 4388 present some unique opportunities.The flux variability timescales in the direct and reprocessed emission from NGC 4388 can be translated into gravitational radii (r g = GM/c 2 ) with unusual precision and can therefore help to resolve the nature and origin of the obscuring torus geometry (e.g., Antonucci 1993).At least one past study has found that the obscuration in NGC 4388 can vary by orders of magnitude on a timescale of just hours (Elvis et al. 2004), suggesting that it extends much closer to the black hole than the parsec-scale geometry that is classically envisioned (also see Risaliti et al. 2002).
The narrow Fe Kα emission line is the strongest atomic feature in the X-ray spectra of Seyfert AGN (for a recent review, see Gallo et al. 2023).Owing to its large geometric covering factor of the torus, at least some Fe Kα line flux must originate on its face, regardless of its distance from the central engine.However, recent studies of some Seyfert-1 AGN suggest that the bulk of the line flux is instead produced within the optical broad-line region (BLR) or at even smaller radii (see, e.g., Miller et al. 2018;Zoghbi et al. 2019;Noda et al. 2023).These studies are particularly sensitive to the innermost extent of the torus.The high flux observed in NGC 4388 has even permitted a tentative detection of Ni Kα emission (Fukazawa et al. 2016).However, sensitive observations of NGC 4388 have not revealed relativistic reflection from the innermost disk (Kamraj et al. 2017;Yaqoob et al. 2023).
It is possible that the line of sight determines which geometry contributes the bulk of the Fe K line flux, and NGC 4388 may provide a key angle on this question.The observed line flux may be dominated by a distant part of the torus along our line of sight, since the BLR is obscured.
Although the innermost extent of the torus may only be a few times larger than the optical BLR (e.g., Minezaki et al. 2019), its radial extent may span orders of magnitude, out to the scale of parsecs (or larger).Complete quenching of the line flux from such a geometry should be nearly impossible in this circumstance, and variations should be slow and muted relative to the central engine and the response of the BLR.If the innermost edge of the obscuring gas instead extends as close to the black hole as prior studies may indicate, and if the innermost edge contributes the bulk of the narrow Fe K line flux, it may also vary on short timescales.We note that weak Fe K emission on kiloparsec scales has been detected in NGC 4388 (Iwasawa et al. 2003;Yi et al. 2021), but this emission is not expected to vary on human timescales.
The spectral resolution and collecting area of NICER (Gendreau et al. 2016) enable spectroscopic monitoring of bright Seyferts, opening new windows on the accretion flow geometry in these AGN.Motivated by this opportunity and the strong variations previously found in NGC 4388, we obtained a series of short monitoring exposures to search for variations in the obscuring column and Fe Kα line flux as diagnostics of the torus.We also summed these new observations to create a very sensitive time-averaged spectrum, similar to another series of prior observations (Miller et al. 2019).
In this paper, we use both the prior and new sets of NICER observations to leverage the variability of spectral features to better understand the inner accretion flow in NGC 4388.We additionally use archival Chandra observations to compare these results with those derived from the detailed Fe K line profile in NGC 4388.In Section 2, we detail the sets of observations and our reduction of the data.Our analysis and results are presented in Section 3. Finally, we discuss the strengths and limitations of our work in Section 4, and we speculate on the potential of future observations.

NICER
NICER observations of NGC 4388 can be grouped into three sets.The first set was obtained near the start of the NICER mission, beginning on 2017 December 3 (MJD 58090) and continuing almost daily until 2018 January 11 (MJD 58129).A series of 34 observations accumulated 109.8 ks of total exposure time.The second set consists of 10 observations from 2018 December 5 (MJD 58457) to 2019 April 4 (MJD 58577) for a total of 9.2 ks.An analysis of the time-averaged spectrum of these two campaigns is reported in Miller et al. (2019).The third set includes 11 daily observations from 2022 June 9 (MJD 59739) to 2022 June 19 (MJD 59739), for a total of 100.0 ks.
In this investigation, we were concerned with separately exploring variability on short and long timescales, on the order of days and years, respectively.The earliest and latest of these three epochs suited these goals, since both have roughly 100 ks of total exposure time and they are separated from each other by over four years.Because the observations in the middle epoch are spaced relatively far apart and have comparatively low exposure times, we chose not to include them in our analysis.While there is some variation, the earliest epoch had a typical exposure time of about 1.5-4 ks per observation.Compared to the latest epoch, with the majority of observations having exposure times greater than 9 ks, the early epoch covers a longer time period and is therefore additionally well suited to investigating variability on a timescale of days.
Table 1 presents the NICER observations used in this analysis.Each of the 45 observations from the first and third epochs (ObsID 1117010101-34 and 5117010101-11) was processed using HEASoft 6.29c.The data reduction was accomplished using NICERDAS 8c, with the xti20210707 CALDB calibration file.The more recent observations all had a considerably larger high-energy background that was not flagged by the standard nicerl2 pipeline.Using the program nicermaketime, we enforced a cap on the magnetic cutoff rigidity (COR_SAX), which is an index describing the level of particle background, as well as the number of overshoot events that are caused by high-energy photons, charged particles, and cosmic rays.In particular, we imposed the nicermaketime parameters cor_range=2.0-*(no default value) and overonly_range=0-1 (default 0-1.5) on all observations.We additionally excised two isolated time intervals (∼100 s each) that were responsible for low-energy background flares in ObsIDs 1117010133 and 1117010134 that were not caught by the previous methods.
After this screening, the net exposure time was 103.3 ks for the first epoch and 96.0 ks for the third epoch.The individual observations from each period were also combined using nicermergeclean, and then all 47 spectra (45 individual observations plus two time averages) were extracted using XSELECT v2.4m.After spectra had been extracted, ARF and RMF response files were generated using nicerarf and nicerrmf, and the spectra were then binned with the optimal method (Kaastra & Bleeker 2016), using the opt parameter in ftgrouppha v1.6.
The background that NICER experiences depends on several factors, including the location of the International Space Station in its orbit and the space weather conditions at the time of observation.The program nibackgen3C50 (Remillard et al. 2022) generates a background estimate by dividing the good time intervals from a given observation into smaller segments and then matching them to corresponding intervals from a library of reference spectra.The program uses different libraries depending on the gain epoch initially used to extract the spectrum, so for each observation we used gaine-poch=2020, which corresponds to the latest CALDB gain file at the time of our analysis.Note.Details of the observations used in this analysis.Our main spectral analysis uses two epochs of NICER observations, separated by roughly 4.5 yr, with each amounting to about 100 ks of total exposure time.We additionally used two Chandra observations, totalling 269 ks, for the supplementary analysis.
Finally, the fully calibrated NICER passband is generally taken to be 0.2-12.0keV, but we limited our energy range to 0.6-9.0keV for the present analysis.The low-energy band was excluded because the detector flux calibration is less certain in this band, and we excluded the high-energy band because most of the individual spectra are dominated by background flux above 10 keV.An energy cutoff at 9 keV retains as many data as possible without background contamination adversely affecting the results of our analysis.

Chandra
Chandra observed NGC 4388 on two occasions, starting on 2008 April 16, 02:03:04 (170.6 ks total exposure) and 2008 April 24, 14:02:59 (98.3 ks total exposure), with the High Energy Transmission Gratings (HETG) in place.We downloaded calibrated spectra and responses using the tgcat facility (http://tgcat.mit.edu).The exposures were obtained sequentially and are part of the same observation, so we coadded the first-order spectra from these observations using the CIAO tool combine_grating_spectra.After the processing applied in the tgcat sequencing, the net exposure time is 269 ks.
Our focus is on the detailed shape of the strong, neutral Fe Kα emission line, so we restrict our analysis to the summed first-order high-energy grating (HEG) spectrum; it has the best spectral resolution and sensitivity in this band.Prior to analysis, the data were binned using the HEASOFT tool ftgrouppha, again using the optimal binning algorithm.Within XSPEC, the data were fit within the 5.0-7.5 keV band to ensure that our results are driven by aspects of the Fe Kα line.A Cash statistic (Cash 1979) was minimized to derive the best overall models and 1σ confidence limits.

Spectral Modeling
All fitting to NICER spectra was performed in XSPEC v12.12.0 (Arnaud 1996) and minimized a χ 2 fit statistic using  Note.Key model parameters from fits to time-averaged NICER spectra for NGC 4388.All reported errors are 1σ uncertainties in the model parameters, and the numbers in parentheses represent symmetric errors in the last digit quoted.The lowend continuum is dominated by three apec plasma components, each with a characteristic temperature kT and normalization.The main power-law continuum cutoffpl has a photon index Γ and normalization, and it is subject to two different absorbers.The neutral TBpcf has a column density N H and partial covering fraction pcf, and the ionized ionabsorb has an ionization parameter ξ and separate column density.The reflection component pexmon shares its photon index and normalization with cutoffpl, so its strength is controlled by the reflection fraction rel_refl.The reflection fraction is presented as a negative number, as it appears when fitting, because pexmon uses a minus sign to indicate that only the reflection has been modeled and not the power law and reflection together.We have chosen to model the continuum and reflection as separate components, hence the minus sign.Also presented here are the total received flux f tot , unobscured power-law flux f PL,0.6-9.0 , and reflection component flux f refl,0.6-9.0 , along with adjusted values in the more standard range 2.0-10.0keV, denoted f PL,2-10 and f refl,2-10 .The 2-10 keV X-ray luminosity L X,2-10 is calculated using f PL,2-10 and a distance of 18.0 Mpc (Sorce et al. 2014), and the Eddington fraction L X,2-10 /L Edd. is presented for reference (see Section 3.3).
standard model weighting.The model used for our primary investigation is as follows: TBabs * (apec 1 + apec 2 + apec 3 + TBpcf * ionabsorb * cutoffpl + pexmon) TBabs (Wilms et al. 2000) models neutral absorption within the Milky Way, which in all fits is fixed at the value 2.57 × 10 20 cm −2 resulting from using the FTOOL nH (Kalberla et al. 2005).APEC models emission from collisionally ionized diffuse gas (Smith et al. 2001), and three components of different temperatures were required to produce a good fit of the low-energy continuum at 3.0 keV.The main continuum from the central engine is modeled by a cutoff power law cutoffpl, which is modified both by an ionized absorber ionabsorb and a partial covering neutral absorber TPpcf.The ionized absorber is an XSTAR table model (see below); its function in these spectra is to model an absorption feature slightly above the Fe Kα line.Initial fits poorly constrained the cutoff energy for cutoffpl, so it was frozen at 100 keV in all subsequent fits.Finally, a neutral reflection spectrum is modeled using pexmon (Nandra et al. 2007), which includes self-consistent Fe and Ni Kα and Kβ lines.The most prominent reflection feature in these spectra, necessitating the inclusion of this component, is the narrow Fe Kα line at 6.4 keV.The parameters in pexmon for power-law photon index, normalization, and cutoff energy are tied to those of cutoffpl.For cutoffpl, pexmon, and all apec components, the model parameter for redshift is frozen at the global redshift of the galaxy (z = 8.4 × 10 −3 , Lu et al. 1993).
The XSTAR table model was constructed using the xstar2xspec functionality within the suite.We assumed solar abundances for all elements based on Grevesse et al. (1996), and a turbulent velocity of v turb = 300 km s −1 .The input spectral energy distribution consisted of a hot, T = 25,000 K blackbody and Γ = 1.7 power-law continuum (artificially bent to zero at low energy) in a 10:1 flux ratio typical of bright, unobscured AGN.The ionizing luminosity was assumed to be L = 1 × 10 44 erg s −1 over the 0.0136-13.6keV band.The grid points were selected to ensure at least 10 points per decade in the critical parameters: 80 points between 2 log 6 x -< < , and 40 points between 1 × 10 20 cm −2 < N H < 6 × 10 23 cm −2 .
After fits were made with this model, we calculated the power-law flux by adding a cflux component to cutoffpl, such that the power-law component of the model reads TBpcf * ionabsorb * cflux * cutoffpl.We also calculated the reflected flux in a similar manner, using cflux * pexmon in separate fits.In both cases, the energy bounds of cflux were set to the full energy range of the present analysis, 0.6-9.0keV, and all model parameters were then frozen, except for the lg10Flux parameter of cflux.To facilitate a better comparison with other measurements in the literature, we also repeated this process with energy bounds at 2-10 keV to calculate fluxes in this more standard energy range.
All uncertainties presented in this analysis reflect 1σ confidence intervals calculated with the error command in XSPEC, unless otherwise specified.These values were corroborated with the results from running steppar as necessary when nonmonotonicity in the fit statistic was detected.

Long-term Variability
Figure 1 and Table 2 summarize the results of the best fits to the time-averaged spectra from each epoch of observations.Comparing the model parameters of the two time-averaged spectra allows us to determine how the overall behavior of NGC 4388 changes on a timescale of approximately 1650 days (4.5 yr). Figure 1 clearly reveals a significant difference in the flux detected above approximately 3 keV compared to the relatively constant flux measured below 3 keV.While our modeling of the emission below 3 keV with apec does result in statistically significant variations in temperature and normalization, there are also a number of factors complicating our confidence in this result.Most noticeable among them are the large residuals between ∼1 and 2 keV, which lie in the range of known NICER detector features.As seen in the left panel of Figure 1, the profiles of corresponding apec components between the two models are generally similar but are slightly shifted in energy, and there is also a significant change between the two power-law components of the larger model; the cutoffpl component contributes more flux in this energy band in the later than in the earlier epoch.Along with the fact that the spectra become locally background dominated around 2 keV, it is possible that detector features, small changes in the background estimation, or differences in the cutoffpl model component contribute to different best-fit values for the various apec components.
The column density of the neutral absorber TBpcf is consistent between the two epochs, with both fits yielding a similar value of around 4 × 10 23 cm −2 .This value is higher than recent attempts at modeling these data (Miller et al. 2019), but it broadly agrees with previous measurements of the column density of this source (Elvis et al. 2004).The associated covering factors for each of these column densities are also very similar, so that together, both pairs of parameters are consistent with a neutral absorber whose average behavior does not change drastically over a period of several years, either in column density or in its covering fraction toward our line of sight.
Despite this apparent constancy in the absorbing column, we still observe changes in the continuum emission on these long timescales.This variability in power-law flux-a change of a factor of two between epochs-therefore likely indicates true changes in the behavior of the central engine, rather than just changes in obscuration in our line of sight.As the flux of this part of the continuum decreases from the first to the second epoch, so too does its spectral slope; the photon index drops significantly from Γ = 1.52 ± 0.02 to 1.32 0.04 0.03 G = -+ .

Short-term Variability
Motivated by the spectral variability seen between these two observational epochs, we next focused on the variability among the 45 individual observations that compose them.The most drastic and most evident change between the two timeaveraged spectra in Section 3.2 is in unobscured power-law flux.Figure 2 additionally plots the source flux in the energy range dominated by power-law emission, showing differences within and between the two epochs on finer timescales.Since power-law flux is an indicator of the behavior of the central engine, the magnitude of this observed variability therefore enables us to investigate how other parameters change with or respond to the central engine.
For sufficiently luminous AGN, we typically observe a softer-when-brighter phenomenon, in which an increasing proportion of the power-law continuum comes from soft X-rays as the total luminosity increases (Yang et al. 2014).Physically, this can be explained by a more luminous central engine-emitting more photons per unit time-that cools the corona more effectively.Compton upscattering by electrons in the corona is theorized to be most responsible for X-ray continuum emission, so a cooler corona would result in less emission in harder X-rays and a softer continuum overall (Fanali et al. 2013).
This phenomenon can be observed by investigating the relation between the power-law flux and the spectral slope, Γ, which is shown in Figure 3. From fitting all individual observations to a line using a Monte Carlo fitting procedure, we find a best fit of Γ = (0.27 ± 0.06) × f PL,2-10 + (1.11 ± 0.09) using the same units as in the figure.Since a higher value of Γ corresponds to a steeper spectral slope and lower flux at higher energies, this positive correlation does indeed indicate softerwhen-brighter behavior in NGC 4388.This is expected for AGN with L X /L Edd.> 10 −3 (Yang et al. 2014;Connolly et al. 2016), and NGC 4388 agrees with this trend in both epochs (see Table 2).
We note that some individual observations reached a hard lower limit of Γ = 1.1 when fitting, which is the lowest value allowed by pexmon.This limit is the point below which calculation of the Fe and Ni Kα and Kβ lines breaks down, indicating more unphysical behavior.The majority of these offending observations have short exposure times and/ or excess background contamination, so these values of Γ and their uncertainties are less reliable and could affect the numerical results of the best fits.However, fitting without these observations still yields the same qualitative result of a positive correlation between the photon index and the power-law flux; a simple Spearman's rank correlation test returns a correlation coefficient of r s = 0.635 (p = 3 × 10 −6 ) for all observations and r s = 0.596 (p = 6 × 10 −5 ) without the Γ = 1.1 observations.We note that joint fits to INTEGRAL and Swift observations of NGC 4388 measured a power-law index of Γ = 1.2 in some flux and hardness windows (Fedorova et al. 2011); it is possible that we have sampled similar phases.

Lag Analysis
After investigating the relation between the flux of the power-law component and its spectral index, we then analyzed the relation between the power law and the reflected continuum.Figure 4 shows the reflected flux from pexmon versus the power-law flux from cutoffpl; there is evidently a wide spread in the data.This relatively weak correlation between the reflected flux and the direct flux could be the result of failing to account for light-travel times between the compact central corona close to the black hole and the primary reflector.Lags of this type are commonly measured between the central engine and BLR in optical studies of AGN (Peterson 1993;Cackett et al. 2021), and were recently measured in X-ray studies as well (Zoghbi et al. 2019;Noda et al. 2023).To test  this possibility, we next undertook a lag analysis using a package that is commonly employed to characterize lags in optical studies of Seyferts and quasars.
JAVELIN is a Python software package that computes lags between continuum and emission-line spectra (Zu et al. 2011).It is fully general and can be applied to X-ray data (see, e.g., Zoghbi et al. 2019).JAVELIN can estimate lags even within relatively sparse data by introducing two key assumptions about the nature of the light curves themselves.First, JAVELIN assumes that a continuum light curve can be modeled as a damped random walk.Second, JAVELIN further assumes that any emission-line light curve can be modeled as a shifted, scaled, and smoothed version of the continuum light curve.
JAVELIN models this damped random walk with two parameters, one for a characteristic timescale (τ), and the other for a characteristic flux amplitude variability (σ).This model is then convolved with a top-hat function to fit the line emission light curve, using three additional parameters: (1) a horizontal shift, describing a lag time between the continuum and line emission; (2) an amplitude, or a scaling factor between the magnitude of continuum and line flux; and (3) a top-hat width, effectively blurring the continuum signal within a given range.The value of the lag time parameter is of primary interest to the present investigation, since a measure of the lag t between a continuum signal and emission line yields a characteristic distance scale d = ct for the line-emitting region.(In principle, the top-hat width parameter could also be used to constrain the extent of such a region, giving an inner and outer characteristic radius.) To construct the light curves necessary for use as inputs to JAVELIN, we used the flux calculated from the cutoffpl and pexmon components, respectively, from the best-fit model of each individual observation in the first epoch of data collection.The pexmon component is dominated by the 6.4 keV iron line, but also encompasses a more complex reflection spectrum.We only use data from the first epoch of observations because initial attempts using both epochs resulted in computed lag times on the order of the time between them.We therefore limited our scope to the available epoch with the longer period of observations.For each of these remaining 34 observations-spanning a period of 39 days-we calculated fluxes for the continuum and reflected components using the method described in Section 3.1.The associated time for each point on the light curve is the median of its observation time interval in MJD.
Figure 5 shows the resulting parameter distributions of a Markov Chain Monte Carlo (MCMC) analysis in JAVELIN.The best-fit models resulting from these distributions are shown in Figure 6, superimposed on the light curves used for the analysis.In JAVELIN, the argument laglimit was set to the range 0-30 days to prevent negative lag times and to reduce the prevalence of calculated lags with a length comparable to the total window of data collection.From these results, we find a lag time between continuum and reflection of 16.37  .This derived lag time can be seen quite clearly by looking at the prominent dips in the modeled light curves of Figure 6; there is a strong dip in the continuum close to MJD 58100 and a similar response by the reflection component some time later.The two individual observations closest to these two dips-ObsIDs 1007010108 and 10007010121-are shown in Figure 7.The time between the medians of both observations is 16.20 days, very similar to the lag found by JAVELIN.In the first observation, the power-law continuum drops to significantly lower values than average, but the Fe K line-dominated reflected continuum does not.The parameter rel_refl in pexmon, which is a scaling factor describing the relative strength of the reflected component compared to the direct continuum, has a magnitude of 0.99 0.28 0.37 -+ (note that the epoch average is 0.141 0.008 0.008 -+ ).In the second observation, the powerlaw continuum returns close to the epoch average, but the Fe K line is very faint and almost appears to be absent; the magnitude of rel_refl can only be constrained with an upper bound of 0.028.Even without relying on spectral modeling, the third panel of Figure 7 shows clear variability in the energy range corresponding to the power-law continuum, as well as significant but distinct variability in the Fe K line.Between these two observations, then, we see strong variability in the direct continuum followed by strong variability in the reflected continuum on a timescale consistent with that found by the JAVELIN analysis.
These two observations happen to have particularly strong changes in power law or reflected flux, and they seem to significantly affect the JAVELIN determination of a tentative lag timescale.Indeed, fitting the light curves from Figure 6 without these two observations yields a much broader distribution of possible lags than those found in Figure 5.Because of their strong effect on our variability analysis, we wished to ascertain the significance of these two dip points within their respective light curves; our flux calculations are inherently model dependent, so we wished to have some way of quantifying a level of confidence in the reliability of this lag determination and results.
In addition to fitting with JAVELIN, we also fit these light curves with a naïve linear model to see by how much the flux values at the dip points vary from the average behavior of the rest of the light curve.In this sense, the power-law flux at ObsID 1117010108 is lower than the model by 7.5σ, and the reflected flux at ObsID 1117010121 is lower by 4.0σ.We then included a Gaussian model component at the locations of the dips to find out by how much the normalization of the Gaussian excludes zero.Using this metric, the two dips are significant at the level of 12σ and 3.6σ.In either case, these tests show that these two instances of variability are indeed significant; the fact that our estimation of a lag depends strongly on these points still means that it is strongly grounded in the data.
One additional check on the reliability of the results found by JAVELIN is understanding how adjusting for a lag of t 16.37 0.38 0.46 = -+ days improves the (lack of) correlation seen in Figure 4. Figure 8 presents fits to the reflected flux versus power-law flux data for the first epoch of observations alongside points that have been adjusted by this tentative time delay.For example, the data point in the bottom left corner in the right panel of the figure corresponds to the two observations in Figure 7, showing the power-law flux from ObsID 1117010108 and the reflected flux from ObsID 1117010121.Similarly, the other points in the panel correspond to any pairs of observations whose separation in time is within the 1σ errors of the JAVELIN lag.While fits to the original data show no correlation, as was evident in Figure 4, fits to these lag-adjusted points have a best-fit slope of m = (2.2 ± 0.9) × 10 −2 .This indicates a positive correlation between power-law flux and reflected flux after a characteristic lag.This positive slope admittedly excludes zero at the level  6), and the increasingly lighter blue contours represent 1σ, 2σ, and 3σ levels of confidence, respectively.The lag time, top-hat width, and τ are all in units of days, and the scale and σ are unitless.The inset shows a more detailed view of the distribution of lag times found by JAVELIN; the vast majority of samples are located around a single peak.The solid black line shows the mean lag value, t = 16.37 days, and the dotted lines bound a 1σ credible interval. of 2.4σ, but this still shows that this tentative lag is an improvement over the noncorrelation of Figure 4 that first motivated our investigation.

Fits to the Fe Kα Line at High Resolution
For fits to the Fe Kα line using Chandra data, we initially considered a simple model consisting of an absorbed cutoff power-law and a simple Gaussian line, ztbabs * cutoffpl +zgauss.In this and all subsequent models, the redshift of the absorption, power-law continuum, and line were fixed to that of the host galaxy, z = 0.008.The power-law cutoff energy was fixed at E cut = 100 keV following Fabian et al. (2015), and we fixed the Gaussian line energy at 6.400 keV.Our results confirm those previously reported by Shu et al. (2011): we measure a line width of FWHM = 2400 ± 600 km s −1 .Assuming (a) that the broadening is due to Keplerian orbital motion and (b) that we observe the accretion flow at an inclination of θ = 60°, the line broadening would imply a radius of r GM c 6.4 10 The assumed inclination angle is broadly appropriate for Seyfert-2 AGN, but also the most probable angle at which a plane is likely to be viewed in 3D, and so a useful point of reference.
The single Gaussian modeling of the Fe K emission line is basic; the feature is actually composed of two lines.Moreover, even two Gaussians would fail to account for the possibility of a Compton shoulder on the line owing to 180°scattering (this feature extends down to 6.25 keV; for a review, see Gallo et al. 2023).Moreover, even if the line is produced at distances of r ; 10 3-4 GM/c 2 , special relativistic distortions to the line can still be important.We therefore considered fits with pexmon, which self-consistently describes Compton-thick reflection from neutral gas, including the Fe Kα and Fe Kβ lines as doublets, and with mytorus (Murphy & Yaqoob 2009).The latter model is a table calculated at very high (2 eV) resolution, wherein the column density of the line-emitting gas is a variable parameter.This parameter is sensitive to the relative strength of the line core and Compton shoulder.We used the mytorus line table calculated for a power law with a cutoff energy of E cut = 100 keV.
Closely following fits to the line profile in Chandra observations of NGC 4151 (Miller et al. 2018), we allowed for broadening by convolving pexmon and mytorus with rdblur (Fabian et al. 1989).This is an analytic blurring function that self-consistently accounts for relativistic Doppler shifts and gravitational redshifts around a black hole.It explicitly assumes a zero-spin black hole, which is unlikely to describe any astrophysical black hole, but the potential far from a spinning black hole is very similar to the zero-spin case.The four parameters of the rdblur model are the line emissivity (we assumed R −3 as for an isotropic point source), the inner line-production radius (allowed to vary freely in our fits), the outer line-production radius (arbitarily fixed at r out = 10 6 GM/c 2 ), and the inclination angle (assumed to be θ = 60°in all fits).The XSPEC formulae for the two models are as follows: TBabs * cutoffpl+rdblur * pexmon and TBabs * cutoffpl+rdblur * mytorus.Difference spectrum of panels (b) and (a), with the range 6.1-6.7 keV highlighted in red.Note that the spectra are similar at 3 keV and differ substantially for >3 keV.The difference in time between the medians of these two observations is 16.20 days, which is comparable to the lag time found by JAVELIN.
The results of our fits with these models are detailed in Table 3.They include one fit with pexmon and a series of four mytorus models that sample Compton-thin and Compton-thick gas.The most important result from these fits is that all measurements suggest a line-production radius of approximately r = 3 × 10 4 GM/c 2 , with approximately 30% uncertainties.The best overall model is the mytorus fit that assumes Comptonthick gas with N H = 3 × 10 24 cm −2 ; it measures a line-production radius of r GM c 2.9 10 0.7 . This model is better than that assuming Compton-thin gas, with N H = 1 × 10 23 cm −2 , by ΔC = −18 for the same number of degrees of freedom.A comparison of these two particular models is shown in Figure 9.While our models effectively capture the details of the atomic features in the Fe K band, it must be noted that the neutral absorption is treated in a simplified manner; the associated column densities should be regarded as fiducial and should not be contrasted with proper fits to the full passband.
Finally, we included the ionized absorber model that was fit to the time-averaged NICER spectra to test whether this has any effect on the properties of the Fe Kα emission-line modeling.We find that the inclusion of an ionized absorber has no impact on the line properties derived using a Gaussian, pexmon, or -.This nominally signals that the ionized absorption in NGC 4388 is indeed an outflow.However, the absorber only improves the fit by Δχ 2 = 12 for ν = 3 degrees of freedom, signaling that it is not highly significant.

Discussion
We have analyzed NICER monitoring observations of the nearby Seyfert-2 AGN, NGC 4388.The two densest monitoring periods were selected, and the summed spectrum from each was analyzed in detail.We find no evidence for strong variations in the column density or covering factor of neutral obscuration along our line of sight to the central engine.However, we measure significant changes in both the direct and reflected flux, and in highly ionized absorption that is evident in the Fe K band.In single monitoring observations, only a weak correlation is observed between the reflected and direct  Note.Spectral fits to the Fe K region of the combined first-order Chandra HEG spectra of NGC 4388.The data were binned using the optimal binning scheme of Kaastra & Bleeker (2016) and fit over the 5.0-7.5 keV range.The line was alternately described using the mytorus model or pexmon.The former allows the column density of the line-emitting gas to vary, while the latter assumes Compton-thick reflection.Both models were blurred using rdblur assuming an isotropic point source emissivity of J ∝ r −3 and an inclination of θ = 60°.The K parameter is the flux normalization of each model; in both cases, this is the flux normalization of the power law, ph cm −2 s −1 keV −1 at 1 keV.The models all agree on an inner line-emission radius that is formally consistent with the lag radius found by JAVELIN, and they prefer a solution with Compton-thick emitting gas.
flux.A search for lags between the direct and reflected flux reveals a potential lag of t 16.37 0.38 0.46

= -+
days; a simple conversion of this lag to distance gives r 1.374 10 0.032 0.039 2

= -
+ pc = 3.4 ± 0.1 × 10 4 GM/c 2 , consistent with the optical BLR in unobscured AGN.In this section, we discuss the strengths and weaknesses of our analysis and results, their implications, and plausible future observations that can build on these results.
The total spectral model that we have fit to the summed spectra from the first and third epochs (see Table 2) and the individual spectra within these monitoring periods is based on fits made to other highly absorbed Seyfert-2 AGN (see, e.g., Kammoun et al. 2020).It is likely that the soft X-ray spectrum is at least partly comprised of photoionized emission, whereas our modeling assumes collisional plasma emission.However, for the purposes of this analysis, the inclusion of the simple apec components allows the flux to be characterized well, facilitating the detection of variations in the spectrum above 3 keV.
Our spectral model also includes a component that measures highly ionized X-ray absorption in the Fe K band, following the detection of ionized Fe K absorption lines with, e.g., Suzaku and NICER (Shirai et al. 2008;Miller et al. 2019).This absorption is only evident in the summed spectra from the first and third epochs; spectra from individual exposures lack the required sensitivity.It is likely to be a low-velocity X-ray warm absorber.When the ionization parameter formalism of ξ = L/nr 2 is used, along with an assumed bolometric correction of 15 (broadly consistent with Vasudevan & Fabian 2007) and a unity filling factor such that nr 2 = Nr, the fits to the first and third epochs imply radii of r 0.92 0.27

= -+
) for the latter.These large radii are broadly consistent with the lack of a significant outflow velocity; however, the gas could originate at much smaller radii if the filling factor is small.We note that the observed spectra have flatter power-law indices than was assumed by our photoionization table models, leading to systematic errors that are likely small compared to the uncertainties on the filling factor.Future fits with the pion photoionization package can deliver more accurate results (Miller et al. 2015;Mehdipour et al. 2016).
Fits to the spectra from individual observations with the models developed for the summed spectra yielded useful results.A significant positive correlation is observed between the power-law photon index, Γ, and the 2-10 keV power-law flux.This correlation is typical of Seyfert-1 AGN, but has not previously been observed in NGC 4388 because its strong nuclear obscuration made it difficult to study the direct continuum in prior monitoring efforts.Only a weak correlation is found between the reflected flux and direct power-law flux in the individual exposures; Figure 4 shows that a significant level of scatter is present in the relevant data.This suggests that light-travel times may not be properly accounted for, leading us to use JAVELIN to systematically search for lags between the direct continuum and reflection components.

= -+
day lag between the direct continuum and reflected flux measured with JAVELIN is partially driven by a sharp drop in the continuum flux around MJD 58099 (see Figure 6).The delayed dip is less pronounced in the reflected flux.However, Figure 7 compares individual spectra obtained on MJD 58099 and MJD 58115, the two observations corresponding to these dips.In the former, the line flux is unusually strong owing to the drop in the direct continuum; in the latter, the Fe K line is not evident at all.This indicates that JAVELIN likely identified a delay time that is strongly grounded in the data.However, this lag must be regarded as tentative because the monitoring period is only twice the lag time.Moreover, we have pragmatically used a simple, neutral reflector and applied it over a limited passband, although some prior observations find evidence of ionized reflection in NGC 4388 (Risaliti et al. 2002).Future NICER observations that monitor over a period that is several times the duration of the tentative lag are required to confirm it.Corresponding observations with NuSTAR will help to constrain the reflection.
In the meantime, current high-resolution X-ray spectroscopy has provided one independent check.Fits to archival highresolution Chandra spectra of NGC 4388 measure line broadening consistent with a radius of r GM c 6.4 10 2.3 4.7 4 2

= -
+ for realistic assumptions.This is nominally a factor of 2 larger than indicated by the tentative lag, although the 1σ error is quite close to the radius indicated by the tentative lag.More physical models for the Fe K region suggest line-production radii that are formally consistent with that indicated by the tentative lag; the best-fit model in Table 3

= -
+ .The physical models also prefer that the line-emitting gas be Compton thick.If this is robust, it indicates that the line is produced through X-ray reflection, as classically envisioned.It is nominally easiest to associate the line-production region with the cold, optically thick accretion disk.The fact that the fits pick out a specific radius suggests that the profile of the disk may not be locally flat, but rather have some vertical extent.It is possible that a warp could provide the extra solid angle to anchor the narrow Fe K line at a specific radius.However, this scenario is likely also consistent with a wind that is infused with cold, optically thick clumps.One model of the optical BLR posits that it is a wind that is initially launched by radiation pressure on iron-rich dust, which has a larger cross section than pure gas (e.g., Czerny et al. 2015;Baskin & Laor 2017).Our results are at least qualitatively consistent with Figure 9.The summed Chandra HEG spectrum of NGC 4388 in the Fe K band, shifted to the frame of the host galaxy.The data represent 269 ks of net exposure and were binned using the optimal scheme of Kaastra & Bleeker (2016); they were then binned additionally here for visual clarity.The models are blurred line functions describing power-law irradiation of cold gas, rdblur * mytorus.A fit with mytorus appropriate for Compton-thin gas (N H = 1 × 10 23 cm −2 ) is shown in blue.A superior model appropriate for Compton-thick gas (N H = 3 × 10 24 cm −2 ) is shown in red; it provides a better fit to the Compton shoulder on the Fe Kα line that extends to 6.25 keV.
this model, with some important caveats.For a wind that is viewed primarily in emission-by force, primarily on the far side of the central engine-not to appear significantly redshifted, its velocity must be low.This effectively requires that the Fe K line flux must originate close to the base of the wind, before it is accelerated (this is broadly consistent with a dusty origin for the BLR; Czerny et al. 2015).
We note that an equally distinctive lag is not identified via JAVELIN when the Fe K line flux alone-via, e.g., a Gaussian model component-is used rather than the entire reflection component.Peaks within the lag spectrum are broad and indistinct from other potential lags.This is likely due to the limited sensitivity of our data.It may also be due to related difficulties in isolating the line from the continuum, which includes both direct power-law flux and reflected flux.In the limit of data with modest sensitivity, it is likely more pragmatic to rely on models that self-consistently account for the Fe K line and reflected continuum.Searches for lags between the direct power-law flux and continuum bands that exclude the Fe K region do not show distinctive lags either.
In the near future, observations of NGC 4388 with XRISM (Tashiro et al. 2018) can effectively test and extend all of our results.The calorimeter spectrometer on board XRISM is expected to deliver a resolution of just 5 eV (E/ΔE ; 1300 at the Fe Kα line), coupled with a collecting area 10 times greater than the Chandra HEG.Moderately deep observations with XRISM will reveal the details of the dynamical content of the Fe K line in NGC 4388, and they may reveal a complex of lines associated with distinct geometries.At this high sensitivity, any Compton shoulder on the Fe K line will be readily detected, enabling a clear measurement of the gas column density.As the mission moves beyond initial observations of numerous targets, it may be possible to monitor NGC 4388 with XRISM to explore the line properties while also testing the tentative lag that we have found.
The maser emission in NGC 4388 previously enabled a precise measurement of its black hole mass (Kuo et al. 2010).However, the situation is not so fortunate in many other Seyfert-2 AGN, wherein there is no maser emission and the optical BLR is not visible.Our results suggest that black hole masses in these systems could plausibly be measured by adapting optical BLR reverberation mapping techniques, replacing the optical continuum and, e.g., Hβ line with the X-ray continuum reflection spectrum and narrow Fe K line.We expect that efforts like this will need to be limited to Comptonthin AGN, as the line emission in Compton-thick AGN may be dominated by large-scale outflows and shaped by related motions (e.g., Kallman et al. 2014).

Figure 1 .
Figure 1.Time-averaged spectra of NGC 4388 from NICER.Both sets of data are individually fit to the model as described in Section 3.1, with the best fit for the early epoch (2017-18) shown in red and the fit for the late epoch (2022) shown in blue.Left: Full spectral fits from 0.6 to 9.0 keV.For each epoch, the total best-fit model is shown with a solid line, and individual additive model components are shown with different line styles, as indicated in the legend.Background estimates generated by nibackgen3C50 are plotted with stars.Right: The power-law-dominated parts of both spectra shown from 4.0 to 9.0 keV, repeated here to highlight the Fe Kα line in each spectrum as well as a small absorption feature captured by ionabsorb, which is more prominent in the early epoch.Both plots have been rebinned for visual clarity.

Figure 2 .
Figure2.NICER light curve of NGC 4388 in the energy range 3-10 keV, indicative of the behavior of the power-law continuum and reflection components of the X-ray spectrum.Each bin represents the average count rate over a period of 1000 s.The two epochs of observations used in this analysis are shown in adjacent panels, with over four years separating them.Note the greater breadth of the first epoch of observations, which made it better suited for our investigation of day-timescale variability, as detailed in Section 3.4.

Figure 3 .
Figure3.Photon index vs.unabsorbed power-law flux, f PL,2-10 , for individual observations in both epochs.The dashed black line is Γ = 1.115, indicating the hard lower bound imposed when fitting due to the limitations of pexmon.The shaded pink region shows sample fits from a Monte Carlo procedure, fitting both epochs jointly.The samples were drawn from the 68% contour around the mean values for the best-fit slope and intercept, so this band represents a 1σ region of uncertainty.
this is treated as light-travel time along with a black hole mass of M BH = 8.4 ± 0.2 × 10 6 M e(Kuo et al. 2010), this corresponds to a radius of r =

Figure 4 .
Figure 4.The reflected flux (modeled with pexmon) vs. the direct power-law continuum flux for NGC 4388.Individual observations of the first (second) epoch are shown in red (blue), and their time-average spectrum is shown in magenta (cyan).Although the late epoch tends toward lower power-law flux values (see Figure 1), there is no strong correlation regardless of epoch.

Figure 5 .
Figure 5. Parameter distributions from a joint fit of the direct continuum and reflected flux light curves using JAVELIN.The light curves used encompass the first epoch of observation (see Figure6), and the increasingly lighter blue contours represent 1σ, 2σ, and 3σ levels of confidence, respectively.The lag time, top-hat width, and τ are all in units of days, and the scale and σ are unitless.The inset shows a more detailed view of the distribution of lag times found by JAVELIN; the vast majority of samples are located around a single peak.The solid black line shows the mean lag value, t = 16.37 days, and the dotted lines bound a 1σ credible interval.

Figure 6 .
Figure 6.Light curves of NGC 4388 based on calculated power-law flux ( f PL,2-10 , top panel) and reflected flux ( f refl,2-10 , bottom panel).All plotted points were derived using the cflux command in XSPEC on either the powerlaw component of the model (cutoffpl) or the neutral reflection component (pexmon), which includes the neutral Fe Kα line.The solid blue and red lines are the best-fit mean light curves based on fitting with JAVELIN, and the shaded blue and red regions represent regions of uncertainty around these bestfit models.The two dashed gray lines highlight the observations ObsID 1117010108 and 1117010121, whose spectra are presented in Figure 7.

Figure 7 .
Figure 7.An example of short-term variability between two individual NICER observations of NGC 4388.(a) Spectrum of ObsID 1117010108 (MJD 58099), with the best-fit model shown in red.(b) Spectrum of ObsID 1117010121 (MJD 58115), with the best-fit model shown in blue.In both panels, the best-fit model of the other observation is shown with a dotted line for reference.(c)Difference spectrum of panels (b) and (a), with the range 6.1-6.7 keV highlighted in red.Note that the spectra are similar at 3 keV and differ substantially for >3 keV.The difference in time between the medians of these two observations is 16.20 days, which is comparable to the lag time found by JAVELIN.
mytorus.The best-fit values and confidence intervals are entirely consistent.The best-fit absorber was measured to have a column density of N 2

Figure 8 .
Figure 8. Fits to plots of reflected flux ( f refl,2-10 ) vs. power-law flux ( f PL,2-10 ) for the first epoch of NICER observations.The shaded pink regions encompass sample fits drawn from the 68% contour from Monte Carlo fitting procedures, representing an approximate 1σ region of uncertainty in the best fit in each panel.Left: Original data, where each data point corresponds to a distinct observation.The best fit to these data has a slope consistent with zero, indicating no correlation between the two quantities.Right: Lag-adjusted data, where each data point corresponds to the power-law flux of one observation and the reflected flux of a different observation taken within the range of t 16.37 0.38 0.46

Table 2
Results of Time-averaged Epoch Spectral Fitting

Table 3
High-resolution Chandra Results