phy-4660/lag/report/report.tex

\documentclass[11pt,letterpaper]{article}

\usepackage{aas_macros}
%\usepackage{natbib}
%\usepackage{cite}
\usepackage{biblatex}
\usepackage{graphicx}
\usepackage[margin=1.in,centering]{geometry}
\usepackage{hyperref}
\usepackage{caption}
\usepackage[export]{adjustbox}
\usepackage{float}


\bibliography{/home/caes/wmu/phy-4660/adv_lab.bib}


\begin{document}



\title{Lab 3: Continuum Reverberation in NGC 5548}
\author{Otho Ulrich, Ed Cackett}
\maketitle

\begin{abstract}
    New data describing the time-variable spectrum from NGC 5548 provide the most detailed look so far at this prototypical Seyfert galaxy. A reverberation mapping of NGC 5548 is attempted using a novel frequency-domain technique. The power spectral densities of lightcurves in 19 wavelength bands are recovered. Time lags are computed from the cross spectra, but results are marred by some computational issues. Empirical fits are obtained for 6 lightcurves, and the time lags are compared with those computed by Fausnaugh et al. in STORM III. Results are mixed but overall protective of the thermal reprocessing by an accretion disk hypothesis.
\end{abstract}

%─────────────
\section{Active Galactic Nuclei}
    \label{sec:agn}
    Active Galactic Nuclei, or AGN, are a class of compact objects that are of great interest to modern astronomers. AGN are superluminous ($\sim10^{44}-10^{50} ergs\ s^{-1})$ systems that surround supermassive ($\sim 10^5 - 10^{10} M_{\odot}$) blackholes. Seyfert galaxies (also called active galaxies) are defined by these objects at their center, and the AGN often outshine all the stars in those galaxies. Quasars are too distant to resolve optically, but their spectra are typical of AGN, leading researchers to believe these are also galaxies with AGN at their core. \cite{2014A&ARv..22...72U}  AGN have a variable spectrum in time across all observed wavelengths, and with no known period. The physics underlying these variable spectra are not well-understood and astronomers are highly motivated to quantify them for many reasons.

    AGN science is a strong candidate for computing more accurate and precise values of the Hubble constant, which is a cosmological parameter that argues the rate of expansion of the universe \cite{1999MNRAS.302L..24C}. Type IA supernovae are the standard used to calculate this value now, because their spectra are well-predicted with known physics, and they are very bright. AGN are much brighter (arguably the brightest objects in the cosmos in terms of intrinsic luminosity), so if their physics were well-understood, we could use them to probe the universe beyond the supernova horizon. Toward the goal of computing a more accurate and precise Hubble constant, we study the fundamental physics at play in AGN, but that is not the only concern.

    Evidence exists that AGN are the engines that power evolution of their host galaxies; this field of study is called AGN feedback. \cite{2012ARA&A..50..455F} They also present an interesting challenge on their own merits, because they are complex systems of gasses and plasmas in a variety of thermodynamic conditions and spatial configurations. Working out those dynamics is particularly difficult because AGN in Seyfert galaxies, and especially in quasars, are too small and too distant to resolve optically. The standard model AGN structure is shown in figure \ref{fig:AGN_standard}. It is believed that an accretion disk surrounds the supermassive blackhole, and thermal emission from this disk is the primary contributer to the continuum emission.

    \begin{figure}[h]
        \center
        \includegraphics[width=3in]{ngc5548.png}
        \caption{NGC 5548: the most well-studied Seyfert galaxy. The physics underlying the time-variable spectrum of its nucleus remain an object of great interest and study. The nucleus at the center of the galaxy has approximately the same luminosity as the rest of the galaxy. \cite{ngc5548photo}}
        \label{fig:ngc5548}
    \end{figure}

    \begin{figure}[h]
        \centering
        \includegraphics[width=3.5in]{AGN_standard.png}
        \caption{Standard model of AGN structure. An accretion disk surrounds the supermassive blackhole at the center, and thermal emission from the accretion disk is believed to be the primary contributor to the continuum component of the observed spectrum. The broad line region extends beyond the accretion disk, as does an obscuring torus, narrow line region, and relativistic jet; these features are not covered by this work. \cite{2016ASSL..439..249B}}
        \label{fig:AGN_standard}
    \end{figure}

    The Seyfert galaxy NGC 5548, shown in figure \ref{fig:ngc5548}, is the most well-studied active galaxy, and the physics of its nucleus are still mysterious. An ongoing study of this object, the Space Telescope and Optical Reverberation Mapping Project (STORM), is being conducted now, and is the most in-depth study of an AGN to date. Data from the third edition of STORM \cite{2016ApJ...821...56F} is used here to quantify the continuum emission variability and time lags seen in NGC 5548.





%─────────────    
\section{Reverberation Mapping}
    \label{sec:reverbmap}

    Reverberation mapping is a popular method for measuring AGN features, and has been employed for about three decades. Since AGN cannot be resolved optically, reverberation mapping uses observed lags between wavelengths in the time-dependent spectrum to infer spatial distributions, thermodynamic properties, and blackhole mass. \cite{2014A&ARv..22...72U} Work is underway to use reveberation mapping to measure blackhole spin, as well. \cite{2016Natur.535..388K}. This study uses reveberation mapping to study the UV/Optical continuum emission in NGC 5548 as an attempt to protect the thermal reprocessing by an accretion disk hypothesis.

    Reverberation mapping is essentially mapping using light echoes. Some light escapes directly to the observer, while other light is reprocessed in the surrounding gas. Assuming that light travels with velocity $c$, additional path length traversed by reprocessed light results in an observed time delay. As an illustrative example, consider a spherically-distributed gas. This geometry results in ellipsoid isodelay surfaces with the source of light at one focus and the observer at the other. Because the distances to AGN are so great, these surfaces can be approximated as paraboloids; see Figure~\ref{fig:isodelay}. The time delay $\tau$ is associated with the geometry of the system by the equation

    \begin{center}
        $\tau = (1+cos(\theta)) \frac{r}{c}$.
    \end{center}

    If the time delay is observable, this model can be used to compute the radius of the gas cloud. AGN are more complicated than the spherical model, but the concept is the same, and the spherical model has been used to compute upper-limit radii for some AGN. A transfer function can capture the arbitrarily complex reverberation observed in AGN. This will be discussed in detail in section \ref{subsec:transferfunc}.

    \begin{figure}
        \hfill
        \includegraphics[width=3in]{isodelay.jpg}[A]
        \includegraphics[width=3in]{isodelay_detail.jpg}[B]
        \caption{[A] Isodelay paraboloids predicted for a gas cloud that is spherically distributed. This is a simple case of reverberation mapping, but AGN systems are more complicated. [B] Radius and angle dependence of a single reprocessed emission. \cite{2001sac..conf....3P}}
        \label{fig:isodelay}
    \end{figure}

    \subsection{Continuum Reverberation in the Accretion Disk of NGC 5548}
        \label{subsec:contreverb}

        The thermal reprocessing hypothesis suggests that a hot accretion disk surrounds and feeds the supermassive blackhole. This disk responds thermally to high-energy EM emission from a corona near the event horizon. The temperature increases along the disk toward the gravitational well of the blackhole, and as the coronal emission varies, the disk responds with an increase in temperature, and in turn an increase in overall emission and Wien wavelength. Collier et al. \cite{1999MNRAS.302L..24C} provide a model for the average time lag $\tau$ as a function of Wien wavelength $\lambda$. This model assumes the accretion disk consists of a distribution of blackbody radiators with temperature distribution $T(R) = T_0 \left(\frac{R}{R_0}\right)^{-\frac{3}{4}}$. The average time lag under these conditions can be predicted numerically, but for these purposes it is sufficient to conclude that $\tau \propto \lambda^{\frac{4}{3}}$. The precise nature of the accretion disk remains a topic of debate, so this should be considered only a rough estimate. A simple picture is presented in Figure~\ref{fig:simple_geometry}

        \begin{figure}
            \center
            \includegraphics[width=4in]{basic_geometry.png}
            \caption{Cross-section of the accretion disk at the center of a standard model AGN. High-energy emission irradiates the accretion disk from a corona just outside the event horizon of the supermassive blackhole. Shorter-wavelength emission emerges from deeper within the system as it responds to that emission. The observer sees longer time delays accompanying longer wavelengths.}
            \label{fig:simple_geometry}
        \end{figure}

        At the very least, we may conclude two things from these ideas. The Wien wavelength decreases with temperature, and the temperature increases toward the center of the disk, so we expect to see an increase in time delay as wavelength increases. Hotter gas emits more power at every wavelength, so we also expect to see decreasing variability as wavelength increases. If these characteristics are observed, it will protect the thermal reprocessing hypothesis. More robust analysis tools may also shed light on more detailed aspects of this hypothesis.


%─────────────
\section{Analysis of Reverberating Lightcurves}
    \label{sec:analysis}

    Observational astronomy in the optical and UV bands suffers from unevenly sampled data due mainly to weather effects and periodic orbital blindness of space observatories. Fausnaugh et al., in STORM III \cite{2016ApJ...821...56F}, published the longest and most detailed to date time-variable lightcurves in the optical and UV spectrum of NGC 5548; see figure \ref{fig:lightcurves}. Reverberation mapping in the optical and UV wavelengths has thus far involved time-domain analyses, and Fausnaugh et al. published one in their paper, included here as figure~\ref{fig:reverb_time}. These analyses have the advantage of weak sensitivity to the unevenness of the data sampling, but the disadvantage of returning only the average time lag in each band. More information than the average time lag is present in the observed data, and a frequency-domain analysis has the potential to deconvolve that information. An important test of our approach will be to confirm the average time lags we compute agree with those reported by Fausnaugh et al.

    Fourier-frequency analysis using the Fast Fourier Transform (FFT) is a common method for analysing the correlations between time-dependent signals. However, even a small deviation in the time interval between samples can create artifacts in the frequency-dependent output. X-ray astronomy suffers from unevenly sampled data in the time domain as well. X-ray astronomer Abdu Zoghbi from the University of Michigan is developing the program $psdlag$ to solve this problem by populating frequency-space bins using a maximum likelihood fitting. We apply this program in an attempt to solve the problem of unevenly sampled data in the optical/UV lightcurves reported by Fausnaugh et al.

    \begin{figure}
        \center
        \includegraphics[width=6.5in]{lightcurves.pdf}
        \caption{Optical and UV lightcurves published in STORM III. \cite{2016ApJ...821...56F}. The variable nature is very clear, and by observing the promiment bumps one can see that time delays exist between the curves. The x axis is the observation day and the y axis indicates the wavelength filter and observed flux.}
        \label{fig:lightcurves}
    \end{figure}

    \begin{figure}
        \center
        \includegraphics[width=4in]{TCCF_fausnaugh.pdf}
        \caption{Time-domain analyses of lightcurves published by Fausnaugh et al. The average time lags are computed from the light curves (left column in legend) and the best fit is compared against competing accretion disk models (right column).}
        \label{fig:reverb_time}
    \end{figure}

    

    \subsection{Time Delay and Transfer Function Theory}
        \label{subsec:transferfunc}

        A transfer function $g(t,\lambda)$ encodes the time-dependent response of each lightcurve relative to a reference lightcurve. By constructing the transfer function from the observed lightcurves, we obtain a description of the reprocessing of the coronal emission by the hypothesized accretion disk. One way to define the transfer function for each wavelength is as the function by which one convolves the reference lightcurve $x(t)$ to return the reprocessed lightcurve $y(t)$, i.e.,

        \begin{equation}
            \label{eq:time_transfunc}
            y(t) = \int_{-\infty}^{\infty} g(\tau) x(t-\tau)  {\rm d}\tau.
        \end{equation}

        The convolution theorem provides that convolution in the time $t$ domain is equivalent to multiplication in the frequency $\nu$ domain, so we can also say 

        \begin{equation}
            \label{eq:freq_transfunc}
            Y(\nu) = G(\nu) X(\nu).
        \end{equation}

        The power spectral density (PSD) provides a measure of the variability of a lightcurve. We define the PSD as $|X(\nu)|^2 = X^*(\nu) X(\nu)$, where $X^*$ is the complex conjugate of $X$. This quantity can be calculated using Fourier transforms, and suggests a useful relationship involving the lightcurves and the transfer function.

       \begin{equation}
            |Y(\nu)|^2 = |G(\nu)|^2 |X(\nu)|^2.
        \end{equation}

        We also define a cross spectrum between two lightcurves as $C(\nu) = X^*(\nu) Y(\nu)$. The argument $\phi$ of the cross spectrum is the phase between those two signals. The time lag can be computed from the phase using equation \ref{eq:timelag}. An empirical fit of the time delay in frequency space can be produced, and then transformed to the time domain. The average time delay can be obtained from this new function, but this function describes more fully the time-dependent response of the reprocessed lightcurve. Composing the many specific time-dependent responses into a single function $g(t,\lambda)$ provides the desired description of the reveberation in the system.

        \begin{equation}
            \tau(\nu) = \frac{\phi(\nu)}{2\pi\nu}.
            \label{eq:timelag}
        \end{equation}

%        Given equation \ref{eq:freq_transfunc}, the cross-spectrum can be written as

%        \begin{equation}
%            C(\nu) = X^*(\nu) G(\nu) X(\nu) =  G(\nu) |X(\nu)|^2
%        \end{equation}

%        Discrete samples of the transfer function can therefore be recovered as 

%        \begin{equation}
%            G(\nu) = \frac{C(\nu)}{|X(\nu)|^2}.
%        \end{equation}

    

    \subsection{Computational Details}
        \label{subsec:computation}

        1365\AA\ is chosen to be the reference wavelength. The $psdlag$ program is used to generate the power spectral densities, cross spectra, and time lags, all in the frequency domain. A least-squares method is used to fit an empirical function to the computed time lags. The inverse FFT transforms that function to the time domain. The most-likely average time lag reported by Fausnaugh et al. is plotted over that function for comparison.

        A tophat function describes one possible time-response mode, and examples of this function in the time domain and their corresponding FFTs to the frequency domain are plotted in figure \ref{fig:tophat_theory}. A tophat function is a rectangular function with mean $t_0$ and width $\sigma$ such that when $(t_0-\sigma) < t < (t_0+\sigma)$, this function returns 1, and otherwise it returns 0. In frequency space, it is known that the function takes the form

        \begin{equation}
            A T sinc\left(\pi \nu T\right) e^{\left(-i 2 \pi \nu t_0\right)},
        \end{equation}

        where $A$ is a normalization amplitude and $T$ is a factor describing the periodicity in frequency space. We fit the real component of this function to the data, i.e.,

        \begin{equation}
            A T sinc\left(\pi \nu T\right) cos\left(2 \pi \nu t_0\right)    
        \end{equation}

        in frequency space, then perform an inverse real FFT to retrieve the response in the time domain. This is perhaps the simplest response function we could use, but it still provides a rough measure for comparison with the time lags reported by Fausnaugh et al. in Storm III.


        

        \begin{figure}
            \centering
            \includegraphics[width=0.5\linewidth]{tophat_timedomain.pdf}[A]
            \includegraphics[width=0.5\linewidth]{tophat_freqdomain.pdf}[B]
            \caption{[A] Tophat transfer functions in the time domain show an average time lag of the reverberating curve and a constant distribution in time over an interval. [B] The frequency-dependent time lags associated with each tophat function. Important features related to the average time lag are observable (maximum, corresponding to the average time delay; value of $\nu$ at steepest change), and complicated relationships with higher frequency waves can be noted.}
            \label{fig:tophat_theory}
        \end{figure}





%─────────────
\section{Results}
\label{sec:results}
    Figure~\ref{fig:PSD_ref} shows the power spectral density computed for the reference curve, illustrating its strong variability. The power spectral densities for all wavelengths were computed as a cohort and compiled into a grid, presented in figure~\ref{fig:grids}. The time lags were also computed for all wavelengths, but observation of the output grid indicates significant computational issues. The errors reported in these grids were sampled from the covariance matrix, which does not account for all correlations between frequency bins. These can therefore only be considered lower-limits estimates of the uncertainty.

    Of the 19 lightcurves for which a time lag analysis was attempted, 10 produced possibly meaningful results through this method. A more rigorously computation was attempted for each of these, by sampling the error from the maximum likelihood function, which does account for error correlations. It was also found that a good fit for 3465\AA\ could be obtained by eliminating the highest-frequency bin from the model. The 6 models that ran successfully are compiled in figures \ref{fig:results1} and \ref{fig:results2}. The $psdlag$-computed time delays are reported as black error-barred dots, and the empirical fits are drawn in blue. The average time lag reported in STORM III is drawn as a black vertical line in the time domain plot.

    \begin{figure}
        \centering
        \includegraphics[width=4in]{psd1367A.pdf}
        \caption{The power spectral density for the chosen reference lightcurve. It is strongly variable across all computed Fourier frequencies.}
        \label{fig:PSD_ref}
    \end{figure}

    \begin{figure}
        \includegraphics[width=.46\textwidth]{psd_atlas.pdf}[A]
        \includegraphics[width=.46\textwidth]{timelag_atlas.pdf}[B]
        \caption{[A] Power spectral densities for all observed light curves, excluding the reference curve. A decrease in variability is observed with increasing wavelength. [B] Time lags as a function of Fourier frequency. Computational issues make many of these results dubious, but ten curves provide at least a lower-limit estimate of the uncertainty. It can nevertheless be noted that time delay increases with wavelength.}
        \label{fig:grids}
    \end{figure}


    \begin{figure}
        \centering
        \includegraphics[width=.45\textwidth]{freq_3465A.png}
        \includegraphics[width=.45\textwidth]{time_3465A.png}[3465\AA]
        \centering
        \includegraphics[width=.45\textwidth]{freq_3471A.png}
        \includegraphics[width=.45\textwidth]{time_3471A.png}[3471\AA]
        \centering
        \includegraphics[width=.45\textwidth]{freq_6175A.png}
        \includegraphics[width=.45\textwidth]{time_6175A.png}[6175\AA]
        \caption{The shorter wavelengths, corresponding to emission from the inner region of the accretion disk, do not show strong agreement with the average time lags reported by Fausnaugh et al. This may be due to difficulty fitting the function because the computer is unable to distinguish between emission from the source and reprocessing.}
        \label{fig:results1}
    \end{figure}
    
    \begin{figure}
        \centering
        \includegraphics[width=.45\textwidth]{freq_6439A.png}
        \includegraphics[width=.45\textwidth]{time_6439A.png}[6439\AA]
        \centering
        \includegraphics[width=.45\textwidth]{freq_7657A.png}
        \includegraphics[width=.45\textwidth]{time_7657A.png}[7657\AA]
        \centering
        \includegraphics[width=.45\textwidth]{freq_9157A.png}
        \includegraphics[width=.45\textwidth]{time_9157A.png}[9157\AA]
        \caption{The empirical fits in these wavelengths show a strong agreement with the average time lags reported by Fausnaugh et al. The strong response trending downward from time 0 may be from the source emission, while the bump is from thermal emission further out in the accretion disk. These components appear to overlap in 6439\AA.}
        \label{fig:results2}
    \end{figure}

\section{Conclusion}
    \label{sec:conclusion}

    Response curves in the time domain were successfully computed for 6 of the 19 lightcurves. All indicate some variable time-dependent response. The functions fit to 3465\AA, 3471\AA, and 6175\AA\ do not show strong agreement with the average time delays reported by Fausnaugh et al. 7657\AA\ and 9157\AA\ show a response bump precisely where Fausnaugh et al. predicted. The response bump in 6439\AA\ overlaps with the driving emission, but still appears to be in agreement with the time lag computed in STORM III. The model used to compute these functions is rudimentary, and a more rigorous analysis may provide better fits and more insight.

    The thermal reprocessing hypothesis predicted that the variability would decrease as the wavelength increases, and this is certainly observed in the PSD grid in figure \ref{fig:grids}. It was also predicted that the time lag would increase with wavelength due to the decreasing temperature along the accretion disk's increasing radius. This feature is observed in a general sense along the grid of time lags in figure \ref{fig:grids}, and also between the fitted curves for 6439\AA\ and 7657\AA/9157\AA, but not in the other fitted curves. This may be due to inaccuracies from the rudimentary fitting function. This information does seem protective of the accretion disk hypothesis. One thing is certain: further study is warranted.

\printbibliography


\end{document}