text:
\section{\label{sec:introduction}Introduction}

The physics of graphene is unusual in that its electrons can be adequately modelled as relativistic massless Dirac fermions, which admit a linear energy dispersion  -- the famous Dirac cones. This property itself is known to induce highly nonlinear dynamics for light \cite{Mikhailov2007}. Since this electronic dispersion is ungapped, with the bands extrema touching at the Dirac points (termed $\mathbf{K}$ and $\mathbf{K'}$), graphene behaves like a zero-gap semiconductor. \\
However, this property is only expected for free-standing, pristine graphene samples. More physically realisable samples are normally deposited on particular dielectric substrates. These intrinsic factors are known to modify the electronic and optical properties of the sample and can be successfully taken into account by simply opening a gap in the two-band spectrum \cite{Zhou2007}.
Using various synthesis and preparation techniques, impurities, local lattice defects and vacancies \cite{Kang2008}, and strain effects \cite{Ni2008} may be physically realised and have also been shown to gap the spectrum. More challenging procedures to achieve this rely on electric biasing of graphene bilayers \cite{Castro2007} and monolayer nanostructuring into nanoribbons \cite{Li2008}. The appearance of a gap can also be conceptualised with a staggered sublattice potential, in which each triangular sublattice of the honeycomb lattice admits opposite non-zero on-site potentials (for instance when graphene is deposited on hexagonal boron nitride). \\

Each process admits a characteristic gap scale. Substrate-induced effects seem to be the most efficient to open a gap which, with the aid of ARPES measurement techniques, has been estimated to be $0.26$ eV for epitaxially-grown graphene on silicon carbide (SiC) \cite{Zhou2007}.
Density Functional Theory calculations estimate monolayer graphene can acquire a gap of $0.35$ eV when deposited on a SiO$_2$ substrate \cite{Shemella2009}.
Note that the extent of such a gap opening is linked to the relative geometrical configurations of the substrate and the sample alongside the dominant chemical bonds in their interaction. For instance, graphene deposited on Si-terminated silica surface with inactive
dangling bonds has been proposed as a configuration to retrieve the linear, gapless dispersion typical of free-stranding graphene \cite{Kang2008bis}.
The transition to a semiconducting regime leads to substantially different optoelectronic features for which devices such as graphene-based transistors and photodetectors rely on \cite{Schwierz2010}. The optical behaviour of the plane-confined carriers is further modified by excitonic effects, in turn caused by screening mechanisms. These may be appreciated through theoretical models of the optical conductivity spectra and phenomenological dependence on disorder and imperfections in Ref. \cite{Pedersen2009}.

Although advancing, the theoretical understanding of these effects on the \emph{ultrafast nonlinear} optical properties of graphene remain broadly uncharted. In this paper, we investigate the role of the energy gap in the ultrafast generation of high harmonic radiation along with related nonlinear processes, within a semiclassical quasirelativistic formalism, by explicitly solving the Dirac equation modelling the carrier dynamics.

As a centrosymmetric material, graphene should not allow the generation of even harmonics. However, intense and ultrashort pulses provide a regime where odd harmonics interfere generating even harmonics, once gapped. Besides, the dynamical centrosymmetry breaking mechanism, a field-driven effect that globally displaces the electronic dispersion and breaks the centrosymmetry  $\mathbf{k} \leftrightarrow - \mathbf{k}$ (and consequently $\mathbf{r} \leftrightarrow - \mathbf{r}$ in direct space), has been shown to predict second-harmonic generation in freestanding, ungapped graphene \cite{Carvalho2016}.
Gapping the spectrum renders the electrons massive and, as will be demonstrated, induces a Berry phase in the carrier dynamics. Such a phase leads to interesting qualitatively different optical behaviour which is studied by extending the Dirac-Bloch equations and their framework, previously applied to massless Dirac fermions (namely monolayer graphene) \cite{Ishikawa2013, Ishikawa2010, Carvalho2016}, to incorporate a gapped Dirac spectrum.

\section{\label{sec:quasirelativisticdynamics}Quasirelativitic Dynamics}

Graphene is a two-dimensional crystal composed of carbon atoms and disposed in a honeycomb lattice. This arrangement stems from particular orbital hybridisation and strong covalent in-plane bonding. Linearisation of the electronic dispersion computed from tight-binding methods yields a linear dependence, which vanishes at two non-equivalent Dirac points in momentum space termed $\mathbf{K}$ and $\mathbf{K'}$. Such a linear dispersion admits a conduction and valence bands which are symmetric and touch at the Dirac points, rendering the monolayer a zero-gap system. However, a gap can be opened at the Dirac points, which are located on the edge of the Brillouin zone and shown in Fig.~\ref{fig1}(a). For low-momentum states around these points, the dispersion attains its extrema (such regions are called valleys) and carriers can be endowed with an effective mass.

Unless an imbalance is physically realised, for instance through an electric field bias or sample inhomogeneities, intervalley scattering is highly unlikely\cite{GrapheneCarbonTwoDimensions}, as it requires exceedingly large phonon momenta, roughly of the order of the separation $|\mathbf{K}-\mathbf{K'}|$ . It can thus be reasonably assumed that the dynamics of both valleys is decoupled of each other and the carriers can be endowed with an additional degree of freedom, the \emph{valley isospin} $\xi$, where $\xi = +1$ $(-1)$ refers to states in the $\mathbf{K} \ (\mathbf{K'})$ valley. A further degree of freedom, the \emph{pseudospin} $\lambda$, with $\lambda = +1$ $(-1)$ denoting conduction (valence) band states, distinguishes between electron and hole states. \\
\noindent In order to understand light-matter interactions in this gapped structure, we proceed by obtaining the wavefunction of an electron of effective mass $m \equiv \Delta/(2 v_{\rm F}^2)$ and momentum $\mathbf{p} = \hbar \mathbf{k}$ in the vicinity of a particular Dirac point in valley $\xi$, which must obey a two-dimensional Dirac equation:
\eq{dirac}{
i \hbar \partial_{t} \ket{\Psi^{
\xi}_{\mathbf{k}}(t)} = H^{\xi}_{\mathbf{k}}(t)  \ket{\Psi^{\xi}_{\mathbf{k}}(t)}.}
$\Delta$ is the energy dispersion gap and $v_{\rm F} \approx c/300$ the electronic Fermi velocity. To obtain the appropriate Hamiltonian for such interactions, the canonical momentum is introduced through the minimal substitution $\mathbf{p} \mapsto \mathbf{p} + (e/c) \mathbf{A}(t) \equiv \bm{\pi}_{\mathbf{k}}(t)$ in the field-free Hamiltonian.  yielding: \begin{equation}
\label{hamiltonian}
H^{\xi}_{\mathbf{k}}(t) = v_{\rm F}~\left(\bm{\sigma}(\xi) \cdot \left(\mathbf{p} + \frac{e}{c} \mathbf{A}(t) \right) \right) + \frac{\Delta}{2} \sigma_{z},\end{equation}
where $\bm{\sigma}(\xi) \equiv (\xi \sigma_{x},  \sigma_{y})$ is a vector comprised of the 2D Pauli matrices, $e > 0$ is the absolute value of the electron charge and $c$ is the speed of light in vacuum. \\
\noindent The pulse is further assumed to be normally-incident and linearly polarised along an arbitrary direction, here taken along $\hat{x}$. Its electromagnetic vector potential $\mathbf{A}$, which is chosen to satisfy the Coulomb gauge $\nabla \cdot \mathbf{A} = 0$, can thus be written as $\mathbf{A}(t) = (A(t),0,0)$. Consequently, the canonical momentum becomes $\bm{\pi}_{\mathbf{k}}(t) = (p_{x} + (e/c)A(t), p_y)$.  \\

\begin{figure}
\includegraphics[width=8cm]{Fig1.png}
\caption{\label{fig1}(colour online) (a) Sketch of the Hamiltonian spectra for both valleys in the low-momentum regime. Each valley admits two bands, gapped by $\Delta$. The relative sign of the field-induced Berry phase is represented by the silver arrows. (b) Depiction of the time-dependent electronic dispersion in momentum space, as given in Eq.~(\ref{dispersion}), for a particular valley. Note that the pulse shifts the dispersion globally by the time-dependent photon momentum $\mathbf{A}(t)$. This field-driven effect is only appreciable for ultrashort and intense pulses.}
\end{figure}

\noindent Inconveniently, general analytical solutions of Eq.~(\ref{dirac}) cannot be obtained due to the time dependence of the Hamiltonian through the external parameter $\mathbf{A}(t)$. To tackle this, a general \emph{ansatz} is constructed through expansion over a basis comprised of the so-called \emph{instantaneous eigenstates} -- two linearly-independent spinors which satisfy $H^{\xi}_{\mathbf{k}}(t)\kkett{u} = \qt{\epsilon} \kkett{u}$. If orthonormalised i.e. $\av{{u^{\xi}_{\lambda}}_{\mathbf{k}}(t)|{u^{\xi}_{\lambda'}}_{\mathbf{k}}(t)} = \delta_{\lambda \lambda'}$, they take the form:
\eq{kets}{
\kkett{u} = \frac{v_{\rm F} |\bm{\pi}_{\mathbf{k}}|}{\sqrt{ \epsilon_{\mathbf{k}} (\lambda \Delta + 2 \epsilon_{\mathbf{k}})}} \begin{pmatrix}
\left(\frac{\lambda \Delta + 2 \epsilon_{\mathbf{k}}}{2 \xi v_{\rm F}|\bm{\pi}_{\mathbf{k}}|} \right)e^{-i \xi \theta_{\mathbf{k}}/2} \\
\lambda e^{i \xi \theta_{\mathbf{k}}/2}
\end{pmatrix},}
where $ \theta_{\mathbf{k}}(t) = \arctan(p_{y}/[p_{x} + \frac{e}{c}A(t)])$ is the \emph{dynamical angle} of the canonical momentum vector and $\epsilon_{\mathbf{k}}(t)$ the positive branch of their \emph{instantaneous} energy:
\eq{dispersion}{ \qt{\epsilon} = \lambda \sqrt{ \left(\frac{\Delta}{2}\right)^2 + ( v_{\rm F} |\bm{\pi}_{\mathbf{k}}(t)|)^2 } \equiv \lambda \epsilon_{\mathbf{k}}(t).}
\noindent As seen in Fig.~\ref{fig1}(a), the spectra of both valley Hamiltonians are \emph{globally} equivalent. These solutions have a straightforward interpretation: for a particular valley $\xi$, electron and hole states exist  respectively in the conduction ($\lambda = +1$) and valence ($\lambda = -1$) bands, which are gapped by $\Delta$. The upper and lower components of the spinor can be construed as amplitudes in each of the triangular sublattices that decompose the honeycomb lattice. \\
\noindent The addition of the gap leads to an inequivalence of these sublattices and consequently the appearance of a Berry phase, whose derivative can be obtained through the instantaneous eigenstates of Eq.~(\ref{kets}) as $ \dot{\gamma}^{\xi}_{\lambda \mathbf{k}}(t) = i \av{\qt{u}|\partial_{t}|\qt{u}} = \xi \lambda \left( \Delta \dot{\theta}_{\mathbf{k}}(t)/(4 \epsilon_{\mathbf{k}}(t)) \right) \equiv \xi \lambda \dot{\gamma}_{\mathbf{k}}(t)$.
The associated wavefunctions of these four states, solutions of the Dirac equation, must evolve in time as:
\eq{bandwavefunction}{
\kkett{\psi} = \kkett{u}e^{-i \lambda \Omega_{\mathbf{k}}(t)}e^{i \xi \lambda \gamma_{\mathbf{k}}(t)}}
and are therefore further phase-shifted by the \emph{dynamical} phase $\Omega_{\mathbf{k}}(t) =(1/\hbar) \int_{-\infty}^{t} \epsilon_{\mathbf{k}}(t')dt'$ and the \emph{Berry} phase $\gamma_{\mathbf{k}}(t) = \int_{-\infty}^{t}\dot{\gamma}_{\mathbf{k}}(t')dt'$, the latter taking the analytical form $\gamma_{\mathbf{k}}(t) = \left(\Lambda_{\mathbf{k}}(t) - \Lambda_{\mathbf{k}}(-\infty) \right)/4$, with:
\begin{equation}
\label{Berry}
\Lambda_{\mathbf{k}}(t) \equiv \arctan \left[ \frac{4  \Delta \epsilon_{\mathbf{k}}(t) \tan \theta_{\mathbf{k}}(t)}{ \Delta^2 - 4 \epsilon^2_{\mathbf{k}}(t) \tan^2 \theta_{\mathbf{k}}(t) } \right],
\end{equation}
where the time-independent term $\Lambda_{\mathbf{k}}(-\infty) \equiv \gamma^{0}_{\mathbf{k}}$ is a global phase factor arbitrarily included to ensure the Berry phase vanishes in the absence of field interaction. The time evolution of the Berry phase for particular states may be appreciated in Fig.~\ref{fig4}.
Finally, the {\em ansatz} $\ket{\Psi^{\xi}_{\mathbf{k}}(t)}$ of Eq.~(\ref{dirac}) is taken through expansion over the band wavefunctions at that valley, i.e.
\begin{equation}
\label{ansatz}
{\ket{\Psi^{\xi}_{\mathbf{k}}(t)} = c^{\xi}_{+1}(t)\ket{\psi^{\xi}_{+1, \mathbf{k}}(t)}} + c^{\xi}_{-1}(t)\ket{\psi^{\xi}_{-1, \mathbf{k}}(t)}.
\end{equation}
A word of caution is in order: although the gapped field-free Hamiltonian of Eq.~(\ref{hamiltonian}) (i.e. with $A(t)=0$)  arises from the breaking of the sublattice inversion symmetry, such Hamiltonian is only a first-order $\mathbf{k} \cdot \mathbf{p}$ approximation of the full tight-binding Hamiltonian and accounts only for its centrosymmetric part. Several works, e.g. in Refs. \cite{Xiao2012,Yao2014}, use our Hamiltonian to model transition metal dichalcogenide (TMD) monolayers, promising two-dimensional relativistic-like semiconductors lacking an inversion centre and hence non-centrosymmetric. Such an approximation for TMDs is adequate only to describe the linear optical properties of such media that are accounted by low momentum states, where this approximation is accurate. It is nonetheless clearly insufficient to accurately capture nonlinear light-matter phenomena of non-centrosymmetric two-dimensional media, for which higher-order terms in the $\mathbf{k} \cdot \mathbf{p}$ expansion explicitly break the centrosymmetry  $\mathbf{k} \leftrightarrow -\mathbf{k}$ (rendering the conduction and valence bands of Fig.~\ref{fig1}(a) asymmetric).

\section{\label{sec:massiveDBEs} The Massive Dirac-Bloch Equations}

The electron dynamics can be more easily understood by obtaining the time derivatives of $c^{\xi}_{\lambda}$ in and introducing new dynamical variables: the `population inversion' $w^{\xi}_{\mathbf{k}}\equiv|c^{\xi}_{+}|^{2}-|c^{\xi}_{-}|^{2}$ and the `microscopic polarisation' $q^{\xi}_{\mathbf{k}} \equiv c^{\xi}_{+} (c^{\xi}_{-})^{*} e^{-i(2 \Omega_{\mathbf{k}} - \omega_{0}t)}$. The full Dirac equation (\ref{dirac}) can be recast in a more transparent set of equations, akin to the Bloch equations of a two-level system -- \emph{the Dirac-Bloch equations} (DBEs) -- which are derived and shown for the case of massless Dirac fermions in Ref. \cite{Carvalho2016}. When generalised to the massive case, they take the form:
\begin{widetext}
\begin{eqnarray}
\label{db1}
\dot{w}^{\xi}_{\mathbf{k}} + \gamma_{1} (w^{\xi}_{\mathbf{k}} - {w_{\mathbf{k}}}_{0}) -  \left( \frac{v_{\rm F} |\bm{\pi}_{\mathbf{k}}|} {\epsilon_{\mathbf{k}}}\right) \left(2 \xi \dot{\theta}_{\mathbf{k}} \text{Im}(q^{\xi}_{\mathbf{k}} e^{i( 2 \xi \gamma_{\mathbf{k}} -  \omega_{0} t)}) + 4 \cot \theta_{\mathbf{k}} \dot{\gamma}_{\mathbf{k}}\text{Re}(q^{\xi}_{\mathbf{k}} e^{i( 2 \xi \gamma_{\mathbf{k}} -  \omega_{0} t)})\right)  = 0 \\
\label{db2}
\dot{q}^{\xi}_{\mathbf{k}} + i \left(2 \dot{\Omega}_{\mathbf{k}} - \omega_{0} -  i \gamma_{2} \right) q^{\xi}_{\mathbf{k}}  + \left( \frac{ v_{\rm F} |\bm{\pi}_{\mathbf{k}}|}{\epsilon_{\mathbf{k}}} \right) \left(\cot \theta_{\mathbf{k}} \dot{\gamma}_{\mathbf{k}} + \frac{i \xi \dot{\theta}_{\mathbf{k}}}{2} \right) e^{-i( 2 \xi \gamma_{\mathbf{k}} -  \omega_{0} t)}   w^{\xi}_{\mathbf{k}}= 0.
\end{eqnarray}
\end{widetext}
Here, $\dot{\Omega}_{\mathbf{k}}(t) = \epsilon_{\mathbf{k}}(t)/\hbar$ and $\dot{\theta}_{\mathbf{k}}(t) \equiv ep_{y}E(t)/|\bm{\pi}_{\mathbf{k}}(t)|^2$. The constants $\gamma_{1 (2)} \equiv 1/T_{1 (2)}$ are phenomenological decay rates of the population inversion (microscopic polarisation) whereas ${w_{\mathbf{k}}}_{0}$ is the equilibrium value of inversion; if the system is undoped, i.e. $\mu = 0$, and at temperature $T = 0$, it has ${w_{\mathbf{k}}}_{0} = -1$, implying that all carriers are initially found in the valence band, regardless of their momentum. \\
Otherwise, for arbitrary doping and temperature, it becomes ${w_{\mathbf{k}}}_{0} = - \sinh(y)/[\cosh(x) + \cosh(y)]$, with $y = \epsilon_{\mathbf{k}}/(k_{\rm B}T)$ and $x = - \mu/(k_{\rm B}T)$.
These two newly-defined fields modelled by Eqs.~\ref{db1}-\ref{db2} depend on a particular valley but are nonetheless connected by precise relations. The real-valued inversions are equal, i.e.  $ w^{\xi}_{\mathbf{k}} = w^{-\xi}_{\mathbf{k}}$, while the complex-valued microscopic polarisations are statically shifted by the momentum-dependent phase $\gamma^{0}_{\mathbf{k}}$, i.e. $ q^{\xi}_{\mathbf{k}} = e^{i \xi \gamma^{0}_{\mathbf{k}}} q^{-\xi}_{\mathbf{k}}$. In the limiting case of a vanishing gap, they satisfy $ q^{\xi}_{\mathbf{k}} = - q^{-\xi}_{\mathbf{k}}$. \\
The massive DBEs [Eqs.~(\ref{db1}-\ref{db2})] contain terms not present in their massless counterparts as derived and shown in Ref. \cite{Carvalho2016}.
In the driving term in Eq.~(\ref{db2}) (the one containing $w^{\xi}_{\mathbf{k}}$), the quantity that multiplies the electric field $E(t)$ may be identified as a valley and time dependent {\em complex-valued} electric dipole moment:
\begin{equation}
\mu^{\xi}_{\mathbf{k}}(t) = e v_{\rm F} \left(\frac{ \xi \sin \theta_{\mathbf{k}}(t)}{2 \epsilon_{\mathbf{k}}(t)} - i \frac{\Delta \cos \theta_{\mathbf{k}}(t)}{4 \epsilon^{2}_{\mathbf{k}}(t)} \right)
\end{equation}

We remark that the singularity found in the DBEs when $\theta_{\mathbf{k}}(t) = 0$ is not problematic since both equations can be identically re-expressed so that no real singularities are present. \\
Finally, we emphasise that the Coulomb interactions amongst the carriers are not included in the massive DBEs. These are known to lead to Fermi velocity and energy band renormalisation \cite{QTheorySemiconductors}. Such effects can in principle be included by coupling the dynamics of two-level systems of all momenta and have been previously implemented for graphene, see for instance Ref. \cite{Malic2011}.

\section{\label{sec:currentanalytics}Current Analytics}

Signatures of nonlinear light-matter interactions can be found and analysed through the electric current generated by the interaction between the monolayer and the pulse. Such a current admits, in general, two components: $\mathbf{J}(t) = (\mathbf{J}_x (t), \mathbf{J}_y (t))^\text{T}$. We investigate the role of the gap (and consequently the Berry phase) in the valley-dependent current contributions. To attain this, we proceed by first determining the $\mu$-component ($\mu = x,y$) of the current contribution of a particular momentum state $\mathbf{p}$ in a valley $\xi$ in time domain, here termed a \emph{microscopic current} $j^{\xi}_{\mu, \mathbf{k}}$, by applying the current density operator $\hat{j}^{\xi}_{\mu, \mathbf{k}}$ to the {\em ansatz} $\ket{\Psi^{\xi}_{\mathbf{k}}}$ of Eq. (\ref{ansatz}):
\eq{currentsandwich}{
 j^{\xi}_{\mu, \mathbf{k}} = \av{\Psi^{\xi}_{\mathbf{k}}|\hat{j}^{\xi}_{\mu, \mathbf{k}}|\Psi^{\xi}_{\mathbf{k}}} - \av{\psi^{\xi}_{-1, \mathbf{k}}|\hat{j}_{\mu, \mathbf{k}}|\psi^{\xi}_{-1, \mathbf{k}}}.}
Energy bands obtained with tight-binding methods must satisfy a sum rule that prevents currents in the valence bands to be produced \cite{Sipe1993}. However, since the dispersion of Eq.~(\ref{dispersion}) and spinors of Eq.~(\ref{kets}) are only applicable over a particular, low-momentum range where these are relativistic, the first current term in Eq.~(\ref{currentsandwich}) is insufficient to describe the actual current generated, as it contains unphysical divergences. The current can nonetheless be regularised through the introduction of the second term, which acts as an ad-hoc subtraction of valence band generated current. \\
By using the definition $\hat{j}^{\xi}_{\mu, \mathbf{k}} \equiv -(e/\hbar) (\partial H^{\xi}_{\mathbf{k}}/ \partial k_{\mu})$, the valley-dependent current density operator is obtained for each component as $\hat{j}^{\xi}_{x, \mathbf{k}} = -(\xi e v_{\rm F}/\hbar) \sigma_{x}$ and $\hat{j}^{\xi}_{y, \mathbf{k}} = - (e v_{\rm F}/ \hbar) \sigma_{y}$. With these, the contribution of both components to the 2D microscopic current $\mathbf{j}^{\xi}_{\mathbf{k}}(t) \equiv (j^{\xi}_{x, \mathbf{k}}(t), j^{\xi}_{y, \mathbf{k}}(t))^{\text{T}}$ as shown in Eq.~(\ref{currentsandwich}) is computed exactly as:

\begin{figure*}
\includegraphics[scale=0.43]{Fig_CurrentTime.png}
\caption{\label{fig2}(colour online) The full photo-generated current and its separate contributions in time domain, rescaled in units of $J_{0} = -e \omega^{2}_{0}/(4 d v_{\rm F})$. (a) The total current, composed of both intraband and interband contributions. Its overall dependence on the mass stems primarily from the interband contribution. (b) The intraband current, generated from electronic transitions within the same band. Its amplitude is monotonically decreasing as the gap increases and maximal when the dispersion is ungapped. (c) The interband contribution, generated from electronic transitions across the bands. It is comprised of two terms, one being exclusively present only for gapped dispersions. The interband current amplitude is maximal when the photon energy is resonant with the gap, rapidly decreasing for larger gaps.}
\end{figure*}

\begin{figure*}
\includegraphics[scale=0.43]{Fig_Spectra.png}
\caption{\label{fig3}(colour online) Current spectra. (a) The total spectrum shows stronger odd harmonic generation than even harmonic generation, which can only be induced by centrosymmetry breaking mechanisms. For gapped systems, even-harmonic peaks, which are plotted in (c), are generated through third harmonic generation (THG) in disguise of higher harmonic generation (HHG). Such peaks are shown with the second and fourth-harmonic enhancements, respectively for gaps satisfying $\tilde{\Delta} = 2,4$. For vanishing gaps, even harmonic generation originates from the centrosymmetry breaking mechanism, which breaks the static centrosymmetry of the lattice and is seen in the intraband spectrum of (b). Generally, the intraband harmonic peaks decrease monotonically as the gap is increased.}
\end{figure*}

\begin{widetext}
\begin{equation}
\label{microcurrent}
\mathbf{j}^{\xi}_{\mathbf{k}}(t) = -e v_{\rm F}
\begin{pmatrix}
\cos\theta_{\mathbf{k}} & \sin\theta_{\mathbf{k}} \\
\sin\theta_{\mathbf{k}} &  -\cos\theta_{\mathbf{k}}
\end{pmatrix}\begin{pmatrix}
\frac{v_{\rm F}|\bm{\pi}_{\mathbf{k}}|}{\epsilon_{\mathbf{k}}}(w^{\xi}_{\mathbf{k}} +1) - \frac{\Delta}{\epsilon_{\mathbf{k}}} \text{Re} \left(q^{\xi}_{\mathbf{k}}e^{i( 2 \xi \gamma_{\mathbf{k}} - \omega_{0}t)} \right )\\
-2 \xi \text{Im} \left(q^{\xi}_{\mathbf{k}}e^{i( 2 \xi \gamma_{\mathbf{k}} - \omega_{0}t)} \right )
\end{pmatrix}
\end{equation}
\end{widetext}
The physical current is finally obtained by appropriately taking all momentum contributions of both valleys into account. In the continuum limit, it is:
 \begin{equation}
\label{macrocurrent}
\mathbf{J}(t) = \frac{g_{s}}{d(2 \pi)^2} \sum_{\xi} {\int \mathbf{j}^{\xi}_{\mathbf{k}}(t)}d \mathbf{k},
\end{equation}
where $d$ is the thickness of the monolayer and $g_{s}=2$ is a spin degeneracy factor, $d \mathbf{k} = k dk d \phi$ is the 2-dimensional differential in momentum space and the sum is performed over both valleys.
The analytical expression in Eq.~(\ref{microcurrent}) encapsulates the \emph{exact} light-matter interactions predicted by the Dirac equation (as no approximations were applied) and displays remarkable physics richness. Two current contributions are present, depending on whether the current is originated from electronic transitions within the same band (\emph{intraband}), or across different bands (\emph{interband}).  These can be identified in Eq.~(\ref{microcurrent}) -- intraband contributions are proportional to $(w^{\xi}_{\mathbf{k}} + 1)$, whereas interband contributions depend on the microscopic polarisation $q^{\xi}_{\mathbf{k}}$, leading to two distinct terms. The one proportional to $(\Delta/\epsilon_{\mathbf{k}})$ is a mass-induced contribution and naturally vanishes for ungapped dispersions. It can be seen that, when taking  $\Delta = 0$, both valleys contribute exactly the same to the current, i.e. $\mathbf{j}^{\xi}_{\mathbf{k}}(t) = \mathbf{j}^{-\xi}_{\mathbf{k}}(t)$, leading to a valley degeneracy factor $g_{\rm v} =2$ in the current of Eq.~(\ref{macrocurrent}) as previously reported in Refs. \cite{Carvalho2016, Ishikawa2013}. \\

\section{\label{sec:simulations}Simulations}

The massive Dirac-Bloch equations encapsulate a breadth of optical phenomena which become highly nontrivial in the nonlinear optical regime, once the electrons are coupled to ultrashort and intense light fields. To probe such behaviour, the graphene monolayer is pumped with a normally-incident pulse of duration $t_{0} = 31.9$ fs, central wavelength $\lambda_{0} = 4 \mu$m and frequency $\omega_{0} = 4.71 \times 10^{14}$ s$^{-1}$, photon energy $\hbar \omega_{0} = 0.31$ eV, intensity $I = 0.45$ GW/cm$^2$ and at temperature $T = 0^{\circ}$K. Additionally, realistic localised zero-averaged fields are assumed: $A(t) = A_{0} \sech{t/t_{0}} \sin{(\omega_{0}t)}$ and $E(t) = -\partial_{t}A/c$. We remark that, in order not to introduce unphysical static fields, these fields satisfy $\int_{-\infty}^{\infty}A(t)dt = \int_{-\infty}^{\infty}E(t)dt = 0 $.  \\
For an ultrashort intense pulse, the dephasing mechanisms that account for decay of the populations and polarisations (and phenomenologically accounted for by the decay rates $\gamma_{1}$ and $\gamma_{2}$) are much longer than the pulse input time $t_{0}$ and can thus be safely neglected. Time and Angle-Resolved Photoemission Spectroscopy (ARPES) techniques estimate these relaxation times as $T_1 \approx 150$ fs and $T_2 \approx 0.8$ ps. These figures are heavily affected by a combination of initial temperature, doping, pump fluence, excitation energy
and substrate type and we refer the reader to Ref.~\cite{Gierz2015} for further information on the preparation and underlying physics of the dephasing mechanisms.
Therefore, $\gamma_{1} = \gamma_2 = 0$ are set throughout all simulations. In this coherent regime, the two-level systems in either valley are conservative, leading to a probability conservation law, namely $\partial_{t} (  4 |q^{\xi}_{\mathbf{k}}|^2 + |w^{\xi}_{\mathbf{k}}|^2) = 0 $, which was used to obtain numerical outputs within a strict tolerance threshold of $10^{-9}$.

\begin{figure*}
\includegraphics[scale=0.35]{Fig_BerryTime.png}
\caption{\label{fig4}(colour online) Berry phase in time domain acquired by carriers in the conduction band in the $\mathbf{K}$ valley. With a fixed angle $\phi_{\mathbf{k}} =  \pi/3$, the dependence of the phase on the rescaled momentum magnitude $\tilde{k}$ reveals highly nontrivial dynamics for momentum states close to the Dirac points, where transitions are resonant and hence strongest, showing a step-like behaviour. For off-resonant states, this phase becomes negligible as its amplitude vanishes. A juxtaposition of the electric field envelope reveals that such phase oscillations are highly asymmetrical.}
\end{figure*}

\subsection{\label{sec:currentsandspectra}Currents and respective spectra}

In order to simulate the microscopic current predicted in Eq.~(\ref{microcurrent}), the massive DBEs [Eqs.(\ref{db1}-\ref{db2})] were solved numerically with an explicit, adaptive, parallelised fourth-order Runge-Kutta algorithm.
The physical generated current is obtained once the microscopic contributions are appropriately integrated in momentum space which, when parametrised in its radial and angular components, is respectively composed of a mesh of $1000 \times 500$ states. Since this model assumes infinitely-extending bands, a radial cutoff was imposed such that all relevant microscopic contributions were accounted for. \\
The role of energy gap in the generated current is now studied with the aid of a dimensionless parameter $\tilde{\Delta} \equiv \Delta/(\hbar \omega_0)$, conveniently rescaled such that a gap satisfying $\tilde{\Delta} =1$ is exactly resonant with the pump photons. As previously mentioned, massless Dirac electrons in either valley contribute equally to the generation of current. The linearly-polarised pulse, along the $\hat{x}$ direction, does not create $\mathbf{J}_{y}$ currents which must therefore vanish identically, once their corresponding microscopic currents are integrated over all momenta and valley contributions; this is indeed observed in our simulations, and is a crucial indicator of the validity of our numerics  \cite{Carvalho2016}. In the massive regime, both components are addressed differently by the valleys, even in this simple polarisation configuration. Both valleys contribute equally to the $\mathbf{J}_{x}$ component of the current. As for the $\mathbf{J}_{y}$ component, both valleys create non-zero currents fully out-of-phase which, upon summation, cancel each other out identically. \\

The effect of the mass and Berry phase on the current may be seen in Fig. \ref{fig2}(a), where the full current in time domain $J_{x}(t)$ is shown. Its amplitude increases as the energy gap is increased, until a maximum is reached when the photon is resonant with the energy gap i.e. when $\tilde{\Delta} = 1$. Subsequently, the current amplitude vanishes for increasingly larger gaps.
 This behaviour is best understood if the intraband and interband currents are plotted separately. Fig.~\ref{fig2}(b) shows the intraband current contribution, where it can be seen that its amplitude is maximal when $\tilde{\Delta} = 0$ and monotonically decreasing with increasing energy gap. Fig.~\ref{fig2}(c) shows the interband current, itself composed of the two polarisation-dependent terms in Eq.~(\ref{microcurrent}), once integrated over momentum and valley isospin. The full current dependence on the mass stems primarily from the interband contributions, as Fig. \ref{fig2}(c) follows the pattern just described.
We remark that both interband current terms are in phase. Figs~\ref{fig2}(b,c) further reveal that the full current emerges from a very complex interplay of the competing, out-of-phase contributions of intraband and interband currents. \\
More optically pertinent information can be obtained by analysing the full current spectrum $S(\omega) = |\omega \mathbf{J}(\omega)|^2$, in dB units, versus the harmonics order $\omega / \omega_{0}$, a dimensionless parameter so that the pump pulse is centred spectrally at $\omega/\omega_{0} = 1$, which is displayed in in Fig.~\ref{fig3}(a). The spectra show strong odd harmonics being generated, commonly expected of a $\chi^{(3)}$ material. The exceedingly small peaks found for $\omega/\omega_{0} = 2, 4, ...$ on this logarithmic scale can be seen as numerical artefacts and suggest that even harmonic generation is generally absent. However, particular gap values can be seen to yield rather enhanced even-harmonic peaks. \\
In order to understand the origin of this behaviour, both the intraband and interband current spectra are respectively shown in Fig.~\ref{fig3}(b,c). For both contributions, odd-order harmonic peaks are predominant over even-order harmonic peaks, given the relatively small efficiency of even-harmonic generation expected from centrosymmetric media. The intraband current reveals second harmonic generation in ungapped samples, as previously reported in Ref. \cite{Carvalho2016}. As for the interband current, clear $n^{\text{th}}$-order harmonic peaks appear when the gap is tuned so that $\tilde{\Delta} = n$, for a positive integer $n \geq 2$. We remark that such peaks are always generated for \emph{any} gap value but will not be contribute to particular harmonic orders unless this tunability condition is met i.e. for integer $\tilde{\Delta}$. In particular, when tuned to even integers, even harmonic peaks are generated in the emission spectrum, as shown in Fig.~\ref{fig3}. Physically, the observed even harmonic peaks do not arise from $\chi^{(2)}$-like processes (occurring only in non-centrosymmetric media) but are rather understood through the coherent interference among odd harmonics, a well-known strong-field effect termed ''\emph{THG in disguise of SHG}''  occurring exclusively at the femtosecond scale. We refer the reader to pages 157-158 of Ref. \cite{Wegener2005}, as well as Ref. \cite{Tritschler2003} for more information on the physics underlying this process.

\begin{figure*}
\includegraphics[scale=0.43]{Fig_BerryONOFF.png}
\caption{\label{fig5}
(colour online) Comparison of current spectra in the presence or absence of the Berry phase for a gap $\Delta = 2 \hbar \omega_{0} = 0.62$ eV. The full current spectrum is shown in (a), displaying a general peak suppression of the dominant, odd harmonics. Rather negligible even harmonic peaks only exist once this phase is considered. These features are caused by the dominant, interband current, plotted in (c). The effect of the harmonic interference on the Berry phase contribution to the interband current for such a gap, here seen through the peak suppression at $\omega/\omega_{0} = 2$. The intraband spectrum, shown in (b), displays a general suppression of odd harmonics, whose extent increases as the harmonic order increases.}
\end{figure*}

\subsection{\label{sec:effectofBerryphase}The effect of the Berry Phase}

The Berry phase in Eq.~(\ref{Berry}) is only present for non-vanishing gaps and induces non-trivial contributions to the current spectra just discussed. Before engaging in determining its role in harmonic generation, the temporal dynamics of the Berry phase for various momentum states is shown in Fig.~(\ref{fig4}), by splitting the radial and angular component of their momentum vector $\mathbf{k}$ respectively as $k$ as $\phi_{\mathbf{k}}$, as well as rescaling $k$ to a dimensionless magnitude $\tilde{k} \equiv (2 v_{\rm F}/\omega_{0})k$. In this fashion, electrons of wavevector $\tilde{k} = 1$ have $v_{\rm F} k = \omega_{0}/2$. By juxtaposing the envelope of the electric field, it can be seen that this field-dependent phase evolves rather nontrivially, oscillates asymmetrically, and is highly dependent on which electronic state is considered.
It is instructive to see how it changes for increasing magnitudes for a fixed angle, here arbitrarily taken as $\phi_{\mathbf{k}} =  \pi/3$. High-momentum states i.e. when $\tilde{k} \geqslant 1 $, extremely detuned from the gap acquire a very small phase, which vanishes monotonically very rapidly, as the magnitude is increased. Understandably, the phase is mostly relevant for quasi-resonant states, found in the vicinity of the Dirac points i.e. $\tilde{k} = 0$, where microscopic polarisations $q_{\mathbf{k}}$ are strongest. At the Dirac points, such contributions are maximal and the Berry phase undergoes continuous step-like transitions between $0$ and $-\pi/2$, independent of the field dynamics. Note that \ref{fig4} shows the phases acquired by electrons in the conduction band in the $\mathbf{K}$ valley. The relative signs acquired for each band and valley, as derived in Sec. \ref{sec:quasirelativisticdynamics}, are depicted in Fig.~\ref{fig1}(a). For instance, valence band carriers acquire a relative negative sign.  \\

The role of the Berry phase on the generation of new harmonics is now discussed. In order to achieve this, the Berry phase and its derivative are neglected by setting $\gamma_{\mathbf{k}}(t) = \dot{\gamma}_{\mathbf{k}}(t) = 0$ in the massive DBEs [Eqs.~(\ref{db1})-(\ref{db2})] and in the microscopic current of Eq.~(\ref{microcurrent}). This procedure is physically consistent since $\mathbf{J}_y (t)$ still vanishes after such terms are disregarded. We proceed by comparing the spectra of the full current and its intraband/interband contributions, obtained by including or excluding such terms. General features can be captured and are exemplified for a particular gap with $\tilde{\Delta} = 2$ (resulting in a realistic energy gap value of $\Delta = 0.62$ eV), whose spectra are shown in Fig.~\ref{fig5}. One can observe that the Berry phase acts on the full current , shown in Fig.~\ref{fig5}(a), and considerably suppresses odd harmonic harmonics and enhances the relevant even harmonic harmonics. The extent of the odd harmonics suppression seems to grow for higher harmonics but is much more prominent in the interband currents, where the peak differences are biggest.\\
The full behaviour can again be seen to originate from the dominant, interband contributions, plotted in Fig.~\ref{fig5}(c). The Berry phase can be identified as the agent that mostly drives odd harmonics interference (and consequently possible even harmonic generation when appropriately tuned) by considering the substantial peak enhancement at $\omega/\omega_{0} = 2$ when the phase is switched on, as previously discussed. These results can again be understood in light of the discussed THG in disguise of SHG. The intraband current spectrum comparison is shown in Fig.~\ref{fig5}(b), where it can be seen that any possibly small even harmonic peak vanishes once the Berry phase is neglected.

\section{\label{conclusions}Conclusions}

In conclusion, we show that relativistic two-dimensional massive fermions acquire a Berry phase when interacting with normally-incident electromagnetic pulses. This phase, whose analytical form is given and plotted for some states, only exists in the presence of an energy gap, a consequence of the inequivalence of both sublattices that decompose the honeycomb lattice of graphene. The spectrum generated by the electronic nonlinear current shows prominent odd-harmonic generation, which is generally suppressed as the energy gap is increased. Although even harmonics are generally absent for gapped dispersions, we show that their generation may be attained at the femtosecond scale through THG in disguise of SHG when the photon energy is appropriately tuned to the energy gap, generating radiation with the desired harmonic order. These processes may be conceptualised as coherent interactions of odd harmonics. Signatures of these interband-driven phenomena can be seen in the enhancement of harmonic peaks. This mechanism is to be contrasted with what has been found in pristine ungapped graphene samples, where the dynamical centrosymmetry breaking mechanism allows for even harmonics to be generated via intraband currents. We also show that the Berry phase plays a major role in the interband current dynamics and hence in the generation of even harmonics. We remark that excitonic effects are absent in the present formalism. These results and methods help establish new techniques to understand and predict the nonlinear optical behaviour of a range of two-dimensional hexagonal relativistic-like semiconductors, and help pave the way to predict quantitatively, in a generalised fashion, the effect of wide range of intrinsic or deliberate properties and phenomena, such as monolayer-substrate interactions, sample imperfections, local defects and strain effects, expected to be found in more realistic samples. \\
D.C. acknowledges insightful discussions with Dr. Lawrence Philips in Prof. Patrik Ohberg's group.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}\label{sec:intro}

Gravitational  lensing  has  become  a tool  of  greater  and  greater
importance  to cosmology  and the  study of  galaxy structure.   Large
scale  weak lensing  surveys are  being used  to measure  cosmological
parameters     and     study     dark     energy\footnote{see     DES:
  http://www.darkenergysurvey.org,                          Pan-Starrs
  :http://pan-starrs.ifa.hawaii.edu}  and larger  ones are  planned in
space       and       with        purpose       built       telescopes
\citep{euclidredbook,anderson01}.  At the same  time strong lensing is
being           used            to           study           cosmology
\citep{2012arXiv1208.6010S,2013arXiv1306.1272T}    and     the    mass
distribution  within galaxies  down  to  unprecedentedly small  scales
\citep{metcalf12,2013arXiv1307.4220X,2012Natur.481..341V}.   Long term
monitoring  of strong  lenses is  being carried  out with  networks of
telescopes    spanning   the    globe    to    obtain   time    delays
\citep{2012A&A...538A..99S,treu2013}.   The   Hubble  Space  Telescope
(HST)     has      spent     many     orbits      observing     lenses
\citep{postman12,coe12,coe13,zitrin13}     and    many     more    are
scheduled.\footnote{HST                Frontier                Fields:
  http://www.stsci.edu/hst/campaigns/frontier-fields/}   Getting   the
most out  of these huge  observational efforts will require  very good
simulations of  the lensing phenomena  for interpretation of  the data
and to control  systematic effects.  These simulations should
include all  the physical and  instrumental effects that we  know will
affect the data.  Currently, weak lensing simulations have limitations
in  resolution and  the contribution  of baryonic  physics is  usually
neglected.  In addition, they often use unperturbed light rays instead
of  fully  tracing their  paths  through  the simulation.   Individual
galaxies and the distortions they experience are not resolved in these
simulations.  Strong  lensing and  cluster lens simulations  often use
highly idealized  lens models and  usually neglect the  environment of
the lens and contributions to the lens that lie at different redshifts
along the  line of sight.  The  {\small GLAMER} code is  an attempt to
improve this situation in a way  that can be applied to many different
lensing situations on different scales.

Typically, simulations of strong  gravitational lensing treat the lens
as  a  single  plane  without  foreground  and  background  structure.
Systematic  studies,   on  clusters  extracted  from   a  cosmological
numerical      simulation,      have       been      conducted      by
\citet{meneghetti10a,meneghetti11}    With   the    aim   of    better
understanding  of   the  properties   of  galaxy  clusters   that  can
potentially act as strong lenses. In these works however the simulated
clusters are analyzed  using the single lens  approximation.  The same
method has  also been  used by \citet{giocoli12,boldrin13},  where the
individual cluster components are analytically modeled. However, it is
probable that  structures along the  line of sight have  a significant
effect on the observed properties of the lens: both in the weak and in
the strong  lensing regime.  For  example, the inconsistencies  in the
flux         ratios         between        four-images         quasars
\citep{metcalf01,dalal02,2002ApJ...567L...5M,2006MNRAS.367.1367A,maccio06b,2009A&A...498...49C,2009MNRAS.398.1235X,2011ApJ...741..117C}
can  be   at  least  partially  explained   by  line-of-sight  objects
\citep{M04b,2009MNRAS.398.1298P,2005ApJ...635L...1W,2010MNRAS.408.1721X}.
Furthermore, in cluster  lensing, the mass and  concentration of halos
derived from lensing  properties can be affected  by the line-of-sight
contribution.  \citet{2009MNRAS.398.1298P} showed that the presence of
uncorrelated large scale  structures can boost the  Einstein radius by
as much as 50\%.

Weak lensing simulations are usually done by creating a shear map on a
fixed resolution  grid via FFT  of the  potential on one  uniform grid
\citep{2000ApJ...530..547J,  2003ApJ...592..699V, 2007A&A...471..731P,
  2013arXiv1309.1161A,        2001MNRAS.327..169H,2009ApJ...701..945S,
  2011ApJ...742...15T} or on  separate grids for long  and short range
\citep{2007MNRAS.382..121H,2009A&A...499...31H}.   The  resolution  of
the grid is much larger than the  image of a typical galaxy.  The shear at the
positions of  individual galaxies is  found by interpolating  off this
grid and is assumed  to be uniform over the image  of the source.  The
lensing  mass in  the light -cone is  projected onto  multiple lensing
planes and rays are traced from  the observer to the source.  Also the
simulation used to represent the  mass distribution typically does not
have high enough resolution to represent the inner regions of halos or
any  of  the  baryonic  physics happening  there.   Because  of  these
limitations no  strong or quasi-strong lensing  of individual galaxies
will  occur  in these  simulations.   Although  the contribution  from
smaller scales  to the weak lensing  might be small it  is a potential
source of  systematic error and it is  important to characterize it  in the
future high precision weak lensing surveys.

To ensure better resolution and less computational demand, other codes
have  used the  tree  force  algorithm \citep{1986Natur.324..446B}  to
calculate     the    deflection     and     distortion    of     rays.
\cite{1998ApJ...494...29W}   used   this   method   particularly   for
simulations of  quasar microlensing.  And  \cite{ 2007MNRAS.376..113A}
use it for lensing by N-body  particles.  Both of these use the single
plane approximation.

Another  set of  codes \citep{1999MNRAS.306..567F,2012MNRAS.420..155K}
use the ray  bundle method, where a collection of  rays is traced from
the image to the observer, integrating the and computing the shear and
convergence  of the  bundle along  the way.    A third  class of  codes
\citep{1999MNRAS.308..180C,2008MNRAS.388.1618C}      calculate     the
distortion to ray-bundles along unperturbed  paths that pass through a
three dimensional  simulated potential  of a  cosmological simulation.
Still others  use  the  unperturbed  paths,   also  know  as  the  the  Born
approximation,  but put  the  mass from  the  simulations on  multiple
planes \citep{2004APh....22...19W,2007MNRAS.382.1494H}     The     Born
approximation  can   cause  systematic  errors  and   is  clearly  not
sufficient in the case of strong lensing.
\begin{figure*}
\centering
\includegraphics[width=0.95\textwidth]{./Figures/plot.png}
\caption[]{A schematic view of the multi-plane formalism, as described
  in eqn.~\ref{eq:recursion_x}. \label{fig:scheme}} \end{figure*}

In  \cite{M&P2013} (here  after MP)  we describe  the {\small  GLAMER}
code's structure, its adaptive mesh refinement capabilities and how it
calculates the deflection of light rays by a single lens plane.  These
capabilities  alleviate  many  of  the  limitations  codes  have  with
resolution and make it possible  to construct fully resolved images of
millions  of galaxies  in  a  single simulation.   In  this paper,  we
describe the code's  ability to simulate lensing  in three dimensions;
that is  with the  lensing mass and  also the  sources distributed
throughout the light-cone in angle and redshift.  We also discuss some
of the ways a light-cone can be generated and input into the code.

There  are many  interesting questions  that can  be addressed  with a
lensing code  like the  one we  have developed.   Among them  are: the
importance of line-of-sight objects to strong lenses, the importance of
strongly  lensed  individual galaxies  to  weak  lensing surveys,  the
statistical properties of strong lenses, the importance of multi-plane
effects on precision weak lensing, and many more.

We  start  this   paper  by  presenting  the  basics   of  lensing  in
Sec.~\ref{sec:basics}.   We  elaborate   on  the  multiple  plane
formalism   in   Sec.~\ref{sec:multip_formalism}.   We   continue   in
section.~\ref{sec:implementation} by describing  the implementation of
this    formalism     in    the    {\small    GLAMER}     code.     In
Sec.~\ref{sec:light_cones} we  describe the  different types  of light-cones implemented  in the code.  In  Sec.~\ref{sec:mill_lc} we discuss
the   use a  simulated halo catalog to populate the  light-cone  and   in
Sec.~\ref{sec:massfunc_lc}, the generation of a light-cone from a mass
function.  Then  in Sec.~\ref{sec:tests}  we perform several  tests of
the  code, determining  the optimal  number  of lensing  planes to  be
used.  We also  show some  future applications  of the  code, such  as
Einstein radius mapping of galaxy clusters and producing galaxy-galaxy
strong   lens    catalogs.    We    conclude   our    discussions   in
Sec.~\ref{sec:conclusions}.   For more  extended  information we  have
also included several appendices.

\section{lensing basics}\label{sec:basics}
We will assume  a weak field approximation to  General Relativity.  In
this case, light  propagates through the universe on  the geodesics of
the  perturbed  Robertson-Walker metric.   The  line  element of  this
metric in the Newtonian gauge is \begin{equation} \mathrm{d}s^2 = a^2\left[ -\left(
  1+\frac{2\Phi}{c^2}\right)\mathrm{d}\tau      ^2       +      \left(
  1-\frac{2\Phi}{c^2}\right)\mathrm{d}l^2 \right],
\label{eqn:metric}
\end{equation}  where  $a$  is  the  cosmological scale  factor,  $\Phi$  is  the
Newtonian  potential,  $c$  is  the  speed of  light,  $\tau$  is  the
conformal time  time, and \begin{eqnarray}  \mathrm{d}l^2 &=& \mathrm{d}\chi  ^2 +
\mathcal{D}^2(\chi)\mathrm{d}\omega    ^2   ,   \\    \mathrm{d}\omega^2   &=&
\mathrm{d}\theta^2+\sin^2 \theta \mathrm{d}\phi^2  .  \end{eqnarray} Here $\chi$
is  the  comoving  distance  and  $\mathcal{D}(\chi)$,  the  comoving  angular
diameter distance, is defined for different curvatures $K$ as follows:
\begin{eqnarray} \mathcal{D}(\chi) =
\begin{cases}
R \sin( \chi/R) & K=+1, \\
\chi & K=0, \\
R \sinh( \chi/R) & K=-1.
\end{cases}
\end{eqnarray}    where    $R$   is    the    curvature    scale   ($R^{-1}    =
H_o\sqrt{1-\Omega_m-\Omega_\Lambda}/c$  in  the standard  cosmological
model before radiation domination).

The geodesic  equation for a  null geodesic  on small enough  scale to
ignore the curvature $R$ is
\begin{equation}
\frac{\mathrm{d}^2 {\bf x}}{\mathrm{d}\lambda^2} = - \frac{2}{c^2} \left[ \nabla \Phi - \left( \frac{\partial {\bf x}}{\partial \lambda } \cdot \nabla \Phi \right) \frac{\partial {\bf x}}{\partial \lambda } \right]
\end{equation}
neglecting higher order terms in  $\Phi$ and $\lambda$ being an affine
distance along the  path.  From the prospective of  the light ray this
can be written more simple as
\begin{equation} \label{eq:geodesic}
\frac{\mathrm{d}^2 {\bf x}_\perp}{\mathrm{d}\lambda^2} = - \frac{2}{c^2} \nabla _\perp \Phi\left( {\bf x} \right)  ~~~,~~~ \frac{\mathrm{d}^2 {\bf x}_\parallel}{\mathrm{d}\lambda^2} = 0
\end{equation}
where  ${\bf x}_\perp$  and  ${\bf x}_\parallel$  are the  coordinates
perpendicular   and  parallel   to   the   direction  of   propagation
respectively.

If we imagine  a thin slab that the ray  enters perpendicularly we can
integrate equation~(\ref{eq:geodesic}) with  a unperturbed path within
the slab.   We can then calculate  the deflection angle of  the ray on
leaving the slab,
\begin{equation}
{\pmb \alpha}\left({\bf x}_\perp\right) = - \frac{2}{c^2} \nabla _\perp \phi({\bf x}_\perp)  ,
\end{equation}
where $\phi({\bf x}_\perp)$ is now the surface potential.

The surface potential is related to the surface density by
\begin{equation}
\nabla _\perp^2 \phi({\bf x}_\perp) = 4 \pi G \,\Sigma\left({\bf x}_\perp\right)
\end{equation}
where $\Sigma({\bf x}_\perp)$ is the surface density within the slab
which can be solved to get
\begin{equation}\label{eq:poissons}
\phi({\bf x}) = \frac{4G}{c^2} \int \mathrm{d}^2x'~ \Sigma({\bf x}') \ln|{\bf x} - {\bf x}'|.
\end{equation}

The position on the sky will be ${\pmb \theta}$.  The position of this
point where  there no  lensing will be  ${\pmb \theta}^S$,  the source
position.  The magnification matrix, or the Jacobian matrix of the map
between these coordinates, is
\begin{equation}\label{eq:A}
{\bf A} \equiv \left[ \frac{\partial \theta^S_i}{\partial \theta_j}\right].
\end{equation}
This matrix  is often  decomposed into the  convergence, $\kappa({\pmb
  \theta})$ and  the components of shear  $\gamma({\pmb \theta})$.  In
the    notation     of    appendix~\ref{sec:ellipt-matr-repr},    this
decomposition is
\begin{equation}\label{eq:Adecomposition}
{\bf A} = (1-\kappa) {\bf I} - \gamma_1{\pmb \sigma}_1 - \gamma_2{\pmb \sigma}_2 - \gamma_3{\pmb \sigma}_3
\end{equation}
We have included a third  component of ``shear'', $\gamma_3$, which is
not usually included, but, as we  will see, is required when there are
more than one lens plane.   The magnification of an infinitesimal point
on the sky is then \begin{equation} \mu = \frac{1}{{\rm det} {\bf A}}.  \end{equation}

\subsection{single lensing plane}
It is instructive to consider the case of a single lens plane where it
is assumed that the light rays follow unperturbed geodesics outside of
the lens plane.  This is the case usually considered in the literature
so  expressing it  here  in our  notation will  serve  to clarify  our
approach and what new phenomena arise from multi-plane lensing.

In this case, the source and image positions are related by
\begin{equation}
{\pmb \theta}^S = {\pmb \theta} - \frac{\mathcal{D}_{ls}}{\mathcal{D}_s}{\pmb \alpha}\left( {\pmb \theta} \right)
\end{equation}
where  $\mathcal{D}_s$ is  the angular  diameter  distance to  the source  and
$\mathcal{D}_{ls}$, between lens plane and  the source (${\bf x}_\perp = {\pmb
  \theta} \mathcal{D}_l$).   In this case,  the convergence can be  written in
terms of the surface density on the plane
\begin{equation}
\kappa({\pmb \theta}) = \Sigma({\pmb \theta})/\Sigma_{\rm crit} ~~~,~~~~
\Sigma_{\rm crit} \equiv \frac{c^2}{4\pi G} \frac{\mathcal{D}_s}{\mathcal{D}_l\mathcal{D}_{ls}}
\end{equation}
so  in   this  case  $\kappa({\pmb  \theta})$  can   be  considered  a
dimensionless surface density.  This is not the case for multiple lens
planes.  The  deflection field  is a potential  field which  mean that
$\gamma^3 =0$, the field ${\pmb \theta}^S - {\pmb \theta}$ has no curl
and the  shear field  has no  B-modes.  Also every  point on  the lens
plane has a  one-to-one correspondence to a point  on the image plane.
As  we  shall see  these  things are  not  guaranteed  when there  are
multiple lens planes.

\begin{figure*}
\centering
\includegraphics[width=0.33\textwidth]{./Figures/massfunctioninFieldPS.pdf}
\includegraphics[width=0.33\textwidth]{./Figures/massfunctioninFieldST.pdf}
\includegraphics[width=0.33\textwidth]{./Figures/massfunctioninFieldPL.pdf}
\caption[]{Halo  mass function,  recreated in  a light -cone of  $10^6
  \mbox{arcsec}^2$, for different redshift bins.   The left panel shows the
  Press-Schechter mass  function, the  middle one --  the Sheth-Tormen
  one, and  the right one --  a power law  with a slope of  $1/6$. The
  distribution of  masses (histograms), for  which we put  together 64
  different light-cone  realizations, is  in very good  agreement with
  the     theoretical    mass     functions     (lines)    for     all
  redshifts.\label{fig:massF}}
\end{figure*}

\section{multi-plane formalism}\label{sec:multip_formalism}
In the  multi-plane formalism, the  mass in the universe  is projected
onto multiple lensing  planes, perpendicular to the line  of sight. As
light travels from  the lensed object to the  observer, it experiences
deflection from  each plane and  passes unperturbed through  the space
between  the  planes.   The  multiple  plane  lens  equation  has  the
following  form \begin{equation}  {\pmb \theta}^S  ={\pmb \theta}  - \sum_{i=0}^{N}
\frac{\mathcal{D}_{is}}{\mathcal{D}_s} {\pmb \alpha}^i({\pmb \theta}^i) \, ,
\label{eqn:basicM}
\end{equation} where  $N$ is  the number  of planes, $\rm  \mathcal{D}_s$ is  the angular
diameter distance  to the source  and $\mathcal{D}_{is}$ -- between  plane $i$
and the source. See figure~\ref{fig:scheme}.

For the  purpose of our  work, from now on  we will use  only comoving
positions  and angular  distances: \begin{eqnarray}  {\bf x}  = (1+z)x  ,\\ {D}  =
(1+z)\mathcal{D}  ,  \end{eqnarray}   where  $z$  is  the  redshift  of   the  plane  in
consideration.

\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{./Figures/pos.pdf}
\caption[]{Position of the four images (different symbols) of a lensed
  quasar  at redshift  $z=3.62$ by  a galaxy  at redshift  $z=0.42$ in
  simulations with different numbers of  lensing planes (6, 8, 12, 14,
  16, 18, 20, 22,  24, 26, 28, 30, 32, 44, 56, 68,  80, 92, 104).  The
  color of the symbols darkens  with increasing number of planes.  The
  positions clearly converge as the number of planes increases.
\label{fig:pos}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.53\textwidth]{./Figures/all.pdf}
\caption[]{Image  positions and  magnification  as a  function of  the
  comoving distance  between the  lensing planes (different  number of
  planes) for the lensing simulation  of a quasar at redshift $z=3.62$
  by a  galaxy at redshift  $z=0.42$. Separation of  approximately 300
  \mbox{comoving Mpc}   \,    corresponds   to    22   lensing   planes    in   this
  case. \label{fig:all}} \end{figure}

It  is important  to note  here that  since we  calculate the  angular
diameter  distances  between   planes  assuming  the  Robertson-Walker
metric, the mass added to those lenses planes will cause our light-cone
to  effectively   contain  more  mass   than  the average  for  the
Friedman-Lemaitre universe used  to calculate the distances  if it is
not compensated  for --  the average angular  size through  the planes
will  the  smaller on  average.   This  is not  a  small  effect if  a
significant fraction of the mass in the universe is represented on the
planes.  When  needed we correct  for this  by add a  uniform negative
mass density to each of the  planes that compensates for the clustered
mass  so that  the  average  surface density  is  zero  on each  plane
\citep{1988ApJ...327..526S}\footnote{This could be  avoided by using a
  Dryer-Roeder  angular diameter  distances.}.  This  will be  further
discussed in section~\ref{sec:light_cones}.

The deflection angles at each plane add to give the angle with respect
to    the    radial   direction    at    the    $i$th   plane.    From
equation~(\ref{eqn:basicM}) it  follows then that the  position on the
$(i+1)$th plane is
\begin{equation}
{\bf x}^{i+1} = {\bf x}^i - D_{i+1,i} \left\{  {\pmb \theta} +
\sum_{j=1}^{i-1} {\pmb \alpha}^j({\bf x}^j) \right\}.
\end{equation}
Applying this to  the $i$th plane in terms of  the $(i-1)$th plane and
subtracting the equations gives a recursion relation for the positions
on the planes
\begin{equation} \label{eq:recursion_x}
{\bf x}^{i+1} = \left( \frac{D_{i+1,i} }{D_{i,i-1}}+ 1 \right){\bf x}^i -
\frac{D_{i+1,i}}{D_{i,i-1}} {\bf x}^{i-1} - D_{i+1,i} {\pmb \alpha}^i({\bf x}^i)
\end{equation}
with the initial conditions
\begin{equation}
{\bf x}^0 = 0~~~~,~~~~{\bf x}^1 = {\pmb \theta} D_1.
\end{equation}
The first  of these is the  requirement that all the  rays converge at
the observer.

The magnification matrix  defined in \eqref{eq:A} can  be expressed on
each plane as
\begin{equation}
{\bf A}^i \equiv \frac{1}{D_i} \frac{\partial {\bf x}^i}{\partial
  {\pmb \theta} }.
\end{equation}
This is what the normal magnification matrix would be if the $(i+1)$th
plane      where      the     source      plane.       Differentiating
equation~\eqref{eq:recursion_x} yields
\begin{align}
{\bf A}^{i+1} =
\frac{D_{i+1,i-1}}{D_{i,i-1}}  \frac{D_{i}}{D_{i+1}} ~{\bf A}^i
- \frac{D_{i+1,i}}{D_{i,i-1}}  \frac{D_{i-1}}{D_{i+1}} ~{\bf A}^{i-1}  \\
- \frac{D_{i+1,i} D_{i} }{D_{i+1}} ~ {\bf G}^i {\bf A}^i \label{eq:1}
\end{align}
 with the initial conditions
\begin{equation}
{\bf A}^0 = 0 ~~~~~~ {\bf A}^1 = {\bf I}.
\end{equation}
The forcing term is
\begin{equation}
{\bf G}^i \equiv \frac{\partial {\pmb \alpha}^i}{\partial
  {\bf x}^i } = g_0^i {\bf I} + g_1^i {\pmb \sigma}_1 + g_2^i {\pmb \sigma}_2 .
\end{equation}
where the notation of Appendix~\ref{sec:ellipt-matr-repr} is used.
The coefficients are related to the Newtonian surface potential on
the planes by
\begin{align}\label{eq:shears}
g_0 = \frac{1}{2} \nabla^2\phi ~~~~~~
g_1 = \frac{1}{2} \left(\phi,_{11}
  - \phi,_{22} \right) ~~~~ g_2 = \phi,_{12}
\end{align}
The term \eqref{eq:1} can be written out as
\begin{align} \label{eq:GA}
{\bf G A} = \left[ g_0 (1-\kappa) - g_1\gamma_1 - g_2\gamma_2 \right]
{\bf I} \nonumber \\
+ \left[ -g_0 \gamma_1+ g_1(1-\kappa) + g_2\gamma_3 \right]
{\pmb \sigma}_1 \nonumber \\
+ \left[ -g_0 \gamma_2+ g_2(1-\kappa) - g_1\gamma_3 \right]
{\pmb \sigma}_2 \nonumber \\
+ \left[ -g_0 \gamma_3- g_1\gamma_2 + g_2\gamma_1 \right]
{\pmb \sigma}_3
\end{align}

It  can be  seen  in  equation~\eqref{eq:GA} that  there  is a  ${\pmb
  \sigma}_3$  component to  the  magnification matrix  which does  not
exist      in      the      case       of      one      lens.       In
appendix~\ref{sec:ellipt-matr-repr}, it  is shown that this  implies a
rotation in  the image  \citep{2006MNRAS.367.1543P}.  The case  of two
lens  planes is  worked  out explicitly  for  pedagogical purposes  in
appendix~\ref{sec:two-plane-lens}.

\section{code implementation}\label{sec:implementation}

Multi-plane lensing is implemented as  an extension to the preexisting
{\small GLAMER} code.  Other  aspects of the {\small GLAMER} code are  described in more
detail  in \cite{M&P2013}  (MP).  The  code is  written in  C++ in  an
object oriented manor so that characteristics of the lenses and sources
can be  chosen by  the user in  a very flexible  fashion.  There  are a
number of options allowed which will be described briefly here.

Rays  are   shot  from  the   observer  to  the  source   plane  using
equations~\eqref{eq:recursion_x} and \eqref{eq:GA}.  In each plane the
deflection,  convergence and  shear are  calculated using  the methods
described  in  MP.   The  mass  distribution  on  each  plane  can  be
represented in several  different ways.  A surface density  map can be
input in  FITS format.  This  option is useful for  representing the
output   of  N-body   simulations   and   semi-analytic  methods   for
constructing   galaxies    and   galaxy   clusters   such    as   {\small MOKA}
\citep{2012MNRAS.421.3343G}.   The   lens  planes  can   also  contain
analytic ``halos''.   The halos can  have a variety of  different mass
profiles -- Navarro-Frenk-White  (NFW) \citep{navarro97}, non-singular
isothermal sphere  (NSIE), power-law,  \cite{hernquist90}, point-mass,
etc.   These  ``halos''  can   represent  the  baryonic  component  of
galaxies, the dark  matter halos or even stars.   The contributions of
the halos  to the  lensing quantities are  calculated using  a modified
tree algorithm when  there are enough of them to  warrant it (see MP).
Finally,  the  mass  on  a  plane can  be  represented  by  simulation
particles  with  an  adaptive  smoothing in  which  case  the  lensing
quantities are calculated by tree  algorithm.  The code allows for any
combination of  these representations.   For example, the  dark matter
halos could be represented by NFW  profiles, the baryonic galaxy by an
NSIE, the mass outside  of halos by a mass map and  the stars by point
masses all at the same time.

Once the  lens has been constructed  rays are shot back  to the source
plane in a uniform grid on the image plane.  This initial gridding can
be used to make shear or  magnification map with uniform resolution if
that is desired.  As described in MP there are a number of options for
refining the  grid of  rays.  The  code can be  made to  find critical
curves and caustics and increase the resolution of them to the desired
level.  It can also be made to refine the grid to resolve a particular
source.  The  ray shooting is  parallelized with POSIX  threads, which
increases the performance significantly  especially when the number of
planes is large.

There are also a variety of options for representing the sources.
The simplest source is circular with uniform surface brightness.  The
source can also have an analytic surface brightness distribution or a
pixelized surface brightness distribution that can be input in FITS
format.  There can be any number of sources and they can be of mixed
type.  Each source can also have a different redshift.  Because of the
adaptive gridding of the rays each source can be resolved by shooting a
relatively small number of rays to the particular redshift of the
source.

\begin{figure}
\centering
\includegraphics[width=0.53\textwidth]{./Figures/mean.pdf}
\caption[]{Relative  error~\eqref{eq:error} of the  lensing properties
  -- inverse  of the magnification (top), magnitude  of the deflection
  angle,  magnitude  of the  shear  and  convergence  (bottom) --  for
  simulations with the same  light-cone, but differing comoving distance between
  the  lensing planes  (different  number of  planes).  Separation  of
  approximately  300 \mbox{comoving Mpc} \,  corresponds to  18 lensing  planes. The
  field of view  is $10^4 \, \mbox{arcsec}^2$, the  mass function goes down
  to $10^9\mbox{M}_\odot$ halos  and up to redshift $z=2.0$.  The true solution
  in this  case is assumed  to be the  one with one plane for each halo
  at precisely its redshift. \label{fig:mean}} \end{figure}

\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{./Figures/kappa_pdf.pdf}
\caption[]{Distribution  of  the  deflection angle,  convergence,  and
  shear of  a set of  32 random and  uncorrelated fields of  view. The
  area is  $10^4 \, \mbox{arcsec}^2$ and the  NFW halos are from  a ST mass
  function down  to $10^9\mbox{M}_\odot$, distributed on 20  planes. The source
  is at redshift $z=2.0$. The results agree very well with simulations
  by \citet{Pace2011} \label{fig:pdf}}
\end{figure}

\section{Generating Light Cones}\label{sec:light_cones}
In this  section, we describe the  different types of light-cones that are
implemented within  or can be  input into  the {\small GLAMER} \,  code.  The
mass within  the light-cone can  be generated in four  different ways:
($i$) particles  can be taken directly  from a N-body simulation  of a
light-cone and  sorted onto planes, ($ii$) halos can  be identified in
the  N-body simulation  and inputted  into {\small GLAMER}  as a  halo catalog,
($iii$) random  halos can be  generated within {\small GLAMER} from  a redshift
dependent  mass  function and  other  analytic  descriptions of  their
internal parameters,  ($iv$) two dimensional maps  of surface density,
which  are  generated  from   a  N-body  simulation  or  semi-analytic
simulation, can be  input into {\small GLAMER}.  A light-cone  can have lensing
planes of any  mixture of these types -- for  example a galaxy cluster
might be represented  by a mass map and the  foreground and background
objects generated randomly from a mass function.  We will discuss each
of these cases below.

\subsection{Light-Cones from Simulation Particles}\label{sec:particle_lc}
{\small GLAMER}  allows the  user to  input a  list of  particle positions  and
smoothing  scales  for a  light -cone.   The  smoothing is  best  done
adaptively according to particle density  in three dimensions -- as in
Smooth Particle  Hydrodynamics (SPH).   The particles are  sorted onto
the  desired number  of  lens  planes and  then  a  tree structure  is
constructed  on  each plane  to  calculate  the deflection  and  other
lensing  quantities  (see  MP   for  details).   For  state-of-the-art
numerical simulations the number of  particles within a complete
light-cone
can be  very large which makes memory  management difficult.  For
this reason, we favor using methods  ($ii$) or ($iv$) or a combination
of  them to  represent  the output  of a  N-body  simulation, but  the
ability to use the particles  directly is maintained.  A simulation of
a single object with relatively  few particles can be implanted within
a light-cone that is represented in a different way.

\subsection{Light-Cones from Halo Catalog}\label{sec:mill_lc}
{\small GLAMER} \,  is able to read in a catalog of halos and their
properties.  These halos can be taken from an N-body simulation or any
other source.   The halos are projected onto the lensing planes
keeping their angular positions fixed and a tree structure is
constructed on each plane.  The force and deflections are calculated
by tree algorithm (see MP).

To represent  an N-body halo  as an  analytic halo, parameters  of the
halo model must be matched  to quantities calculated from the particle
distribution  --  virial  mass,  half mass  radius,  maximum  circular
velocity,  moments  of  inertia,  etc.   This  procedure  is  somewhat
dependent on the  halo catalog used and what  information it provides.
One can also attempt to represent the baryons that might not be in the
simulation with  an additional ``halo''  within the dark  matter halo.
This requires further  modeling of the relation between  halo mass and
stellar mass  and the  distribution of baryonic  mass within  a N-body
halo of a  given type.  These are interesting subjects  that are being
actively investigated.  Lensing has  great potential for investigating
these relationships.   We will not  discuss them in detail  here since
our purpose is  to introduce a tool for investigating  them and not to
advocate any particular model.

We have implemented one full pipeline  for N-body to sky images using
the Millennium Simulation  catalog of halos and with  sources from the
Millennium Run Observatory  \citep{2013MNRAS.428..778O}.  This project
will be discussed in a separate paper.

\subsection{Light-Cones from Generated Halos}\label{sec:massfunc_lc}

The halos catalog can also be created internally and then sorted onto
lens planes as in the case a catalog from a N-body simulation.   No
attempt has been made yet to cluster these halos or to create subhalos
within halos for the full light-cone.  The halos masses are drawn
randomly from a mass function.   The mass
functions implemented in our code so far are:
\begin{itemize}

\item \citet{press74} (PS);

\item \citet{sheth99b} (ST);

\item Power-Law (PL), given by the following equation:
\begin{eqnarray}
\nonumber \dfrac{\mathrm{d}N}{\mathrm{d}\ln M} &=& \dfrac{A}{M^2} \sqrt{\dfrac{1}{2 \pi}}
 \left[ \dfrac{M}{M^*(z)} \right]^{\alpha}  \\
&\times& \exp \left\{ -\dfrac{1}{2} \left[\dfrac{M}{M^*(z)} (1+z)^{3/2} \right]^{\zeta \alpha}\right\}\,;
\end{eqnarray}
where  $\alpha$  represents  the  slope  of  the  power-law,  $A=0.2$,
$\zeta=1.3$ and  $M^*(z)$ the  non-linear mass  at redshift  $z$ given
such  that  $\delta_c(z)=S(M)$,  with  $\delta_c(z)$  the  overdensity
required for spherical collapse at $z$,  and $S(M)$ is the variance in
the linear fluctuation field when smoothed with a top-hat filter of scale
$R = (3M/4  \pi \bar{\rho})^{1/3}$ where $\bar{\rho}$  is the comoving
density  of the   background.   The \citet{sheth99b} mass function in
$\mathrm{\Lambda}$CD on cluster mass scales is well approximated by
$\alpha=1/6$.
\end{itemize}

The three different mass  function are plotted in Fig.~\ref{fig:massF}
.  For a test, we have  populated light-cones of $10^6 \mbox{arcsec}^2$ and
measured the mass function in  different redshift bins. The histograms
of  the  halo  masses,  for  which  we  put  together  $64$  different
realizations    of   the    light-cone,    are    also   plotted    in
Fig.~\ref{fig:massF}.  The agreement is very good.

In order  to model the  halo lensing properties  we need to  make some
assumption about the dark matter distribution inside the virial radius
of each system - the virial  radius defines the scale inside which the
overdensity  exceeds   the  critical  value  predicted   for  collapse
\citep{gunn77,bond91,lacey93}.  Numerical simulations predict that the
dark  matter follows  an  NFW  profile \citep{navarro97,navarro04}  or
something very  close to it.   Using this model  we must also  fix the
concentration (the ratio between the scale and the virial radius: $c =
r_s/R$ ) for  each halo.  In {\small GLAMER} \, we  have implemented four model
for the mass-concentration relation (we  refer the reader to the cited
papers  for  more  details):  \citet{zhao09},  \citet{munozcuartas11},
\citet{giocoli12b} and a power-law relation:
\begin{equation}
c =  10 \, \left( \dfrac{M}{10^{12}} \right)^{\alpha} \left[ \dfrac{H_0}{H(z)} \right]^{2/3}
\end{equation}
where the redshift evolution factor is  motivated by the change of the
normalization as seen in cosmological numerical simulations.

This  method of  populating  the  light-cone  with  halos has  several
limitations  with respect  to using  a halos  catalog derived  from an
simulation.  Since the positioning of  the halos is random, we discard
any clustering or particular environmental effects.  Halos do not have
subhalos within  them.  The advantage  is that there is  no resolution
limit and no limit to how many light-cones can be created.

As  in the  halos  from N-body  case, the  baryonic  component can  be
modeled by adding  an additional halo to the center  of each NFW halo.
So far we  have implemented one simple way of  doing this.  The galaxy
is modeled as a NSIE with a  random orientation.  The mass of the NSIE
is   calculated  using   the  stellar   mass/halo  mass   relation  of
\cite{2010ApJ...710..903M}:
\begin{equation}
M_{\rm star} = m_o \frac{ (M_{\rm halo}/M_1)^{\gamma_1}}{\left[ 1 +
    (M_{\rm halo}/M_1)^\beta\right]^{(\gamma_1-\gamma_2)/\beta}}
\end{equation}
with                                 $m_o=7.3113\times10^{10}\,\mbox{M}_\odot$,
$M_1=2.8575\times10^{10}\,\mbox{M}_\odot$,  $\gamma_1=7.17$ ,  $\gamma_2=0.201$
and $\beta=0.557$.   This is  derived from  matching the  abundance of
observed galaxies to the number of  predicted halos in mass bins.  The
NFW halos mass  is reduced accordingly.  Much  more sophisticated ways
of doing this are possible and will be implemented in the future.

The centers  of some halos will  lie outside the light-cone, but they
will  intersect the  cone's  boundary.  This  can  result in  boundary
effects  because  the  density  out  the  boundary  of  the  field  is
artificially low.   A buffering option  is implemented to  avoid this.
All the halos are created in the region outside of the cone but within
a fixed proper distance ($\sim 1 \mbox{Mpc}$) from the boundary.  The result
is a funnel shape instead of a cone.

\subsection{Mass Maps}\label{sec:mass_map}
Maps of  surface density in  FITS format can  be read into  {\small GLAMER} and
placed on  a lens plane at  the desired redshift.   The deflection and
shear are  found by Fourier transforming the  surface density, solving
equation    \eqref{eq:poissons}   for   the    potential,   evaluating
equations~\eqref{eq:shears} and then Fourier transforming back to real
space.  Maps of  shear and deflection are stored so  that this need be
done  only  once.   The  lensing  quantities  are  then  evaluated  by
interpolating from these maps during the ray shooting procedure.

When using a light-cone taken from a large-scale simulation it is very
advantageous to project the particles in shells around the observer onto
planes and then construct two  dimensional density maps by CIC (clouds
in cell) or some other method.  Because the shells can be quite thick,
as we will  see in section~\ref{sec:tests}, this can  save an enormous
amount  of memory  and  computational time  over  using the  particles
directly without significant lose in accuracy.  We have also used this
capability to implant simulated galaxy  clusters into light-cones that
are otherwise populated  with halos.  Fig.~\ref{fig:cluster_map} shows
and example of this.

\subsection{Sources}
Sources can  be distributed  throughout the light-cone.  They  can be
represented by a number of analytic forms including a circular uniform
surface brightness  profile, a \cite{1963BAAA....6...41S}  profile and
an inclined exponential disk profile.  A pixelized image of the source
can  also be  used.  A  catalog  of sources  can  be read  in from  an
external source or individual sources can be created randomly.  {\small GLAMER}
can adaptively  find and refine  around the  image even when  they are
highly distorted  and there  are multiple images  of the  source.  The
rays are shot only back to the redshift of this source.

A numerical effect can occur if  the source is placed in the center of
a halo and  then the halo is projected to the  nearest lens plane that
is at  lower redshift than  the source.  This will  produce alignments
with significant differences in  radial distance and spurious lensing.
The distance between planes should be small enough to make the lensing
by  the  nearest plane  small,  but  the  perfect alignment  can  have
unwanted effects.  To avoid this, the deflection is turned off for the
lens plane that contains the  source's halo.

\begin{figure*}
\centering
\includegraphics[width=1\textwidth]{./Figures/CluterandClusterLSS.pdf}
\caption[]{Convergence  map  of  a  cluster without  (left)  and  with
  (right)  structure along  the line  of sight.   The cluster  lens is
  located at redshift $z=0.3$, it posses a mass of $9.2 \times 10^{14}
  \mbox{M}_\odot$ and  a concentration of  $7$. The  green contours in  the two
  figures  show  the  location   of  radial  and  tangential  critical
  lines. The cluster on the left  possess an Einstein radius of $18.9$
  arcsec while in the right panel,  where the structure along the line
  of   sight   are   included,   the   Einstein   radius   is   $21.3$
  arcsec.
\label{fig:cluster_map}}
 \end{figure*}

\begin{figure}
\centering
\includegraphics[width=0.54\textwidth]{./Figures/np_z.pdf}
\caption[]{The number of planes needed in a simulation versus source redshift. The values are chosen  for separation  distance between  lensing planes  of 300
  \mbox{comoving Mpc}.
\label{fig:np_z}}
\end{figure}

\section{Tests}\label{sec:tests}
We perform several tests to check the accuracy and convergence of our
multi-plane gravitational lensing implementation in the code {\small GLAMER}.
In our first test we determine the
convergence of our code by comparing the  image
positions and magnification of a
four-imaged quasar, as a function of the number of lensing planes.
In Fig.~\ref{fig:pos}  we plot  the position  of the  four images  of the
quasar  at  redshift   $z=3.62$,  lensed  by  a   galaxy  at  redshift
$z=0.42$. The  field of  view of  $10^4 \mbox{arcsec}^2$  with a  buffer of
1~\mbox{Mpc},  filled  with  NFW  halos  from a  ST  mass  function  down  to
$10^9\mbox{M}_\odot$. We compute the lensing properties and change the number
of planes in two plane interval as the main galaxy lens stays fixed on a
plane at redshift $z=0.42$.

The  changes in the  image positions are no  greater than
0.4  \mbox{arcsec}  \, in  one  x-direction  and  1.2  \mbox{arcsec} \,  in  the
y-direction. This  can also  be observed in  Fig.~\ref{fig:all}, where
plot the image  positions and inverse magnifications as  a function of
the distance between  the planes.  The positions show  very small ($<$
0.2 \mbox{arcsec})  scatter for  300  \mbox{comoving Mpc}  \, or  less,  equivalent to  22
lensing planes. The same holds for the inverse magnification.

Our next test aims to examine how other lensing properties, such  as
shear  and  convergence, converge
with increasing the number of lensing planes. For the study  we
generate a light-cone with a field of view of
$10^4 \, \mbox{arcsec}^2$ that extends  to redshift $z=2.0$. We populate it
with NFW  halos from  a ST mass function down  to $10^9\mbox{M}_\odot$.  As a
``true''  result we  consider the  case  where an  infinite number  of
planes  is used,  i.e. each  lensing  halo occupies  a single  lensing
plane. We will denote this results as the ``true'' solution.

We compare the relative error of the lensing properties of simulations
with different  number of planes to  the true solution. To do this we
compare the image maps pixel by pixel and we formulate the
the relative  error  $\epsilon$ in the following way:
\begin{equation}  \label{eq:error}\epsilon =  \frac{1}{N_{\rm
    pixels}}  \sum_{i=0}^{N_{\rm pixels}-1}\frac{x_i  - x_i^{\rm true}}{x^{\rm
    true}_i},  \end{equation}  where  \begin{equation}  N_{\rm  pixels}  =  \sum_{i=0}^{N_{\rm
    pixels}-1} x_i , \end{equation}
and $x_i$ is  the value of the quantity at the
respective pixel $i$.

In Fig.~\ref{fig:mean} we plot the relative error for the convergence,
shear, deflection  angle, and inverse  magnification as a  function of
the comoving separation between the lensing planes. The relative error
for the deflection  angle and the inverse magnification  is very small
--  for  less than  300  \mbox{comoving Mpc}  \, between  the  planes  it is  around
0.1\%. The  error for  the convergence and  shear is ten  times larger
 and drops  down to less  than 5\%
for  less than  300 \mbox{comoving Mpc},  which  supports the  conclusions from  our
previous test. This separation corresponds  to 18 lensing planes for a
source at redshift $z=2$.

In the next test, we  investigate the agreement of our distributions of
the  lensing quantities  with already  published results.   We compute
histograms of the lensing properties of  a field of view with the same
parameters as the  one from the previous test (field  of view of $10^4
\,  \mbox{arcsec}^2$, extending  to  redshift $z=2.0$,  populated with  NFW
halos from  a ST mass function  down to $10^9\mbox{M}_\odot$),  but keeping the
number of planes  fixed to 20. We change  randomly the distribution of
the NFW  halos in the light-cone and  obtain a set of  32 different and
uncorrelated fields of view. In Fig.~\ref{fig:pdf} we show the average
of  the   convergence,  shear,  and   deflection  angle  for   the  32
realizations  of  the  image   plane  and  we  overplot  the  standard
deviation.   The  results  agree   very  well   with  the   ones  from
\citet{Pace2011}, despite  the much  smaller field of  view.

Finally,  we show  how structure  along  the line  of sight  may change  the
Einstein radius of  a galaxy cluster. In our simulation  we define the
Einstein  radius as  the  median distance  of  the tangential  central
critical points from the cluster center \citet{meneghetti10a}.  In the
left panel of Fig.~\ref{fig:cluster_map} we  show convergence map of a
halo  created   with  {\small  MOKA}  \citep{giocoli12}   at  redshift
$z=0.3$.  The cluster  possess a  virial mass  of $9.2  \times 10^{14}
\mbox{M}_\odot$ and a  concentration of $7$.  In the right  panel we include in
the map the effect of uncorrelated structures along the line of sight,
assuming an  empty light-cone. The  green lines in the  maps represent
the radial and tangential critical  curves, where the magnification of
the background  sources, located at  $z_z=2$, is infinite.   The green
connected  points in  the center  of  the maps  indicate the  Einstein
radius  of  the  cluster,  measured  as the  median  distance  of  the
tangential  central  critical points  to  the  halo center.  From  the
analysis we notice  that while the Einstein radius in  the left figure
is $18.8$ arcsec,  in the right increases to $21.3$  arcsec, where LSS
are  included.  We  will  discuss motivate  and  better quantify  this
effect in a future paper Giocoli et al. in preparation.

We have  also performed further small  test to ensure  the accuracy of
our multi-plane lensing  scheme.  For example, we have  tested that the
average  lensing mass  on each  plane is  zero and that  an empty
light-cone with an arbitrary  number of lensing planes always produces
zero deflection.  We also checked that an  analytic lens placed  in an empty  light-cone
with an arbitrary  number of planes agrees fully  with the single lens
plane calculation.

From our tests we  conclude that a separation of 300  \mbox{comoving Mpc} \, or less
is sufficient  for accurate  representation of  the lensing  field. In
Fig.~\ref{fig:np_z}  we  show  the correspondence  between  number  of
planes and  source redshift,  that gives  approximately 300  \mbox{comoving Mpc} \,
plane separation.   For redshift $z=2$,  18 planes are  sufficient. For
much higher redshifts, $z=10$, the number increases to 32.

\section{Some applications}
\label{sec:some-applications}

In this section,  we present some example  applications to demonstrate
{\small GLAMER}'s applicability  to interface with  other lensing codes  and to
produce simulated  observations.  We  do this without  fully describing
the parts  of the  simulations that  are not  directly related  to the
functionality of {\small GLAMER}.

As   previously   discussed,   Fig.~\ref{fig:cluster_map}   shows   the
convergence  map for  a simulated  galaxy cluster.   The cluster  lens
located    at    redshift    $z=0.3$   was    created    using    {\small MOKA}
\citep{2012MNRAS.421.3343G},  and has  as input  parameters a  mass of
$9.2\times 10^{14}M_{\odot}/h$ and a concentration of $7$, and sources
are  located at  redshift $z_s=2$.   The axial  ratios of  the cluster
ellipsoid are  randomly drown  from the  \citet{jing02} distributions,
they are $a/b=0.62$ and $a/c=0.3$ where $a$, $b$ and $c$ represent the
smallest,  intermediate  and  largest  axis  of  the  halo  ellipsoid,
respectively; the random  orientation of the systems,  with respect to
the line of sight, gives in the plane of the sky an axial ratio of the
ellipse that is  $a'/b'=0.32$. As a result the cluster  on the left is
characterized by  an Einstein radius $\theta_E=18.8$  arcsec. In {\small GLAMER}
we  have  created  an  interface   that  allows  {\small MOKA}  outputs  to  be
reprocessed   including   background   and   foreground   objects   by
constructing    a    random    light-cone    as    described    in
section~\ref{sec:massfunc_lc}.  The {\small MOKA} map  was in this way inserted
into the  light-cone.  As in  all these examples, the  ray tracing was
done without  the Born  approximation.  Similar map  of the  shear and
deflection are created.  The critical lines and caustics can be mapped
and the images of lensed background  objects can be found. The created
light-cone realization  modify the  size  of the  Einstein radius  to
$\theta_E=21.3$ arcsec.

Fig.~\ref{fig:invmag_map}  shows   a  map   of  the  inverse   of  the
magnification for a random 8.4 by 8.4 arcminute field.  The light-cone
was  constructed from  a halo  catalog extracted  from the  Millennium
simulation.  Each  halo is represented by  a NFW profile for  the dark
matter and  a NSIE profile for  the inner regions.  The  grid in which
the rays have been shot was  refined to give higher angular resolution
where the magnification is high.

Fig.~\ref{fig:12lens_images}  shows   twelve  examples   of  simulated
galaxy-galaxy  strong lenses.   The ray  shooting grid  is refined  to
adapt to  the surface brightness  of the  images.  Each source  is ray
traced back  to its own  redshift and has  a NFW+NSIE halo.   Here the
caustics in a random light-cone were found and an extra source was put
close to  the caustic.  Many thousands  of such lenses can  be created
easily.  Such  images are being  used to  test lens finding  robots, lens
modeling algorithms and predict the statistics of detectable lenses.

\begin{figure*}
\centering
\includegraphics[width=0.8\textwidth]{./Figures/invmag.png}
\caption[]{Map of the  inverse of the magnification,  $1/|\mu|$, for a
  8.4 by 8.4 arcmin random field with the sources plane at z=2.5.  The
  image  is on  a regular  grid,  but the  rays have  been refined  to
  resolve the critical curves so the high magnification ( and the $\mu
  <  0$)  regions   had  higher  resolution  before   this  image  was
  interpolated to a  regular grid.  In some of the  regions where $\mu
  \simeq 1$  the resolution  of the  adaptive grid  is lower  then the
  final image resolution. \label{fig:invmag_map}}\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=0.8\textwidth]{./Figures/12lenses.png}
\caption[]{Twelve simulated galaxy-galaxy lenses. The mass
  distributions of the galaxies and their dark matter halos are taken
  to be NSIE within NFW profiles.  The visible galaxies are modeled
  as \cite{1963BAAA....6...41S} profiles at the centers of the
  halos.  Noise has been added and a Gaussian point spread function
  (psf) applied.  Two of these are the same lens with different
  sources behind it.  The output images are in standard fits files.
\label{fig:12lens_images}}\end{figure*}

\section{discussion and conclusions}\label{sec:conclusions}
We  have  presented an  extension  of  the ray-shooting  lensing  code
{\small GLAMER} to multiple planes  and its ability to perform lensing
simulations in  a light-cone,  rather than  using a single  lens. This
allows  us  to perform  more  accurate  statistical strong  and  weak
lensing simulations and study properties of the strong lensing galaxies
and galaxy clusters, as well as the structure along the line of sight
outside of the primary lens.

We   have  presented   in   detail  the   methodology  and   numerical
implementation that we have use. The mass in the universe, in the form
of halos,  is projected onto a  discrete number of lensing  planes and
the  lensing  properties,  such  as deflection  and  convergence,  are
computed on each  of these planes and summed up  according to the lens
equations. The summation of the  contributions from the different halos
is performed via a spacial tree algorithm in  order to improve
accuracy and speed.  Details of  the performance and limitations of the
tree can be found in MP.

The halos,  populating a light-cone,  can be either taken  from a dark
matter  simulation  or randomly  generated  from  a mass  function  at
different  redshifts.  For the
second  one  we  have  adopted   three  different  mass  functions  --
Press-Schechter, Sheth- Tormen and a simple power law. In this context
more general  and useful relations  can be implemented and  inputed as
desired.

We have presented various tests and  applications of our code. We have
shown that the optimal number of planes per light-cone is achieved for
a separation  of circa 300  \mbox{comoving Mpc}.  We  have presented that  the image
position of  a lensed  quasar by  a galaxy  converges for  this value.
Moreover, the relative uncertainty between the lensing properties in a
simulation with that  given plane separation and  an infinitesimal one
is  less  than   3\%.   We  have  also  shown   that  the  statistical
distribution for the lensing properties  in our simulations agree very
well with previously published results.
Finally, we have demonstrated some of the applications that the code is
already being applied to.

\section*{Acknowledgments}
This research is part of the project GLENCO, funded
under the Seventh Framework Programme, Ideas, Grant
Agreement n. 259349.

\bibliographystyle{mn2e}
\input{paper2.bbl}