Modeling optical and UV polarization of AGNs V. Dilution by interstellar polarization and the host galaxy

One of the main challenges for polarimetric observations of active galactic nuclei (AGN) is to properly estimate the amount of parasitic light that contaminates the polarization signal. Removing this unpolarized flux is a complex task that has only been achieved in a couple of objects. In this fifth paper of the series, we present a new version of the Monte Carlo code STOKES that accounts for dilution by interstellar polarization and host starlight in radiative transfer modeling. We upgrade our code by including spectral energy distribution (SED) templates for a set of representative host galaxies. The unpolarized light emitted by those hosts alters the observer polarization while being coherently radiatively coupled to the AGN structure. We also include in our analysis tool a routine that may add, depending on the user's objectives, an interstellar component. Using a generic AGN model, we illustrate how interstellar polarization and starlight dilution impact the observed polarimetric signal of AGN. We apply our code to NGC-1068, an archetypal edge-on AGN and demonstrate that STOKES can reproduce its SED, the expected wavelength-dependent polarimetric signatures, and the observed high-angular resolution polarimetric maps. Using the flexibility of the code, we derive several intrinsic parameters such as the system inclination and the torus opening angle. The new version of our publicly available code now allows observers to better prepare their observations, interpret their data and simulate the three-dimensional geometry and physics of AGN in order to probe unresolved structures. Additionally, the radiative interaction between the host and the AGN can be used to probe the co-evolution of the system.


Introduction
Active galactic nuclei (AGN) are very compact and highluminosity engines located at the center of their host galaxies. Since the innermost structures of AGN cannot be resolved by current UV, optical, or IR observational techniques, morphological and kinematic models have to be derived indirectly, most often from broadband spectroscopy and variability studies.
In so-called thermal radio-quiet AGN, the emission is due to accretion onto a supermassive black hole (SMBH),with typical masses M SMBH >10 6 M (Salpeter 1964;Lynden-Bell 1969). The accretion disk produces the optical and UV continuum radiation (see Shields 1978;Shakura & Sunyaev 1973;Pringle et al. 1973). The disk axis points out a preferred direction that for the remainder of this paper we call the polar direction of the system, to be distinguished from the equatorial directions that enclose the plane of the accretion disk. Observationally, thermal AGN can then be divided into two classes: type 1 sources are assumed to be viewed at a polar line of sight and show Doppler-broadened Balmer emission lines produced in the so-called broad line region (BLR), whereas type 2 sources are viewed at equatorial lines of sight, at which the optical view of the BLR is obscured by a dusty region commonly called the dusty torus (Antonucci 1993). It extends out to a spatial scale of a few parsecs for a 10 7 M SMBH. Only in very few objects do the outer parts of the circumnuclear dust start to be resolvable by mid-IR interferometry techniques (Tristram et al. 2007(Tristram et al. , 2014Raban et al. 2009). Gaskell (2009) introduced the so-called bird's nest paradigm for the central region of thermal AGN, adopting the idea that its overall structure is clumpy. Clumpy models successfully reproduce the dust emission seen in the IR (Nenkova et al. 2008;Siebenmorgen et al. 2015;Hönig & Kishimoto 2010). The BLR is expected to be fragmented and turbulent, but radially stratified, as can be inferred from the velocity structure of low-and highionization emission lines (Gaskell 1982;Gaskell et al. 2008;Gaskell & Goosmann 2013). The BLR is expected to form inside the sublimation radius, possibly from inflowing torus material. It is therefore embedded in the torus funnel. Since the outer parts of the accretion disk around the SMBH must be self-gravitating, it is straightforward to identify it with the inner parts of the BLR. A coherent picture for the circumnuclear matter of AGN can be found with a morphology that is non-spherical, but symmetric Article number, page 1 of 8 arXiv:1712.01147v1 [astro-ph.HE] 4 Dec 2017 A&A proofs: manuscript no. 31331_P_A_Rojas_Lobos with respect to the polar axis, while being fragmented on small spatial scales.
Not all matter that is pulled in by the black hole is finally accreted. Depending on the accretion state of a given thermal radio-quiet AGN, a significant fraction of matter may be ejected along the polar directions as moderate to ultrafast winds (Pounds & Page 2006;King & Pounds 2003;Tombesi et al. 2012). Depending on the distance from the SMBH, these winds are mildly to strongly ionized and scatter much radiation that is emitted by the central black hole. The common picture of the inner structure of AGN was recently described by Marin (2016), and the reproduction in Fig 1 depicts the regions inside the AGN (not to scale): the black hole sits at the center and is surrounded by the accretion disk, the BLR, and the circumnuclear dust. Along the polar directions, ionized winds and the narrow line region (NLR) are located, which in radio-loud objects may be pierced by collimated jets.
Polarimetry is a valuable tool to probe the innermost structures of AGN. In the past, optical spectropolarimetry gave proof of the unified theory of AGN when Antonucci & Miller (1985) discovered polarized broad Balmer lines in Seyfert 2 galaxies and thus pointed out a relation between the observer's inclination with respect to the polar axis and the position angle of optical continuum polarization in type 1 and type 2 AGN. It was shown in that the typical polarization properties in type 1 objects can be produced by scattering in a flattened distribution of material situated around the accretion disk (Antonucci 1984;Smith et al. 2004;Goosmann & Gaskell 2007).
More recently, AGN polarimetry has been brought to a different level by introducing the technique of polarization reverberation mapping  for the Seyfert 1 nucleus NGC 4151). In this method, the time-dependent total (polar-ized+unpolarized) continuum emission is cross-correlated with the polarized emission of a given AGN. A characteristic time delay is found between the two components. When we assume that the continuum polarization is induced by scattering of the continuum radiation by surrounding structures, the reverberation places constraints on the light travel distance between the emitting source and the reprocessing mirror(s).
In this work, we present initial modeling to compare with current and future campaigns of polarization reverberation mapping. We conduct radiative transfer simulations for the continuum radiation in different AGN geometries, taking the time dependence of the polarized radiation in the optical band into account. The goal of this research is to numerically explore the new method and indirectly resolve the innermost region in AGN by polarization reverberation. To do this, we provide simulated results for polarization time lags from different AGN constituents as a function of inclination. The paper is organized as follows: in Sect. 2 we present the details of our model and specific parametrization. In Sect. 3 we present and analyze our results. Finally, in Sect. 4 we discuss our results and follow-up projects.

Modeling the circumnuclear dust scattering
Following the observational example of Gaskell et al. (2012), we here attempt to theoretically investigate the details of the inner AGN structure by studying the expected time lag between the total and the polarized continuum flux. The presence of the obscuring torus on the line of sight should be derived from the time lag, together with its geometrical shape and composition. Figure 2 shows a sketch of our different modeling options. In the following, we call "dusty-torus" a toroidal region with an elliptical cross section, and "flared-disk" a region with  Marin (2016). The cartoon shows that AGN with narrow and broad Balmer emission lines (type 1, i = 0 − 60 • ) and those with only narrow lines (type 2, i = 60 − 90 • ) derive from the same morphology, but are viewed from different angles with respect to the AGN axis. This is the original basis of the unified model. . a wedge-shaped cross section. We investigate polarization timing for four versions of a circumnuclear geometry: a uniformdensity dusty torus, a uniform-density flared disk, a clumpy torus, and a clumpy flared disk. Every geometry is tested for two types of optically thick dust: AGN-dust (Gaskell et al. 2004) and Milky Way dust (so-called MRN dust, according to Mathis et al. 1977) at different optical depths. We increase the complexity of the model by adding ionized outflows in the polar direction in Sect. 3.2 and a scattering ring, responsible for the specific polarization position angle seen in type 1 AGN (Smith et al. 2004), in Sect. 3.3. The latter region is fully ionized and is located between the source and the circumnuclear dust region. As the model setup becomes more sophisticated, we follow the evolution of the polarization time lags.
We applied an extended version of the Monte Carlo radiative transfer code stokes. This version includes the timing information of the registered photons. The code was initially written by Goosmann & Gaskell (2007) and further developed by Marin et al. (2012) and Marin et al. (2015); an intermediate version of the code is publicly available 1 . It simulates radiative transfer in different geometries for emitting and scattering structures Mie Mie and determines the net Stokes parameters as a function of wavelength, time, and the observer's line of sight. Subsequently, it is possible to obtain the polarization fraction, position angle, polarized flux, and total flux from the Stokes parameters. For compatibility with paper I (Goosmann & Gaskell 2007), we define the polarization position angle to be 0 • when the polarization angle is perpendicular to the projected symmetry axis of the model. When the polarization angle is equal to 90 • , it is said to be parallel to the projected symmetry axis of the model. To compute the time lags, the code compares the path between the direct flux from the source and the trajectory of the photons that have been scattered from the dust and/or electron regions. For each viewing direction, the time lag is measured relative to the plane of sky, which goes through the origin of the central source and is oriented perpendicularly to the line of sight. The lag is expressed in pc/[c] and represents the light traveling distance. We modeled the flared-disk and toroidal geometries shown in Fig 2 assuming uniform-density and clumpy dust distributions. Two types of dust were tested: AGN-dust composed of 85% of silicate and 15% of graphite, and Milky Way dust composed of 62.5% of silicate and 37.5% of graphite. The grain radii lay between 0.005µ (0.005µ) and 0.25µ (0.2µ), with a distribution n(a) ∝ a s with s = −3.5 (s = −2.05) for Milky Way (AGN) dust. The filling factor of the clumpy regions was set to 25% (Marin et al. 2015). At the center we defined a point source with unpolarized flux that emitted a power law spectrum F ν ∝ ν −α with α = 1. For the ionized regions we assumed an inner radius of 8 light days such as found by Gaskell et al. (2012) in the case of NGC 4151, but without detailed modeling. In the case of the circumnuclear dusty regions, we adopted the observational constraints derived by near-IR observation of NGC 4151, namely 0.061 pc, according to Kishimoto et al. (2013).
The outer radius of the torus and its half opening angle were based on the findings by Ruiz et al. (2003), Table 2, who conducted near-IR polarimetry observations of the same object. We defined an external radius of 15.061 pc and a half opening angle of 60 • from the disk symmetry axis. The dust is opaque, with a radial optical depth of ∼ 150 inside the equatorial plane. The modeling results can be evaluated over the optical wavelength range and as a function of the observer's viewing angle.
In a second step of this study, we added polar winds to this baseline model. They have the geometrical shape of doublecones and are assumed to be ionized. In the model they are filled with electrons with a radial Thomson optical depth of 0.03 and  0.3. The winds have a common interface with the circumnuclear dust at the inclination of its horizon. Finally, we added a third region, an equatorial scattering ring located between the source and the circumnuclear dusty region. This area, modeled as a flared disk, was found to be necessary to reproduce the parallel polarization position angle observed in the majority of Seyfert 1s, see Smith et al. (2004). Following the composition constraints from Paper II, we opted for a radial optical depth of unity. The inner radius of the scattering ring was set to 0.0067 pc (i.e., eight lights days, Gaskell et al. 2012) and the outer radius to 0.0577 pc, so that the disk and the torus did not mix (see Table 1 for a summary of the parameters).

Results and analysis
All modeling cases were tested for the two different dust prescriptions: Milky Way dust (Mathis et al. 1977), and AGN-dust (Gaskell et al. 2004). A detailed comparison of the effect of the dust prescription on the polarization results was given in Goosmann & Gaskell (2007). The new simulated observable added in this work, the time lag of the polarized emission, is not much af-Article number, page 3 of 8 A&A proofs: manuscript no. 31331_P_A_Rojas_Lobos fected when the dust prescription is modified (Fig 3). The small differences between the two cases are due to the specific albedo and scattering phase functions of the two dust types. Since there is no major effect on the results, we only present modeling for the case of standard Milky Way dust in this paper. The response in terms of flux, polarization, and time lag only mildly changes with wavelength (Fig. 4). Only at type 2 viewing angles does the uniform dusty torus show a diminution in flux and an increase in polarization around the 2175Å absorption feature. Overall, the results shown in Fig. 4 confirm our previous work: for a type 1 inclination, we obtain a higher flux level and lower polarization than at type 2 viewing angles. The new simulated observable in this work is the average time lag of the radiation. We note that this time lag is taken relative to the direct emission coming from the central source. The units are set to pc/c, which allows us to interpret the time lag in a straightforward manner in terms of light traveling distances. This time lag is always higher at type 2 viewing directions than for type 1, which is expected because below the torus horizon, we only detect (multiply) scattered radiation that has taken a non-direct path from the central source to the observer.
Except for these features, the results in terms of spectral flux, polarization, and time lag do not vary much with wavelength. For the remainder of this paper we therefore show our results for the V band centered on 5500Å and as a function of viewing angle. All models show significant differences between face-on and edge-on views.

Baseline model for circumnuclear dust
Our baseline model contains the axisymmetric dust region with a half opening angle of 60 • as measured from the axis (see Fig. 2). The modeling results are plotted in Fig. 5. At face-on inclinations, all four models show roughly the same spectral shape in total flux, polarization degree, and time lag. The normalization of the spectral flux is almost the same for all four morphological setups because at these viewing angles the radiation predomi- nantly comes from the central source. Still, a face-on observer finds a higher time lag and lower polarization degree for the clumpy dust distribution than for the uniform density cases. This is due to radiation scattered into the type 1 lines of sight. In a clumpy circumnuclear environment, the escaping photons experience on average a higher number of scatterings before they end up on a type 1 line of sight (see Marin et al. 2015). For optically thick uniform-density regions, the scattered radiation at face-on viewing angles is mostly dominated by the first-scattering component that is more strongly polarized and accumulates less of a time lag. At edge-on inclination, the observer's line of sight is obscured by dust, and the flux strongly decreases. The polarization degree slowly increases with increasing viewing angle for all four geometries. At all type 2 viewing angles, the time lag accumulated for a clumpy or uniform torus geometry is much shorter than for the case of a flared dust distribution.
This difference is related to the shape of the inner surfaces. For a torus, the inner surface is convex, while for the flared disk it is parametrized as a concave shape. As discussed in Goosmann & Gaskell (2007), when comparing the obscuration efficiency between compact and extended dusty tori, we find again here that a more convex shape of the inner surfaces makes it more likely for a photon to escape after only a few scattering events. Therefore, the accumulated time lag for the case of a toroidal dust distribution remains short even when the distribution is clumpy. For the case of our flared geometry, the inner surface is parametrized as a sphere segment against which the photons are injected from the inside. In general, multiple backscattering inside this sphere segment is necessary for the photons to finally escape, and if they do, the escape cone is strictly limited by the half opening angle of the flared disk as measured from its symmetry axis. Figure 5 shows that the rise in time lag and polarization degree is steep around this inclination for the flared dusty disk, while for the toroidal geometry, the rise is much more shallow. For clumpy dust distributions, the horizon of the torus and of the flared disk is more blurry, which decreases the slope of the polarization degree and the time lag even further. The polarization position angle remains at 90 • (i.e., parallel to the projected symmetry axis, following the usual convention in this paper series) over all viewing angles for models without polar winds, with the exception of a clumpy flared geometry, for which the polarization position angle stays at 90 • at face-on viewing angles, then shifts to 0 • (perpendicular polarization) around the transition between type 1 and type 2 to return to parallel polarization at more extreme type 2 viewing angles. We argue that this behavior is related to the actual size of the clumps and to their distribution along the line of sight.

Adding polar outflows to the baseline model
Next, we added ionized winds stretched along the symmetry axis to our model. As in previous work of this series, the wind geometry was realized by a double-cone filled with electrons of a radial optical depth of 0.03 or 0.3 (see Fig. 6). When adding a low-density wind, the type 1 cases exhibit almost the same spectral flux level and polarization degree as for the baseline models. The additional scattering in the polar direction has a more significant effect on the time lag recorded at face-on viewing angles. It increases by a factor of roughly 2 when compared to the case without polar winds. When we add low-density polar winds, the polarization angle changes from 90 • (parallel) at face-on view to 0 • (perpendicular) in face-on for uniform models, while for a clumpy distribution, the polarization anlge remains at 0 • for both geometries at all viewing angles.
The situation is slightly more spectacular for the edge-on case. The presence of even a low-density wind scatters a significant part of the radiation "around" the circumnuclear dust and thereby diminishes the strong absorption efficiency of the flared dusty disk.
When the winds are added, we find longer time lags with a more shallow dependence on viewing angle than for the baseline model. This is related to the fact that before escaping, the photons cross a second scattering region that adds more scattering events and increases the light travel path. The polarization position angle remains at 0 • over all viewing angles for models with higher-density polar winds. This behavior, linked to the optical depth of the polar outflows, is in agreement with the polarization state observed for six Seyfert 1 AGN by Batcheldor et al. (2011). Since Smith et al. (2002), it has been a well known fact that a small fraction of type 1s shows optical polarization spectra similar to those of Seyfert 2 galaxies in which polarized broad lines are detected. In the compendium of AGN presented in Marin (2014), it was shown that these peculiar objects are associated with polarizations of a few percent and represent a rather small subclass of Seyfert 1s. Similarly to our previous conclusion on the dependence of optical depth versus inclination, we note that a clumpy dust distribution is not very sensitive to inclination angles. Across type-1 and moderate type 2 viewing angles, the time lag distribution remains rather flat for both geometrical clumpy configurations; the time lag only increases for extreme type 2 lines of sight. This ambiguity is stronger for the wind with higher optical depth. In contrast to this, the time lag remains a discriminator between type 1 and type 2 inclination for the flared-disk geometry and is far higher with a uniform dust distribution.

Adding a scattering ring to the baseline model
We now include a third scattering region, an equatorial, ionized disk located between the source and the circumnuclear region. Compared to a model without equatorial electron disk (Fig. 5), our new modeling shows that only the polarization properties of light and its time lags are affected (see Fig. 7). An additional scattering component that does not absorb radiation has almost no effect on the flux level observed at various inclinations. However, the electron disk stabilizes the degree of polarization observed in type 1 viewing angles, independently of the circumnuclear dust structure. The disk also sets the polarization of radi-Article number, page 5 of 8  ation to 90 • at type 1 orientations, as has been observed (Smith et al. 2004). More importantly, the proximity of this region to the central source means that the observed time lags are shorter for type 1 inclinations. A fraction of the input radiation, depending on the solid angle covered by the opening angle of the electron ring, is directly scattered toward the observer, inducing a shorter time lag. For higher inclinations, radiation passes through the optically thin medium and scatters inside the optically thick dust structure, so that the effect of the equatorial ring becomes almost undetectable.
When we include polar outflows in addition to this scattering ring, see Fig. 8, we reach similar conclusions. The addition of an equatorial electron ring remains marginal in terms of total fluxes, but strongly helps the polarization degree and angle to stabilize at typical values seen for local AGN (see Marin 2014). The degree of polarization remains lower than 1% for type 1 inclinations and sharply rises for type 2 orientations. The polarization angle rotates from parallel to perpendicular as soon as the observer's line of sight is obscured by the circumnuclear dusty region. The effect of the disk on the time lags at type 1 views is even stronger: the time lags are 5 -10 times shorter, reflecting the presence of the ring. The effect of different wind optical depths on time lags is the same as what we described in Sect. 3.2.

Summary and conclusion
In this initial study, we extended polarization modeling of AGN into the timing domain. Investigating the time lag of the polarized emission as a function of viewing angle, we showed that owing to multiple scattering effects, the polarized signal in type 2 objects lags farther behind the total intensity than for type-1 sources. Furthermore, the geometrical shape of the dust distribution such as the shape of the inner surfaces and also their optical depth has an effect on the resulting time delay at type 2 lines of sight, and much less for type 1. A significant increase in time delay is recorded when the dust distribution is clumpy.
When adding polar winds to the model, scattered radiation further increases the polarization time lag. The presence of a scattering equatorial electron disk naturally reduces the observed time lags at type 1 inclinations since a fraction of the photon flux is reprocessed onto the electrons and escapes along the polar directions.
The overall picture thus remains complex and certainly still includes degeneracy. From the modeling we present here, it is not yet straightforward, for instance, to disthinguish between the density structure of polar outflows and the level of clumpiness in the torus. Nonetheless, treating polarization as a function of time offers another independent way to constrain the model parame-  ters. For this, it is important to combine the polarization time lag, polarization degree, and position angle and compare these observables to the type of simulation results that we present here.
This study is inspired by the discovery of polarization reverberation mapping in NGC 4151 . In this paper we also refer to the initial modeling of polarization time lags with stokes. The various uniform-density tori that were tested at that time revealed that the polarization time lag at type 1 viewing angles is only slightly longer than the light-crossing time of the inner torus radius. In the following the torus was ruled out as a reverberation source because the observed time lags were shorter by a factor of five. They instead corresponded to the light-crossing time of the BLR. Polarization reverberation contains important seismological information about AGN, but as for polarization in general, careful modeling is required to correctly interpret it. Next, more central structures like the BLR and the outer parts of the accretion disk need to be included to recover shorter reverberation times. This work is currently underway and will be presented in a forthcoming publication, including depen-dencies on the density of the circumnuclear region and the cross correlation between polarized and unpolarized flux.
In the meantime, we would like to advocate systematic polarization monitoring campaigns of AGN such as those by Afanasiev et al. (2014Afanasiev et al. ( , 2015. Although these are long-term projects, they have a high potential to reveal more details of the emission and scattering geometry in AGN.