\section{Introduction}

The $\sim$2.4 sr field of view (FoV) and
broad energy coverage from 20 MeV to $>$300 GeV of the Large Area Telescope (LAT, \citealt{2009ApJ...697.1071A}) on the {\it Fermi Gamma-ray Space Telescope} mission make it a powerful instrument to monitor the sky for high-energy transients.  As the LAT surveys the entire sky every 3 hours, it is sensitive to transient emission from a variety of sources, including stellar-mass compact objects (Neutron Stars -- NS, and Black Holes -- BH) over timescales from milliseconds to years, including those predicted to be associated with gravitational waves (GWs). Current GW detectors are sensitive to signals from the merging of compact objects. Some of these mergers, like the putative progenitors of short Gamma-Ray Bursts (sGRBs, \citealt{1986ApJ...308L..43P,1989Natur.340..126E,1992ApJ...395L..83N,2011ApJ...732L...6R,2012ApJ...746...48M}), emit both a short-lived $\gamma$-ray signal ($\lesssim 2$ s) immediately after the merger (``prompt'' emission), and a long-lived and broadband ``afterglow'' signal lasting minutes to hours. If the GRB happens to be in the FoV at the time of the trigger, the LAT can detect the short-lived prompt emission phase.
If the GRB is outside the FoV, because of its survey capability the LAT can still detect the GRB by measuring its temporally-extended afterglow emission as soon as it enters the FoV. The LAT has detected high-energy $\gamma$-ray emission from $>$130 GRBs to-date \citep{2016AAS...22741601V}, including $\sim$10 sGRBs. Given the uncertainty in theoretical predictions for counterparts to GW sources and the demonstrated emission of high-energy $\gamma$-rays from systems containing compact objects, searching the LAT data for $\gamma$-ray counterparts to GW events.

The era of gravitational wave astronomy has begun with the first science run (`O1') of the recently upgraded
LIGO \citep{1992Sci...256..325A,2009PhRvD..80j2001A} from 2015 September to 2016 January. The Virgo Observatory \citep{1999APh....10..369C, 2009CQGra..26h5009A} will soon be added to the network for the second science run in late 2016. The GW frequency range that LIGO and Virgo are sensitive to is expected to be dominated by mergers of compact stellar-mass objects that are most likely remnants of
stellar evolution: two neutron stars (NS-NS), two black holes (BH-BH), or a NS and a BH. The sensitivity and horizon distance of the GW network to these mergers scales with the masses of the systems; therefore the accessible volume of the Universe for NS mergers is significantly smaller than that of BH mergers.
When those mergers include at least one NS, an electromagnetic counterpart is predicted to accompany the merger signal in the form of a sGRB. The electromagnetic outcome of a BH-BH merger is less well understood. Finding the counterpart of a GW event is important for understanding the nature of the source. It also has an additional yet less-evident benefit: it improves significantly the accuracy with which all parameters (distance, mass, spin, inclination, etc.) can be estimated. This is obtained by better constraining the localization of the event, which is normally poorly estimated using only the GW signal \citep{2009LRR....12....2S}. \emph{Fermi}\xspace-LAT can localize a transient source with sub-degree accuracy, a very big improvement with respect to a typical localization region from a GW detector, which will typically cover hundreds of square degrees.

On 2015 September 14 at 09:50:45 UTC the LIGO Hanford and Livingston installations  detected a coincident signal, within 10 ms, from \GWNAME, a high-significance trigger in the engineering run just prior to the start of O1. The trigger was determined to be consistent with a waveform predicted by General Relativity from the inspiral and merger of a stellar-mass binary BH system, with constituent BHs of masses around 29 M$_\odot$ and 36 M$_\odot$ \citep{LIGO_GW150914}. The GW luminosity expected theoretically for a massive BH-BH merger leads to an estimate of around 400 Mpc (i.e., z$\sim 0.09$, \citealt{2041-8205-818-2-L22}) for the distance of the source. This observation provides evidence for the existence of isolated and binary stellar-mass black holes, and the first observation of such a system merging. Two days later, on 2015 September 16, LIGO notified the electromagnetic (EM) partner observatories operating within a Memorandum of Understanding (MOU). The EM partner observatories executed follow-up programs \citep{LVEM} to search for a counterpart within the 750 deg$^2$ localization region ($\sim 90\%$ confidence), which was later refined to 601 deg$^2$ \citep{LVEM}. \emph{Fermi}\xspace was operating in normal survey mode at the time of the trigger. Hence, the LAT autonomously observed the entire LIGO localization region within $\sim 70$ minutes of the GW trigger, independently of any notification from LIGO, in the high-energy $\gamma$-ray band.
The LAT Collaboration reported a preliminary search throughout the LIGO localization area that did not reveal any new $\gamma$-ray sources \citep{2015GCN..18709...1O}. The results of a search of the data of the other instrument on board \emph{Fermi}\xspace, the Gamma-Ray Burst Monitor \citep[GBM,][]{GBMinstrument} and the evidence for a weak counterpart are discussed separately in \cite{GBMLIGO}.

In this paper, we describe LAT observations of the localization area of \GWNAME around the time of the trigger, including a dedicated search for an EM $\gamma$-ray counterpart.  No candidate counterparts were found. We describe the details of the data analysis in  \S\ref{sec:data_analysis}, discuss the implications of the counterpart search and prospects for future GW triggers in \S\ref{sec:discussion}, and conclude in \S\ref{sec:conc}.

\section{Data analysis} \label{sec:data_analysis}

We performed two complementary sets of searches for transient high-energy $\gamma$-ray emission:  automated searches (\S\ref{sec:auto_search}) that are performed routinely on all LAT data, and targeted searches in the LIGO localization region (\S\ref{sec:ligo_cont}) on short and long time baselines that exploit the full sensitivity of the standard LAT analysis chain. In Appendix~\ref{apex:likelihood}, and \ref{sec_tsmaps} we provide more details on the various analysis steps.

\subsection{Automated Searches}\label{sec:auto_search}

Since the launch of \emph{Fermi}\xspace in 2008, automated on-board and on-ground analyses of GBM and LAT data have been in place to search for new transients at various time- and energy-scales. As our understanding of the instruments, data, and the variable and transient $\gamma$-ray sky has improved, so have our automated analyses.
Three LAT pipelines were relevant to the counterpart search for \GWNAME:
\begin{itemize}
\item The Burst Advocate (BA) Tool and LAT Transient Factory (LTF) are automated pipelines that search for excess emission in the LAT data at the positions of triggers from GBM, \emph{Swift}\xspace, INTEGRAL and MAXI at the time of the trigger and intervals in the hours afterwards.  As there were no on-board triggers by any of these instruments coincident with \GWNAME, the BA Tool and LTF were not initiated.  However, in the event of a LAT on-board trigger or a trigger from these observatories coincident with a GW trigger, the pipelines would perform an automated search once the LAT data were available on the ground ($\sim 8$ hours after trigger), with results monitored by on-shift personnel.
\item Automated Science Processing (ASP; \citealt{ASP}) is the standard LAT search for transient or flaring sources on 6-hour, 24-hour, and 7-day timescales. The ASP pipeline performs a detection step via a blind search for sources on all-sky counts maps constructed from the event data acquired at each timescale; then, a standard likelihood analysis is run on those data using a source model that includes the candidate sources found in the detection step as well as the already known sources that have been designated for regular monitoring.
LAT Flare Advocates (a.k.a. Gamma-ray Sky Watchers) offer a prompt human verification service to the automatic pipelines and review the results daily, providing an internal report to the LAT Collaboration.
Relevant information on monitored, flaring and transient sources is released to the astrophysical community using the LAT multiwavelength mailing-list\footnote{To sign up for the LAT multiwavelength list visit \url{http://fermi.gsfc.nasa.gov/ssc/library/newsletter/}}, Astronomer's Telegrams\footnote{\url{https://www-glast.stanford.edu/cgi-bin/pub_rapid}} and Gamma-ray Coordinates Network notices\footnote{\url{http://gcn.gsfc.nasa.gov/fermi_lat_mon_trans.html}}. Weekly summary digests are made available through the \emph{Fermi}\xspace Sky Blog\footnote{\url{http://fermisky.blogspot.com}}.
The LAT Flare Advocate service has been very effective in identifying potential candidates for quick follow-up and coordinated observations at other wavelengths \citep{2013arXiv1303.4054C}.
ASP discovers an average of 8 previously unknown $\gamma$-ray transients per year, and has also detected bright GRB afterglows \citep[e.g., GRB 130427A,][]{2014Sci...343...42A}.
\item \emph{Fermi}\xspace All-sky Variability Analysis (FAVA) is a photometric technique that searches for new transients and variable sources on 7-day timescales \citep{2013ApJ...771...57A}. This method compares the number of detected $\gamma$ rays with the average number of expected  $\gamma$ rays based on the observed long-term average in a given region of the sky. In this way, FAVA provides a computationally inexpensive blind search of flux variations over the entire sky that is independent of both an assumed spectral shape of the flaring source and any model for the diffuse $\gamma$-ray background. The FAVA pipeline detects an average of 16 flares per week; about 10\% of these are not associated with $\gamma$-ray catalog sources \citep[e.g.,][]{2014ATel.6098....1K}. Seven-day FAVA lightcurves for any position on the sky are publicly available at the FAVA Data Portal\footnote{\url{http://fermi.gsfc.nasa.gov/ssc/data/access/lat/FAVA/}} hosted at NASA's \emph{Fermi}\xspace Science Support Center (FSSC).

\end{itemize}
During the 6-hour interval\footnote{The ASP 6-hour interval containing the LIGO trigger time includes LAT data between 2015 September 14, 06:11:33$-$12:00:00 UTC.} containing the LIGO trigger \GWNAME, ASP detected ($>$3$\sigma$) twelve known $\gamma$-ray sources and three low-significance ($>$1$\sigma$) unidentified transients, none consistent with the LIGO event localization. None of the LAT pipelines found a possible counterpart to \GWNAME.

\clearpage

\subsection{Search in the LIGO localization region}\label{sec:ligo_cont}
The LIGO Scientific Collaboration reported results from a Bayesian parameter estimation analysis of \GWNAME under the assumption that the signal arises from a compact binary coalescence (CBC) using the latest offline calibration of the GW strain data.
The most accurate localization map for this event (\texttt{LALInference}) is
based on a Bayesian Markov-Chain Monte Carlo and nested sampling to forward model the full GW signal including spin precession and regression of systematic calibration errors.
The localization probability is primarily in the southern portion of the annulus determined by the arrival time difference between LIGO Hanford and LIGO Livingston of $\sim 7$ ms.

Given the uncertainty on EM signals from the merging of two BHs,
we searched different time windows by carrying out two customized analyses of the LAT data. Both analyses are based on the standard maximum likelihood analysis technique used for LAT data, and summarized in Appendix~\ref{apex:likelihood}. In all our searches we included in the likelihood model all sources (point-like and extended) from the LAT source catalog \citep[3FGL,][]{2015ApJS..218...23A}, as well as the Galactic and isotropic diffuse templates provided by the \emph{Fermi}\xspace-LAT collaboration\footnote{\url{http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html}}. We used the Pass~8 \texttt{P8\_TRANSIENTR010E\_V6} event class and the corresponding instrument response functions. These searches are described in the following sections.\\

\subsubsection{Short-baseline search}\label{sec:short_search}

This search focuses on the hours immediately after the GW trigger $t_{GW}$ (2015 September 14, 09:50:45 UTC). The LAT can detect long and short GRB afterglows up to thousands of seconds after the end of the prompt emission \citep{2010ApJ...709L.146D,2013ApJS..209...11A, 2015arXiv150203122V}. Thus a search in this time window is the most likely to find a counterpart to \GWNAME if it is similar to a sGRB.
\emph{Fermi}\xspace was in normal survey mode operations around $t_{GW}$, rocked 50$^\circ$ North from the orbital plane. We consider a point in the sky observable by LAT if it is within the 65$^\circ$ radius FoV and has an angle with respect to the local zenith smaller than 100$^\circ$. The latter requirement is used to exclude contamination from $\gamma$-ray emission from the Earth's atmosphere. The {\it coverage} is the integral of the probability densities of all points in the LIGO localization probability map observable by LAT at a given time, and it is shown in Fig.~\ref{fig:coverage} as a function of time.
\begin{figure*}[t]
\begin{center}
  \includegraphics[width=\textwidth]{coverage.png}
\caption{\emph{Fermi}\xspace-LAT coverage (see text) of the LIGO localization map, as a function of time, before $t_{gw}$ (dashed curve) and after $t_{gw}$ (solid curve).
The green dashed vertical lines denote the boundary of the interval $T_{1}$ used in the analysis (see text for details), while the horizontal dotted lines mark respectively 0\% and 100\% coverage.}
\label{fig:coverage}
\end{center}
\end{figure*}
While the coverage was between 50\% and 90\% in the hours before the trigger (dashed line), at $t_{GW}$ the LAT was unfortunately viewing  the opposite side of the sky from the LIGO localization region. The coverage was zero until $\sim t_{GW} + 4200$. The time interval $t_{GW}+4442$--$4867$\,s ($T_1$) had coverage $>$90\%, while during $t_{GW}+4867$--$10000$\,s coverage varied between 50 and 98\%, decreasing back to zero at around $t_{GW} + 10$ ks. We searched for a transient source in the time interval having more than 90\% coverage ($T_1$), and we did not find any significant excess.

We then derived upper limits for the $\gamma$-ray flux of \GWNAME.
Since the sensitivity of the LAT depends strongly on the angle from the source to the boresight of the instrument, the continuous variation of the LAT viewing direction in survey mode makes any flux limit for a particular source position time-dependent.  Flux limits are also sensitive to astrophysical backgrounds, particularly in the Galactic plane, so that positions along the LIGO arc will have different flux limits, even for the same observing conditions.  These effects mean that flux limits vary according to both time of observation and position in the localization region.
We show a map of the derived upper limits (95\% confidence level) for the $\gamma$-ray flux of GW150914 in the band 100 MeV $-$ 1 GeV in the left-hand panel of Fig.~\ref{fig:ul_map}, and a histogram of the upper limits in the right-hand panel, both for interval $T_{1}$. Assuming a power-law spectrum for the source with a photon
index of $\alpha = -2$, which is typical for GRB afterglows at LAT
energies, the upper limits we find have a median of $1.7 \times
10^{-9}$ erg cm$^{-2}$s$^{-1}$, and 5\% and 95\% percentiles of $0.9 \times
10^{-9}$ and $3.7\times 10^{-9}$ erg cm$^{-2}$s$^{-1}$, respectively.
These upper limits are only weakly dependent on the choice of
$\alpha$ as shown in the right-hand panel of Fig. \ref{fig:ul_map}.
We now describe the upper limits analysis in more detail.\\
\begin{figure*}[tb]
\begin{center}
  \includegraphics[width=0.45\textwidth]{ul_map.png}
  \includegraphics[width=0.45\textwidth]{ul_hist.png}
\caption{Flux upper limits (95\% c.l.) in the energy range 100 MeV--1 GeV for \GWNAME during the interval $T_{1}$ (4442--4867 s from $t_{GW}$). Left panel: upper limits map covering the 90\% region of the LIGO probability map. Right panel: histogram of the upper limits in the map. We assumed a power-law spectrum for the source, with a photon index $\alpha$ of $-2$ (typical of afterglows of GRBs, green histogram), $-1$ (blue histogram) and $-3$ (red histogram). While the distributions are slightly different for the three cases, the ranges spanned by the upper limits are largely independent of the photon index.
\label{fig:ul_map}
\end{center}
\end{figure*}

We considered all $\gamma$-rays with reconstructed energies between 100 MeV and 100 GeV. We then considered the \texttt{LALInference} probability map provided by LIGO, which is a \texttt{HEALPix} map with \texttt{NSIDE = 512} \citep{2005ApJ...622..759G}, corresponding to a typical pixel size of $\sim$$0\fdg11$. The point-spread function (PSF) of the LAT has a 68\% containment radius at 1 GeV of $\sim$1$^\circ$. To save computation time, we resampled the map to \texttt{NSIDE = 256}, which corresponds to a typical pixel size of $\sim$0$\fdg2$. We considered the 11463 pixels in the resampled map contained within the 90\% containment region provided by LIGO. Let us denote $\vec{v}_{i}$ the direction of the center of the i-th pixel. For each i-th pixel we performed an independent likelihood analysis (see Appendix~\ref{sec:auto_search}), considering a Region Of Interest (ROI) of 8$^\circ$ centered in $\vec{v}_{i}$ and testing for a new source at the position $\vec{v}_{i}$. For each likelihood analysis we required the zenith angle of the events to be no more than 100$^\circ$. We did not find any new source with a test statistic (TS) above our adopted threshold of 25, corresponding to $\sim 5 \sigma$ (pre-trials). We then computed the 95\% confidence level upper limit for the flux of a source at each $\vec{v}_{i}$. In order to obtain upper limits reasonably independent of the photon index $\alpha$ adopted in the analysis we need to choose an energy range small enough. We have chosen the range 100 MeV--1 GeV for the upper limits, which provides the largest photon statistic and the maximum sensitivity for sources similar to GRBs. The right-hand panel of Figure~\ref{fig:ul_map} shows that our measurement is indeed largely independent of the choice of $\alpha$.\\

\cite{GBMLIGO} reported the weak transient $\gamma$-ray source \GWNAME-GBM lasting $\sim 1$ s, 0.4 s after the LIGO trigger on \GWNAME. \GWNAME-GBM is consistent with being due to a low-fluence sGRB at an unfavorable viewing geometry to the GBM detectors, although this is not expected from a BH-BH merger.
Assuming the two signals have a common origin, the
combined LIGO and GBM observations reduce the 90\% confidence region from 601 deg$^2$ to 199 deg$^2$. Within the combined LIGO/GBM localization, and in the time interval $T_{1}$, the most significant excess in LAT data has TS = 18. We estimate for this excess a p-value of $\sim 0.05$, taking into account the number of trials of our analysis, which corresponds to a significance that is well below our threshold of $5 \sigma$. The excess has a spectrum well modeled with a power law with a soft photon index $\alpha = -3.2 \pm 0.8$, and it is located close the limb of the Earth (which has indeed a soft spectrum). Therefore we consider this excess in the LAT data very likely to be either a statistical fluctuation or due to Earth limb contamination, and therefore unrelated to \GWNAME or \GWNAME-GBM.

\subsubsection{Long-baseline search}\label{sec:long_search}

In this second search we considered data gathered during a 2-month interval centered on $t_{GW}$. In order to increase the number of $\gamma$-rays we included all photons with energies between 60 MeV and 100 GeV. Since the PSF at 60 MeV is broad, we applied a zenith cut of 95$^\circ$ to further limit Earth limb contamination.
We looked both for a long-duration signal of the order of one day, as well as for a short-duration signal but not necessarily in strict temporal coincidence with the LIGO trigger.
To this end, we covered the entire 90\% probability region provided by LIGO with a set of nine partially overlapping ROIs, each with a radius of 10$^\circ$. Figure \ref{fig:contours} shows the locations of these ROIs and the confidence contours obtained from the \texttt{LALInference} probability map. They are overlaid on a sky map of the $\gamma$-rays detected by the LAT over the interval $t_{GW}$ to $t_{GW} + 10000$ s. In Table~\ref{table:roi} we provide the location of the center of each ROI, listing all the LAT sources from the 3FGL catalog \citep{2015ApJS..218...23A} within each ROI, along with their associations and their classes, using the same notation as the 3FGL catalog.

\begin{figure*}[h!t!b!]
\begin{center}
  \includegraphics[width=1.0\textwidth,trim={1cm 0 1cm 0},clip]{figure_map_ait_10000.png}

\caption{LAT-detected $\gamma$-rays in celestial coordinates (J2000) in the interval from $t_{GW}$ to $t_{GW}$+$10000$s, and in the energy range 60 MeV--100 GeV. The dots are colored according to energies of the $\gamma$ rays as indicated in the color bar. The highest-energy photon recorded in this time interval has an energy of $\sim 60$ GeV, corresponding to the maximum of the color bar. The 90\% and 50\% contour levels from \texttt{LALInference} are overlaid on the image, in red. The map is in \texttt{Hammer-Aitoff} projection and is centered at R.A., Dec$=90^{\circ}$, $-40^{\circ}$. The Galactic plane is highlighted in green. The nine LAT ROIs where we perform the searches described in \S\ref{sec:long_search} are shown as yellow circles.}

\label{fig:contours}
\end{center}
\end{figure*}

For the first analysis of the second search we divided the data in 10~ks time bins. For each time bin and for each ROI we calculated a TS map (see Appendix~\ref{sec_tsmaps}) and  determined the location of the grid position with the maximum TS (TS$_{max}$). We considered the position of TS$_{max}$ as the location of a possible counterpart, and we ran an unbinned likelihood analysis adding a point source at the position of TS$_{max}$. This gave a value of TS$_{src}$ (which is normally similar to TS$_{max}$).
In these maps derived from low-statistics data single high-energy $\gamma$-rays can cause a high value of TS$_{max}$.
To reduce the number of false positives from random coincidences, we required that the number of photons N$_{\gamma}$ that have a probability larger than 0.9 to be associated with the candidate counterpart to be greater than 2.  No excesses met this requirement. In Appendix~\ref{sec_tsmaps} we present the TS maps for the 9 ROIs for the time interval 0--10~ks since $t_{GW}$. We repeated the same analysis considering time bins of 1 day, and again did not find any significant excess.

We also considered the possibility of excesses over shorter timescales ($<$1 h), similar to the typical duration of high-energy emission from GRBs \citep{2013ApJS..209...11A,2015arXiv150203122V}, but not in temporal coincidence with the GW trigger and hence not covered by the ``short-baseline'' search described in \S\ref{sec:short_search}.
We calculated the entry and exit times for each ROI in the FoV of the LAT (a "FoV passage"), requiring that the distance between the LAT boresight and the center of the ROI be $<$60$^\circ$.
In standard survey mode, the duration of a FoV passage varies from a few hundred seconds to nearly one hour. Since we do not know if an EM signal would be in temporal coincidence with the GW signal, we searched for possible excesses in every passage, corresponding to a total of 6615 passages for each ROI. We did not detect any significant excess in any of the passages before or after $t_{GW}$ for any ROI.

\begin{figure}[!ht]
\begin{center}
  \includegraphics[width=0.45\textwidth]{tsmap_150_-6_1513509.png}
  \includegraphics[width=0.45\textwidth]{tsmap_131_-68_-1785773.png}\\
  \includegraphics[width=0.45\textwidth]{tsmap_63_-66_-2117671.png}
  \includegraphics[width=0.45\textwidth]{tsmap_52_-57_270539.png}
\caption{TS maps of the most significant excesses detected over one-orbit timescales ($\sim$ minutes) within 30 days before or after the LIGO trigger. Yellow stars are sources from the 3FGL catalog, and circles are individual LAT $\gamma$ rays with their size proportional to the reconstructed energy. The yellow dots trace the LIGO 50\% and 90\% confidence contour.}
\label{fig:tsmap}
\end{center}
\end{figure}

Figure~\ref{fig:tsmap} displays the four TS maps with the highest TS$_{src}$ and N$_{\gamma}>2$.
The first map (TS$_{src}=21$) is from the first ROI and corresponds to the passage from 22:15:53 to 22:20:30 UTC on 2015 October 10 (17.5 days after the trigger time $t_{GW}$).
The second map is from the third ROI in the passage from 17:47:41 to 17:56:41 UTC on 2015 August 21 (20.7 days before $t_{GW}$), and has also a TS$_{src}=21$.
The third map (TS$_{src}=22$) is from the seventh ROI in the interval from 21:35:48 to  21:50:01 UTC on 2015 August 20 (24.5 days before $t_{GW}$).
Finally, the last map is from the eighth ROI, during 12:59:40 to 13:39:41 UTC on 2015 September 17 (3.1 days after the trigger LIGO trigger time $t_{GW}$) has the highest value of TS with TS$_{src}=23$.
The peak TS values correspond to low-energy $\gamma$ rays in random coincidence with high-energy $\gamma$ rays, and the highest values obtained in this analysis should not be considered indicative of a possible EM counterpart.
Moreover, the above search involves a large number of trials. To estimate their impact on the peak TS values we performed a Monte Carlo study described in the next section.

\subsection{Comparison with Monte Carlo simulations}

To validate our interpretations of TS values we performed a detailed Monte Carlo simulation of 2 months of data (the same interval used in our analysis).
The actual pointing history of the satellite was used; therefore the correct exposure of the sky was automatically taken into account.
All the sources from the 3FGL catalog were kept fixed at their 3FGL catalog fluxes.
As a result, the simulation is suitable for computing the distribution of TS under the null hypothesis that no transient signal is present.
With the simulated data we repeated exactly the same analysis used on real data and described in the previous section.

In Figure \ref{fig:tsdistribution} we compare the distribution of TS$_{src}$ obtained from flight data (filled lines) and Monte Carlo (dashed lines).
We note that the Monte Carlo distributions are a good match to the distributions of the TS$_{src}$ values obtained from the flight data, and the good absolute agreement is consistent with no statistically-significant transient counterpart being present in the flight data.
Also, given the large number of trials, relatively high values of TS can be obtained in Monte Carlo simulations even if no transient signal was added.
In the flight data we found 4 cases with TS$>$20 and with N$_\gamma > 2$, and this must be compared with the 9 cases we obtain when we analyzed the simulated data. In other words, we expected 9 false positives with TS$_{src}>20$ (and N$_{\gamma}>2$) in 2 months of flight data, and we observe 4.

\begin{figure}[ht]
\begin{center}
  \includegraphics[width=0.7\textwidth]{ts_distribution.png}
\caption{Data-Monte Carlo comparison of the distributions of the TS$_{src}$ value. Solid lines correspond to TS distributions obtained from data, while dashed lines correspond to distributions obtained from Monte Carlo simulations. Blue lines are the distributions of the TS values when no additional selections are applied, while red corresponds to the final distribution when we require that at least three $\gamma$-rays have at least a 90\% probability to be associated to the additional point source.}
\label{fig:tsdistribution}
\end{center}
\end{figure}

\section{Discussion} \label{sec:discussion}

The most promising astrophysical GW sources in the frequency range of LIGO/Virgo are the mergers of compact object binaries with NS and/or stellar-mass BH constituents.  The detection of \GWNAME is undoubtedly a major breakthrough in astronomy, being the first detection of GWs and the first detection of a merging binary BH system.

Maximizing the science return from GW detections requires the identification and study of coincident EM counterparts, which
would help resolve degeneracies associated with the inferred binary parameters. For example, a strong $\gamma$-ray signal may only be seen at small binary inclination relative to the sky, such that the jet (along the direction of the total angular momentum) is pointed toward us. A high-significance counterpart association would also decrease the significance threshold
necessary for a confident GW detection, thereby effectively increasing the distance to which the GW signal can be detected and the searchable volume as the cube of the distance.  It would help identify the host galaxy and thereby constrain or measure the merger redshift,
setting the luminosity scale and allowing an independent measurement of the Hubble constant or other cosmological parameters \citep{2012ApJ...746...48M,2012ApJ...748..136C}.  The complementary information encoded in the EM signal (spectral and temporal evolution, energetics, inferred environment) is likewise essential to unravel the astrophysical context of the coalescence event.

The most commonly hypothesized EM counterpart of an NS-NS/NS-BH merger is a sGRB, powered by accretion onto one of the two central compact objects
\citep{1986ApJ...308L..43P,1989Natur.340..126E,1992ApJ...395L..83N,2011ApJ...732L...6R,2012ApJ...746...48M},
which launches relativistic jets that produce a short ($<$2 s) bright flash of keV-MeV-peak $\gamma$ rays followed by a broadband longer-lasting afterglow. The LAT detects approximately two sGRBs per year, with $>$100 MeV afterglows lasting up to hundreds of seconds after the trigger \citep{2010ApJ...716.1178A,2013ApJS..209...11A}, seeded by sGRB triggers from the GBM and {\it Swift}-BAT.  The most prominent example of a LAT-detected sGRB is GRB 090510, which simultaneously triggered \emph{Swift}\xspace and GBM, as well as causing a LAT on-board transient source trigger, resulting in a LAT localization being circulated within seconds. The LAT detected both a spike during the prompt $\gamma$-ray emission and an extended afterglow lasting hundreds of seconds, consistent with the observed optical and X-ray emission \citep{2010ApJ...709L.146D}.
Figure \ref{fig:grb090510} shows the $>$100 MeV $\gamma$-ray lightcurve of GRB 090510 scaled to z=0.09, the nominal redshift of \GWNAME inferred from the GW observations, in comparison to the LAT upper limits described in \S\ref{sec:short_search}.  If \GWNAME had had a high-energy $\gamma$-ray lightcurve similar to GRB 090510 and had been more favorably placed relative to the LAT boresight at the trigger time, it would have been easily detectable by the LAT during observations similar to those described in this paper. A GRB at the redshift of the LIGO GW source entering the field of view of \emph{Fermi}\xspace-LAT within the first 100 s would have been detected if it were more than an order of magnitude fainter than GRB 090510.

LAT-detected long GRBs have also shown similar behavior, with a prompt spike contemporaneous with sub-MeV emission followed by long-lived emission lasting from minutes to hours with similar time dependence as radio-to-X-ray afterglows \citep{2013ApJ...763...71A,2014Sci...343...42A}.
Although the detection rate of LAT sGRBs is low, the low redshifts of potentially detectable GW sources and therefore their potentially bright EM counterpart emission, the uncertainty in the observational signatures of GW sources, and the continuous observations of the entire sky during survey operations, make searching the LAT data for counterpart sources worthwhile.

\begin{figure}[t]
\centering
\includegraphics[width=0.7\textwidth]{GRB090510_lc_scaled.png}
\caption{
GRB~090510 is the only LAT-detected sGRB with a measured redshift \citep[z = 0.903,][]{2010ApJ...716.1178A}. We compare it here to LAT observations of \GWNAME.  The 100 MeV--10 GeV lightcurve of GRB~090510 has been scaled to z=0.09, the redshift corresponding to the best-fit distance from the GW observation \citep{2041-8205-818-2-L22}. The red arrow indicates the 95\% confidence upper limit measured during $T_1$ across the LIGO localization region and is similar to the upper limits from the measured lightcurve. Even though GRB~090510 had exceptionally bright prompt and afterglow emission, this figure demonstrates that the LAT would detect transients more than an order of magnitude fainter than GRB 090510, provided a more favorable placement of the GW candidate so that it enters the LAT FoV within 100~s of the trigger time.
}
\label{fig:grb090510}
\end{figure}

\subsection{LAT Searches for sGRBs}

Two strategies are useful for associating EM detections of transients with GW sources. The first is to search GW data for counterparts to EM events.  This was done archivally using GRB triggers for the years 2006--2011 \citep{2014PhRvD..89l2004A} and initial LIGO data runs, yielding no credible candidate sources (with a $<20$ Mpc horizon distance). This approach is not very efficient since most of the EM events are outside the current volume sampled by the initial LIGO runs, which increases the trials factors and diminishes the GW sensitivity.

The second approach is to search for EM bursts (e.g., sGRB) related to GW events.  The most promising way to find a counterpart is via a prompt trigger from the burst itself using wide FoV instruments like GBM or \emph{Swift}\xspace-BAT.  Given its large FoV and sensitivity to energies from 8 keV to 40 MeV, GBM is the most prolific detector of sGRBs currently operating, with $\sim$40--80 detections per year, and therefore it is very well suited for EM counterpart searches \citep{GBMLIGO}. Yet the LAT offers several capabilities that make it a unique resource.

In the normal mode of \emph{Fermi}\xspace GRB detection and measurement, GBM issues triggers on-board for both short and long GRBs (as well as solar flares, terrestrial $\gamma$-ray flashes, magnetar flares, and other short, hard X-ray/soft $\gamma$-ray transients), and approximately half of the GBM GRBs occur in the LAT FoV.  However, in cases of high peak brightness, GBM initiates an automated repoint of the \emph{Fermi}\xspace spacecraft to optimize LAT observations over the next 2.5 hours.
Regardless of whether an automated pointing is triggered, even during normal sky survey operations, the LAT will observe the entire sky within 2 orbits ($\sim 3$ hours), fully covering 100\% of any GW localization contour, making follow-up observations of GW triggers both automatic and routine.  Neither instrument on \emph{Fermi}\xspace has observing constraints limiting the viewing near the Sun or Moon.

If the LAT were to detect a counterpart to a GW trigger, it has the distinct advantage of providing smaller localization uncertainties than GBM.  Typical GBM localizations have radii of several degrees plus systematic uncertainties of 3--14$^\circ$ \citep{2015ApJS..216...32C}, while LAT localizations have radii $\sim 0.1-1^\circ$.  A LAT detection would substantially reduce the sky area for follow-up observations at other wavelengths, some of which require extensive tiling campaigns over days after a trigger to cover a significant fraction of the GW localization region \citep{2016MNRAS.455.1522E,LVEM}.  Although the LAT obtained an on-board localization for short GRB 090510 \citep{2010ApJ...716.1178A} within 15 seconds, LAT localizations typically have an $\sim$8-hr latency for data transmission and ground processing.

The LAT has also recently (as of 2015 June) benefited from a major upgrade to the event reconstruction analysis pipeline, known as Pass 8 \citep{Pass8}.  Pass 8 improves the LAT sensitivity due to increased photon acceptance, especially at low ($<$100 MeV) and high energies ($>$10 GeV), and reduces localization radii by $\sim 30\%$.  Combined with new analysis pipelines that search for transients on all timescales \citep{2015arXiv150203122V}, the LAT is now better suited to discover counterparts to GW sources in both automated pipelines and the specialized searches described in \S\ref{sec:data_analysis}.

The LIGO EM follow-up partners conducted large follow-up campaigns in the optical, radio, and X-ray to search for a counterpart to \GWNAME \citep{LVEM}.  Similar campaigns will occur for other GW triggers no matter the likely progenitor and regardless of whether it is observed along the putative jet axis.
Although no orphan GRB afterglows (which are expected to be associated with off-axis observing angles) have conclusively been detected, models \citep{2002ApJ...570L..61G,2010ApJ...722..235V} predict delayed, fainter transient emission from the afterglow itself, an isotropic optical-near-infrared kilonova signature \citep{1998ApJ...507L..59L}.  Other potential EM counterparts have been theorized, including late radio emission from mildly relativistic material that is dynamically ejected during the merger and drives a shock into the external medium \citep{2011Natur.478...82N}, and radio-to-X-ray (or $\gamma$-ray) emission over seconds to days from relativistic ejecta created as the NSs collide during their merger \citep{2014MNRAS.437L...6K}.
Given the capability of the LAT to detect transients on timescales from milliseconds to years, careful searches for possible afterglow emission within the LIGO localization area over relevant timescales after the GW trigger are essential and require no change to the \emph{Fermi}\xspace observing strategy.

High-energy $\gamma$-ray emission from GRBs observed by the LAT above 100 MeV lasts much longer (minutes to hours) than the prompt emission observed by the GBM \citep{2013ApJS..209...11A}, for both long and short GRBs \citep{2010ApJ...712..558A,2010ApJ...709L.146D}.  Therefore, even if a localization probability region from LIGO/Virgo is outside the FoV of the LAT during the time of the trigger, the instrument could detect temporally-extended emission in the minutes following the prompt signal, when the region re-enters its FoV. This could result in a much more precise localization of the EM counterpart.

\subsection{EM Counterparts to BH Binary Mergers}

As discussed in \cite{LIGO_GW150914}, the \GWNAME waveform is consistent with the expectation for the merger of two stellar-mass BHs.  The comparatively clean waveform in the ``chirp'' phase before the merger and the ``ring-down'' stage after the merger would not naturally be expected from coalescing NSs, which possess matter distributed outside an event horizon. The expected progenitors of sGRBs are NS-NS or NS-BH binaries \citep{2013ApJ...762L..18G,1986ApJ...308L..43P,1989Natur.340..126E,1992ApJ...395L..83N,2011ApJ...732L...6R,2012ApJ...746...48M}.  Therefore, a classical sGRB counterpart to \GWNAME is not expected.

Prior to the discovery of \GWNAME, studies focusing on EM counterparts to stellar-mass BH mergers were few in number.  Most of the numerical simulations have focused on supermassive BH mergers, where circumbinary disk formation is expected with ample gas supply available to power an EM counterpart \citep{2007Sci...316.1874M}.  Stellar-mass BH mergers should require a substantial quantity of nearby gas to form the disk-jet system that is expected to be necessary for an EM counterpart to be detectable at tens-to-hundreds of Mpc distances.  The weak counterpart candidate detected by the GBM described by \cite{GBMLIGO} poses an interesting dilemma for theoretical models, if it is connected to \GWNAME.  Unfortunately, since the GW localization region was not in the FoV of the LAT at the time of the GW trigger, the prompt signal from the GBM candidate could not be addressed with the LAT observations.  Should such $\gamma$-ray
associations arise for future LIGO events, a new paradigm for these mergers will be indicated.

\section{Conclusions} \label{sec:conc}

The \emph{Fermi}\xspace-LAT is uniquely capable of searching for high-energy $\gamma$-ray counterparts to sources detected by GW observatories.  We use this capability to undertake a detailed search in the regular LAT survey data for a counterpart temporally and spatially coincident with the LIGO trigger on \GWNAME. Although \GWNAME was not in the LAT FoV at the trigger time, the LAT observed the entire region within $\sim 70$ minutes of the GW trigger.
We searched on short and long timescales for evidence of a transient $\gamma$-ray source contemporaneous with \GWNAME.  No LAT counterpart is detected, and upper limits have been set on the GeV $\gamma$-ray flux within the LIGO localization.

As the sensitivity of LIGO and Virgo improve over the next few years, their detection horizon for NS-NS and NS-BH binary mergers, and thus the likelihood of sGRB coincidence, will increase greatly.  The merger rates are highly uncertain and depend on the populations and evolution of binaries and opening angles of sGRB jets.  The coincident detection of prompt $\gamma$-ray signals, with on-axis or potentially off-axis afterglow signals will teach us about the physics of binary mergers.
With the discovery of \GWNAME, the search for EM-GW coincidences enters a new phase. It is important to test the strong expectation that stellar BH mergers do not radiate much light; it is equally important to refine the techniques that will be needed to associate sGRBs with GW events from NS mergers. The approach described in this communication is well suited to achieve these twin goals when LIGO detects more GW events.

\acknowledgments{
  The \emph{Fermi}\xspace LAT Collaboration acknowledges generous ongoing support
from a number of agencies and institutes that have supported both the
development and the operation of the LAT as well as scientific data analysis.
These include the National Aeronautics and Space Administration and the
Department of Energy in the United States, the Commissariat \`a l'Energie Atomique
and the Centre National de la Recherche Scientifique / Institut National de Physique
Nucl\'eaire et de Physique des Particules in France, the Agenzia Spaziale Italiana
and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education,
Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research
Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and
the K.~A.~Wallenberg Foundation, the Swedish Research Council and the
Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'\'Etudes Spatiales in France.
}

\bibliographystyle{yahapj}


\section{Introduction}
Fueled primarily by the promise \cite{Pendry2000} of optical imaging beyond the diffraction limit,  metamaterials, \cite{ref:metamaterials} artificial composites with the unit cell size well below the optical wavelength, now for more than a decade remain one of the primary foci of modern research. However, while  metamaterial devices, from the originally promised superlens \cite{Pendry2000} based on the negative index media, to its younger challenger the hyperlens \cite{OE2006,Engheta2006} that uses hyperbolic metamaterials,  have been demonstrated in experiment \cite{ZhangHyperlens, SmolyaninovHyperlens}  in broad range of frequencies and even extended to other wave phenomena such as e.g. ultrasound imaging, \cite{UltrasoundHyperlens} they have yet to find their way to broad acceptance in real-world applications.

The inherent curvature of the metamaterial device that is necessary for optical magnification (either in the general formfactor of the device \cite{OE2006,Engheta2006,MagnifyingSuperlens} or in its unit cell topology \cite{PlanarHyperlens}), leading to the requirement on precise positioning of the object at subwavelength scale, remains a major practical obstacle to a wider acceptance of metamaterials-based imaging systems. More important, however, is the result of  substantial loss in metamaterial structures, that is generally a direct and inevitable consequence of using metallic (or other negative permittivity) components that remain necessary for either the negative-index metamaterials or hyperbolic media. Faced with  an order of magnitude intensity loss \cite{ZhangHyperlens} as the cost of producing a magnified  optical image of a subwavelength pattern, one must either accept a prohibitive cost in the signal-to-noise ratio, or substantially increase the intensity of the field that illuminates the target. The latter option however does not fit well with the constraints of biomedical imaging, where high optical fields may severely damage biological tissue.

While directly related to the fundamental restrictions on optical metamaterials \cite{Stockman}  via the Kramers-Kronig relations, this problem however allows a natural and straightforward solution: the metamaterial component of a super-resolution imaging setup must be located ``in front'', rather that ``after'', the object -- so that it's probed by the already attenuated field, thus avoiding damage. This is precisely the case in the structured illumination \cite{Gustaffson} approach to optical imaging, where  a grid pattern, usually generated through the interference of incident light, is superimposed on the specimen.\cite{Gustaffson,si} The imaging resolution  is then determined by the periodicity of the illuminating field, and for illumination from free space this limit
 reduces to one quarter of the wavelength, $\Delta = \lambda_0/4$.
 However, as we demonstrate in the present paper, hyperbolic metamaterials that are free from the conventional diffraction limit and can support optical patterns with the periodicity that is orders of magnitude below the free-space wavelength $\lambda_0$, can dramatically improve the imaging resolution in structured illumination. The resulting   {\it hyper-structured illumination} approach to super-resolution imaging is the main objective of this work.

\section{Structured Illumination in Hyperbolic Media}

In a conventional transparent medium with the refractive index $n$, the frequency $\omega$ generally limits the magnitude of the wavenumber $k$ to the value of $n \omega/c$. However, the situation is  dramatically different in the world of hyperbolic metamaterials, \cite{OE2006} where
the opposite signs of the dielectric permittivity components in two orthogonal directions ($\epsilon_\tau \epsilon_n < 0$ -- see the inset in Fig. 1(a))  lead  to a hyperbolic dispersion of $p$-polarized propagating waves
\begin{equation}
\frac{k_\tau^2}{\epsilon_n} + \frac{k_n^2}{\epsilon_\tau}  =  \frac{\omega^2}{c^2},
\label{eq:hyperbola}
\end{equation}
with the wave numbers unlimited by the frequency. For $k \gg \omega/c$, the hyperbolic iso-frequency surface of (\ref{eq:hyperbola}) asymptotically approaches a cone, with
$
k_n \simeq \sqrt{\frac{\epsilon_\tau}{\epsilon_n}} k_\tau.
$

As the direction of the group velocity ${\bf v}_g$ corresponding to a given wavevector ${\bf{k}}$ coincides with the normal to the iso-frequency surface, at large wavenumbers in a uniaxial hyperbolic medium  the angle between the group velocity and the material
symmetry axis $\hat{{\bf z}} \equiv \hat{\bf{n}}$ (see the inset to Fig. 1(a)) approaches the value of $\theta_c = \arctan\sqrt{{\epsilon_\tau}/{\epsilon_n}}$. As a result, for  a point source placed at the edge of such medium, the emission intensity diagram forms a conical pattern (see Fig. 1 (b)-(d), calculated for a practical realization of hyperbolic media using silver-glass layered metamaterial), with the ``thickness'' of the expanding conical ``sheet'' of electromagnetic radiation well below the free-space wavelength. This phenomenon is now well-established, and has recently been successfully used for subwavelength optical lithography. \cite{hyper-lithorgraphy}

As strong material dispersion is inherent to hyperbolic materials, a change in the frequency $\omega$ will result in a substantial variation of the emission cone angle $\theta_c$, as illustrated in Fig. 1(a) for three different wavelengths. While generally considered detrimental to metamaterial applications (since this would immediately push the system away from whatever resonant condition is was designed for), it is precisely this property that is essential to the concept of the hyper-structured illumination. While a single source (or an opening) at the bottom of the hyperbolic ``substrate'' only illuminates a small part of the object plane at its ``top'', a change in the electromagnetic frequency would allow to sweep the entire object plane (see Fig. 1 (b) - (d)). As a result, a complete subwavelength image can be obtained in a  single hyper-spectral measurement.

This behavior is further illustrated in Fig. 2, which shows the object formed by four silicon nanowires on top of the silver-glass metamaterial ($10$ silver and $11$ glass layers, each $5$ nm  thick  -- see Fig. 2 (a)), illuminated at the wavelengths of $425$ nm (panel (b)) and $575$ nm (panel (c)). Note the subwavelength localization of the electromagnentic field near the object, and a  dramatic difference in the intensity distribution for different wavelengths that allows to clearly ``resolve'' distinct elements of the target despite their close spacing.

A finite size of the unit cell of the metamaterfial, $a$, introduces an upper limit on the (Bloch) wavenumber $k_n \equiv k_z$, equal to $k_n^{\rm max} = \pi/a$, as well as quantitative corrections to the dispersion $\omega\left({\bf k}\right)$ and to the geometry of the iso-frequency surface,  leading to its deviations from the hyperbolic shape. However,  the overall topology of the iso-frequency surface, together with  its dependence on the wavelength, are generally preserved -- see Fig. 1 (a).  As a result, one will still observe the general profile of the conical illumination pattern,  as well as its evolution with  the wavelength that are consistent with the predictions of the effective medium theory, as clearly seen in Fig. 1 (b)-(d).  However, as the size of the unit cell limits the maximum wavenumber that can be supported by the propagating waves in the hyperbolic metamaterial, it defines the ultimate resolution limit $\Delta$ of the hyper-structured illumination:
\begin{equation}
\Delta \simeq a
\label{eq:hsi_limit}
\end{equation}
With the scale of the unit cell in high-quality hyperbolic metamaterials that have already been demonstrated, on the order of a few nanometers, \cite{Liu} Eqn. (\ref{eq:hsi_limit}) allows for optical imaging of deeply subwavelength objects, to the point of resolving large individual molecules.

\begin{figure}[htbp]
\centering
\fbox{\includegraphics[width=3 in]{fig1}}
\caption{Panel (a): iso-frequency curves in a silver-glass hyperbolic metamaterial, for three different free-space wavelengths $\lambda_0$ : 425 nm (blue curves), 550 nm (green) and 750 nm (red). Solid curves are obtained from the exact calculation taking into account the finite width and the actual material loss of each layer ($a = 5 \ {\rm nm}$), while the dashed lines correspond to the effective medium approximation (see Eqn. (\ref{eq:hyperbola})). The arrows normal to the iso-frequency curves, indicate the directions of the group velocity. The inset shows the schematics of the silver-glass hyperbolic metamaterial, and introduces the coordinate system.
 Panels (b)-(d) : the magnetic field intensity for light of different wavelengths ($\lambda_0 = 425 \ {\rm nm}$ (b), $ 550 \ {\rm nm}$ (c) and $ 750 \ {\rm nm}$ (d)), incident from $z < 0$ half-space onto the hyperbolic metamaterial formed on top of a thick silver film with a narrow slit centered at $x = 0$. The slit is  parallel to the $y$-axis and is $5 \ {\rm nm}$ wide in the $x$-direction.  Note the conical emission pattern typical for any point source in  a hyperbolic medium, and strong sensitivity of its directionality to the wavelength.
 }
\label{fig:1}
\end{figure}

\begin{figure}[htbp]
\centering
\fbox{\includegraphics[width=3.25 in]{fig2}}
\caption{Hyper-structured illumination, with silver-glass hyperbolic metamaterial substrate.
 Panel  (a) shows the geometry of the imaging system, with dark and light-gray regions corresponding to silver and glass layers (each $5$ nm thick) respectively. Depending on the illumination wavelength, light incident from the slit at the bottom of the metamaterial, illuminates different areas of the target -- see panels (b) and (c), with a sub-diffraction resolution. The imaging target consists of four long silicon nanowires with rectangular cross-section of $20 \ {\rm nm} \times 10$ nm, aligned along the $y$-direction. The false-color representation shows the magnitude of the magnetic field, calculated taking full account of the material absorption in the system.}
\label{fig:2}
\end{figure}

\section{Depth of field in hyper-structured illumination and the Motti projection}

While a hyperbolic metamaterial supports propagating electromagnetic waves with the characteristic length scale up to its unit cell size, these waves will evanescently decay outsize this medium. As a result, the imaging resolution of the proposed hyper-structured illumination setup will rapidly deteriorate with  $\Delta z$ (see Fig. 2(a)), approaching the diffraction-limited value of $\lambda_0 / 4$ for $\Delta z \gg \lambda_0$. For a 3D object structured at a deeply subwavelength scale, the proposed imaging system based on hyper-structures illumination, will resolve the Motti projection \cite{MorseFeshchbach}  $\Delta_M\epsilon\left(x,y\right)$ of its tree-dimensional permittivity distribution $\epsilon\left({\bf r}\right)$:
\begin{eqnarray}
\Delta_M\epsilon\left(x,y\right) & = & \int \frac{d^3{\bf r'} \cdot \left( z' / 2 \pi\right) \cdot \left(\epsilon\left({\bf r'}\right) - 1\right)  }{\left(\left(x - x'\right)^2 + \left(y - y'\right)^2 + z'^2\right)^{3/2}}.
\label{eq:Motti}
\end{eqnarray}
Note that the integral transformation in Eqn. (\ref{eq:Motti}) is independent of the illumination wavelength, and does not introduce an additional length scale: the effective averaging in (\ref{eq:Motti}) is defined by the distance $z'$ to the object plane. As a result, Eqn. (\ref{eq:Motti}) reduces the resulting resolution $\Delta$ to the corresponding depth of fields $\Delta z$.
This is illustrated in Fig. 3, which shows the Motti projection (panel (b)) for several solid spheres near the object plane $\Delta z = 0$.

\begin{figure}[htbp]
\centering
 \fbox{\includegraphics[width=2.5 in]{fig3}}
\caption{A group of five identical solid spheres positioned near the object plane (panel (a)) and their Motti projection (panel (b)), shown in gray scale. Note the reduced
contribution to $\Delta_M\epsilon\left(x,y\right)$ at longer separation from the object plane $\Delta z = 0$.}
\label{fig:3}
\end{figure}

\section{ Theoretical Framework}

When illuminated from the slits at the ``bottom'' of the hyperbolic substrate (see Fig. \ref{fig:2}(a)),
an object with the dielectric permittivity $\epsilon\left({\bf r}\right)$ placed near the image plane $\Delta z = 0$, will
induce scattered light with the far-field amplitudes for the $s$ and $p$ polarizations
\begin{eqnarray}
E_s\left({\bf k}; \omega\right) & =& \frac{4 \pi^2 \omega^2}{c^2 k_n k_\tau} \left\{ \left( {\bf p}_{\bf k}  \cdot  \left[\hat{{\bf n }} \times {\bf k}\right]
\right) \right. \nonumber \\
& + & \left. r_s\left({\bf k'}\right)  \left( {\bf p}_{\bf k'}  \cdot  \left[\hat{{\bf n }} \times {\bf k'}\right]
\right)   \right\}, \label{eq:mtr_s}\\
E_p\left({\bf k}; \omega\right) & =& \frac{4 \pi^2 \omega}{c k_n k_\tau} \left\{ \left( \left[{\bf p}_{\bf k} \times {\bf k}\right] \cdot  \left[\hat{{\bf n }} \times {\bf k}\right]
\right) \right. \nonumber \\
& - & \left. r_p\left({\bf k'}\right) \left( \left[{\bf p}_{\bf k'} \times {\bf k'}\right] \cdot  \left[\hat{{\bf n }}\times {\bf k'}\right] \right)  \right\},
\label{eq:mtr_p}
\end{eqnarray}
where $k_\tau$ and $k_n$ correspond to the tangential and normal to the surface of the hyperbolic substrate components of the far-field wavevector
${\bf k} \equiv \left({\bf k}_\tau, {\bf k_n}\right)$,  $\hat{\bf n}$ is the  normal to the surface unitary vector,  $r_s$ and $r_p$ are the reflection coefficients at the air-substrate interface for the $s$- and $p$-polarizations,   ${\bf k'} \equiv  \left({\bf k}_\tau, - {\bf k_n}\right)$, and ${\bf p}_{\bf k}$ is the spatial Fourier transform of the polarization of the object, ${\bf p}\left({\bf r}\right) \equiv (1/4\pi) \left(\epsilon\left({\bf r}\right) - 1\right){\bf E}_i\left({\bf r},\omega\right)$, when it is illuminated by the incident field ${\bf E}_i\left({\bf r},\omega\right)$. For an object placed close to the image plane and with the vertical dimension well below the free-space wavelength ($\Delta z \ll \lambda_0$), the polarization of the object
${\bf p}_{\bf k}$ can be related to the spatial Fourier transform $\Delta_M\epsilon\left({\bf q}\right)$ of its Motti projection (\ref{eq:Motti}) :
\begin{eqnarray}
{\bf p}_{\bf k} & = & \frac{1}{8 \pi^2} \int d^2{\bf q} \ {\bf E}_i\left({\bf k + q};\omega\right) \ \Delta_M\epsilon\left({\bf q}\right).
\end{eqnarray}
Given the hyper-spectral measurement of the far field  ${\bf E}\left({\bf k}, \omega\right)$, the system of coupled linear equations (\ref{eq:mtr_s}),(\ref{eq:mtr_p}) can be used to  calculate $\Delta_M \epsilon$.
In particular, when the hyper-structured illumination pattern  ${\bf E}_i\left(x,y\right)$ at the image plane $\Delta z = 0$ originates from a periodic
array of the slits at the bottom of the substrate, separated by the distances of $d_x$ and $d_y$ in the $x$ and $y$ directions,
\begin{eqnarray}
{\bf E}_i\left(x,y; \omega\right) =   \sum_{m_x, m_y = - \infty}^{\infty}
{\bf E}_0\left(x - m_x d_x, y - m_y d_y; \omega \right),  \
\end{eqnarray}
then each spatial Fourier component of the object  profile $\Delta \epsilon\left({\bf q}\right)$ is only coupled to a discrete set of  other spatial components, that are different by an integer multiple of the reciprocal lattice vectors $2\pi/d_x \ \hat{{\bf x}}$ and
$2\pi/d_y \  \hat{{\bf y}}$:
\begin{eqnarray}
\sum_{n_x, n_y = - \infty}^{\infty} G^{(s,p)}_{n_x,n_y}\left( {\bf k}; \omega \right) \ \delta\epsilon_{n_x,n_y}
 =
\frac{2 c k_n k_\tau} { \omega}  E_{(s,p)}\left( {\bf k}; \omega \right) \
\label{eq:InvMtr}
\end{eqnarray}
where
\begin{eqnarray}
\delta\epsilon_{n_x,n_y} \equiv \Delta_M\epsilon\left({\bf k}_\tau + \hat{\bf x} \ 2\pi  \ n_x / d_x +  \hat{\bf y} \ 2 \pi n_y / d_y   \right),
\end{eqnarray}
 and
\begin{eqnarray}
G^{(s)}_{n_x,n_y}\left( {\bf k}; \omega \right) & = & \left(1 + r_s\right) \left( \left[\hat{{\bf n }} \times {\bf k}\right]   \cdot  { \bf  e}_{n_x,n_y}\left(\omega\right) \right), \nonumber \\
G^{(p)}_{n_x,n_y}\left( {\bf k}; \omega \right) & = & \frac{c}{\omega} \  \left[\hat{{\bf n }} \times {\bf k}\right]  \cdot
 \left[ { \bf   e}_{n_x,n_y}\left(\omega\right) \times  \left({\bf k} - r_p \ {\bf k'}  \right)
\right], \nonumber
\end{eqnarray}
with
\begin{eqnarray}
 { \bf   e}_{n_x,n_y}\left(\omega\right) \equiv \frac{1}{d_x d_y} \int d^2{\bf r} \ {\bf E}_0\left({\bf r};\omega\right)
 \ e^{ \ 2\pi  i \left( x n_x / d_x +   y n_y / d_y   \right) } .
 \nonumber
\end{eqnarray}
Given the desired resolution $\Delta$, this implies a summation by $n_{\rm max} \sim d/\Delta$ terms in either direction. Even with $\Delta \sim 10$ nm and diffraction-limited slit-to-slit spacing $d \sim \lambda_0/2$, we find  $n_{\rm max} \sim 10$ -- thus reducing the numerical complexity of solving Eqn. (\ref{eq:InvMtr}) to that of inverting a $100 \ {\rm x} \ 100$ matrix. On a modern CPU, such computational task can be completed within a 0.1 ms time frame. Thus, despite the need for data post-processing that is inherent to structured illumination methods, the proposed approach can be
successfully used for real-time imaging of the dynamical processes at sub-ms time scales.

\section{Hyper-structured Illumination in a noisy environment}

The  linear inversion procedure of the previous  section is robust to the presence of the noise and material absorption -- which is illustrated in Fig. 4, where panels (a) and (b) present the gray-scale images of the original
object and its image obtained in a noisy environment. For a quantitative comparison, panel (c) shows the dimensions of the original object (left scale) and the  recovered Motti projection (right scale). Note that the $30$ nm separation, which is less than  $1/10$ of the shortest illumination wavelength, is clearly resolved -- despite a substantial amount of noise in the system.

\begin{figure}[htbp]
\centering
\fbox{\includegraphics[width=2.75 in ]{fig4v2}}
\caption{Super-resolution imaging with hyper-structured illumination. Panel (a) : the ``top view'' of the object  (four identical silicon
nanowires with 10 x 5 nm cross-section) on the silver-glass hyperbolic substrate (see Fig. 2(a)). Panel (b) : the reconstructed image,
in gray-scale. The calculation includes the effects of the actual material loss in the silver-glass hyperbolic substrate. Panel (c): the object profile (red) and the  (normalized to the silicon permittivity
$\epsilon_{\rm Si}$)  reconstructed Motti projection $\Delta_M\epsilon /( \epsilon_{\rm Si} - 1)$ (blue).}
\label{fig:4}
\end{figure}

\section{Discusson and Conclusions}

The hyper-structures illumination approach introduced in the present paper, combines deeply subwavelength resolution
with the relative simplicity of {\it planar} geometry that would facilitate surface functionalization for cell targeting and imaging, essential for biomedical applications. Compared to the existing super-resolving extensions of structured illumination, such as the
plasmonic structured illumination \cite{pSI}, the proposed approach is not limited to particular resonant wavelength (such as e.g.
the surface plasmon wavelength), and is not sensitive to the material absorption.

In terms of the actual fabrication of the required metamaterial substrate, it's well within the limits of already demonstrated
capabilities, \cite{Liu} and experimental demonstration of the proposed imaging method will be straightforward.

$ { \ } $

This work was partially supported by NSF  (grant 1120923) and  Gordon and Betty Moore Foundation.


\section{Introduction}

Tungsten (W), being an important constituent of tokamak reactor walls, is perhaps the most important element for studying  fusion plasmas. It immensely radiates at almost all ionisation stages, and therefore to assess  radiation loss and for modelling plasmas, atomic data (including energy levels and  oscillator strengths or radiative decay rates) are required for many of its ions. The developing  ITER project has raised urgency for the data requirements, and as a result several groups of people are engaged in producing atomic data for its ions. However, the most desirable requirement for (any) atomic data is its accuracy$^{\cite{fst}}$ without which the modelling of plasmas will not be reliable.  Recently, S. Aggarwal {\em et al.}$^{\cite{ajm}}$, henceforth to be referred to as AJM,  have reported results for energy levels, oscillator strengths, radiative rates,  and lifetimes for Al-like W. For their calculations, they adopted the modified version of the {\sc grasp} (general-purpose relativistic atomic structure package) code, available at the website {\tt http://web.am.qub.ac.uk/DARC/}. It is a fully relativistic code,  based on the $jj$ coupling scheme, and hence ideal for generating atomic data for heavier ions, such as of W. Further relativistic corrections arising from the Breit interaction and QED (quantum electrodynamics) effects are also included.

AJM have reported energies and lifetimes ($\tau$) for the lowest 148 levels of the 3s$^2$3p, 3s3p$^2$, 3s$^2$3d, 3s3p3d, 3p$^3$, 3p$^2$3d, 3s3d$^2$, 3p3d$^2$, and 3d$^3$ (nine) configurations. They have also listed radiative rates (A-values), oscillator strengths (f-values) and  line strengths (S-values) for four types of transitions, namely electric dipole (E1), electric quadrupole (E2), magnetic dipole (M1), and magnetic quadrupole (M2), but only from the ground to higher excited levels. However, unfortunately there are several discrepancies and a few serious errors in their paper, and here we focus only on three, discussed below.

Firstly, their use of the nomenclature of the ion (being wrong) in the title and throughout the text of the paper is confusing and misleading. This is because Al-like tungsten is W LXII, and not W XLVII.   Secondly, to calculate atomic data they have included configuration interaction (CI) among 894 levels of  35 configurations, the additional 26 are: 3s3p4$\ell$, 3s3d4$\ell$, 3p3d4$\ell$, 3s$^2$4$\ell$, 3p$^2$4$\ell$ (except 3p$^2$4d), 3p4$\ell^2$ (except 3p4p$^2$), and 3d4$\ell^2$. However, these 35 configurations generate 1007 levels in total (see section 2), and hence there is a discrepancy of 113 levels. Both of these discrepancies may, at best, be attributed to oversight and do not affect the calculated results. However, the third and the final one is the gross error in their calculated lifetimes, of up to 14 orders of magnitude for over 90\% of the levels -- see section 3. Therefore, the purpose of this comment is to report the correct $\tau$ values and to explain the reasons for these large errors.

\section{Calculations}

For our calculations we adopt the same version of the  {\sc grasp}  code as employed by AJM$^{\cite{ajm}}$. Similarly,  for compatibility we use the same  option of {\em extended average level} (EAL), as adopted by them. To make comparisons, we have performed two sets of calculations, one with 148 levels of 9 configurations (GRASP1), i.e. 3s$^2$3p, 3s3p$^2$, 3s$^2$3d, 3s3p3d, 3p$^3$, 3p$^2$3d, 3s3d$^2$, 3p3d$^2$, and 3d$^3$; and another with 37 (GRASP2), the additional 28  being 3s3p4$\ell$, 3s3d4$\ell$, 3p3d4$\ell$, 3s$^2$4$\ell$, 3p$^2$4$\ell$, 3p4$\ell^2$, and 3d4$\ell^2$. All these configurations, along with the number of levels generated by each, are listed in Table 1. In total, these configurations generate 1056 levels. AJM included the same configurations (except two) in their calculations (GRASP3), but it is not clear why they ignored  3p$^2$4d and 3p4p$^2$. These two configurations (number 28 and 31 in Table 1) generate 49 levels and therefore the number of levels included by AJM should have been 1007 (1056-49), and not 894 as stated by them. However, as far as the calculations of $\tau$ (our parameter of interest) are concerned, the inclusion/exclusion of these two configurations is of no consequence as we will see in the next section.

The energies obtained in our calculations  for the lowest 148  levels of W LXII are comparable to those reported by AJM in their Table 1, and hence are not repeated here. Similarly, there is no (significant) discrepancy with their A-values, listed in their Table 3. Therefore, we now focus on the $\tau$ values for which we find large discrepancies.

\section{Lifetimes}

The lifetime $\tau$ of a level $j$ is determined as  ${\tau}_j =1/{\sum_{i} A_{ji}}$, i.e. the summation is over all transitions from lower levels $i$ to higher $j$. Generally, A-values for E1 transitions dominate for the determination of $\tau$, but for other type of transitions become important, particularly when the E1 does not exist, such as for 1--3/5/9 -- see Table~2 for level definitions. For this reason,  A-values from all types of transitions are included in a calculation, as was also done by  AJM$^{\cite{ajm}}$. In Table~2, we list our calculated $\tau$ from both GRASP1 and GRASP2 calculations, and also include the results of AJM (GRASP3) for a ready comparison.

For all 148 levels listed in Table~2, there is no (appreciable) discrepancy between the GRASP1 and GRASP2 results. Therefore, as stated earlier, the inclusion or exclusion of the 3p$^2$4d and 3p4p$^2$ configurations, omitted by AJM, is of no consequence, and hence the  values of $\tau$ calculated by AJM should have been comparable to those of ours. Unfortunately, there are large discrepancies for almost all levels, over 90\% to be precise. For all levels, the reported values of $\tau$ by AJM are invariably {\em higher},  by up to 14 orders of magnitude, see for example levels 11, 12, 25, 26, 37, and 82. Below we explain the possible reasons for these large discrepancies.

For level 3 (3s$^2$3p $^2$P$^o_{3/2}$), the $\tau$ by AJM is higher by a factor of two. For this level, the AJM calculation appears to include only the A-value (2.18$\times$10$^8$ s$^{-1}$) for the 1--3 M1 transition  -- see their Table~3, whereas the A-value for the 2--3 E1 transition (2.64$\times$10$^8$ s$^{-1}$, not listed by AJM) has been completely ignored. Similarly for level 11 (3s3p3d $^4$F$^o_{3/2}$), their calculation is based on the A-value (2.95$\times$10$^6$ s$^{-1}$) of 1--11 M1 transition alone. Interestingly, they have also listed the A-value for the 1--11 E2 transition, i.e. 9.38$\times$10$^6$ s$^{-1}$. The inclusion of these two A-values alone will give $\tau$ = 8.1$\times$10$^{-8}$~s, and {\em not} 3.39$\times$10$^{-7}$~s, as listed by them in their Table~1. However, the correct value for this level is $\sim$2.1$\times$10$^{-12}$~s, as listed in our Table~2. This is because the main contributing transitions are:  2--11 E1 (A=3.72$\times$10$^{11}$ s$^{-1}$), 5--11 E1 (A=5.26$\times$10$^{10}$ s$^{-1}$) and 6--11 E1 (A=3.29$\times$10$^{10}$ s$^{-1}$). Finally, the largest discrepancy of 14 orders of magnitude is for level 82, i.e. 3p3d$^2$ $^4$F$^o_{5/2}$. For this level, AJM have listed A=3.66$\times$10$^{5}$ s$^{-1}$ for the 1--82 E2 transition, which alone gives $\tau$ = 2.73$\times10^{-6}$~s, and {\em not} 5.87$\times$10$^{1}$~s, as listed by them in their Table~1. For this level the dominant contributing transition is 23--82 E1 for which A=2.24$\times$10$^{12}$ s$^{-1}$, giving $\tau$ = 4.46$\times10^{-13}$~s, a value within $\sim$10\% of that listed in present Table~2. Therefore, it is abundantly clear that the $\tau$ values reported by AJM$^{\cite{ajm}}$ are simply wrong for most of the levels.

\section{Conclusions}

In this work, we have calculated energy levels, radiative rates and lifetimes for the lowest 148 levels/transitions of Al-like tungsten, i.e. W LXII. For the calculations, the well known {\sc grasp} code has been adopted, as by S. Aggarwal {\em et al.}$^{\cite{ajm}}$ who have recently reported similar results. However, their results for energy levels and A-values (although for limited transitions) are correct, the corresponding values of $\tau$ are in huge errors, of up to 14 orders of magnitude, for almost all levels. In fact, their reported $\tau$ appear to have no relationship with the A-values for transitions of this ion. Therefore, for the benefit of users as well as for future workers, we have listed the correct values of $\tau$, and have also explained the reasons for discrepancies. Apart from the errors in $\tau$, other discrepancies in the paper of AJM have been noted. However, we stress that the listed anomalies are not the only ones, there are a few more. As an example, the 101--129 levels listed in their Table~1 have no (odd) parity indication.

Finally, the A-values listed by S. Aggarwal {\em et al.}$^{\cite{ajm}}$ for transitions from the ground level alone are not sufficient for applications, because a {\em complete} set of data for {\em all} transitions are required for any modelling application. Besides this, there is scope for improvement in their work, because in their CI calculations they have ignored the inclusion of some of the important configurations, such as 3p$^2$4d, 3p4p$^2$, 3d$^2$4$\ell$,  and 3$\ell$4$\ell\ell'$, apart from those of $n$ = 5. Therefore, an improved set of complete data for all transitions of W LXII are reported in a separate paper$^{\cite{w62}}$.

\newpage