\section{Introduction}
Gamma-ray bursts (GRBs) were discovered almost fifty years ago \cite{bib:GRB_review}, but despite thorough study, details of their particle acceleration mechanisms still remain in the dark. A recent observation of the Large Area Telescope (LAT) on board the \emph{Fermi Gamma-Ray Space Telescope} (\emph{Fermi}-LAT) has established that GRBs produce photons in the very-high-energy (VHE, $>100$~GeV) regime, when a $95.3~\rm{GeV}$ photon (or $128~\rm{GeV}$ when corrected for redshift) was detected in GRB~130427A \cite{bib:GRB130427A_Fermi_Science}. The effective area of satellite experiments is very limited at these energies and since the photon flux decreases steeply with energy, there is little hope of detecting more than just single photons. The highest energies and temporal evolution of GRB spectra carry important information about GRB physics, e.g. would the observation of VHE photons at late or early times challenge scenarios, in which those photons are produced by synchrotron emission. Furthermore, VHE observations can also provide insights into the extragalactic background light (EBL) and Lorentz invariance violation.

At VHEs, ground-based experiments might be able to supply those important measurements. Imaging Atmospheric Cherenkov Telescopes (IACTs) are pointed instruments that need to slew to the GRB position and will therefore in general miss the prompt and early afterglow phase. In contrast, the High Altitude Water Cherenkov (HAWC) observatory has a large instantaneous field of view ($\sim2~\rm{sr}$ or 16\% of the sky), $>95$\% duty cycle and the lack of observational delays will allow observing the GRB prompt phase. This makes HAWC an ideal detector for studying transient sources like GRBs.

\section{High Altitude Water Cherenkov Observatory}
HAWC is a new, ground-based VHE gamma-ray air-shower detector that recently finished construction at Sierra Negra, Mexico, at an altitude of 4,100~m above sea level \cite{bib:HAWC_highlights}. It applies the water Cherenkov technique, where VHE photons are detected by measuring Cherenkov light from secondary particles in an extensive air shower. The full detector consist of 300 steel tanks, containing light-tight bladders of 7.3~m diameter and 4.5~m depth, holding $\sim200,000$~litres of filtered water. The Cherenkov light is measured by three $8^{\prime\prime}$ photomultiplier tubes (PMTs) and one $10^{\prime\prime}$ PMT at the bottom of each tank. The PMT pulses are shaped and discriminated at two thresholds and the threshold crossing times are digitised. All signals are recorded to memory and triggering is entirely done in software by applying a threshold on the number of coincident PMT signals in a 150~ns time window (trigger threshold).

HAWC features two data acquisition (DAQ) systems. One is designed to read out full air-shower events (main DAQ), recording the time and charge of individual PMT pulses and the signal arrival time in different tanks. This makes it possible to reconstruct the incident direction of the shower. The dominant background when observing VHE gamma-rays comes from the abundant population of hadronic cosmic rays. A set of cuts can be used to remove hadronic events (gamma-hadron separation). The other DAQ (scaler DAQ) counts the signals in each PMT, allowing to detect GRBs by a statistical excess over the noise rate (``single particle technique'' \cite{bib:scaler_method}). Both DAQs have sensitivity to GRBs at different energies and therefore complement each other \cite{bib:HAWC_GRBs}. In this contribution, the analysis of the main DAQ data is described.

Science operations of the HAWC array started with a partially built array in August 2013 \cite{bib:Andy-ICRC-2015}. The stage of the detector between August 2, 2013 and July 8, 2014 is called HAWC-111. During this time, the uptime fraction was 83\%, which includes shutdowns to accommodate ongoing construction. The detector operated with a trigger threshold between 15 and 25 PMTs. The data presented in this contribution are from the HAWC-111 period.

Data was collected before the HAWC-111 period with a smaller detector. Into that period falls GRB~130427A, the most powerful burst ever detected with a redshift $z\lesssim0.5$. HAWC provided the first limits on the prompt emission of GRB~130427A in the VHE range, but due to the high zenith angle of the GRB and the incomplete HAWC detector, the limits are about two orders of magnitude higher than a simple extrapolation of the \emph{Fermi} data \cite{bib:HAWC_GRB130427A}.

\section{GRB Selection}
GRBs in the field of view during HAWC-111 operations were found using the LAT\footnote{http://swift.gsfc.nasa.gov/archive/grb\_table/}, \emph{Fermi} Gamma-ray Burst Monitor (\emph{Fermi}-GBM)\footnote{https://heasarc.gsfc.nasa.gov/W3Browse/all/fermigbrst.html} and \emph{Swift}\footnote{http://swift.gsfc.nasa.gov/archive/grb\_table/} online tables. The GRBOX: Gamma-Ray Burst Online Index\footnote{http://www.astro.caltech.edu/grbox/grbox.php} was searched without a match for additional GRBs from other experiments.

Considering GRBs down to a zenith angle of 51 degrees, there was only one LAT detected GRB in the HAWC field of view (GRB~130907A) at a zenith angle of $27\arcdeg$. Unfortunately, this burst occurred during one of the detector downtimes and no data is available. About 40 GBM bursts happened during HAWC-111 and for six of them there is no data available. Since October 2013, only one GBM burst has no data available, which demonstrates the excellent uptime fraction of HAWC. GBM-detected bursts are often poorly localised and require a different analysis technique than what is presented here, so this contribution focuses on the well localised \emph{Swift} bursts.

During HAWC-111, 22 \emph{Swift}-detected GRBs occurred in the HAWC field of view. Three of them (GRB~130907A, GRB~130912A and GRB~131018A) were missed due to detector downtime. GRB~140226A, which is the first burst discovered optically and independent of a high-energy trigger by the Intermediate Palomar Transient Factory \cite{bib:iPTF_GRB140226A}, peaked in Konus-\emph{Wind} at 10:02:57~UTC. Unfortunately, no data is available for the prompt phases of this GRB due to a run transition. Table~\ref{tab1} lists the 18 bursts analysed in this contribution.

\begin{table}
\footnotesize
\begin{tabular}{l|ccccccc}
\hline
GRB & \linebreakcell[t]{Trigger \\ Number} & \linebreakcell[t]{Time \\ UTC} & RA J2000 & DEC J2000 & \linebreakcell[t]{Zenith Angle \\ deg} & \linebreakcell[t]{BAT T90 \\ s} & \linebreakcell[t]{Significance \\ $\sigma$} \\
\hline
140628A & 602803 & 13:35:37 & 02h42m39.88s & -0d23m05.7s & 26.0 & 10.5 & -0.74 \\
140622A & 602278 & 09:36:04 & 21h08m41.56s & -14d25m09.5s & 33.4 & 0.13 & -0.93 \\
140607A & 601051 & 17:13:31 & 05h45m29.52s & 18d54m14.4s & 27.9 & 109.9 & 3.42 \\
140518A & 599287 & 09:17:46 & 15h09m00.60s & 42d25m05.6s & 48.6 & 60.5 & -0.61 \\
140430A & 597722 & 20:33:36 & 06h51m44.61s & 23d01m25.2s & 31.3 & 173.6 & -1.75 \\
140423A & 596901 & 08:31:53 & 13h09m08.54s & 49d50m29.4s & 46.9 & 134 & 0.21 \\
140419A & 596426 & 04:06:51 & 08h27m57.56s & 46d14m25.3s & 45.3 & 94.7 & 1.35 \\
140414A & GA     & 06:06:29 & 13h01m14.40s & 56d54m07.2s & 37.8 & 0.7 & -0.18 \\
140408A & 595141 & 13:15:54 & 19h22m51.83s & -12d35m42.5s & 32.4 & 4.00 & -0.02 \\
140331A & 594081 & 05:49:48 & 08h59m27.46s & 02d43m02.3s & 45.7 & 209 & -2.18 \\
140215A & 586680 & 04:07:10 & 06h56m35.81s & 41d47m11.7s & 23.2 & 84.2 & 0.30 \\
140206A & 585834 & 07:17:20 & 09h41m20.26s & 66d45m38.6s & 47.7 & 93.6 & -1.86 \\
140129A & 585128 & 03:23:59 & 02h31m33.78s & -01d35m43.4s & 47.8 & 2.99 & 1.65 \\
140114A & 583861 & 11:57:40 & 12h34m05.16s & 27d57m02.6s & 11.1 & 139.7 & 0.29 \\
131229A & 582374 & 06:39:24 & 05h40m55.61s & -04d23m46.7s & 27.7 & 13.86 & 1.23 \\
131227A & 582184 & 04:44:51 & 04h29m30.78s & 28d52m58.9s & 10.1 & 18.0 & -0.48 \\
131117A & 577968 & 00:34:04 & 22h09m19.36s & -31d45m44.3s & 50.9 & 11.00 & 0.27 \\
131001A & GA     & 05:37:24 & 00h33m12.96s & 25d33m25.2s & 12.4 & 4.9 & 0.96 \\
\end{tabular}
\caption{The 18 \emph{Swift}-detected GRBs that occurred in the HAWC field of view during HAWC-111 operations for which data is available. GA means the burst was found in a ground analysis.}
\label{tab1}
\end{table}

\section{Analysis Method}
The analysis described here was originally developed to run at the HAWC site (online), but is currently also used for offline data analysis. No gamma-hadron separation is applied, because the most efficient gamma-hadron separation were not available online. Since the EBL cuts of the GRB spectra at TeV energies, most of the signal is expected at lower energies, were the efficiency of the gamma-hadron separation is lowest. Therefore, at the lowest energies the gamma-hadron separation will only modestly increase the sensitivity.

The analysis consist in defining an angular bin (the search bin) around the position of the GRB, determining an estimated number of background events that will appear in the search bin and measuring the number of events above background in the bin. Besides the angular bin, the only other cut applied is that a successful angular fit of the event is required.

The size of the angular bin ($\Delta\Psi$) is optimised using Monte Carlo simulations of gamma and cosmic-ray air showers and a detailed simulation of the detector. The optimisation was done for a GRB with an $E^{-2}$ spectrum, which is the typical spectrum observed by \emph{Fermi}-LAT at high energies \cite{bib:Fermi_LAT_GRB_catalogue}. The burst was simulated at various zenith angles and redshifts, modelling the absorption on the EBL according to the fiducial model in \cite{bib:EBL}. The different trigger thresholds for the HAWC-111 operation were also simulated. The optimal $\Delta\Psi$ does not depend strongly on zenith angle or trigger threshold. It does however depend on the assumed redshift, because the angular resolution of the detector depends on the number of PMTs with signal in the array, which correlates with the energy of the primary gamma ray. As the redshift increases, the signal is dominated by lower energetic gamma rays that have poorer resolution. For most GRBs, the redshift remains unknown, therefore the angular bin is conservatively chosen to be $3\arcdeg$. The optimisation of the search circle for two redshifts is shown in Fig.~\ref{angular_cut_optimisation}.

\begin{figure}[htbp]
\centering
\includegraphics[width=.5\textwidth]{icrc2015-0237-01.pdf}
\caption{The figure shows the signal (S) over the square root of the background (B) versus the size of the angular bin for a GRB with an $E^{-2}$ spectrum in HAWC-111. The solid (dotted) curve corresponds to a GRB at redshift 0.3 (2.0), assuming EBL attenuation according to the fiducial model in \cite{bib:EBL}. The maximum of both curves is scaled to 1 for clarity.}
\label{angular_cut_optimisation}
\end{figure}

Currently, the search duration is $T_{90}$, the central time interval, in which 90\% of the prompt \emph{Swift} Burst Alert Telescope (BAT) flux is detected, starting at the BAT trigger time. Figure~\ref{allsky_event_rate_GRB140423A} shows the allsky event rate around the time of GRB~140423A in bins of width $T_{90}$. The rate is fit with a constant plus a sinusoidal oscillation with a 12~h period on top. The sinusoidal variation of the order of $\sim0.5$\% comes from the atmospheric tide that affects the air shower propagation. The magnitude of this effect is very stable over the last two years. The distribution of excess between the data and the fit is only slightly broader than expected for a perfect Poissonian counting experiment.

\begin{figure}[htbp]
\centering
\includegraphics[width=.7\textwidth]{icrc2015-0237-02.pdf}
\caption{Allsky event rate around the trigger time of GRB~140423A at zero. The data are well fitted by a constant plus a sinusoidal modulation on top. The fitted amplitude is relative to the constant.}
\label{allsky_event_rate_GRB140423A}
\end{figure}

The background estimate is obtained using the classical ON/OFF method (for example used in the Whipple observatory \cite{bib:ON_OFF}). In Fig.~\ref{on_source_event_rate_GRB140423A}, when the time equals zero, the search bin is centred on the GRB position and follows it for the duration of $T_{90}$ (on-observation). For all other times (off-observations), the search bin is offset in Right Ascension by multiples of $T_{90}$, thus observing an empty field of the sky that covers equal zenith angles as the on-observation. Due to the continuous data taking of HAWC, off-observations are available before and after the GRB. A background estimate at each point is obtained by taking the fit of the allsky event rate and scale it to match the summed counts in Fig.~\ref{on_source_event_rate_GRB140423A}. Counts between -380~s and 680~s are not used in the sum, because due to the size of the HAWC point spread function it takes some time before the GRB position has moved out of the search bin. It can be seen that the sinusoidal oscillation is a negligible correction in a $3\arcdeg$ search circle and fitting the off-observations with a constant yields consistent results. The counts in the on-observation are compared to the estimate from the background model obtained using the off-observations and the significance is calculated assuming a Poisson distribution.

\begin{figure}[htbp]
\centering
\includegraphics[width=.7\textwidth]{icrc2015-0237-03.pdf}
\caption{Event rate for on-observation (at zero) and off-observations before and after the trigger time of GRB~140423A. The red line shows the background model that is created by scaling the allsky event rate fit. The horizontal dashed lines show the region where data is not used for the scaling of the background model.}
\label{on_source_event_rate_GRB140423A}
\end{figure}

\section{Results}
Fig.~\ref{significance_distribution} shows the significance distribution that is obtained by applying the analysis method described above to 25,188 fake GRB positions in HAWC-111 data, which are placed at zenith angles between 0\arcdeg and 60\arcdeg. The distribution is very well fitted by a Gaussian, showing that the background is well understood to the 4~$\sigma$ level.

\begin{figure}[htbp]
\centering
\includegraphics[width=.7\textwidth]{icrc2015-0237-04.pdf}
\caption{The blue histogram shows the significance distribution obtained from analysing the 18 \emph{Swift}-detected GRBs in Table~\protect\ref{tab1}. The blue points show the significance distribution that is obtained by applying the analysis method to 25,188 fake GRB positions. The results of the Gaussian fit (red line) are given in the inset box.}
\label{significance_distribution}
\end{figure}

Applying the analysis method to the 18 \emph{Swift}-detected GRBs from Table~\ref{tab1} gives the blue histogram of significances in Fig.~\ref{significance_distribution} (the individual significance for each burst are listed in Table~\ref{tab1}). GRB~140607A is the only burst with a significance above 3~$\sigma$. Accounting for the trial factor of looking at 18 bursts at once, the significance drops to 2.5~$\sigma$. GRB~140607A was detected by the BAT only, because due to a Sun observing constraint, \emph{Swift} couldn't slew to the BAT position. No other observations were reported for this GRB.

\section{Outlook}
In this contribution initial results for GRBs in the HAWC field of view during HAWC-111 were shown. None of the GRBs is significant above 3~$\sigma$ when accounting for trial factors. An alternative self-triggered search for GRBs is also planned for HAWC \cite{bib:Josh-ICRC-2015}. This search has a worse detection threshold than the one presented here, but it may find bursts not reported by satellites.

\emph{Fermi}-LAT found that the $>100\nobreakspace\rm{MeV}$ emission of GRBs is temporally extended and starts later than the keV--MeV emission \cite{bib:Fermi_LAT_GRB_catalogue}. This implies that time windows other than the BAT $T_{90}$ should be searched for VHE emission as well. The angular bin size has been optimised for the next stage of the detector, called HAWC-250, and the application of gamma-hadron separation cuts for both online and offline analysis are being explored. Gamma-hadron separation is most efficient at $>$TeV energies, but the application of cuts on HAWC-250 can still improve the detection significance by 50\% to 100\%. Results of offline and online analysis of HAWC-250 and the full HAWC array will be presented elsewhere.

It has been show elsewhere that HAWC is sensitive enough to detect several historical GRBs (like GRB~090510 and 090902B) if their emission extends only slightly beyond the highest energy observed by \emph{Fermi}-LAT \cite{bib:HAWC_GRBs}. Furthermore, a GRB similar to 130427A would be detected even if the spectrum doesn't extend beyond what was observed by the LAT \cite{bib:HAWC_GRB130427A}. In addition to exceptional bursts, HAWC might detect other GRBs with a rate as high as 1--2 GRBs per year \cite{bib:HAWC_GRB_rate}. Even the absence of detections will be useful in constraining the highest energy of GRB spectra.

\section*{Acknowledgments}
\footnotesize{
We acknowledge the support from: the US National Science Foundation (NSF);
the US Department of Energy Office of High-Energy Physics;
the Laboratory Directed Research and Development (LDRD) program of
Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnolog\'{\i}a (CONACyT),
Mexico (grants 260378, 55155, 105666, 122331, 132197, 167281, 167733);
Red de F\'{\i}sica de Altas Energ\'{\i}as, Mexico;
DGAPA-UNAM (grants IG100414-3, IN108713,  IN121309, IN115409, IN111315);
VIEP-BUAP (grant 161-EXC-2011);
the University of Wisconsin Alumni Research Foundation;
the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory;
the Luc Binette Foundation UNAM Postdoctoral Fellowship program.
}


\section{Introduction}

Numerous scientific, medical, industrial, and verification applications demand high-resolution spectral measurements, which can be accomplished using Fabry-Perot (FP) interferometry, among other methods~\cite{VanDerWal1979,Zhou2012,Vogus2015}. Recently, FP interferometry has been adopted in pulsed measurements via laser-induced breakdown spectroscopy (LIBS)~\cite{Effenberger2012}. Spectra from laser-produced plasmas are commonly measured by use of angular dispersion, such as in Czerny-Turner (CT) or echelle spectrometers~\cite{Aragon2008}. In LIBS, high-resolution measurements are especially valuable for measurements of isotope ratios, where the typical isotopic shifts are on order a few to several tens of picometers. Currently, rapid and economical techniques with isotopic specificity are highly sought for nuclear safeguards and non-proliferation applications. Development of LIBS to detect and monitor facilities for illicit activities (\textit{e.g.}, isotope enrichment to undeclared levels)~\cite{Hanson2014,Ko2013} is critical to risk preparedness and treaty verification. However, high-resolution isotope measurements presently require very large (few-meter scale) CT spectrographs~\cite{Smith2002} or echelle spectrographs with complex instrument functions and limited dynamic range~\cite{Cremers2012}. More portable, simple, and cost-effective options are desirable in field applications. A practical alternative is an inexpensive, compact FP etalon coupled to a much lower resolution CT imaging spectrometer~\cite{Effenberger2012}.

Accurate analysis of the 2-D intensity pattern measured by the instrument is important for interpreting profiles of spectral lines. A variety of algorithms have been developed to deconvolve the input spectrum from the interferometric measurement, but they rely on assumptions made about the shape (\textit{e.g.}, Gaussian, Lorentzian, Voigt functions) of the source spectrum~\cite{Larson1967,Koo1987,Naylor1991}. Other techniques include analysis in the Fourier domain, which also assumes a known source distribution~\cite{Cooper1971}. Direct least squares fitting of 2-D ring images from a single spectral line has been proposed, but is mainly aimed at calibrating and characterizing visible laser diodes~\cite{OHora2005}. Another alternative is data inversion~\cite{Zipoy1979,Heeg2007,Abbiss2008}, which can be effective in cases where the spectrum is unknown. Abbiss and Heeg measured spectra by varying the input angle of a tunable filter~\cite{Heeg2007}. While such scanning methods are not suitable for pulsed, imaging LIBS applications, similar mathematical principles were used in this case. The reconstruction method employed to interpret the high-resolution spectrum captured by a hybrid interferometric/dispersive spectrometer takes into account the intensity distribution of fringes observed in a 2-D image, and can be applied to spectra of a broadband source without presuming the shape of the input spectrum. In order to combine the benefits of 2-D algorithms with the speed necessary for rapid and high-resolution fringe measurement, this paper presents a simplified and robust fitting algorithm. The technique was also experimentally verified with an actual FP fringe imaging system.

The hybrid instrument, when used in conjunction with the reconstruction technique discussed here, represents an effective tool for measuring the high-resolution spectrum with accurate peak intensity ratios. Rapid, high-resolution spectral reconstruction capability from this type of hybrid measurement is ideal in certain applications of LIBS. The resolution of the hybrid interferometric-dispersive approach is at least an order of magnitude greater than the resolution of the CT spectrometer alone. This instrument can achieve the desired high resolution for measurement of the atomic isotopic shift of uranium, while reducing the size and cost of the instrument by approximately an order of magnitude.

\section{Instrument model and spectral analysis}

The fringe pattern from an etalon exhibits angular and spectral dependence arising from constructive and destructive interference. A broad and complex spectrum incident on the FP etalon may produce a transmission pattern that cannot be uniquely reconstructed. This is due to degeneracy inherent in the interference pattern produced by the etalon, which exhibits a characteristic periodicity and wavelength-angle coupling, whereby an identical phase shift can be produced by different combinations of wavelength and transmission angle. This degeneracy can be broken by coupling the etalon to a dispersive spectrometer, which provides an additional angular dispersion, as described by Effenberger~\textit{et~al.}~\cite{Effenberger2012}. A diagram illustrating the approach to image sampling and dispersion is shown in Fig.~\ref{figure:Setup} and the image incident on the 2-D detector array is shown in Fig~\ref{figure:Normalize}(a). The high-resolution spectral information is encoded on the axis parallel to the spectrometer slit ({\it i.e.}, orthogonal to the dispersion plane of the spectrometer). The production of the interference ring pattern by the FP etalon can be described as
\begin{equation}
\label{equation:MatEq}
 B(\theta) = \int_0^\infty{T(\theta,\lambda)A(\theta,\lambda)d\lambda},
\end{equation}
where $A(\theta,\lambda)$ is the spectrally and angularly dependent intensity of the light incident on the etalon and $B(\theta)$ is the spectrally integrated intensity of the transmitted light. The FP etalon transmission function is described by the Airy equation,
\begin{equation}
\label{equation:Airy}
T(\theta, \lambda) = \left(1 + F\sin^2\left(
2 \pi d\,n/\lambda  \cos \theta
\right) \right)^{-1},
\end{equation}
where $F$ is the coefficient of finesse, given by $F=4R/\left(1-R\right)^2$ where $R$ is the etalon reflectivity, $d$ is the separation of the etalon mirrors, $n$ is the refractive index in the space between the mirrors, $\lambda$ is the wavelength, and $\theta$ is the angle measured with respect to the etalon axis. In a typical FP optical arrangement, the transmission pattern is measured in the focal plane of a lens (far field). The fringe pattern measured in this way is comprised of concentric rings of light corresponding to the angularly and spectrally dependent interference maxima arising from propagation through the FP etalon. The maxima appear when the condition $\cos \left(\theta\right)/\lambda = k \pi /\left(2 \pi d n\right)$ is met (where $k$ is an integer). Degeneracy arising from the coupling of angle and wavelength can be broken by introduction of angular dispersion, which leads to the measurement of
\begin{equation}
\label{equation:intpatt2D}
B(\theta,\lambda) = T(\theta,\lambda)A(\theta,\lambda).
\end{equation}
Here, $B(\theta,\lambda)$ represents the recorded 2-D interference pattern. The measurement of this pattern is typically limited by the available spectral resolution, which is dictated by the properties of the arrangement used to produce dispersion and the size of the detector pixels. At any wavelength $\lambda_0$, the quantity measured is
\begin{equation}
B(\theta,\lambda)|_{\lambda=\lambda_0} = \int_{\lambda_0-\epsilon/2}^{\lambda_0+\epsilon/2}T(\theta,\lambda)A(\theta,\lambda)d\lambda,
\end{equation}
where $\epsilon$ is the spectral resolution at wavelength $\lambda_0$. It is assumed for the purpose of this analysis that the resolution in the angular direction ($\theta$) is sufficient to accurately measure the shape of $B(\theta,\lambda)|_{\lambda=\lambda_0}$. It is further assumed that light incident on the etalon does not exhibit angular dependence of the spectrum, so that the angular dependence of $B(\theta,\lambda)$ arises solely from the nonuniform angular intensity profile and $B(\theta)=T(\theta,\lambda)A(\lambda)$. The light is incident on the etalon from the output of an optical fiber, and the Gaussian model chosen for the angular dependence of the intensity distribution is
\begin{equation}
\label{equation:fit}
B(\theta) = B_0 \exp \left(-(\theta-\theta_0)^2/w^2 \right) +C,
\end{equation}
where $B_0$ is the peak intensity of the Gaussian envelope, $\theta_0$ is the angle corresponding to the etalon axis, $C$ is a constant, and $w$ is the width of the distribution envelope~\cite{Scholl2004,Hirschberger2011}. The envelope function can be used to correct the measured FP pattern for its angular intensity dependence by fitting it to the Gaussian envelope $B(\theta)$.

This model of the measurement can be used for reconstructing the high-resolution features in the input spectrum directly from the measured pattern produced from the hybrid FP-dispersive instrument, without making additional assumptions regarding the shape of the spectrum. In a real measurement, the discretization and normalization (to the best-fit envelope function $A(\theta)$) of the measured transmission pattern in the neighborhood of the wavelength $\lambda_0$ yields the vector
\begin{equation}
\mathbf{B}_\theta = \mathbf{T}_{\theta,\lambda} \mathbf{A}_{\lambda}.
\end{equation}
Here, $\mathbf{T}_{\theta,\lambda}$ is the transmission matrix constructed by discretization of Eq.~(\ref{equation:Airy}), and $\mathbf{A}_{\lambda}$ is the vector representing the spectrum in the vicinity of the wavelength $\lambda_0$. The high-resolution spectrum can then be reconstructed as
\begin{equation}
\label{equation:inverse}
\mathbf{A}_\lambda = \mathbf{T}^{-1}_{\theta,\lambda} \mathbf{B}_{\theta},
\end{equation}
which is, in practice, reduced to finding the least-squares solution to this inverse problem.

\section{Experimental results and discussion}

The reconstruction method was applied to data collected from a hybrid FP instrument that used a 550 mm focal length imaging spectrometer (Horiba Jobin Yvon) equipped with a 1800~mm$^{-1}$ grating and a 1024$\times$1024 pixel array \mbox{iStar} ICCD (Andor) composed of 13~\textmu m pixels. Coupled to the entrance of the spectrometer was a FP etalon (SLS Optics). Light was transported to the etalon by a multimode 400~\textmu m core diameter optical fiber (Ocean Optics). Fiber-optic connector, etalon, and lens were affixed in a cage setup, with multi-axis adjustments to enable accurate alignment of the circular transmission pattern to the spectrometer slit. A 25~mm diameter lens with a focal length of 100~mm was used to focus the output of the etalon onto the entrance slit of the spectrometer. The light source used was a Hg(Ar) calibration pen lamp (Oriel).

The FP etalon was coupled to the entrance slit of the spectrometer to produce fringes, as shown in Fig.~\ref{figure:Normalize}(a). The section of the spectrum corresponding to the 313.1555/313.1844~nm doublet was summed along the horizontal axis to obtain the 1-D profile of the transmitted intensity, shown in Fig.~\ref{figure:Normalize}(b). The symmetric pattern was collapsed to summing about the etalon axis to improve the signal-to-noise ratio (SNR). A Gaussian envelope function was fit to the fringes exhibiting a characteristic double-peaked shape to obtain the $A(\theta)$ correction, as shown in Fig.~\ref{figure:Normalize}(c). The normalized transmitted intensity pattern is shown in Fig.~\ref{figure:Normalize}(d).

It is necessary to calibrate the angular axis before formulating the instrument function. This calibration was performed using the neighboring 312.5674~nm line (Fig.~\ref{figure:Normalize}(a)). The angular position of highest intensity (peak) in the fringe profile was used to perform the calibration in Fig.~\ref{figure:FitR}. Because the paraxial approximation is valid, the calibration can be described by a linear function, also shown in Fig.~\ref{figure:FitR}.
The region of the fringes selected for this fit did not include the area near the optical axis or the outermost fringes, which may exhibit distortions due to their low intensity and possible optical aberrations.

Next, the instrument function for the 313~nm doublet was determined. The ranges of $\theta$ and $\lambda$ used to construct the instrument function were 20 mrad $<\theta<$ 32 mrad and 313.106~nm $<\lambda<$ 313.234~nm, respectively. The etalon had air-spaced ($n=1$) mirrors separated by a distance of $d=0.44$~mm. The etalon reflectivity ($R=0.73$) was determined by fitting the measured intensity profile to the known FP etalon transmission function, Eq. (\ref{equation:Airy}), as shown in Fig.~\ref{figure:FitR}.

Interpolation between points and increased sampling of the experimental spectrum was performed prior to reconstruction of the source spectrum, $\mathbf{A}_\lambda$. The least-squares solution to the system of linear equations associated with Eq.~(\ref{equation:inverse}) could have been used to reconstruct $\mathbf{A}_\lambda$. However, Eq.~(\ref{equation:inverse}) is an example of a problem where $\mathbf{T}_{\theta, \lambda}$ is an ill-conditioned matrix~\cite{Abbiss2008}. Hence, the pseudoinverse $\mathbf{T}_{\theta, \lambda}^{+}$ of the matrix was used to calculate a least-squares solution to the system of equations using
\begin{equation}
\mathbf{T}_{\theta,\lambda}^{+} \mathbf{B}_\theta = \mathbf{A}_\lambda.
\end{equation}
The pseudoinverse solution was implemented in Mathematica (Wolfram Research) using the function \texttt{pseudoinverse}. Because of the ill-condition of the matrix $\mathbf{T}_{\theta,\lambda}$, it was necessary to specify the tolerance parameter. The singular values matrix are the non-negative square roots of the eigenvalues of $\mathbf{T}_{\theta,\lambda}^* \mathbf{T}_{\theta,\lambda}$, where $\mathbf{T}_{\theta,\lambda}^*$ is the conjugated transposed matrix. Singular values smaller than the product of the tolerance and the largest singular value of $\mathbf{T}_{\theta,\lambda}$ were omitted~\cite{Golub1970}. In practice, the choice of tolerance controls the influence of smaller singular values in the composite error of the reconstruction, and is a trade-off between reducing round-off error at the expense of a larger residual. A tolerance of 0.1 was used throughout the analysis to yield the optimal approximation to the solution.

Minimization of the residual error of the reconstruction was the next step. The widths of Hg lines from a low-pressure discharge calibration lamp are approximately 1~pm~\cite{Cremers2012}. The width of the reconstructed spectrum is dictated primarily by the reflectivity $R$. Because the widths are narrow compared to the instrument function, an accurate reconstruction assuming this experimentally determined reflectivity was not possible. By assuming a broader source spectrum, it was possible to minimize residual error in the spectral reconstruction. A reflectivity value of $R_\alpha = \alpha R$, where $\alpha$ is a multiplicative constant of $R$ that minimizes the residual error, $\parallel \mathbf{B}_{\theta} - \mathbf{T}_{\lambda, \theta} \mathbf{A}_{\lambda} \parallel$ was assumed. An optimal value $\alpha = 1.16$, obtained by scanning a range of values 1.00$<\alpha<$1.30, was utilized to minimize the residual error of the reconstruction.

A spectrum of Hg measured using the CT spectrometer without the FP etalon is shown in Fig.~4(a) for comparison. In this case the resolution ($\sim$70~pm) was not sufficient for resolving the 313~nm doublet. The high-resolution reconstruction provided by the hybrid instrument is shown in Fig.~\ref{figure:Reconstruct}(b) and summarized in Table~1. The resolution of the 313.1844~nm line was $\sim$4~pm, which exceeds the resolution of the CT spectrometer by more than one order of magnitude. These results demonstrate that, while the resolution results from the design of the entire hybrid instrument, the resolution enhancement is a consequence of the introduction of the FP etalon. An accurate reconstruction method is critical to fully realize the utility of a hybrid instrument. For comparison, the 313.1844~nm line in Fig.~\ref{figure:Normalize}(d) has a resolution of $\sim$8 pm using more traditional methods for converting the fringe pattern to a spectrum. Additionally, Effenberger~\textit{et~al.}~estimated the resolution to be $\sim$11~pm from the measured peak separation in spectra acquired with an Hg(Ar) lamp coupled into a FP attached to a 500~mm CT spectrometer with a 512$\times$512, 24~\textmu m/pixel detector array.

Correcting for the angular drop in intensity of the measured FP pattern was critical for peak ratio determination; the reconstructed peaks of non-normalized data yielded a ratio of 0.629, which is lower due to the decreasing intensity as a function of distance from the interferogram center. Validation of the reconstructed hybrid instrument spectrum was performed by comparing the results with those acquired with a high-resolution echelle spectrometer (DEMON, Lasertechnik Berlin) using the same Hg(Ar) lamp. The DEMON results are also provided in Fig.~\ref{figure:Reconstruct} and Table~\ref{table:Ratios}, along with literature values~\cite{Reader1996}.

Averages reported in Table~\ref{table:Ratios} are based on 10, 5, and 20 spectra for sources 1, 2, and 3, respectively. The reconstructed peak wavelengths were 313.156~nm and 313.185~nm for the Hg doublet, and agree well with actual values. While the resolution is not as high for the hybrid instrument, the spectra are similar (Fig.~\ref{figure:Reconstruct}(b)) and give a very similar ratio accuracy (Table~\ref{table:Ratios}).
\begin{table}[h!]
  \caption{Peak ratios of the 313~nm doublet of a Hg(Ar) calibration lamp measured by different sources. The same lamp was measured by method 1 and 2.}
  \begin{center}
    \begin{tabular}{ c c c | c c c | c c c | c c c }
    & Source & & & Ratio & & & Uncert.(\%) & & & $\lambda/\Delta\lambda$ & \\
    \hline
    1. & Hybrid device & & & 0.682 & & & 1.8 & & & $\sim$37,000 & \\
    2. & DEMON & & & 0.681 & & & 2.4 & & & $\sim$75,000 & \\
    3. & Ref.~\cite{Reader1996} & & & 0.687 & & & 7.2 & & & N/A & \\
    \hline
    \end{tabular}
  \end{center}
  \label{table:Ratios}
\end{table}

To further explore the utility of this spectral reconstruction method, it was applied to the analysis of LIBS data from the mineral cinnabar (HgS) acquired with a different hybrid spectrometer~\cite{Effenberger2012}. In this case, plasma was generated by 1064 nm laser pulses at 25 mJ per pulse. The spectrum was collected from the accumulation of 600 laser pulses at a rate of 10 Hz and under 10 Torr He atmosphere, which is the optimum pressure for the Hg 313~nm doublet~\cite{Effenberger2012}. The reconstructed doublet, shown in Fig.~\ref{figure:Reconstruct}(b), exhibits a ratio of 0.682 and is in good agreement with values reported in Table~\ref{table:Ratios}. Peak broadening in the LIBS spectrum is due to plasma dynamics~\cite{Scott2014}.

\section{Conclusion}

This spectral reconstruction approach for hybrid FP/dispersive spectrometers demonstrates that simple, low-cost instruments with high resolution and accurate relative peak intensities can be realized. Accurate wavelength and ratio determination of Hg lines using this instrument has been presented, where the accuracy depends on proper normalization and calibration of the measured interferogram. Because only a fraction of light is coupled to the CT spectrometer, a fraction of the signal is lost in this method. Summation of signal from different orders produced by the FP etalon may be able to increase the SNR, but was not implemented in the reconstruction discussed here. Experimentally, a higher repetition rate laser may be used to speed collection time and further increase the SNR. While such instruments would be invaluable for LIBS measurements requiring accurate isotope ratio measurements, the technique could enable practical instruments for a variety of high resolution field applications.\\

\newpage

The authors would like to thank Andrew J.~Effenberger, Jr.~and Elizabeth J.~Judge for their assistance. Research was performed under appointment to the Nuclear Nonproliferation International Safeguards Graduate Fellowship Program sponsored by the National Nuclear Security Administration's Next Generation Safeguards Initiative (NGSI). Material is based upon work supported by the U.S. Department of Homeland Security under Grant Award Number, 2012-DN-130-NF0001-02 and in part funded by the Consortium for Verification Technology under Department of Energy National Nuclear Security Administration award number DE-NA0002534. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Department of Homeland Security. Research was also sponsored by the U.S. Department of Energy under DOE Idaho Operations Office Contract DE-AC07-05ID14517.

\newpage


\section{Introduction}

\label{sec:1}

In 1930 Sergey L.\ Sobolev, one of the most distinguished mathematicians of
the 20$^{\text{st}}$ century, has published two papers \cite{Sob1,Sob2}
where he has constructed the solution of the hyperbolic equation of the
second order in the 3-d case with variable coefficients in the principle
part of the operator. It was assumed that these coefficient depend on
spatial variables. This result is readily available in the textbook of
Smirnov \cite{Smirnov}. Sobolev did not find the fundamental solution. The
main reason of this was that the notion of the fundamental solution was
unknown in 1930. Still, in his books about inverse problems \cite{R2,R3}
Romanov has modified the method of Sobolev to construct the fundamental
solution of that equation. However, the constructions of both Sobolev and
Romanov impose a quite restrictive assumption on the variable coefficients
in the principal part of the hyperbolic operator. Let the time variable
t\in \left( 0,T\right) .$ It is assumed in above cited publications that
geodesic lines generated by these coefficient are regular in a domain
Q\left( T\right) \subset \mathbb{R}^{3}.$ The larger $T$ is, the larger
Q\left( T\right) $ is. So, $Q\left( \infty \right) =\mathbb{R}^{3}.$ The
regularity of geodesic lines in the domain $Q\left( T\right) $ means that
for any two points $x,y\in Q\left( T\right) $ there exists a single geodesic
line connecting them.

Another construction of the fundamental solution of that equation can be
found in the book of Vainberg \cite{V}. The construction of \cite{V} imposes
the non-trapping condition on variable coefficients in the principal part of
the hyperbolic operator. In addition, the technique of \cite{V} relies on
the canonical Maslov operator, which is not easy to obtain explicitly. Thus,
even though the structure of the fundamental solution in any of above
constructions can be seen, formally at least, still many elements of this
structure cannot be expressed via explicit formulas. For example, neither
the solution of the eikonal equation of the method of Sobolev, nor the
Maslov construction cannot be expressed via explicit formulas.

The numerical factor is important nowadays.\ However, it is not immediately
clear how to compute numerically fundamental solutions obtained in above
references. Therefore, it is also unclear how to compute solutions of Cauchy
problems for heterogeneous hyperbolic equations if using those fundamental
solutions.

In this paper we propose a simple method of the construction of a function,
which is associated in a clear way with the fundamental solution for the
acoustic equation in the 3-d case with the variable coefficient. We impose
almost minimal assumptions on this coefficient. Geodesic lines are not used.
It is important that our function can be both accurately and effectively
approximated numerically via the Galerkin method as well as via the Finite
Difference Method. Thus, we show that this function can be straightforwardly
used for computations of the Cauchy problem for the heterogeneous acoustic
equation.

The author was prompted to work on this paper while he was working on
publications \cite{KR1,KR2} about reconstruction procedures for phaseless
inverse scattering problems. Indeed, in \cite{KR1,KR2} the structure of the
fundamental solution of the acoustic equation in time domain with the
variable coefficient in its principal part is substantially used.

\section{Construction}

Below $x\in \mathbb{R}^{3}$ and $\Omega \subset \mathbb{R}^{3}$ is a bounded
domain. Let $c_{0}$ and $c_{1}$ be two constants such that $0<c_{0}\leq
c_{1}.$ We assume that the function $c\left( x\right) $ satisfies the
following condition
\begin{equation}
c\in C^{1}\left( \mathbb{R}^{3}\right) ,c\in \left[ c_{0},c_{1}\right] ,
\label{2.1}
\end{equation
\begin{equation}
c\left( x\right) =1\text{ for }x\in \mathbb{R}^{3}\diagdown \Omega .
\label{2.3}
\end{equation}

Consider the following Cauchy proble
\begin{equation}
c\left( x\right) y_{tt}-\Delta _{x}y=f\left( x,t\right) ,\left( x,t\right)
\in \mathbb{R}^{3}\times \left[ 0,\infty \right) ,  \label{2.15}
\end{equation
\begin{equation}
y\mid _{t=0}=y_{t}\mid _{t=0}=0.  \label{2.16}
\end{equation
It is assumed here that $f$ is an appropriate function such that there is a
guarantee of the existence of the unique solution $y\in H^{2}\left( \mathbb{
}^{3}\times \left( 0,T\right) \right) ,\forall T>0$ of this problem, see,
e.g. Theorem 4.1 in \S 4 of Chapter 4 of the book of Ladyzhenskaya \cite{Lad}
for sufficient conditions for the latter. We will specify $f$ later. The
fundamental solution for the operator $c\left( x\right) \partial
_{t}^{2}-\Delta $ is such a function $P\left( x,\xi ,t,\tau \right) $ that
the solution of the problem (\ref{2.15}), (\ref{2.16}) can be represented in
the for
\begin{equation}
y\left( x,t\right) =\dint\limits_{0}^{\infty }\dint\limits_{\mathbb{R
^{3}}P\left( x,\xi ,t,\tau \right) f\left( \xi ,\tau \right) d\xi d\tau .
\label{2.30}
\end{equation
Below we modify formula (\ref{2.30}).

Let $\xi \in \mathbb{R}^{3}$ and $\tau \geq 0$ be parameters. Consider the
following Cauchy proble
\begin{equation}
c\left( x\right) u_{tt}=\Delta _{x}u+\delta \left( x-\xi \right) \delta
\left( t-\tau \right) ,  \label{2.4}
\end{equation
\begin{equation}
u\mid _{t=0}=u_{t}\mid _{t=0}=0.  \label{2.5}
\end{equation}

\subsection{Heuristic part of the construction}

It is convenient for us to work in this subsection with a purely heuristic
derivation. Represent the solution of the problem (\ref{2.4}), (\ref{2.5})
as $u=u_{0}+v$, where $u_{0}$ is the fundamental solution of the wave
equation
\begin{equation}
u_{0}=\frac{\delta \left( t-\tau -\left\vert x-\xi \right\vert \right) }
4\pi \left\vert x-\xi \right\vert }.  \label{2.6}
\end{equation
Thus
\begin{equation}
\partial _{t}^{2}u_{0}=\Delta _{x}u_{0}+\delta \left( x-\xi \right) \delta
\left( t-\tau \right) ,  \label{2.60}
\end{equation
\begin{equation}
u_{0}\mid _{t=0}=u_{0t}\mid _{t=0}=0.  \label{2.61}
\end{equation
Hence, the function $v$ satisfies the following condition
\begin{equation}
c\left( x\right) v_{tt}=\Delta _{x}v-\left( c\left( x\right) -1\right) \frac
\delta ^{\prime \prime }\left( t-\tau -\left\vert x-\xi \right\vert \right)
}{4\pi \left\vert x-\xi \right\vert },  \label{2.7}
\end{equation
\begin{equation}
v\mid _{t=0}=v_{t}\mid _{t=0}=0.  \label{2.8}
\end{equation}

Consider the operator $A$,
\begin{equation*}
A\left( f\right) =\dint\limits_{0}^{t}f\left( y\right) dy
\end{equation*
for appropriate functions $f$. Purely heuristically again apply the operator
$A^{4}$ to both sides of equation (\ref{2.7}).\ Denote $\widetilde{w}\left(
x,\xi ,t,\tau \right) =A^{4}\left( v\right) .$ Then (\ref{2.7}) and (\re
{2.8}) imply that
\begin{equation}
c\left( x\right) \widetilde{w}_{tt}-\Delta _{x}\widetilde{w}=-\left( c\left(
x\right) -1\right) \frac{\left( t-\tau -\left\vert x-\xi \right\vert \right)
}{4\pi \left\vert x-\xi \right\vert }H\left( t-\tau -\left\vert x-\xi
\right\vert \right) ,  \label{2.9}
\end{equation
\begin{equation}
\widetilde{w}\mid _{t=0}=\widetilde{w}_{t}\mid _{t=0}=0,  \label{2.10}
\end{equation
where $H\left( z\right) $ is the Heavyside function
\begin{equation*}
H\left( z\right) =\left\{
\begin{array}{c}
1,z>0, \\
0,z<0
\end{array
\right.
\end{equation*}

\subsection{Rigorous part of the construction}

Starting from this point, we are not acting heuristically anymore. To the
contrary, we work rigorously everywhere below.

Consider the Cauchy problem (\ref{2.9}), (\ref{2.10}). Denot
\begin{equation*}
g\left( x,\xi ,t,\tau \right) =-\left( c\left( x\right) -1\right) \frac
\left( t-\tau -\left\vert x-\xi \right\vert \right) }{4\pi \left\vert x-\xi
\right\vert }H\left( t-\tau -\left\vert x-\xi \right\vert \right) .
\end{equation*
It follows from (\ref{2.3}) that for any fixed pair $\left( \xi ,\tau
\right) \in \mathbb{R}^{3}\times \left[ 0,\infty \right) $ and for any $T>0$
functions $g,g_{t}\in L_{2}\left( \mathbb{R}^{3}\times \left( 0,T\right)
\right) .$ Hence, Theorem 4.1 of \S 4 of Chapter 4 of the book of
Ladyzhenskaya \cite{Lad}, (\ref{2.1}), (\ref{2.3}) as well as other results
of that chapter imply that for any fixed pair $\left( \xi ,\tau \right) \in
\mathbb{R}^{3}\times \left[ 0,\infty \right) $ and for any $T>0$ there
exists unique solution $\widetilde{w}\in H^{2}\left( \mathbb{R}^{3}\times
\left( 0,T\right) \right) $ of the Cauchy problem (\ref{2.9}), (\ref{2.10}).
Furthermore, it was shown in the proof of that theorem of \cite{Lad} that
this solution can be effectively constructed numerically via the Galerkin
method. It is also well known that it can be numerically constructed via the
Finite Difference Method. In addition, the energy estimate implies that
\begin{equation}
\widetilde{w}\left( x,\xi ,t,\tau \right) =0\text{ for }t\leq \tau .
\label{2.11}
\end{equation
Also, it follows from (\ref{2.9}) that $\widetilde{w}=\widetilde{w}\left(
x,\xi ,t-\tau \right) $

Consider now the function $w$ defined a
\begin{equation*}
w\left( x,\xi ,t-\tau \right) =\widetilde{w}\left( x,\xi ,t-\tau \right)
+A^{4}\left( u_{0}\right) .
\end{equation*
Hence
\begin{equation}
w\left( x,\xi ,t-\tau \right) =\widetilde{w}\left( x,\xi ,t-\tau \right)
\frac{\left( t-\tau -\left\vert x-\xi \right\vert \right) ^{3}}{6\cdot 4\pi
\left\vert x-\xi \right\vert }H\left( t-\tau -\left\vert x-\xi \right\vert
\right) .  \label{2.12}
\end{equation
Using (\ref{2.6})-(\ref{2.61}), (\ref{2.9}), (\ref{2.10}) and (\ref{2.12})
and applying direct calculations, we obtain that the function $w$ satisfies
the following condition
\begin{equation}
c\left( x\right) w_{tt}-\Delta _{x}w=\delta \left( x-\xi \right) \frac
\left( t-\tau \right) ^{3}}{6}H\left( t-\tau \right) ,  \label{2.13}
\end{equation
\begin{equation}
w\mid _{t=0}=w_{t}\mid _{t=0}=0.  \label{2.14}
\end{equation}

Consider now an arbitrary function $f\left( x,t\right) $ such that
\begin{equation*}
f\left( x,t\right) \in C^{4}\left( \mathbb{R}^{3}\times \left[ 0,\infty
\right) \right) ,
\end{equation*
\begin{equation*}
\partial _{t}^{k}f\left( x,0\right) =0,k=0,1,2,3,4,
\end{equation*
\begin{equation*}
f\left( x,t\right) =0,\forall x\in \mathbb{R}^{3}\diagdown G_{f},
\end{equation*
where $G_{f}$ is a bounded domain depending on the function $f$. Consider
the function $p_{f}\left( x,t\right) $ defined a
\begin{equation}
p_{f}\left( x,t\right) =\dint\limits_{0}^{t}\dint\limits_{G_{f}}w\left(
x,\xi ,t-\tau \right) f\left( \xi ,\tau \right) d\xi d\tau .  \label{2.17}
\end{equation
Using the integration by parts, (\ref{2.11}) and (\ref{2.12}), we obtain
that there exist four derivatives of the function $p_{f}\left( x,t\right) $
with respect to $t$. In particular,
\begin{equation}
\partial _{t}^{4}p_{f}\left( x,t\right)
=\dint\limits_{0}^{t}\dint\limits_{G_{f}}w\left( x,\xi ,t-\tau \right)
\partial _{\tau }^{4}f\left( \xi ,\tau \right) d\xi d\tau .  \label{2.18}
\end{equation
Next, applying the operator $c\left( x\right) \partial _{t}^{2}-\Delta $ to
both sides of (\ref{2.18}) and using (\ref{2.11}) and (\ref{2.13}), we obtai
\begin{equation}
\left( c\left( x\right) \partial _{t}^{2}-\Delta \right) \partial
_{t}^{4}p_{f}\left( x,t\right) =\dint\limits_{0}^{t}\frac{\left( t-\tau
\right) ^{3}}{6}\partial _{\tau }^{4}f\left( x,\tau \right) d\tau .
\label{2.19}
\end{equation
Using integration by parts in (\ref{2.19}) as well as (\ref{2.141}), we
obtai
\begin{equation}
\left( c\left( x\right) \partial _{t}^{2}-\Delta \right) \partial
_{t}^{4}p_{f}\left( x,t\right) =f\left( x,t\right) .  \label{2.20}
\end{equation
In addition, it follows from (\ref{2.11}), (\ref{2.12}) and (\ref{2.18})
that
\begin{equation*}
\partial _{t}^{4}p_{f}\mid _{t=0}=\partial _{t}^{5}p_{f}\mid _{t=0}=0.
\end{equation*
Since the function $\partial _{t}^{4}p_{f}\in H^{2}\left( \mathbb{R
^{3}\times \left( 0,T\right) \right) ,\forall T>0,$ then the uniqueness
theorem for the problem (\ref{2.15}), (\ref{2.16}) and (\ref{2.18}) imply
that
\begin{equation}
y\left( x,t\right) =\dint\limits_{0}^{t}\dint\limits_{G_{f}}w\left( x,\xi
,t-\tau \right) \partial _{\tau }^{4}f\left( \xi ,\tau \right) d\xi d\tau .
\label{2.200}
\end{equation
Thus, the solution of the problem (\ref{2.15}), (\ref{2.16}) is constructed
in (\ref{2.200}) using the function $w$. The formula (\ref{2.200}) is our
first analog of the formula (\ref{2.30}). The integration by parts in (\re
{2.200}) and the use of (\ref{2.11}) and (\ref{2.12}) leads to the second
analog of the formula (\ref{2.30}),
\begin{equation}
y\left( x,t\right) =\dint\limits_{0}^{t}\dint\limits_{G_{f}}\partial
_{t}^{2}w\left( x,\xi ,t-\tau \right) \partial _{\tau }^{2}f\left( \xi ,\tau
\right) d\xi d\tau .  \label{2.26}
\end{equation}

\subsection{Numerical comments}

Recall that the function $w$ can be accurately numerically approximated via
either the Galerkin method or the Finite Difference Method. The derivatives
of the function $f$ can be found analytically, if $f$ is given by an
explicit formula.\ However, if $f$ is given with a noise, then a
regularization method should be applied, see, e.g. the book of Tikhonov and
Arsenin \cite{TA}. In particular, it is explained in this book how to stably
differentiate noisy functions using the regularization.\ Furthermore,
numerical examples of stable computations of first and second derivatives
are presented in \cite{TA}. The form (\ref{2.26}) might be sometimes more
convenient than the form (\ref{2.200}) since second derivatives of noisy
functions are obviously more stable to calculate than fourth derivatives.

In summary, formulas (\ref{2.200}) and (\ref{2.26}) imply that the solution
y\left( x,t\right) $ of the problem (\ref{2.15}), (\ref{2.16}) can be
effectively numerically calculated via three steps:

\begin{enumerate}
\item \emph{Step 1}. Compute the function $w$ either via the Galerkin method
or via the Finite Difference method.

\item \emph{Step 2}. Compute corresponding $t-$derivatives of the function
f $ should be computed either analytically, if $f$ is given by an explicit
formula, or numerically using the regularization, if $f$ is given with a
noise.

\item \emph{Step 3}. Apply one of formulas (\ref{2.200}) or (\ref{2.26}).
\end{enumerate}

\begin{center}
\textbf{Acknowledgments}
\end{center}

This work was supported by the US Army Research Laboratory and US Army
Research Office grant W911NF-15-1-0233 as well as by the Office of Naval
Research grant N00014-15-1-2330.


\section{Motivations}
\label{sec1}

The more observational precision increases, the more inhomogeneous our Universe looks. Supernov\ae\ \cite{Conley:2011ku}, lensing \cite{schneider2006gravitational,Amendola:2012ys,2009arXiv0912.0201L}, and other light-propagated observables will soon encounter associated complications
\cite{Clarkson:2011br,Fleury:2013sna,Lavinto:2015iba,AppEt1933a,Buchert:2011sx,Buchert:2011yu,
Clarkson:2011zq,Kolb:2011zz,Coley:2005ei}.
We present coordinates first developed to simplify averages of scalars on our past light-cone \cite{P1}, next used to estimate the effect of inhomogeneities on luminosity distance \cite{P2,P4,Marozzi:2014kua,ErratumMarozzi2014}, Hubble diagram \cite{P3,P5}, and recently applied to lensing quantities \cite{P6} (and illustrated in a Lema\^itre-Tolman-Bondi model), following this historical order.

\vspace{-0.1cm}

\section{The geodesic light-cone coordinates}
\label{sec2}

We define a light-cone adapted metric (close to ``\emph{observational coordinates}''\cite{Maartens1,1985PhR...124..315E}, but different \cite{NugierThesis}) composed of 6 arbitrary functions ($\Upsilon$, $U^a$, $\gamma_{ab}$) and totally gauged fixed\,:
\begin{equation}
\mathrm{d} s_{\rm GLC}^2 = \Upsilon^2 \mathrm{d} w^2 - 2\Upsilon \mathrm{d} w \mathrm{d} \tau + \gamma_{ab} (\mathrm{d} \widetilde{\theta}^a - U^a \mathrm{d} w) (\mathrm{d} \widetilde{\theta}^b - U^b \mathrm{d} w) ~~.
\end{equation}

\vspace{-0.8cm}

\begin{figure}[ht!]
\centering
\includegraphics[width=5cm]{LCinhom.jpg} ~~
\includegraphics[width=5cm]{LChomo.jpg}
\caption{\label{Fig1} \emph{Left\,:} Inhomogeneous light-cone parametrized by GLC coordinates. \emph{Right\,:} Homogeneous light-cone or an adapted system (like GLC coordinates) in an inhomogeneous geometry.}
\end{figure}

\vspace{-0.4cm}

This metric uses a null coordinate $w$ defining past light cones, the proper time of a geodesic observer $\tau$, and angles $\tilde{\theta}^a$ that photons keep along their path orthogonal to a 2-spheres $\Sigma(w,\tau)$ of constant time in our past light-cone (see Fig. \ref{Fig1}). In the FLRW limit\,: $w = \eta + r$ (conformal time, radius), $\tau=t$ (cosmic time), $(\tilde{\theta}^1,\tilde{\theta}^2) = (\theta,\phi)$, $\Upsilon = a(t)$, $U^a = 0$, $\gamma_{ab} = a^2 r^2 (1,\sin^2\theta)$. In general \cite{NugierThesis} $\Upsilon$ is like an inhomogeneous scale factor (lapse function), $U^a$ like a shift-vector and $\gamma_{ab}$ is the metric inside $\Sigma(w,\tau)$. We can notice two direct simplifications in GLC which, combined together, give the distance-redshift relation\,:
\vspace{-0.3cm}
\begin{eqnarray}
\label{redshift}
& \mbox{Redshift\,:} ~~~ (1+z_s) = \Upsilon(w_o, \tau_o, \widetilde \theta^a) / \Upsilon(w_o, \tau_s, \widetilde \theta^a) ~~, \\
\label{AngDist}
& \mbox{Angular distance\,:} ~~~ d_A = \gamma^{1/4} (\sin\tilde{\theta}^1)^{-1/2} ~~\mbox{with}~ \gamma \equiv \det(\gamma_{ab}) = \frac{|\det (g_{\rm GLC})|}{\Upsilon^2} ~.~~
\end{eqnarray}
These coordinates share similarities with historical ones such as ``observational coordinates" \cite{saunders_observations_1968,saunders_observations_1969,Maartens1,1985PhR...124..315E}, (see elements of comparison in Ref.\cite{NugierThesis}) or the ``optical coordinates" \cite{1938RSPSA.168..122T}.

\section{Simplification of light-cone averages}
\label{sec3}

The light-cone average\cite{P1} of a scalar $S$ (e.g. $d_L ~,~ d_L^{-2}$) is in general given by $\left\langle S \right\rangle (V_0,A_0) = I(S;V_0,A_0) / I(1;V_0,A_0)$ and we define the average integral to be gauge invariant, invariant under $(A,V) \rightarrow (\tilde{A}(A),\tilde{V}(V))$ and general coordinate transformations\,:
\begin{equation}
I(S;V_0,A_0) = \int_{{\cal M}_4} \mathrm{d}^4x \sqrt{-g} ~ D(V_0-V) D(A-A_0) ~ {\mathcal N}(V,A,\partial_\mu) ~ S(x) ~~,
\end{equation}
where ${\mathcal N}(V,A,\partial_\mu)$ is a normalization, $D(X) = \delta_D(X)$ or $\Theta(X)$ (Heaviside function)\,:

\vspace{-0.5cm}

\begin{table}[ht!]
\begin{center}
\begin{tabular}{| c | m{3cm} | m{3cm} | m{3cm} |}
\hline
$\mbox{Average}_{~\delta_D}^{~\Theta}$ & \centering $\left\langle S \right\rangle_{V_0}^{A_0}$ & \centering $\left\langle S \right\rangle_{A_0}^{V_0}$ & ~~~~~ ~~ $\left\langle S \right\rangle_{V_0,A_0}$ \\
Illustration & ~~ \includegraphics[width=2.3cm]{Image1.jpg} & ~~ \includegraphics[width=2.3cm]{Image2.jpg} & ~~ \includegraphics[width=2.3cm]{Image3.jpg} \\
${\mathcal N}(V,A,\partial_\mu)$ & \centering $\frac{|\partial_\mu V \partial^\mu A|}{\sqrt{-\partial_\nu A \partial^\nu A}}$ & \centering $\sqrt{-\partial_\mu A \partial^\mu A}$ & ~~~~~ ~ $|\partial_\mu V \partial^\mu A|$ \\
\hline
\end{tabular}
\end{center}
\end{table}

\vspace{-0.5cm}

Among these 3 types of averages, $\left\langle S \right\rangle_{V_0,A_0}$ is closer to physical observables as it averages over the deformed 2-sphere embedded in the light-cone $V=V_0$ and a spatial hypersurface $A = A_0$. In GLC coordinates (where $V \rightarrow w$, $A \rightarrow \tau$) we can simplify the average and use Eq. \rref{redshift} to get (with $\tau_z \equiv \tau (z_s, w_o, \widetilde \theta^a))$\,:
\begin{equation}
\langle S \rangle_{w_o, z_s} =  \left( \int \mathrm{d}^2 \widetilde \theta \sqrt{\gamma (w_o, \tau_z, \widetilde \theta^b)} \, S(w_o, \tau_z, \widetilde \theta^b) \right) \Big / \left( \int \mathrm{d}^2 \widetilde \theta \sqrt{\gamma (w_o, \tau_z, \widetilde \theta^b)} \right) ~~,
\end{equation}
allowing us to average scalars on the sky, at a certain redshift.

\section{Distance-redshift relation at $\mathcal O(2)$}
\label{sec4}

The GLC metric enables the computation of $d_L(z)$ at $\mathcal O(2)$ in the Newtonian gauge (NG)\,:
\begin{equation}
\mathrm{d} s_{\rm NG}^2 = a^2(\eta) \left[ -(1+ 2 \Phi) \mathrm{d} \eta^2  + (1- 2 \Psi) ( \mathrm{d} r^2 + r^2 [\mathrm{d} \theta^2 + \sin^2 \theta \, \mathrm{d} \phi^2]) \right]
\end{equation}
with $\Phi = \psi + \frac12 \phi^{(2)}$ , $\Psi = \psi + \frac12 \psi^{(2)}$ the gauge invariant Bardeen potentials (no matter shear at $\mathcal O(1)$, or see Ref.\cite{Marozzi:2014kua}). Establishing the full transformation between GLC and NG coordinates at $\mathcal O(2)$ in (scalar) perturbations, we get $(\tau, w,\widetilde{\theta}^1,\widetilde{\theta}^2) = {\rm func}(\eta, r, \theta, \phi)$ and $(\Upsilon, U^a, \gamma^{ab}) = {\rm func}(\psi,\psi^{(2)},\phi^{(2)})$ which allows us to compute $d_L$ to $\mathcal O(2)$ (using Eq. \rref{AngDist})\,:
\begin{equation}
d_L(z_s, \theta^a) = d_{L}^{FLRW}(z_s) \left( 1 + \delta_S^{(1)}(z_s, \theta^a) + \delta_S^{(2)}(z_s, \theta^a) \right) ~~.
\end{equation}

The first order in scalar perturbations is given by\,:
\begin{equation}
\delta_S^{(1)}(z_s, \theta^a) = \left( 1 - \frac{1}{\mathcal H_s \Delta \eta} \right) J - \frac{Q_s}{\Delta \eta} - \psi_s^{(1)} - {\rm Lensing}^{(1)} ~,
\end{equation}
with $\mathcal H_s = a'(\eta_s)/a(\eta_s)$, $\Delta \eta = \eta_o - \eta_s^{(0)}$, containing (Integrated) Sachs-Wolfe ([I]SW), Doppler, and lensing (convergence) effects \,:
\begin{eqnarray}
& J  = ([\partial_+ Q]_s - [\partial_+ Q]_o) - ([\partial_r P]_s - [\partial_r P]_o) \sim {\rm SW} + {\rm ISW} + {\rm Doppler} ~~, \\
\medskip
& {\rm Lensing}^{(1)} = \frac12 \nabla_a \tilde{\theta}^{a (1)} = \int_{\eta_s^{(0)}}^{\eta_o} \frac{\mathrm{d} \eta}{\Delta\eta} ~ \frac{\eta - \eta_s^{(0)}}{\eta_o - \eta} \Delta_2 \psi (\eta, \eta_o-\eta, \theta^a) ~~, \nonumber
\end{eqnarray}
where `$o$' (`$s$') denote quantities evaluated at the observer (source), and we defined\,:
\begin{equation}
Q(\eta_+, \eta_-, \theta^a) = \int_{\eta_+}^{\eta_-} \mathrm{d} x~ \hat{\psi}(\eta_+,x,\theta^a) ~~,~~ P(\eta, r, \theta^a) = \int_{\eta_{in}}^\eta \mathrm{d} \eta'~ \frac{a(\eta')}{a(\eta)} \psi(\eta',r,\theta^a) ~~.
\end{equation}

Similarly obtained $\mathcal O(2)$ corrections contain (see Refs.\cite{P4,NugierThesis})\,:
\begin{itemize}
  \setlength{\itemsep}{-5pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
\item Two dominant terms\,: $(\mbox{Doppler})^2$, $(\mbox{Lensing})^2$, \\
\item Combinations of $\mathcal O(1)$-terms\,: $\psi_s^2$, $(\mbox{[I]SW})^2$, $\mbox{[I]SW} \times \mbox{Doppler}$, $(\psi_s, Q_s) \times $ (Lensing, [I]SW, Doppler), \\
\item Genuine $\mathcal O(2)$-terms\,: $\psi_s^{(2)}$, $Q_s^{(2)}$, $\mbox{Lensing}^{(2)} \! = \! \frac12 \nabla_a \tilde{\theta}^{a (2)}$, \\
\item New integrated effects\,: $\frac{1}{4 \Delta\eta} \int_{\eta_s^{(0)+}}^{\eta_s^{(0)-}} \mathrm{d} x~ \left[4 \hat{\psi} ~ \partial_+ Q + \hat{\gamma}_{0}^{ab} ~ \partial_a Q ~ \partial_b Q \right] (\eta_s^{(0)+},x,\theta^a)$, \\
\item Angle deformations\,: $(\gamma_0)_{ab} \partial_+ \tilde{\theta}^{a (1)}\partial_- \tilde{\theta}^{b (1)} ~~,~~ \partial_a \tilde{\theta}^{b (1)}\partial_b \tilde{\theta}^{a (1)}$ ~, \\
\item Redshift perturbations from Eq. \rref{redshift}, involving transverse peculiar velocities\,: \\ $\gamma_0^{ab} \left( \partial_a P ~ \partial_b P ~,~ \partial_a Q ~ \partial_b Q ~,~ \partial_a Q ~ \partial_b P \right)$ ,
$\partial_+ \int_{\eta_+}^{\eta_-} \mathrm{d} x  \left[ 4 \hat{\psi} ~ \partial_+ Q + \hat{\gamma}_0^{ab} \partial_a Q ~ \partial_b Q \right]$ ,
$\int_{\eta_{in}}^\eta \mathrm{d} \eta'~ \frac{a(\eta')}{a(\eta)} \partial_r \left[\left( \partial_r P \right)^2 + \gamma_0^{ab} \partial_a P ~ \partial_b P \right]$ , \\
\item Other important terms\,: Lens-Lens coupling, corrections to Born approximation.
\end{itemize}

These $\mathcal O(2)$ results were confirmed recently\cite{Fanizza:2015swa} working directly in terms of GLC coordinates rather than NG. An independent derivation \cite{2014CQGra..31t2001U1,2014CQGra..31t5001U2} lead to similar results, but a rigorous comparison with Ref.\cite{P4} is still lacking. Results were also given for vector/tensor perturbations\cite{P4} (Poisson gauge) and extended to the case with $\mathcal O(1)$ anisotropic stress \cite{Marozzi:2014kua,ErratumMarozzi2014}.

\section{Effects of large-scale structure on the Hubble diagram}
\label{sec5}

In Sec. \ref{sec4} we expressed $d_L$ in terms of $(\psi,\psi^{(2)},\phi^{(2)})$, hence it needs a description of the Bardeen potentials at $\mathcal O(1,2)$. We decompose the first order gravitational potential $\psi$ in Fourier modes and denote by $\overline{(...)}$ the ensemble (or stochastic) average over perturbations. The $\mathcal O(2)$ potentials can be related to $\psi$ by\cite{Bartolo:2005kv}\,: $\psi^{(2)}, \phi^{(2)} \propto \nabla^{-2} (\partial_i \psi \partial^i \psi) ~,~ \partial_i \psi \partial^i \psi$. Hence we can combine light-cone and stochastic averages in order to study the effect of statistical perturbations on the whole sky. For example, the (trivial) average of $\psi^2$ gives\,: $\overline{\left\langle \psi_s^2 \right\rangle} = \int_{0}^{\infty} \frac{\mathrm{d} k}{k} \mathcal P_{\psi}(k)$ where $\mathcal P_\psi (k) \equiv (k^3 / 2 \pi^2) |\psi_k(\eta)|^2 = (3 / 5)^2 A (k/k_0)^{n_s-1} T^2(k) \, g^2(z)$ is the power spectrum describing perturbations. At linear order, with $A$, $n_s$, $k_0$ taken from WMAP, $T(k)$ is a transfer function\cite{Eisenstein:1997ik} including a baryonic component (Silk damping), and $g(z)$ is the growth factor describing the recent time evolution of perturbations.
In CDM we get exactly the spectral coefficients coming from each correction of $d_L(z_s, \theta^a)$ described in Sec. \ref{sec4}\,:
$\overline{\left\langle d_L \right\rangle} = \int_0^\infty \frac{\mathrm{d} k}{k} \mathcal P_{\psi}(k) \, C(k \Delta\eta)$. We do the same in $\Lambda$CDM, with reasonable assumptions to simplify integrations \cite{P5}, and also with a non-linear power spectrum \cite{Smith:2002dz,Takahashi:2012em}.

Like $d_L$ one can also average the flux $\Phi = L / (4 \pi d_L^2) \simeq \Phi_0 + \Phi_1 + \Phi_2$. We get\,:
\begin{itemize}
  \setlength{\itemsep}{2pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
\item $\overline{\langle d_L^{-2} \rangle} \equiv (d_L^{FLRW})^{-2} \left[1 + f_\Phi(z) \right]$ where $f_\Phi(z) \simeq f(z) \int_0^{\infty} \frac{\mathrm{d} k}{k} \left(\frac{k}{  {\cal H}_0}\right)^2 \mathcal P_{\psi}(k)$ ~~, \\
\item $\overline{\langle d_L \rangle}(z) = d_L^{FLRW} \left[1 + f_d(z) \right]$ with $f_d = -(1/2) f_{\Phi} + (3/8) \overline{\langle \left(\Phi_1/\Phi_0\right)^2 \rangle}$~~.
\end{itemize}

\medskip

Corrections to $d_L$ involve a flux variance dominated by peculiar velocity and lensing\,:
\scriptsize
\begin{equation}
\overline{\left\langle \left(\phi_1 / \phi_0 \right)^2 \right\rangle}/4 = \overline{\left\langle (\delta_S^{(1)})^2 \right\rangle} \simeq \left( 1 - \frac{1}{\mathcal H_s \Delta\eta} \right)^2 \left\{ \left[\overline{\langle  ([\partial_r P]_s)^2 \rangle} + \overline{\langle  ([\partial_r P]_o)^2  \rangle}\right] + \overline{\left\langle \left({\rm Lensing}^{(1)}\right)^2 \right\rangle} \right\} ~~, \nonumber
\end{equation}
\small
as shown in Fig. \ref{Fig2} for realistic (non-)linear power spectra \cite{P5}. It turns out that the luminosity flux is minimally affected by lensing w.r.t. other scalar observables at large redshift. This calculation can be seen as a check at $\mathcal O(2)$ of Weinberg's argument of flux conservation \cite{1976ApJ...208L...1W} and has also been confirmed by recent papers through different approaches \cite{Kaiser:2015iia,2005ApJ...632..718K,2015JCAP...07..040B}.

\begin{figure}[ht!]
\centering
\includegraphics[width=5cm]{f3b.pdf} ~~~~
\includegraphics[width=5.1cm]{f6.pdf}
\caption{\label{Fig2} Corrections $f_\Phi(z)$ and $f_d(z)$ for a linear (\emph{Left}) and non-linear (\emph{Right}) power spectrum.}
\centering
\end{figure}

Similarly, we can get the average/dispersion of the distance modulus\,:
\begin{equation}
\overline{\langle \mu \rangle} = \mu^{\rm FLRW} - 1.25 (\log_{10}e) \Big[ 2f_\Phi - \overline{\langle \left(\Phi_1/\Phi_0\right)^2 \rangle} \Big] ~~,~~ \sigma_\mu = 2.5 (\log_{10} e) \sqrt{\overline{\langle \left(\Phi_1/\Phi_0\right)^2 \rangle}} ~~.
\end{equation}
Compared to the Union 2 data and using a non-linear power spectrum in $\Lambda$CDM (Fig. \ref{Fig3}, Left), we find that peculiar velocities explain well the scatter at small $z$ and that lensing explains only part of the scatter at large $z$. Finally, we can compare our dispersion on the Hubble diagram with the experimental estimations coming from lensing\cite{Jonsson:2010wx,Kronborg:2010uj} (Fig. \ref{Fig3}, Right). We find that the total effect is well fitted by Doppler ($z \leq 0.2$) + lensing ($z > 0.3$) effects and that the lensing prediction is in great agreement with experiments so far.

\begin{figure}[ht!]
\centering
\includegraphics[width=5cm]{Paper7DataAddingFinalPlusNLPS.pdf} ~~
\includegraphics[width=5.3cm]{SigmaInt.jpeg}
\caption{\label{Fig3} \emph{Left\,:} Theoretical dispersion from power spectrum vs Union 2 data. \emph{Right\,:} Dispersion $\sigma_{\mu}$ from peculiar velocities and lensing (solid line) compared to experimental estimates (grey areas).}
\centering
\end{figure}

\section{Jacobi map and weak lensing}
\label{sec6}

We now consider lensing\cite{P6}, motivated by Sec. \ref{sec5} and recent work\cite{DiDio:2014lka} on galaxy number counts in GLC. The relative separation of two neighbour light rays simultaneously emitted from a source $S$ and converging to an observer $O$ follows the geodesic deviation eq. (GDE)\,: $\nabla_\lambda^2 \xi^\mu = R^{\mu}_{~\,\alpha\beta\nu} k^\alpha k^\nu \xi^\beta$ with $\nabla_\lambda\equiv {\rm D}/\mathrm{d} \lambda \equiv k^\mu \nabla_{\mu}$, $k^\mu$ the photon momentum, $\lambda$ an affine parameter along the photon path, and $\xi^\mu$ an orthogonal displacement w.r.t. to the rays.
We project the GDE on the Sachs basis $\{ s^\mu_A \}_{A=1,2}$ (two zweibeins with flat index $A=1,2$)\,: $g_{\mu\nu} s^\mu_A s^\nu_B = \delta_{AB} ~,~ s^\mu_A u_\mu = 0 ~,~ s^\mu_A k_\mu = 0 ~,~ \Pi^\mu_\nu\nabla_\lambda s^\nu_A = 0$ ; with $u_\mu$ the peculiar velocity of comoving fluid ($S$, $O$ comoving too), $\Pi^\mu_\nu$ a ``screen" projector orthogonal to $u_\mu$ and $u_\mu+\left( u^\alpha k_\alpha \right)^{-1} k_{\mu}$.
We define the Jacobi map $J^A_B$, from the observed sky angle $\bar{\theta}_o^A$ to $\xi^A \equiv \xi^\mu s^A_\mu$, by $\xi^A(\lambda) = J^A_B(\lambda,\lambda_o) \, \bar{\theta}_o^A$. Projected quantities $\xi^A$ and $R^A_B \equiv R_{\alpha\beta\nu\mu}k^\alpha k^\nu s^\beta_B s^\mu_A$ (optical tidal matrix) bring us the Jacobi equation (see e.g. Refs.\cite{Fleury:2013sna,Fanizza:2013doa})\,:
\begin{eqnarray}
\label{JacobiEq}
& \frac{\mathrm{d}^2}{\mathrm{d} \lambda^2} J^A_B (\lambda , \lambda_o) = R^A_C (\lambda) \, J^C_B (\lambda, \lambda_o) ~~, \\
\label{JacobiEqIC}
& \mbox{with I.C.\,:} ~~~ J^A_B(\lambda_o, \lambda_o) = 0 ~~~,~~~ \frac{\mathrm{d}}{\mathrm{d} \lambda} J^A_B (\lambda_o,\lambda_o) = (k^\mu u_\mu)_o \, \delta^A_B ~~.
\end{eqnarray}
A direct resolution of Eq. \rref{JacobiEq} gives the angular distance of the source $d_A(\lambda_s) \equiv \sqrt{d S_s / d^2 \Omega_o} = \sqrt{\det J^A_B(\lambda_s,\lambda_o)}$.
Also, the (unlensed) angular position of the source $\bar\theta^A_s$ and the observed lensed position $\bar\theta^A_o$ (of the image) are given by\,: $\bar\theta^A_s=\left(\xi^A / \bar d_A\right)_s$, $\bar\theta^A_o=\left( k^\mu\partial_\mu\xi^A / k^\mu u_\mu \right)_o$, where $\bar{d}_{A}$ is the homogeneous and isotropic background our model refers to.
This allows us to define the so-called amplification matrix as\,:
\begin{equation}
\mathcal A^A_B \equiv \frac{\mathrm{d} \bar{\theta}^A_s}{\mathrm{d} \bar{\theta}^B_o} = \frac{J^A_B (\lambda_s,\lambda_o)}{\bar d _A(\lambda_s)} = \left( \begin{array}{cc} 1 - \kappa - \hat{\gamma}_1 & - \hat{\gamma}_2 + \hat{\omega} \\ - \hat{\gamma}_2 - \hat{\omega} & 1 - \kappa + \hat{\gamma}_1 \end{array} \right) ~~,
\end{equation}
which defines the lensing quantities\,: $\kappa = 1-\text{tr}J^A_B / 2\bar d_A$ (convergence), $\hat\omega=|J^1_2-J^2_1| / 2\bar d_A$ (vorticity), $|\hat\gamma| \equiv \sqrt{(\hat\gamma_1)^2 + (\hat\gamma_2)^2} = \sqrt{(1-\kappa)^2 + \hat{\omega}^2 - \mu^{-1}}$ (shear) and $\mu \equiv 1/(\det\mathcal A)=\bar d_A^2 / \det J^A_B$ (magnification).

Let us now turn to the GLC coordinates and express these lensing quantities in it. First, the zweibeins are written as $s^\mu_A=(s^\tau_A,0,s^a_A)$ and $k^\mu \equiv \omega \Upsilon^{-1}\delta^\mu_\tau$ (with $\omega$ a pure constant). Second, the solution to Eqs. \rref{JacobiEq} and \rref{JacobiEqIC} is\,:
\begin{equation}
J^A_B(\lambda,\lambda_o) = s_a^A(\lambda)\,\left[ 2 u_{\tau} (\dot{\gamma}_{ab})^{-1} \right]_os^B_b(\lambda_o)
\end{equation}
where $(\ldots)^{\mbox{\Large $\cdot$}} \equiv \partial_\tau (\ldots)$\,.
The angular distance and the magnification become\,:
\begin{equation}
d_A = 2 u_{\tau_o} (\gamma \gamma_o)^{1/4} \big / \sqrt{(\det \dot{\gamma}_{ab})_o} ~~~,~~~ \mu = \left(\bar d_A / d_A\right)^2 = \Phi / \bar{\Phi} ~~,
\end{equation}
involving $\bar{d}_A = a^2(\tau) r^2$ with $r = w - \int \mathrm{d} \tau / a(\tau)$ measured from the observer and $\Phi$ ($\bar{\Phi}$) the flux in the in(homogeneous) geometry. Expressions for the zweibeins can be obtained in the GLC coordinates\cite{P6}, but it is more convenient to compute the squared lensing quantities, combined with $s^A_as^A_b=\gamma_{ab}$ and $\epsilon_{AB}\,s^A_as^B_b=\sqrt{\gamma}\,\epsilon_{ab}$ ($\epsilon$ the anti-sym. symbol), to get\,:
\begin{equation}
\left( \begin{array}{c} \left( 1-\kappa \right)^2+\hat\omega^2 \\ \hat\gamma_1^2+\hat\gamma_2^2 \end{array} \right) = \left( \frac{u_{\tau_o}}{\bar d_A} \right)^2 \left( \left[ \frac{\gamma\,\dot\gamma_{ab}\gamma^{bc}\dot\gamma_{cd}}{\left(\det^{ab}\dot\gamma_{ab}\right)^2} \right]_o\gamma\,\gamma^{ad} \pm  \frac{2 \sqrt{\gamma\,\gamma_o}}{\left( \det^{ab}\dot\gamma_{ab} \right)_o} \right) ~~,
\end{equation}
We thus have general lensing quantities expressed with only 3 metric functions (of $\gamma_{ab}$), showing the great advantage of working in GLC coordinates.

The Jacobi Eq. \rref{JacobiEq} can be rewritten as a first order differential equation for the so-called deformation matrix\,:
\begin{equation}
\mathcal S^A_B\equiv \frac{\mathrm{d} J^A_C}{\mathrm{d} \lambda} (J^{-1})^C_B =\hat\theta\,\delta^A_B+\left(\begin{array}{cc}
\hat\sigma_1&\hat\sigma_2\\
\hat\sigma_2&-\hat\sigma_1
\end{array}\right) ~~,
\end{equation}
involving the optical scalars, $\hat\theta$ (expansion scalar) and $\hat\sigma \equiv \hat\sigma_1 + i \hat\sigma_2$ (shear scalar), and known as the Sachs equations\,:
\begin{equation}
\frac{\mathrm{d} \hat\theta}{\mathrm{d} \lambda}+|\hat\sigma|^2+\hat\theta^2 = \Phi_{00} ~~~~~,~~~~~ \frac{\mathrm{d} \hat\sigma}{\mathrm{d} \lambda}+2\hat\theta \hat\sigma = \Psi_0 ~~.
\end{equation}
The RHS terms are the Ricci and Weyl focusing and are defined as follows\,:
\begin{equation}
\Phi_{00}=-\frac{1}{2}R_{\alpha\beta} k^\alpha k^\beta ~~~,~~~ \Psi_0=\frac{1}{2}C_{\alpha\beta\mu\nu}k^\alpha k^\mu \Sigma^\beta \Sigma^\nu ~~,
\end{equation}
where $R_{\alpha\beta}$ is the Ricci tensor, $C_{\alpha\beta\mu\nu}$ the Weyl tensor, and $\Sigma^\mu \equiv s^\mu_1+is^\mu_2$.
As well as the amplification matrix, the deformation matrix simplifies in the GLC coordinates as\,:
\begin{equation}
\mathcal S^A_B = \frac{\mathrm{d} s^A_a}{\mathrm{d} \lambda} s^a_B = \frac{\omega}{2\Upsilon} s^a_A s^b_B\, \dot \gamma_{ab} ~~.
\end{equation}
Using $s^a_A s^b_A=\gamma^{ab}$ we get the optical scalars\,:
\begin{equation}
\hat\theta = \omega\,\frac{\gamma^{ab}\dot\gamma_{ab}}{4\Upsilon}=\frac{\omega}{4\Upsilon} \frac{\dot\gamma}{\gamma} ~~~~~,~~~~~ |\hat\sigma|^2 = \left(\frac{\omega}{4 \Upsilon} \frac{\dot\gamma}{\gamma}\right)^2 - \frac{\omega^2}{4\Upsilon^2}\,\frac{\det\dot\gamma_{ab}}{\gamma} ~~.
\end{equation}
We also get the Ricci and Weyl focusing in GLC coordinates\,:
\begin{equation}
\Phi_{00} = \frac{\omega^2}{4\Upsilon^2} \gamma^{ab} \, Y_{ab} ~~,~~ |\Psi_0| = \frac{\omega^2}{4\Upsilon^2} \sqrt{\left( \gamma^{ac}\gamma^{bd}+\gamma^{ad}\gamma^{bc}-\gamma^{ab}\gamma^{cd} \right) Y_{ab} Y_{cd}}~,
\end{equation}
where
$Y_{ab} \equiv \ddot \gamma_{ab} - (\dot \Upsilon / \Upsilon) \dot\gamma_{ab} - (1/2) \dot\gamma_{ac}\gamma^{cd}\dot\gamma_{db}$ depends only on $\gamma_{ab}$, $\Upsilon$ and their time derivatives. This proves again the usefulness of GLC coordinates for lensing and it was illustrated\cite{P6} by the computation of lensing quantities in the case of an off-center observer in a Lema\^itre-Tolman-Bondi model (considering only the decaying mode).

\section{Conclusions}
\label{sec7}

We have shown that there are many advantages in using the GLC coordinates. They are indeed adapted to calculations involving light-propagation, they can also be used for weak lensing (where $\gamma_{ab}$ acts as a screen), and may help to get new predictions on cosmology or to study other aspects of lensing (e.g. lensing statistics).

\section*{Acknowledgments}
I want to thank the \textsc{Fourteenth Marcel Grossmann Meeting} for giving me the opportunity to share these properties of GLC coordinates in two parallel sessions and for giving me the occasion to write this document. My researches are supported by the project GLENCO lead by B. Metcalf, funded under the FP7, Ideas, Grant Agreement n. 259349.

\ssmall

\providecommand{\href}[2]{#2}\begingroup\raggedright