\section{INTRODUCTION}

Many novel materials find themselves close to the metal-insulator transition (MIT), where strong electron-electron interactions and disorder are believed to give rise to a variety of complex phenomena.\cite{Vlad-bookintro}  Near the MIT, these two effects are usually comparable in magnitude, so that their competition is expected to lead to glassy behavior of electrons, in analogy with other frustrated systems.\cite{Eduardo-Vlad_review,Vlad-bookchapter}  A common denominator for all glasses is the existence of a complex or ``rugged'' free energy landscape (FEL), consisting of a large number of metastable states, separated by barriers of different heights.  This results in phenomena such as slow, nonexponential relaxations, divergence of the equilibration time, and breaking of ergodicity, i.e. the inability of the system to equilibrate on experimental time scales.  Therefore, such out-of-equilibrium systems also exhibit aging effects,\cite{Struik,aging} where the response to an external excitation (i.e. relaxation) depends on the system history in addition to the time $t$.  A detailed analysis of temporal fluctuations (noise) of the relevant observables yields complementary information on configurational rearrangements or transitions between metastable states.  Non-Gaussian distributions of various observables in glassy systems have been reported,\cite{hetero-review} reflecting the presence of large, \textit{collective} rearrangements.

Most experimental studies of charge or Coulomb glasses have focused on situations where electrons are strongly localized due to disorder, far from the MIT.\cite{Amir-eglass_review}  In recent years, however, a two-dimensional electron system (2DES) in Si metal-oxide-semiconductor field-effect transistors (MOSFETs) has emerged as an excellent model system for studying glassy or out-of-equilibrium charge dynamics near the MIT.\cite{DP-bookchapter}  In particular, studies of both relaxations and noise on very disordered samples have established that the 2DES in Si exhibits all the main manifestations of glassiness: slow, correlated dynamics (non-Gaussian noise),\cite{Snezana2002,Jan2002} nonexponential relaxations,\cite{Jan2006} diverging equilibration time (as temperature $T\rightarrow 0$),\cite{Jan2006} aging and memory.\cite{Jan2007-1,Jan2007-2,Jan2009}  Glassiness is observed for all carrier densities $n_s<n_g$, where $n_g>n_c$ ($n_c$ -- the critical density for the MIT), thus giving rise to an intermediate glassy, ``bad metal'' phase ($k_{F}l <1$; $k_{F}$ -- Fermi wave vector, $l$ -- mean free path).  In a 2DES with a relatively low amount of disorder, noise studies show that the intermediate phase becomes vanishingly small, i.e. $n_g\gtrsim n_c$.\cite{Jan2002,Jan2004}  These observations are consistent with predictions of the theoretical models that describe the MIT as a Mott transition with disorder.\cite{Darko}  The results presented here were obtained on samples with a relatively large amount of disorder, where the intermediate ($n_c<n_s<n_g$) glassy phase is more pronounced, as mentioned above.  We note that aging studies have found\cite{Jan2007-2,Jan2009} that an abrupt change in the nature of the glassy phase occurs at $n_c$, i.e. at the 2D MIT, signifying the intrinsic glassiness of the 2DES.

Glassy dynamics in a 2DES in Si has been studied using several different experimental protocols, all of which involved rapid changes of $n_s$ (``$\Delta n_s$ protocol'').  For example, a large, rapid (within 1~s) change of $n_s$ at a low enough $T$ results in a nonexponential relaxation of conductivity $\sigma(t)$ towards a new equilibrium value.\cite{Jan2006}  Small changes of $n_s$ at low $T$ do not produce any observable relaxations,\cite{Jan2007-2} but have led to the observations of non-Gaussian fluctuations in $\sigma(t)$.\cite{Snezana2002,Jan2002,Jan2004}  On the other hand, in other types of glasses (e.g. structural, spin), studies of aging processes and other glassy properties are typically carried out following a thermal quench from an equilibrium state at a high $T$ to a nonequilibrium state at a low $T$.  The analysis of the relaxations following such a thermal quench, however, is greatly complicated by their nontrivial dependence on the cooling time and, hence, a fairly fast cooling process is essential (see, e.g., Ref.~\onlinecite{SG-fullaging}).  It is of fundamental interest to establish whether the charge glass realized in a 2DES in Si also falls out of equilibrium as $T$ is reduced, in common to other types of glasses.

Here we investigate the charge dynamics in a 2DES in Si at low $n_s$, near the MIT, after cooling from an equilibrium state at a high $T$ with a fixed $n_s$ (``$\Delta T$ protocol'').  Since it is not possible to achieve rapid cooling in our experimental set-up (the shortest cooling time is about 30 minutes), we focus on the resulting fluctuations (noise) in $\sigma(t)$.  After describing the samples and measurement techniques in more detail (Sec.~\ref{expt}), we present the analysis of the noise power spectra obtained using a $\Delta T$ protocol (Sec.~\ref{noisepower}).  The results demonstrate the glassy arrest of the charge dynamics for $n_s<n_g$, where the glass transition density $n_g$ is the same as that found in noise studies that employed $\Delta n_s$ protocols.\cite{Snezana2002,Jan2002}

The power spectrum, being the Fourier transform of a temporal correlation function, only provides information about the second moment of the fluctuating quantity.  Additional information can be extracted from the full probability distribution of the fluctuations.  Therefore, in Sec.~\ref{pdf}, we compare the probability density functions (PDFs) obtained in the two protocols, $\Delta n_s$ and $\Delta T$.  While in both cases the PDFs are non-Gaussian for $n_s<n_g$, the change of the carrier density has a qualitatively different and more dramatic effect than $\Delta T$.  In particular, the results suggest that $\Delta n_s$ induces strong changes in the free energy landscape of the system, which should have important implications for theoretical modeling of the glassy dynamics in a 2DES.

For a $\Delta T$ protocol, the power spectra and the second spectra\cite{Weissman88,Weissman92,Weissman93} of the noise in the glassy regime have also been investigated for different values of the excitation voltage $V_{\mathrm{exc}}$ used in the measurements of $\sigma$.  In indium oxide, another realization of a Coulomb glass, the non-Gaussian noise was attributed to a non-linear mechanism caused by the high applied bias.\cite{Zvi-nonlinear}  As shown in Sec.~\ref{exc}, in the 2DES in Si the large, non-Gaussian noise does not depend on the applied bias over more than two orders of magnitude variation of $V_{\mathrm{exc}}$, providing additional evidence that it reflects the intrinsic glassiness of the 2DES.  The main results are summarized and discussed further in Sec.~\ref{concl}.

\section{EXPERIMENT}
\label{expt}

\subsection{Samples and measurement techniques}

Measurements were carried out on rectangular n-channel (100)-Si MOSFETs with a relatively large amount of disorder.  Similar to previous studies,\cite{Snezana2002,Jan2006,Jan2007-1,Jan2007-2,Jan2009} the back-gate bias of $-2$~V was applied to maximize the 4.2~K peak mobility to $\sim$0.06~m$^2$/Vs.  The sample dimensions $L \times W$ ($L$--length, $W$--width) were $1 \times 90$~$\mu$m$^2$ for samples A and A1, and $2 \times 50$~$\mu$m$^2$ for samples B and B1.  The devices were fabricated with poly-Si gates, self-aligned ion-implanted contacts, substrate doping $N_a \sim 2 \times 10^{17}$~cm$^{-3}$, oxide charge $N_{\mathrm{ox}} \sim 1 \times 10^{11}$~cm$^{-2}$, and oxide thickness $d_{\mathrm{ox}}=50$~nm.  All four samples exhibited the same behavior, and some of them were, in fact, used in previous studies (e.g. A and B in Refs.~\onlinecite{Jan2006,Jan2007-1,Jan2007-2,Jan2009}, A1 in Ref.~\onlinecite{Snezana2002}).

The conductivity $\sigma$ was measured using a standard two-probe ac technique (typically at $\sim 11$ or 13~Hz) with an ITHACO 1211 current preamplifier, and either a PAR124A or a SR7265 lock-in amplifier.  The contact resistances and the contact noise were determined to be negligible in these samples.\cite{Snezana2002}
Unless noted otherwise (see Sec.~\ref{exc}), the excitation voltage $V_{\mathrm{exc}}$ was kept constant and low enough ($2-10~\mu$V) to ensure that the conduction was Ohmic.  The low-frequency of $V_{\mathrm{exc}}$ was limited by the low cut-off frequency of RC filters, which were used to reduce external electromagnetic noise, as well as by high resistance of the samples.  A precision dc voltage standard (EDC MV116J) was used to apply the gate voltage $V_g$, which controls $n_s$: $n_s(10^{11}$cm$^{-2})=4.31(V_g[$V$]-6.3)$.  Measurements of current ($I$) fluctuations as low as $10^{-13}$~A were performed both in a dilution refrigerator and in a $^3$He system (base $T=0.24$~K).  Relatively small fluctuations of $T$, $V_g$, and $V_{\mathrm{exc}}$ were ruled out as possible sources of the measured noise, since no correlation was found between them and the fluctuations of $I$.  The background noise was measured for all $n_s$ and $T$ by setting $V_{\mathrm{exc}}$ to zero.  It was always independent of frequency (i.e. white) and usually several orders of magnitude smaller than the sample noise.

\subsection{Temperature dependence of the time-averaged conductivity}

Figure~\ref{fig:RvT} shows the time-averaged (see Sec.~\ref{noisepower}) conductivity $\langle \sigma \rangle$ as a function of $T$ for
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig1.eps}
\caption {(Color online) Sample B.  Conductivity $\langle\sigma\rangle$ \textit{vs.} $T$ for different $n_s$, as shown.  $n_s$ was varied at high $T\approx 10-12$~K.  Solid lines guide the eye. }\label{fig:RvT}
\end{figure}
different $n_s$ for one of the samples.  The density was varied at $T\approx 10-12$~K, where the 2DES is in a thermal equilibrium.\cite{Jan2006}  The critical density $n_c$ is determined from the fits to $\langle\sigma(n_s, T)\rangle$ on both metallic and insulating sides of the MIT, as discussed in detail elsewhere.\cite{Snezana2002,Jan2002,Jan2004,Jan2009}  In particular, in the insulating regime, the transport is thermally activated, whereas in the intermediate, $n_c<n_s<n_g$ glassy phase, $\langle\sigma(n_s, T)\rangle$ acquires a very specific, non-Fermi liquid $T^{3/2}$ correction.  In all samples, $n_c=(4.5\pm 0.4)\times 10^{11}$~cm$^{-2}$, and the glass transition density determined from $\Delta n_s$ protocols\cite{Snezana2002,Jan2006,Jan2007-1,Jan2007-2,Jan2009} was $n_g=(7.5\pm 0.3)\times 10^{11}$~cm$^{-2}$.  The error bars account for sample-to-sample variations.

\section{Conductance noise from a $\Delta T$ protocol}
\label{noisepower}

In a $\Delta T$ protocol, each noise measurement starts by fixing the carrier density $n_s$ at $T\approx 10-12$~K.  The sample is then cooled to the measurement temperature.  The cooling times varied from 30 minutes to 10 hours, with no effect on the noise results.  In addition, such slow cooling did not result in any visible relaxations, at least on our experimental time scales. Typical data are presented in Fig.~\ref{fig:fig2}(a), which shows the relative fluctuations $\Delta\sigma(t)/\langle\sigma\rangle$ vs time for several $n_s$ measured at $T=0.24$~K.
\begin{figure}
\includegraphics[width=7.1cm]{Lin-noise-resub_fig2a.eps}\\
\includegraphics[width=7.1cm]{Lin-noise-resub_fig2b.eps}
\caption {(Color online)  Sample B.  (a) Relative fluctuations $\Delta\sigma(t)/\langle\sigma\rangle=(\sigma(t)-\langle\sigma\rangle)/\langle\sigma\rangle$ vs time for several $n_s$ shown on the graph (top to bottom) at $T=0.24$~K. The averaging of $\sigma(t)$ was done over a time interval shown on the plot.  Traces are shifted for clarity.  $n_s$ was changed at $T\approx 10-12$~K and the cooling time was about $10^4$~s.  (b) In order to make the change in the character of the noise with $n_s$ more apparent, the same fluctuations $\Delta\sigma(t)=(\sigma(t)-\langle\sigma\rangle)$ are shown normalized by $\delta\sigma$, where $\delta\sigma=\langle(\sigma-\langle\sigma\rangle)^{2}\rangle^{1/2}$.
}\label{fig:fig2}
\end{figure}
It is obvious that, by decreasing $n_s$, the amplitude of the fluctuations increases dramatically.  Even more striking is the change in the character of the noise: at low $n_s$, fluctuations occur on many different time scales, including slow changes over periods of several hours; at high $n_s$, on the other hand, the variance of the noise no longer varies with time.  These features are more apparent in Fig.~\ref{fig:fig2}(b), where we choose to plot the same $\Delta\sigma(t)=(\sigma(t)-\langle\sigma\rangle)$ normalized by $\delta\sigma$, the square root of the variance of the fluctuations.  The increase in the noise magnitude with decreasing $n_s$ is clearly accompanied by the onset of non-Gaussian behavior, similar to the noise results obtained from a $\Delta n_s$ protocol.\cite{Snezana2002,Jan2002,Jan2004}
\begin{figure}
\includegraphics[width=7.1cm]{Lin-noise_fig3a.eps}
\includegraphics[width=7.1cm]{Lin-noise_fig3b.eps}
\caption {(Color online)  Open symbols: sample A1; solid symbols: sample B.  Different symbols correspond to different $T$, as shown, where $n_s$ was changed before cooling down to 0.24~K to measure noise.  (a) The normalized noise power $S_{I}(f=1$~mHz) vs $n_s$.  (b) The exponent $\alpha$ vs $n_s$.}\label{fig:Salpha}
\end{figure}

The normalized power spectra of the relative changes in conductivity $\Delta\sigma(t)/\langle\sigma\rangle$ are found to obey the empirical law $S_I \propto 1/f^{\alpha}$ ($f$ -- frequency).  In order to compare the noise magnitudes under different conditions, the spectra were averaged over octaves, and the resulting fraction of power $S_I(f=1~\mathrm {mHz})$ is taken as the measure of noise.  Figure~\ref{fig:Salpha}(a) shows that $S_I(f=1~\mathrm {mHz})$ increases by several orders of magnitude as $n_s$ is reduced below $\sim 7\times 10^{11}$~cm$^{-2}$, while $\alpha$ rises rapidly to a much higher value [Fig.~\ref{fig:Salpha}(b)] in the same range of $n_s$.  Both effects indicate that the electron dynamics suddenly and dramatically slows down for $n_s<n_g\sim 7\times 10^{11}$~cm$^{-2}$, which is precisely the regime where non-Gaussian behavior of the noise is also observed (e.g. Fig.~\ref{fig:fig2}).  Therefore, the results of the noise measurements obtained using a $\Delta T$ protocol (black triangles in Fig.~\ref{fig:Salpha}) also provide evidence for glassy freezing of the 2DES for $n_s<n_g$, with the value of $n_g$ in agreement with that obtained from $\Delta n_s$ protocols.

Furthermore, the values of $S_I(f=1~\mathrm {mHz},n_s)$ and $\alpha(n_s)$ have been compared to
those from a $\Delta n_s$ protocol (see, e.g., Fig.~\ref{fig:protocols}).  In the latter, $n_s$ was changed
\begin{figure}
\includegraphics[width=7.1cm]{Lin-noise_fig4a.eps}
\includegraphics[width=7.1cm]{Lin-noise_fig4b.eps}
\includegraphics[width=7.1cm]{Lin-noise_fig4c.eps}
\caption {(Color online) Sample B.  Conductivity $\sigma$ vs time at $T=0.24$~K for several $n_s$.  Traces are shifted for clarity.  The data were obtained (a) after changing $n_s$ at 10~K, followed by a slow decrease of $T$
for about 12 hours ($\Delta T$ protocol), and (b) after a small $n_s$ change ($\Delta n_s\approx 0.22\times 10^{11}$cm$^{-2}$) at low $T=0.24$~K ($\Delta n_s$ protocol).  In another run, (c) $n_s$ was changed at 1~K ($\Delta n_s$ protocol); after staying at 1~K for two hours, $T$ was lowered to 0.24~K in one hour, and $\sigma(t)$ was then measured.}\label{fig:protocols}
\end{figure}
in small steps (typically $\sim 0.2\times 10^{11}$~cm$^{-2}$) at low $T$, where the equilibration times exceed the experimental time scales.\cite{Jan2006}  Figure~\ref{fig:Salpha} includes the results for $S_I$ and $\alpha$ obtained by changing $n_s$ at 1~K, 0.8~K (as in Ref.~\onlinecite{Snezana2002}) and 0.24~K.  Both the noise magnitude and the exponent $\alpha$ demonstrate that $\Delta T$ and \emph{small} $\Delta n_s$ protocols produce the same results.

While both protocols, cooling with a fixed $n_s$ and a small change of $n_s$ at low $T$, provide evidence for the onset of glassy dynamics from the fluctuations of $\sigma(t)$, we note that there are no visible relaxations of $\sigma$ following either process.  In both cases, though, the size of the perturbation ($k_{B}\Delta T\sim 10$~K or $\Delta E_{F}\sim 1-2$~K, respectively, where $E_F$ is the Fermi energy) is much smaller than $E_F$ ($\sim 20-60$~K),\cite{Ef-comment} i.e. probably too small to lead to any observable relaxations.  The relaxations were indeed observed for much larger changes of $n_s$,\cite{Jan2006,Jan2007-1,Jan2007-2} such that $\Delta E_{F}\sim 30-150$~K, i.e. for $\Delta E_{F}> E_F\gg k_{B}T$.

Since $\Delta E_{F}>k_{B}T$ in all $\Delta n_s$ protocols on a 2DES in Si MOSFETs, the perturbative treatments of the effect of the gate voltage,\cite{Markus2005-1,Markus2005-2} which may be relevant for other materials, are not suitable for this system.  In particular, here the applied $\Delta V_g$ (i.e. $\Delta E_{F}$) are expected to trigger major, \emph{collective} rearrangements of the electron configuration.  As shown below (Sec.~\ref{pdf-small}), a comparison of the full distribution of the fluctuations obtained from the two protocols demonstrates that even a small $\Delta n_s$ ($\Delta E_F\ll E_F$) has a qualitatively different effect than $\Delta T$, even though the results for the noise power spectra (Fig.~\ref{fig:Salpha}) and $\langle\sigma(n_s,T)\rangle$ are the same within experimental error.

\section{Distribution functions of the fluctuations}
\label{pdf}

In systems that are in equilibrium and away from criticality, the probability density functions (PDFs) of fluctuations in global quantities are Gaussian, according to the central limit theorem, which is based on the hypothesis that the system may be decomposed into many uncorrelated elements.  When this condition is not satisfied, e.g. due to the divergence of the correlation length at criticality or slow, nonexponential decay of correlation functions in glassy systems, non-Gaussian PDFs have been observed in a number of disparate systems.

\subsection{PDFs from $\Delta n_s$ and $\Delta T$ protocols}
\label{pdf-small}

Figure~\ref{fig:pdf} shows some typical PDFs of the conductance ($G$) fluctuations for $n_s<n_g$ measured in a $\Delta n_s$ protocol.  The density $n_s$, lying in the metallic glass regime ($n_c<n_{s}<n_g$), was fixed at
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig5a.eps}
\includegraphics[width=8.1cm]{Lin-noise_fig5b.eps}
\caption {(Color online) Sample A1.  Probability density function (PDF) of the conductance ($G$) fluctuations for $n_c<n_{s}(10^{11}$cm$^{-2})=5.58<n_g$ at $T=0.13$~K.  $n_s$ was fixed at $T\approx 0.8$~K after a small change $\Delta n_s=0.43\times 10^{11}$cm$^{-2}$, corresponding to $\Delta E_{F}\sim 3$~K$\ll E_{F}\approx 41$~K.  The measurement time was divided into sixteen intervals of about 35 minutes each.  PDFs for intervals 1, 4, and 12 are shown in (a), demonstrating the non-Gaussian behavior of the noise. (b)  The color map of the PDF of the fluctuations in $G$ as a function of all sixteen sequential time intervals.  The averaging $\langle\ldots\rangle$ was performed over the total measurement time ($\sim 9$~hours).}\label{fig:pdf}
\end{figure}
$T\approx 0.8$~K after a small change $\Delta n_s=0.43\times 10^{11}$cm$^{-2}$ ($\Delta E_{F}\sim 3$~K$\ll E_{F}\approx 41$~K).  The sample was then cooled down to $T=0.13$~K, and $G(t)$ was recorded over sixteen sequential intervals of about 35 minutes each.  The PDFs, shown in Fig.~\ref{fig:pdf}(a) for three out of sixteen intervals, are clearly not Gaussian, but rather have a complicated, multi-peaked structure that changes with time.

These results indicate that, on experimental time scales, the 2DES ``visits'' only a small number of states, but that new states also become available with time.  This is also reflected in the correlated ``wandering'' of the PDFs with time [Fig.~\ref{fig:pdf}(b)].  For $n_s>n_g$, on the other hand, the PDFs become Gaussian (not shown) on much shorter time scales.

The PDFs obtained from a $\Delta T$ protocol (Fig.~\ref{fig:pdfT}, left column) exhibit the same general behavior:
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig6.eps}
\caption {(Color online)  Sample A.
PDF vs $[\sigma(t)-\langle\sigma(t)\rangle]/\langle\sigma_0\rangle$ at $T=0.24$~K for several $n_0$, as shown.  $\langle\sigma_0\rangle$ is the time-averaged conductivity corresponding to $n_0$ at the measurement $T$, obtained from a $\Delta T$ protocol.  The red curves are fits to a Gaussian distribution.  Left column: PDFs of the noise after cooling from 10~K to 0.24~K for each given $n_0$ ($\Delta T$ protocol).  The cooling time was at least one hour long. In this protocol, $\langle\sigma(t)\rangle\equiv\langle\sigma_0\rangle$, i.e. there are no observable relaxations after cooling.  Right column: PDFs of the noise measured after a subsequent change of $n_s$ from $n_0$ to a much higher value $n_1=20.26\times 10^{11}$cm$^{-2}$ during $t_w=1000$~s.  In this protocol, $\langle\sigma(t)\rangle$ describes the slowly relaxing background.}\label {fig:pdfT}
\end{figure}
they are non-Gaussian for $n_s<n_g$ and Gaussian for $n_s>n_g$.  (Even though the non-Gaussian shape of the PDFs for $n_s<n_g$ is more obvious on a semi-log scale, they are presented on a linear scale for easier comparison to the effects of $\Delta n_s$.  See also Sec.~\ref{exc}.) Here, however, even the non-Gaussian PDFs are smooth, single-peaked functions, reminiscent of PDFs in a variety of systems displaying critical,\cite{criticalPDF1,criticalPDF2} glassy,\cite{hetero-review} or other out-of-equilibrium behavior (e.g. the Danube water level\cite{Danube}).  In all these systems, the PDFs are skewed, resembling a zero-centered Gaussian, which describes pseudo-equilibrium fluctuations, with one (exponential) tail that is due to large, rare events.
In that context, the totally irregular structure of the PDFs obtained after a small $\Delta n_s$ (Fig.~\ref{fig:pdf}) is striking.

The above differences can be understood in the following way.  After cooling, the 2D charge glass settles in a deep energy valley, perhaps a global minimum of the FEL, corresponding to a given $n_s$.  As a function of time, the system will explore the subvalleys of the energy landscape, but it will in general remain in the ``main'', deep valley.  Such pseudo-equilibrium fluctuations will produce PDFs of the type shown in the left column of Fig.~\ref{fig:pdfT}.  In the $\Delta n_s$ protocol, in contrast, a change in the number of charge carriers reshuffles all energies because of the Coulomb interactions, so that the 2DES in general no longer finds itself in the minimum of the FEL.  This will result in a more dramatic wandering of the system through the FEL, producing the PDFs of the type shown in Fig.~\ref{fig:pdf}.

\subsection{PDFs from a waiting-time protocol}
\label{pdf-large}

Large $\Delta n_s$ ($\Delta E_F > E_F$) lead to observable relaxations in $\sigma (t)$.  This effect has been used to investigate aging and memory in a 2DES in Si by employing the so-called waiting-time ($t_w$) protocol,\cite{Jan2007-1,Jan2007-2,Jan2009} which consists of the following.  It starts with a $\Delta T$ protocol, i.e. cooling at a fixed density $n_0$ and measuring noise (see left column in Fig.~\ref{fig:pdfT} for the corresponding PDFs).  The density is then switched rapidly (within 1 s) to a different value $n_1$, where it is kept for a time $t_w$.  During that time, the system relaxes away from its initial (pseudo-)equilibrium state determined by $n_0$ and towards a new equilibrium state determined by $n_1$.  Finally, the density is changed back (within 1 s) to $n_0$ at $t=0$, and the slowly evolving $\sigma(t)$ is measured.  We note that here $t=0$ is defined as the time when the charge carrier density reattains its original value $n_0$.  If $t_w$ is shorter than the equilibration time, then $\sigma (t)$ will depend on $t_w$, i.e. on the measurement history (aging efect).  It is also said that the system has a memory of the time it spent with $n_1$.  Here we focus on that situation, i.e. on noise during aging.

Figure~\ref{fig:pdfT} presents a comparison of the PDFs obtained from $\Delta T$ (left column) and $t_w$ (right column) protocols for different $n_0$ and the same $n_1=20.26\times 10^{11}$cm$^{-2}$ and $t_w=1000$~s.  For the highest $n_0\gtrsim n_g$, the PDFs are Gaussian in both protocols, with roughly the same standard deviation.  For lower $n_0$, however, there are some remarkable differences.  First, the noise after a temporary density change is much larger than the noise before the change, that is, after a $\Delta T$ protocol.  At the lowest $n_0$, this difference amounts to more than one order of magnitude, as first noted in Ref.~\onlinecite{Jan2007-2}.  Second, the change of density during $t_w$ results in complex, multi-peaked, totally random-looking PDFs, analogous to those obtained from $\Delta n_s$ protocols (Sec.~\ref{pdf-small}).  These differences are even more striking considering that, for a given $n_0$, the PDFs in Fig.~\ref{fig:pdfT} were obtained from the measurements that were carried out under \emph{exactly the same} experimental conditions; the only difference was in the sample history.

These results imply the existence of a rugged free energy landscape that is modified in a non-trivial way by all $\Delta n_s$ of practical interest (corresponding to $\Delta E_F >k_{B}T$).  The rearrangements of the electron configuration and the reshuffling of the FEL caused by Coulomb interactions will clearly need to be taken into account in theoretical models of the glassy dynamics of the 2DES in Si.  This complicated problem is further aggravated by the proximity to the MIT, where substantial changes in the screening of the 2DES may be expected.\cite{Pastor,Markus2005-PRL,Pankov2005}

The large noise that accompanies the aging process reflects the collective wandering of the 2DES through the FEL, as the system relaxes towards a global minimum that corresponds to the equilibrium state.  In other out-of-equilibrium systems (e.g. polymers and colloidal glasses,\cite{Ciliberto-PDF}, conventional spin glasses\cite{SG-PDF}), the PDFs have been observed to evolve towards a Gaussian shape as the system ages.  Here the evolution of the PDFs has been investigated from the data obtained in several consecutive time intervals, during the relaxation of $\sigma(t)$ induced using a $t_w$ protocol, as described above.  A typical behavior is presented in Fig.~\ref{fig:pdftw}, which
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig7a.eps}
\includegraphics[width=8.1cm]{Lin-noise_fig7b.eps}
\includegraphics[width=8.1cm]{Lin-noise_fig7c.eps}
\caption {(Color online)  Sample A.  PDF \textit{vs.} $[\sigma(t)-\langle\sigma(t)\rangle]/\langle\sigma_0\rangle$ at $T=1$~K for $n_0=3.02\times 10^{11}$cm$^{-2}$, following a change of $n_s$ to $n_1=20.26\times 10^{11}$cm$^{-2}$ during (a) $t_w=10$~s, (b) $t_w=1000$~s, and (c) $t_w=31600$~s. For each $t_w$, $n_0$ was initially set at 10~K.  Different curves correspond to PDFs for four consecutive 3-hour time intervals, as shown.  The PDFs become more Gaussian with time, as the system ages.}\label{fig:pdftw}
\end{figure}
shows the PDFs for the same density $n_0=3.02\times 10^{11}$cm$^{-2}$ in the insulating regime ($n_0<n_c$) and for three different $t_w$.  In each case, the PDFs for four consecutive 3-hour time intervals demonstrate that the noise indeed becomes more Gaussian with time, as the system slowly approaches equilibrium.

Understanding of the age dependence of the PDFs
is a non-trivial problem in general.  Even in situations where PDFs feature a ``simple'' zero-centered Gaussian with an exponential tail, different models have been proposed.  In addition, numerical simulations of aging models have been performed, including those based on energy landscapes with hierarchically organized metastable states (see, e.g., Ref. \onlinecite{Sibani-hier}).  While it is unclear whether any of that work is relevant to our system, we note that resistance noise measurements on both spin glasses\cite{Weissman93} and the Coulomb glass in a 2DES in Si\cite{Jan2002} have provided support for the hierarchical picture of glassiness.\cite{Binder}  A number of other, remarkable similarities in the phenomenological behavior of these two types of glassy systems have been also reported.\cite{Jan2006,Jan2007-1,Jan2007-2}

\section{Dependence of conductance noise on applied electric field}
\label{exc}

The emergence of the non-Gaussian conductance noise for $n_s<n_g$, i.e. in precisely the same regime where relaxations exhibit the out-of-equilibrium behavior,\cite{Jan2006,Jan2007-1,Jan2007-2,Jan2009} as well as the detailed study and comparison of the noise obtained using different protocols (Secs.~\ref{noisepower} and \ref{pdf}), clearly demonstrate that it reflects the intrinsic glassiness of the 2DES.  For example, the qualitatively different behavior of the PDFs measured under the same conditions, but with different sample histories (Fig.~\ref{fig:pdfT}), rules out any spurious effects as possible sources of non-Gaussian noise from a $t_w$ protocol.  Likewise, the agreement of the results for the noise power spectra obtained using $\Delta T$ and $\Delta n_s$ protocols (Fig.~\ref{fig:Salpha}) implies that they have the same origin.  Experiments on a disordered indium oxide, an electron glass deep in the insulating regime, have shown that,\cite{Zvi-nonlinear} in that system, the degree of non-Gaussian behavior may be affected by the applied electric bias $V_{\mathrm{exc}}$.  In particular, the non-Gaussian nature of the noise was found to be non-monotonic with $V_{\mathrm{exc}}$, vanishing at both low and high bias, and peaking at an intermediate value of $V_{\mathrm{exc}}$.  This behavior was attributed to a nonlinear effect inherent to variable-range hopping transport.  While the indium oxide samples were indeed very deep in the hopping regime, this is in general not true of the experiments on a 2DES, where glassiness sets in already on the metallic side of the MIT ($n_g>n_c$).  Nevertheless, it is interesting to investigate whether similar effects might be observed also in the glassy regime of the 2DES, in particular in the range of parameters where non-Gaussian conductance noise has been studied (Secs.~\ref{noisepower}, \ref{pdf} and Refs.~\onlinecite{Snezana2002,Jan2002,Jan2004}).

In order to explore the effects of $V_{\mathrm{exc}}$, conductance noise $\Delta\sigma(t)/\langle\sigma\rangle$ was measured using a $\Delta T$ protocol.  After changing $n_s$ at a high $T$ (here at $\approx 20$~K), the 2DES was cooled to a desired $T$, and $\sigma (t)$ measured for constant values of $V_{\mathrm{exc}}$ ranging from 0.5~$\mu$V to 1~mV.  Typical data are presented in Fig.~\ref{fig:rawVexc}, which shows $\Delta\sigma/\langle\sigma\rangle$ vs
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig8.eps}
\caption {(Color online) Sample B1. Relative fluctuations of conductivity $(\sigma(t)-\langle\sigma\rangle)/\langle\sigma\rangle$ vs time for $n_s=3.24 \times 10^{11}$~cm$^{-2}<n_c$ and different $V_{\mathrm{exc}}$ at a) $T=0.100$~K and b) $T=0.300$~K.  Traces are shifted for clarity. The data were obtained after changing $n_s$ at $T \approx 20$~K, cooling down to the desired $T$, and then measuring $\sigma(t)$ for different values of $V_{\mathrm{exc}}$.}\label{fig:rawVexc}
\end{figure}
time for $n_s=3.24 \times 10^{11}$~cm$^{-2}$ in the insulating regime ($n_s<n_c$) and for several $V_{\mathrm{exc}}$ at two different $T$.  The noise is obviously non-Gaussian and large, of the order of 100$\%$ at low $T$.  The amplitude of the relative fluctuations decreases with increasing $T$, as established in earlier studies,\cite{Snezana2002,Jan2002} as well as for high enough excitations.

It is indeed obvious from the corresponding normalized power spectra $S_I(f)=C/f^{\alpha}$ (Fig.~\ref{fig:powerexc})
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig9.eps}
\caption {(Color online) Sample B1. The octave-averaged noise power spectra $S_{I}=C/f^{\alpha}$ vs $f$ for several $V_{\mathrm{exc}}$ at a) $T=0.100$~K and b) $T=0.300$~K.  Solid lines are linear least-squares fits with the slopes equal to $\alpha$.}\label{fig:powerexc}
\end{figure}
that the noise does not depend on the applied excitation at low enough $V_{\mathrm{exc}}$, but the power spectra $S_I$ become suppressed as $V_{\mathrm{exc}}$ is increased beyond a certain value.  A detailed comparison of the dependence of the exponent $\alpha$, the parameter $C$, and the time-averaged conductivity $\langle\sigma\rangle$ on $V_{\mathrm{exc}}$ is presented in Fig.~\ref{fig:sigmaexc}.  It is apparent that $\langle \sigma \rangle$ remains
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig10.eps}
\caption {(Color online) Sample B1.  (a) $\langle\sigma\rangle$, (b) exponent $\alpha$, and (c) parameter $C$ (in $S_I=C/f^{\alpha}$) vs $V_{\mathrm{exc}}$ at $T=0.100$~K (red squares) and $T=0.300$~K (blue dots) for $n_s=3.24 \times 10^{11}$~cm$^{-2}<n_c$.  In (a) and (c), the arrows show the onset of the non-Ohmic regime, and dotted lines guide the eye.  In (b), dotted lines represent the average values of $\alpha$ in the Ohmic regime: $\alpha=1.5\pm 0.2$ at $T=0.100$~K and $\alpha=1.7\pm 0.1$ at $T=0.300$~K.}\label{fig:sigmaexc}
\end{figure}
Ohmic over a wide range of $V_{\mathrm{exc}}$ [Fig.~\ref{fig:sigmaexc}(a)], up to an electric field $\sim 50$~V/m ($V_{\mathrm{exc}}\sim 100~\mu$V) at $T\sim 0.1$~K (where $eV_{\mathrm{exc}}\sim k_{B}T$).  In the same range of $V_{\mathrm{exc}}$, both the exponent $\alpha\approx 1.6$ [Fig.~\ref{fig:sigmaexc}(b)] and the noise magnitude, i.e. the parameter $C$ [Fig.~\ref{fig:sigmaexc}(c)], also do not depend on $V_{\mathrm{exc}}$ within experimental error.  The latter result is expected if the applied voltage is sufficiently small to be in the Ohmic regime and if the noise is in fact due to conductance (or resistance) fluctuations.  At higher $V_{\mathrm{exc}}$, where $\langle\sigma\rangle$ deviates from the Ohmic behavior, $\alpha$ does not seem to change with $V_{\mathrm{exc}}$ or $T$ within the scatter of data, but the parameter $C$ clearly decreases rapidly with increasing $V_{\mathrm{exc}}$, similar to the effect of raising $T$.  Therefore, the results strongly suggest that the suppression of the noise by high $V_{\mathrm{exc}}$ is most likely caused by simple heating.

Similar to previous studies of noise on the 2DES in Si\cite{Jan2002,Jan2004} and those in indium oxide\cite{Zvi-nonlinear}, the effect of $V_{\mathrm{exc}}$ on the non-Gaussian character of the noise has been investigated by analyzing the second spectrum $S_2(f_2,f)$, a fourth-order noise statistic.\cite{Weissman88,Weissman92,Weissman93}  In particular, $S_2(f_2,f)$ represents the power spectrum of the fluctuations of $S_{I}(f)$ with time.  For Gaussian noise, $S_2(f_2)$ is white (independent of $f_2$).  On the other hand, for non-Gaussian noise, $S_2 \propto 1/f_2^{1-\beta}$ with the exponent $(1-\beta)\neq 0$. Therefore, the value of the exponent $(1-\beta)$ is a convenient measure of the non-Gaussian character of the noise.  $S_2$ was analyzed by using digital filtering \cite{Seidler-Solin,Abkemeier} in a given frequency range $f=(f_L, f_H)$ (usually $f_H = 2f_L$).  Some examples of the normalized $S_2(f_2)$, with the Gaussian background subtracted, are presented in Fig.~\ref{fig:S2}(a) for the data in Fig.~\ref{fig:rawVexc}(a).  The second spectrum is clearly non-white, in
\begin{figure}
\includegraphics[width=8.1cm]{Lin-noise_fig11.eps}
\caption {(Color online) Sample B1, $T=0.100$~K, $n_s=3.24 \times 10^{11}$~cm$^{-2}$.  (a) Second spectral density $S_2 \propto 1/f_{2}^{1-\beta}$ vs $f_2$ for several values of $V_{\mathrm{exc}}$, as shown.  $S_2$ was measured in the frequency octave $f=(f_L,2f_L)$, where $f_L = 16$~mHz. The solid lines are linear least-squares fits with the slopes equal to $1-\beta$.  (b) Exponent ($1-\beta$) vs $V_{\mathrm{exc}}$ for different values of $f_L$, as shown. Dotted line guides the eye.}\label{fig:S2}
\end{figure}
agreement with earlier work.\cite{Jan2002}  Here we find that the exponent $(1-\beta)\approx 0.8$ does not depend on $V_{\mathrm{exc}}$ [Fig.~\ref{fig:S2}(b)] in the regime studied.  Therefore, non-Gaussian noise is observed over a wide range of electric fields in the regime of linear or Ohmic response.

Detailed studies of non-Gaussian noise in the 2DES in Si (Secs.~~\ref{noisepower}, \ref{pdf} and Refs.~\onlinecite{Snezana2002,Jan2002,Jan2004}) have been performed with the applied electric fields $\sim 5$~V/m, i.e. an order of magnitude lower than the field where the non-Ohmic behavior sets in (see Fig.~\ref{fig:sigmaexc}).  This is in contrast to the experiment on indium oxide,\cite{Zvi-nonlinear}, where only the non-Ohmic regime was explored, and where the non-monotonic dependence of the non-Gaussian behavior on $V_{\mathrm{exc}}$ was reported.  A study of such a deep non-linear regime, however, is beyond the scope of this work, so the existence of a similar effect in our system for much higher $V_{\mathrm{exc}}$ cannot be ruled out.  Our study does demonstrate, on the other hand, that the non-Gaussian noise in a 2DES is not due to the same non-linear mechanism caused by the high bias.  Instead, it provides further evidence that the correlated noise reflects the intrinsic out-of-equilibrium dynamics of the 2DES in Si.

\section{CONCLUSIONS}
\label{concl}

We have presented a study of the conductance noise in a 2DES in Si at low $T$, measured by using three different experimental protocols.  In the so-called $\Delta T$ protocol, the carrier density $n_s$ is fixed at a high $T$ where the system is in an equilibrium state, and the noise is measured after cooling with a fixed $n_s$.  The results provide evidence that, similar to other types of glasses, the 2DES falls out of equilibrium as $T$ is reduced.  This glassy freezing is observed for all densities $n_s<n_g$, as $T\rightarrow 0$.

In a $\Delta n_s$ protocol, $n_s$ is changed by a small amount at such low $T$ that the system is unable to equilibrate on experimental time scales.  While the perturbations introduced by both $\Delta T$ and $\Delta n_s$ protocols are too small to lead to observable relaxations of $\sigma$ (i.e. $k_{B}\Delta T, \Delta E_F\ll E_F$), they produce the same results for the power spectra of the non-Gaussian noise observed for $n_s<n_g$, where the glass transition density $n_g$ is found to be the same in both cases.  The analysis of the probability density functions (PDFs) of the fluctuations, however, reveals that $\Delta n_s$ has a strikingly different effect on the system than $\Delta T$ (both $\Delta E_F, k_{B}\Delta T>k_{B}T$).  In particular, the results strongly suggest that the density change reshuffles all energies, because of the Coulomb interactions, thus modifying the free energy landscape of the 2DES.  For this reason, theoretical modeling of the glassy dynamics in this system might be considerably more difficult than in some other types of glassy materials.

The effect of the density change was studied further by using a waiting-time ($t_w$) protocol, in which $n_s$ is changed only temporarily during $t_w$, but $\Delta E_F>E_F$.  As a result of such a large perturbation, the system exhibits visible relaxations of $\sigma$, i.e. aging.\cite{Jan2006,Jan2007-1,Jan2007-2}  Interestingly, the non-Gaussian PDFs exhibit both history dependence and an evolution towards a Gaussian shape as the system ages and slowly approaches equilibrium, similar to the behavior of a great variety of out-of-equilibrium systems.

The power spectra and the higher order noise statistics of the non-Gaussian noise have been also investigated over a wide range of the applied bias voltage.  The data demonstrate that the non-Gaussian noise is observed in the regime of linear response, i.e. that it is not caused by the application of a high bias, but rather that it reflects the intrinsic out-of-equilibrium behavior of the 2DES.

In summary, we have established several new characteristics that the 2DES in Si has in common with a large class of both 2D and 3D out-of-equilibrium systems.  Thus this work provides additional strong evidence that many such universal features are robust manifestations of glassiness, regardless of the dimensionality of the system.  In addition, we have revealed some effects that are unique to Coulomb glasses.  Therefore, our findings should be helpful in the understanding of the complex behavior near the MIT in a variety of strongly correlated materials.

\acknowledgments
This work was supported by NSF DMR-0905843, the National High Magnetic Field Laboratory through NSF Cooperative Agreement DMR-0654118, and the State of Florida.


\section{Introduction}
The PALFA survey is an on-going, large-scale, survey for radio pulsars tuned to
find millisecond pulsars (MSPs). PALFA is similar to, but more sensitive than
other on-going projects (e.g. HTRU-N/S, GBNCC, SPAN512, described by Ng et al.,
Keith et al., Lynch et al., and Desvignes et al., respectively, in these
proceedings). MSPs are particularly useful for searches for the direct
detection of nano-Hz gravitational waves using Pulsar Timing Arrays
\citep[e.g.][]{vlj+11}, strong-field tests of relativistic theories of gravity
\citep[e.g.][]{ksm+06}, studying the equation of state of ultra-dense matter
\citep[e.g.][]{dpr+10}.  Another main goal of the survey is to expand the known
population of radio pulsars, to gain a better understanding of their overall
numbers, the distribution of their properties, determine birth-rates, and
estimate birth properties. The survey is sensitive to transient pulsars, such
as nullers, rotating radio transients (RRATs), intermittent pulsars, and radio
magnetars.  Discoveries may also be young pulsars associated with high-energy
emission, supernova remnants, or pulsar wind nebulae. PALFA has the potential
for discovering even more exotic systems such as pulsar-black hole binaries,
which could have profound ramifications for the understanding of gravitation
\citep[e.g.][]{wk99}.

The PALFA Survey makes use of the 7-beam L-band (1.4 GHz /
21-cm) receiver installed on the Arecibo Observatory's William E. Gordon
Telescope. Observing has focused on the plane of the Galaxy ($|b| <
5^{\circ}$), in the two regions visible from Arecibo, due to its transit
design, the ``inner Galaxy'' ($34^{\circ} < l < 77^{\circ})$ and the
``anti-centre'' ($168^{\circ} < l < 214^{\circ}$) regions.

The Wideband Arecibo Pulsar Processor (WAPP) spectrometers used by the PALFA
survey since its beginning in 2004 were replaced with the more sensitive Mock
spectrometer system in 2009. The Mock set-up allows for ALFA's full 320 MHz
bandwidth to be recorded, more than tripling the bandwidth available compared
to the older WAPP set-up. In 2009, before completely phasing out the WAPP
spectrometers for use in the survey, both backends were run in parallel to
ensure the quality of the data from the new system. The ALFA receiver's system
temperature is $\sim$24 K, and the gain is $\sim$10.4 K/Jy for the central beam
and $\sim$ 8.2 K/Jy for the outer beams \citep{cfl+06}. The survey's observing
parameters, summarized in Table \ref{tab:params}, are characterized by high
time and frequency resolution. Integration times are 4.5 min for
the inner Galactic region, and 2.25 min for the anti-centre region.

\begin{table}[t]
\footnotesize
\centering
\caption{PALFA Survey observing system parameters.}
\label{tab:params}
\begin{tabular}{rcc}
    \hline
    & WAPP Spectrometers & Mock Spectrometers \\
    \hline
    Centre Frequency & 1420/1440 MHz$^{a}$ & 1375.432 MHz \\
    Sample time & 64 $\mu$s & 65.476 $\mu$s \\
    Bandwidth & 100 MHz & 322.617 MHz \\
    Number of Channels & 256 & 960 \\
    \hline
\end{tabular}
\footnotesize
\\$^{a}$ The centre frequency of the WAPPs was increased in 2005 Nov. to reduce interference.
\end{table}

The high resolution employed by the PALFA survey allows for the dispersive
smearing caused by free electrons in the interstellar medium along the
line-of-sight to be almost completely removed for pulsars at any depth into the
Galactic plane using incoherent dedispersion techniques. Scattering due to
multi-path propagation still limits the depth to which the fastest spinning
pulsars can be detected within the Arecibo-visible sky (see Figure
\ref{fig:broad+timeline}, left).

The large collecting area and hence high gain of the Arecibo Observatory make
the PALFA survey the most sensitive large-scale survey of the Galactic plane to
date. The high sensitivity of the survey enables integration times far shorter
than would be possible when using smaller telescopes, and reduces the total
data volume of the survey. It is important to note that shorter observations
are more computationally efficient to search, especially for pulsars in binary
systems.

\section{Processing \label{sec:processing}}
Processing of PALFA survey data is currently done with two complementary
pipelines. One is based on the PRESTO suite of pulsar search code
\footnote{http://www.cv.nrao.edu/$\sim$sransom/presto/}. The other pipeline
makes use of the BOINC-based\footnote{http://boinc.berkeley.edu/index.php}
Einstein@Home distributed global volunteer computing
platform\footnote{http://einstein.phys.uwm.edu/}.

Before data are distributed to processing sites the data are converted to
4-bits-per-sample PSRFITS format, and archived at Cornell University's Center
for Advanced Computing. All data are tracked by a database at Cornell. Data are
automatically downloaded from the archive by processing sites.

\subsection{PRESTO}
The PRESTO pipeline first searches for, and removes, bright un-dispersed
narrow-band impulsive, and periodic, radio frequency interference (RFI) signals
in the data.

The pipeline then dedisperses the data with an assumed dispersion measure (DM)
value, removing the corresponding frequency-dependent delay. Since the DM of
unknown pulsars are not known \textit{a priori}, 4188 (1140) dedispersed time series,
each with a different trial DM value ranging from 0 to $\sim$1000 pc~cm$^{-3}$,
are generated when searching Mock (WAPP) observations. Each resulting time
series is then searched in the time domain for bright individual pulses, and
in the Fourier domain for periodic signals.

Each dedispersed time series is searched for single pulses. To increase
sensitivity to pulses of width larger than the observation's sampling time a
series of boxcar matched filters of different widths are convolved with the
data. A summary plot of all significant candidate pulses identified at all
trial DMs is produced, and inspected.

In order to increase sensitivity to pulsars in binary systems a series of trial
Fourier-domain templates, each corresponding to a constant acceleration, are
used when searching the Fourier Transform of the dedispersed time series
\citep[see][for more details]{rem02}. Given the length of PALFA observations,
and the number and range of trial acceleration values, the PRESTO pipeline
maintains sensitivity to pulsars in binary orbits with periods as small as
$\sim$1~hour.

The PRESTO pipeline then identifies the most likely pulsar candidates from the
list of periodicity signals identified. For each periodicity candidates the raw
data are folded using its period and DM, and a diagnostic plot is produced. Also,
a series of scoring heuristics are computed for each folded candidates. These
are later used as input to Artificial Intelligence (AI) systems trained to identify
pulsar-like candidates, and ignore RFI-like and noise-like signals that were
folded by the pipeline. All periodic signals identified are summarized in
overview diagnostic plots.

The PRESTO pipeline is run on computing resources at the University of British
Columbia, and on the Guillimin
supercomputer\footnote{http://www.clumeq.ca/index.php/en/about/computers/guillimin}
managed by McGill University. Combined, 3000-6000 observations can be reduced
per month, depending on contention for resources on Guillimin.

\subsection{Einstein@Home}
Prior to distributing data to volunteers' computers, the Einstein@Home-based
pipeline identifies, and masks narrow-band impulsive RFI in the time-domain.
The pipeline also identifies periodic RFI, which is replaced by random noise.
The pipeline then generates 3808 (628) dedispersed time series for Mock (WAPP)
data, each with a different trial DM, up to $\sim$1000 pc~cm$^{-3}$. These
dedispersed time series that are transferred to volunteers' computers to be
analyzed.

Sensitivity to binary pulsars is improved by using 6662 circular orbital
templates to demodulate the time series in the time domain. This technique is
capable of detecting isolated pulsars \citep[e.g.][]{kac+10}, and still
maintains sensitivity to pulsars in binaries with orbital periods as short as
$\sim$11 min. The Fourier transform of each demodulated time series is searched
for periodic signals. The top 100 signals are returned to the Einstein@Home
servers where they are collated, diagnostics plots are created, and heuristic
pulsar-identification algorithms are run.

No single pulse search is performed in the Einstein@Home analysis.

The Einstein@Home pipeline can make use of volunteers' GPUs, if permitted.
Processing on GPUs accounts for $\sim$75\% of the pipeline's processing power.
The performance of the pipeline, and more importantly the GPU code, is
constantly being improved, reducing the time required for analysis with each
iteration. Currently, the Einstein@Home pipeline searches $\sim$3000
observations per month.

\subsection{Candidates}
The PALFA collaboration uses a centralized database to store metadata and
diagnostics for each observation, and details for candidates identified by
processing. The database holds approximately six million candidates.

Candidates are primarily inspected using an online viewer available on the
PALFA Survey's portal on the CyberSKA website.\footnote{www.cyberska.org/palfa.
Access is restricted to members of the PALFA Consortium.} The viewer allows
users to select candidates based on information in the database, including
scoring heuristics, AI recommendations, and the usual observation/candidate
information. The user interface provides one-click access diagnostic
information and plots about the candidate being viewed, and the observation in
which it was found. The candidate viewer also provides direct access to another
web-application designed to track the most promising candidates, and their
follow-up. Both applications, and others not mentioned here, are hosted at
McGill University.

\section{Results}
By the end of August 2012, the PALFA survey has discovered 100 pulsars. The
discovery rate increase significantly when the analysis of the Mock data began
(see Figure \ref{fig:broad+timeline}, right). This is partly due to the
increase of computing resources available. However, the number of discoveries
per sq. deg. is higher for Mock data, owing to their higher sensitivity. The WAPP
spectrometers have been used to discover 56 pulsars in 164 sq. deg. of PALFA
data (1 pulsar per 2.9 sq. deg.) whereas 44 pulsars have been discovered in
only 70 sq. deg. of Mock data (1 pulsar per 1.6 sq. deg.).

\begin{figure}[t]
    \centerline{
    \includegraphics[width=2.5in]{trimmed_palfa_depths.eps}
    ~
    \includegraphics[width=2.75in]{trimmed_palfa_discoveries_timeline.eps}
}
\caption{{\it Left}.  A top-down view of the Galactic plane showing Galactic
MSP overlaid with the depths corresponding to various pulse broadening times
as predicted by the NE2001, for the PALFA survey's observing parameters. \quad
{\it Right}. Number of PALFA discoveries per year since the start of the
survey. The projected number of discoveries for 2012 is shown with a dashed
extension. \label{fig:broad+timeline}}
\end{figure}

\subsection{Highlights}
The PALFA survey has made several exciting discoveries since its start,
including:

\begin{description}
\item[PSR~J1906+0746] \quad (\textit{P=144~ms; DM=212~pc~cm$^{-3}$})
This binary pulsar has a characteristic age of
$\tau_{c}=\frac{P}{2\dot{P}}$=112~kyr, making it the youngest known binary
pulsar \citep{lsf+06}. Also, it is the second most relativistic binary pulsar
system known. \citep[See also,][]{kas12}.

\item[PSR~J1903+0327] \quad (\textit{P=2.15~ms; DM=298~pc~cm$^{-3}$})
This Galactic MSP is in a binary orbit whose eccentricity of $e=0.44$
challenges the standard MSP formation picture \citep{crl+08}.
Pulsar timing has measured the Shapiro delay, and the relativistic advance of
periastron, allowing the determination of the pulsar's mass
($M_p=1.67(2)M_\odot$) and the companion's mass ($M_c=1.029(8)M_\odot$). Near-IR
spectroscopic observations have identified the companion as a main-sequence
star. Simulations suggest the current binary was formed from a triple-star
system \citep{fbw+11,pvvn11}.

\item[PSR~J1856+0245] \quad (\textit{P=80.9~ms; DM=650~pc~cm$^{-3}$})
This young ($\tau_{c}=$ 27~kyr), highly energetic ($\dot{E}=4\pi^{2}I\dot{P}P^{-3}=4.6 \times
10^{36}$~erg~s$^{-1}$) pulsar is coincident with the TeV $\gamma$-ray source
HESS~1857+026, as well as a faint X-ray source AX~J185651+0245 \citep{hng+08,rgv+12}.

\item[PSR~J1949+3106] \quad (\textit{P=13.1~ms; DM=164~pc~cm$^{-3}$})
This binary MSP is in a circular ($e=4 \times 10^{-5}$), whose inclination
($i=80^{\circ}$) makes the detection of Shapiro delay possible in the timing
data \citep{dfc+12}. The measurement of Shapiro delay provides mass
measurements of the pulsar ($M_p=1.47^{+0.43}_{-0.31}~M_{\odot}$) and of the
companion ($M_c=0.85^{+0.14}_{-0.11}~M_{\odot}$). Detection of relativistic
advance of periastron, which could be possible within a few years, will greatly
reduce the uncertainties on the masses. Also, the absence of timing noise, the
decent timing precision ($RMS_{TOA} = 3.96\mu$s), and the as-of-yet unresolved
pulse profile, coupled with Arecibo's new PUPPI backend could make this pulsar
useful for Pulsar Timing Array experiments.
\end{description}

\subsubsection{Global distributed volunteer computing discoveries}
The Einstein@Home platform's first discovery was found in 2010 in a WAPP PALFA
data \citep[J2007+2722,][]{kac+10}. Since then, the Einstein@Home
pipeline has discovered 22 additional pulsars, including four MSPs, and three
binaries.

\subsubsection{Distant MSPs}
The PALFA survey has discovered 17 MSPs (defined here as P $<$ 30~ms). These
pulsars are amongst the highest DM MSPs in the Galactic plane. The ratio DM/P
can be used as a rough measure of the difficulty of finding a pulsar, since the
smearing in time caused by DM is more detrimental to shorter period pulsars.
The MSPs found in the PALFA survey make up a large proportion of the highest
DM/P pulsars known (see Figure \ref{fig:dmp+topdown}, left).

Also, by assuming the NE2001 model for the Galactic distribution of
free-electron density \citep{cl02}, it is possible to infer distances given the
sky-position and DM of a pulsar. Figure \ref{fig:dmp+topdown} (right) shows the
inferred positions of PALFA-discovered MSPs projected onto the plane of the
Galaxy, compared to all other Galactic plane MSPs. The median distance to PALFA
MSPs is 5.8~kpc, double that of other known Galactic MSPs ($D_{median} =
2.9$~kpc), highlighting PALFA's ability to find distant MSPs.

\begin{figure}[t]
    \centerline{
    \includegraphics[width=2.75in]{trimmed_dm_vs_p.eps}
    ~
    \includegraphics[width=2.45in]{trimmed_palfa_discoveries_topdown_msp.eps}
}
\caption{{\it Left}. DM vs. Period for Galactic MSPs. Notice that red Xs,
PALFA-discovered MSPs tend towards the top of the plot. These pulsars are
generally thought to be harder to find. \quad {\it Right}. A top-down view of
the Galactic plane showing Galactic MSP. The PALFA-discovered MSPs (red Xs) are
generally located deeper into the plane of the Galaxy. Distances are inferred
from the DM using the NE2001 model. \label{fig:dmp+topdown}}
\end{figure}

\section{Future Outlook}
The PALFA survey has observed only 10\% of its nominal survey region with the
Mock spectrometers, which are proving to be superior to the WAPPs due to their
tripled bandwidth. The Mock spectrometers will be used as the survey continues,
and will eventually be used to re-observe regions covered with the WAPPs.

In addition, an updated version of the PRESTO pipeline, including enhanced RFI
mitigation, improved candidate optimization, additional diagnostic information,
and a few bug fixes, is being beta-tested on the Guillimin supercomputer. The
updated pipeline is already finding pulsars missed by the first analysis.
Improved single pulse analyses based on the work of \cite{sccs12} are also
being applied to the PALFA data taken thus far, and will be integrated into the
survey's standard data reduction, as will a modified version of the algorithm
developed by Karako-Argaman et al. (in these proceedings).

With the 100 pulsars discovered to date, the improved data quality, additional
computing resources available, and enhancement to the data reduction, the
PALFA survey is well on its way to discovering hundreds of pulsars, including
expanding the population of faint, distant MSPs it has already been successful
at finding.


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

Observations in the GeV range of gamma-ray bursts (GRBs) by the Fermi LAT
\citep{Atwood+09LAT} instrument have uncovered new peculiarities of these
objects, which were unknown previously. Among the  minority of bursts which
show significant emission above 100 MeV \citep{Kocevski+12LATpaucity}, we
concentrate here on a handful of bursts with substantial detected GeV emission
to construct models aimed at interpreting the radiation and the basic
parameters of the bursts.

The bursts with observed GeV emission have rather diverse spectra. So far,
because of relatively low number of such GRBs, we restrict ourselves only to
bursts which have a large enough data coverage to make detailed model fits
meaningful.  In general terms, the GeV emission in some bursts can be described
through an extrapolation of the lower energy spectrum, without strong evidence
for departures from this, while in other bursts an extra power law component
emerges at GeV energies, in addition to the spectral component describing the
MeV radiation.  An important new feature present in both short and long bursts
is the presence of a delay between the MeV and the GeV radiation. Furthermore,
some bursts also exhibit strong spectral evolution.

The origin of the GeV component has been attributed to leptonic mechanisms in
the external forward shock, e.g. in the radiative regime in a pair enriched
interstellar medium \citep{Ghisellini+09-gevrad} or in the adiabatic regime
\citep{Kumar+09gevforsho}, or alternatively to hadronic cascades
\citep[e.g.][]{Asano+09-090510}, among others.  In the external forward shock,
inverse Compton scattering and Klein-Nishina effects were considered by
\citet{He+12grbnu}, and it has been argued that at least the first few seconds
of the early high energy emission cannot be attributed to the forward shock
{\citep{He+11-090510,Gao+09HE}}.  Thus, it is important to consider both the prompt
(largely MeV) emission and the GeV emission jointly, as they interact with each
other. Two recent developments which are relevant for the prompt emission are
the increasing consideration of magnetically dominated models for the GRB
outflow, and of photospheric models for the formation of the prompt radiation
spectrum \citep[e.g.][]{Meszaros+12raa}.

Previously, we have shown  that some of the  brightest {\it LAT } detected bursts can
be modelled by a magnetic two component model \citep{Veres+12magnetic}
\citepalias[][hereafter]{Veres+12magnetic}.  We showed that within this model,
without going into detailed model fits, it is possible to incorporate the
seemingly different behaviors of the high-energy spectrum (which is either a
single Band function or a Band function with an extra high energy component)
into a single framework.  We also addressed the non-detection problem of bursts
by LAT, namely some of the GRBs should have been detected by LAT if the
spectrum continues smoothly from the MeV regime, or generally there is lack of
LAT detected bursts compared to initial estimates, based on the simple
extrapolation of the MeV emission. We have also presented a generalized
formalism for the dynamics of GRBs valid for a generic acceleration law
$\Gamma\propto r^\mu$ with $1/3\le\mu\le1$ \citep{Veres+12-110721a}, where
$\mu$ can be understood as a parameter accounting for the bulk Lorentz factor
$\Gamma(r)$ averaged across the effective cross section of the jet at that
radius.

Here we present detailed numerical fits of the MeV and GeV data to specific
physical models of GRBs based a photospheric origin of the prompt emission in
two dynamic regimes, i.e. the extreme magnetically dominated regime and the
usual baryonically dominated regime, with inclusion of the inverse scattering
effects (both self and external) incurred in the external forward and reverse
shocks.  This is one of the few instances of fitting GRB spectra with
physically detailed  models (for other examples see \citet{Ryde+10phot090902,
Burgess+11synch}). We construct the models based on the generalized dynamics
described above, and convolve them with the response function of the {\it Fermi }
detectors, showing that they give a good statistical description of the data.

In this paper, we analyze four representative bursts, which capture the
diversity observed in LAT: GRBs 080916C, 090510, 090902B and 090926A. We model
their spectrum with the two component magnetic ($\mu=1/3$) model and
demonstrate that it is a viable model to interpret the spectra from
approximately 10 keV to 10 GeV. We also show that the baryonic model ($\mu=1$)
gives similarly good fits to the data, as quantified by specific goodness of
fit criteria.

We present our basic  model in \S 2. The Fermi observational data for the four
GRBs and the fitting method are described in \S 3. The fitting results are
presented in \S 4. We summarize the results and compare the baryonic with the
magnetic model, with a discussion on the physical implications, in \S 5.

\section{Photospheric model with arbitrary acceleration parameter}
\label{sec:model}

Our basic  model consists of a dissipative photosphere, whose location depends
on whether the dynamics there is kinetic energy dominated or magnetically
dominated, and an external shock region, where reverse shocks may be present as
well as a forward blast wave. The bulk of the prompt MeV radiation arises from
the neighborhood of the photosphere in the form of a nonthermal
component\footnote{Depending on model details, the spectrum may start forming
below the photosphere, and further contributions can arise from outside of it
\citep{Peer+06phot,Lazzati+10phot,Beloborodov10phot, Giannios11peak}; for our purposes here, we can assume
that the bulk of the radiation arises at the photosphere.}, and depending on
the model, a separate, subdominant thermal component can also be evident.  At
the deceleration radius, shock accelerated electrons upscatter the prompt
photons to produce a component peaking in the GeV range.

In the present paper we generalize the model presented in \citetalias{Veres+12magnetic}
to cover a wider range of dynamic behavior stretching from the extreme kinetic energy
dominated to the extreme magnetic dominated cases.  The acceleration of the ejecta is
described by the expression (\ref{eq:accel}) for the bulk Lorentz factor
\begin{equation}
{\Gamma(r)}\propto \left\{
\begin{array}{lll}
 r^{\mu}		&	{\rm if}	& r<r_{\rm sat}\\
 {\rm const.}	&	{\rm if}	& r_{\rm sat}<r<r_{\rm dec}.
\end{array}
\right.
\label{eq:accel}
\end{equation}
Here the exponent $\mu$ is taken as $\mu=1/3$ when the dynamics is
governed by dissipation of the energy stored in magnetic fields in the presence of
a sub-dominant inertia of the baryonic component \citep[e.g.][]{Drenkhahn02,
Meszaros+11gevmag}, while if the baryonic component dominates, we have	$\mu=1$.
Intermediate values of	$\mu$ can arise in magnetic outflows with certain symmetries{, and the exact value of $\mu$ also depends on the magnetic field configuration}
\citep[e.g.][]{Tchekhovskoy+10grb,McKinney+11switch}.   $r_{\rm sat}$
and $r_{\rm dec}$ are the saturation (coasting) and the deceleration radii, and
with this parametrization the radius where the ejecta reaches the coasting Lorentz
factor is $r_{\rm sat}= r_0 \eta^{\frac{1}{\mu}}$. For simplicity we consider the
transition from the acceleration phase to coasting as a broken power law function,
although this transition occurs more gradually \citep[e.g.][]{kobayashi99}.

In what follows is it is useful to give the dependence of the physical parameters
in terms of a general $\mu$. The main difference between the models is whether
the photosphere occurs below or above the saturation radius. This is governed by
the relation of $\eta$ and a quantity $\eta_T$ given by eq.(\ref{eq:etat}).
Here $\eta_T$ is a limiting value for the Lorentz factor, which separates the two
cases of $\eta>\eta_T$ (photosphere below saturation) and $\eta<\eta_T$ (photosphere
above saturation, where $\eta_T$ is a function of the dynamics, depending on $\mu$.
This limiting Lorentz factor is obtained from setting $r_{\rm sat}=r_{\rm ph}$ for
a general $\mu$, generalizing the related value introduced in \citet{Meszaros+93gasdyn},
\begin{equation}
\eta_{T}=\left(\frac{ L\sigma_T}{8\pi m_p c^3 r_0}\right)^{{\mu}/{(1+3\mu)}}.
\label{eq:etat}
\end{equation}
The photosphere occurs in the acceleration phase if $\eta>\eta_T$, typical for
a magnetically dominated ($\mu=1/3$) case, where $\eta_T\simeq 150~L_{53}^{1/6}
r_{0,7}^{-1/6}$.  The photosphere is in the coasting phase for $\eta<\eta_T$ which it
is typical for baryonic cases ($\mu=1$) for which $\eta_T\simeq 1900~L_{53}^{1/4}
r_{0,7}^{-1/4}$.  The two radii are equal for $\mu\approx0.6$ if other nominal
parameters are left unchanged.	The photosphere will occur at
\begin{eqnarray}
\frac{r_{\rm ph}}{r_0}&=&\left(\frac{ L\sigma_T}{8\pi m_p
c^3r_0}\right)\frac{1}{\eta \Gamma_{\rm ph}^2} =\\
&=&\left\{
\begin{array}{ll}
 \eta_T^{1/\mu}(\eta_T/\eta)^{1/(1+2\mu)}	&	{\rm if~ }	 \eta>\eta_T \\
 \eta_T^{1/\mu}(\eta_T/\eta)^3			&	{\rm if~ }	 \eta<\eta_T
\end{array}
\right.
\end{eqnarray}
noting that $\Gamma_{\rm ph}$ carries an $r_{\rm ph}$ dependence for
$\eta>\eta_T$. The Lorentz factor at the photosphere is $\Gamma_{\rm
ph}=(r_{\rm ph}/r_0)^{\mu}\propto L^{\mu/(2\mu+1)} \eta^{-\mu/(2\mu+1)}
r_0^{-\mu/(2\mu+1)} $ for $\eta>\eta_T$ and $\Gamma_{\rm ph}\approx\eta$ in
other cases.
The shell has initially a constant width, $\Delta r = c t_w\approx 3 \times
10^{11}  t_{1} \hbox{~cm}$, where $t_w$ is the time during which the central engine is
active. At a later stage the shell starts to spread radially when the
when the relative width $\delta r/r\sim \Delta v/c\sim 1/2\Gamma^2$ becomes of the
same order as the initial shell width $\Delta r /r$.  Depending on the regime, this
can happen at $r_{\rm spread}=
 \left(2 c t_w/r_0^{2\mu}\right)^{1/1-2\mu}	{\rm if } ~r<r_{\rm sat} $ and $
 2 \eta^2 c t_w 	 {\rm if } ~r>r_{\rm sat}$.
For most of the parameters used here the spreading occurs well beyond the
saturation radius, $r_{\rm spread}=6\times 10^{17} \eta_{600}^2 t_1\hbox{~cm} $, and
it can be neglected close to the photosphere.
The scaling of the baryon number density in the comoving frame is
\begin{equation}
n'_b(r)=L/4\pi r^2 m_p c^3 \Gamma(r) \eta.
\label{eq:nb}
\end{equation}
This results in a dependence:
$n'_b(r)\propto V'^{-1}(r)\propto r^{-2-\mu} $ if $ ~r<r_{\rm sat}$ or	$
r^{-2}$ if $ r>r_{\rm sat}$.

{Prompt spectrum:} Observationally the prompt emission spectrum is described by
a Band function \citep{Band+93}, and a spectrum with broadly similar attributes is
shown to arise in photospheres \citep{thompson94,Beloborodov10phot,Vurm+11phot,Giannios11peak,
Peer+12thermal,Beloborodov12peak}.  The location of the peak of the spectrum may be
due to effects related to electron heating by nuclear collisions, reconnection, or
semirelativistic shocks, all of which give the peak energy at the right order of
magnitude. A detailed radiative transfer would be model-dependent and too
cumbersome for performing numerical fits to the data. For the purposes of this study
we choose a simple prescription where the peak energy is determined by the synchrotron
peak from mildly relativistic shocks around the photosphere\footnote{Shocks are
expected as result of magnetic reconnection, as well as possibly from outflow
Lorentz factor fluctuations, the particle acceleration time being of the same order
for reconnection and Fermi acceleration \cite{Meszaros+11gevmag}, and the spectral peaks
are in both cases in the MeV range, as for baryonic thermal photosphere models.}.

It should be noted however, that the above is just a simple heuristic
prescription for the peak, the shape being taken to be the Band shape.  The
Band spectra and peaks deduced in the previously mentioned numerical
simulations arise from distortions of the original black-body component formed
around the launching radius through Compton scattering and synchrotron
radiation, which escapes at the photosphere. In these cases no additional
black-body component is expected below the peak. In our model considered here
the peak is ascribed to synchrotron radiation from electrons accelerated near
the photosphere, which introduces a dependence on the magnetic field; and the
advected blackbody component is also released, but at energies typically below
the peak, for magnetically dominated dynamics. There are possible shortfalls
for a pure synchrotron interpretation, such as the low energy spectral slope
observations which show inconsistencies with the synchrotron model when electrons are
in the fast cooling regime \citep[e.g.][]{Preece+98syn}. While recent analyses
\citep{Sakamoto+11batcat} show that the fraction of bursts showing problems is smaller
then previously thought, there is still a non-negligible fraction of bursts
where other mechanisms may need to be invoked. This could be the Compton
upscattering mentioned above, or modifications of the simple synchrotron
mechanism, e.g. such as discussed by \citet{Medvedev06jitter,Daigne+11syn}. A first
principles discussion of such cases for our fitting exercise undertaken here on
individual bursts would be too extensive, and for this reason we have adopted
the simple prescription mentioned above.

The peak of the prompt spectrum is then given by $\varepsilon_b=({3q_e B'_{\rm ph}}/
{4 \pi m_e c}) \gamma_{\rm e,ph}^2 [{\Gamma_{\rm ph}}/({1+z})]$, where $B'\propto
{n'}_b^{1/2} \Gamma_r$ {(taking $n'_b$ from Eq. \ref{eq:nb})} and $\gamma_{e,ph}\propto \Gamma_r$
\begin{eqnarray}
\varepsilon_b \propto \left\{
\begin{array}{ll}
  L^{\frac{3\mu-1}{4\mu+2}} \eta^{-\frac{3\mu-1}{4\mu+2}}
  r_0^{\frac{-5\mu}{4\mu+2}} \Gamma_r^3/(1+z)		&	{\rm if~ }  \eta>\eta_T\\
  L^{-1/2} \eta^{3}
    \Gamma_r^3/(1+z)		&	{\rm if~ }  \eta<\eta_T.
\end{array}
\right.
\label{eq:ebreak}
\end{eqnarray}
{Here we have assumed an efficiency of 0.5 \citep{Giannios11peak}.}
The two cases imply widely different values for the same set of parameters.
This is mainly caused by the strong dependence on $\eta$  when the
photosphere arises in the coasting phase ($\eta<\eta_T$).  Note furthermore,
that the photosphere loses any memory of the acceleration in the same case and
no $\mu$ dependence is present. Also in this case, the peak energy is
independent of the launching radius.  We take the low and high energy photon
number indices from observations ($\alpha\approx 1$, $\beta\approx2.2$), and in
the fitting sections we directly fit for these.

{Thermal component: } Even though the main emission component is nonthermal, one
expects a thermal component to arise from the photosphere. There are hints of
this separate thermal component emerging below the Band peak
\citep{Page+11Xthermal, Guiriec+10-100724, mcglynn+12phot}. The comoving
temperature goes as ${T'}\propto r^{-\left(\frac{\mu+2}{3}\right)}	$ if $
r<r_{\rm sat} $ and $ r^{-2/3}	$ if $	 r_{\rm sat}<r<r_{\rm spread}.$
This results in a dependence of the observed temperature as:
\begin{eqnarray}
T_{\rm obs}(r_{\rm ph}) \propto \left\{
\begin{array}{ll}
  L^{\frac{14\mu-5}{12(2\mu+1)}} \eta^{\frac{2-2\mu}{6\mu+3}}
  r_0^{-\frac{10\mu-1}{6(2\mu+1)}} /(1+z)	& {\rm if~ } \eta>\eta_T \\
  L^{-5/12} \eta^{8/3} r_0^{1/6} /(1+z) 	&	{\rm if~ } \eta<\eta_T .
\end{array}
\right.
\label{eq:temp}
\end{eqnarray}

{External shock-related components:}
The ejecta starts to decelerate when it has plowed up approximately $1/\eta$
times its own mass from the interstellar matter. This occurs at the
deceleration radius $ r_{\rm dec}= 4.8 \times 10^{16}	L_{t,53}^{1/3}
(1-\zeta_r)^{1/3} t_{1.3}^{1/3} n_0^{-1/3} \eta_{600}^{-2/3} \hbox{~cm}$. Here
$\zeta_r$ is the fraction of the luminosity in radiation, which for simplicity
is fixed to $0.5$ throughout the rest of this article.	At the	deceleration
radius a forward shock develops (and possibly a reverse shock as well). The
prompt emission is up-scattered on the shock-accelerated electrons, giving
rise, to a $\gtrsim 1 \hbox{~GeV}$ peaked radiation. This component dominates at
sufficiently high energies, but radiation from synchrotron self Compton (SSC)
from the forward shock (FS) may also be relevant here.	The deceleration radius
is larger than the photospheric and the saturation radius for every realistic
set of parameters.

\section{Data and fitting method}
\label{sec:datafit}
\subsection{Data}
\label{sec:data}

Currently there are at least 40  LAT-detected GRBs at energies $\gtrsim 100 $
MeV \citep{Omodei+11fermigrb,Zheng+12extraLAT}.  Out of these, there are only a
handful of bursts with	enough data at $\hbox{~GeV}$ energies to derive meaningful
constraints on the model parameters \citep{Zhang+11-latgrb}. In this paper we
choose four of such bright LAT GRBs, namely, GRB 080916C, 090510, 090902B and
090926A as our sample

We processed the Fermi GBM and LAT data using the same method as presented in
\citet{Zhang+11-latgrb}.  For each burst, we choose the 1-2 brightest GBM
spectra and 1 or 2 BGO spectra together with LAT spectrum to compare with our
model.  We have fixed the BGO normalization factor, accounting for differences
between detectors' effective area, to the usual value of $0.9$
\citep{Abdo+09-090902B}.

By varying the free parameters,
we get the best models which can reproduce the observed count spectrum.  For the total
duration of the bursts we managed to fit GRBs 080916C and 090510 with our model.
For GRBs 090926A and 090902B, our model did not give a satisfactory fit when
considering the entire duration of the burst, but time resolved spectra were well
fitted by our model. We put this on the account of the strong spectral evolution
present in these two GRBs as published in their initial discovery papers.

\subsection{Fitting method}
\label{sec:fitmet}
We have taken the theoretical model described above with 9 free parameters.  We
fit the data using the simplex method implemented in the {\it amoeba} routine
\citep[see ][]{nr} modified, in some of the cases with simulated annealing
(Erik Rosolowsky, private communication). This addition makes the method less
sensitive to local minima.  The parameters describe the spectrum  of a
photospheric model (either magnetic or baryonic).

The procedure is to convolve the model with the detectors' response functions
and compare this result with the observed counts.  For this purpose, we have
written our own IDL codes to handle both the observed count spectrum, the
instrumental response matrix (DRM) and read in the multi-parameter theoretical
model to carry out the fitting.

Because the number of fitting parameters is large, we introduce some additional
constraints, in accordance with the model. We require the deceleration time to
be of order of few seconds, which corresponds to the delay between the MeV and
GeV components. This is achieved by introducing a penalty in the goodness of
fit for $t_{\rm dec} \gtrsim $ few $s$.  We vary one by one each of the
parameters around the best fit to calculate their errors.  When the $\chi^2$
statistic increases by the amount characteristic for the one-sigma limit for
given degrees of freedom, we fix the range and thus the errors for that
parameter.

\begin{table}
\begin{center}
  \begin{tabular}{lcccc} \hline
  & 080916C(m) & bar &  090510(m) & bar    \\
\hline \hline
 $\epsilon_e$/0.01	&   $3.1^{+0.6}_{-0.7}$ 	&$2.0^{{+}1.7}_{{-}1.0}$	& $0.50^{+0.50}_{-0.25}$			& $8.3^{+8.3}_{-4.2}$		 \\
 $\eta$ 			&  $1277^{+111}_{-150}$ 	&$312^{+11}_{-8}$		& $318^{+77}_{-159}$			& $235^{+18}_{-14}$		 \\
 $ \epsilon_{\rm B,FS}/10^{-3}$
					& $56^{+54}_{-32}$			&$0.022^{+0.27}_{-0.017}$	& $4.4*$							& $0.49 *$			\\
 $ L_t [10^{53}\frac{\rm erg}{\rm s}]$
					&  $76\pm0.6$				&$20^{+1}_{-1.1}$		& $0.57{\pm 0.11}$				& $0.60^{+0.08}_{-0.09}$	  \\

 $ \alpha $		& $0.47^{+0.13}_{-0.14}$	&$0.95^{+0.06}_{-0.08}$ & $0.77^{+0.15}_{-0.20}$		& $0.55^{+0.20}_{-0.31}$	  \\
 $ \beta $			& $1.92^{+0.01}_{-0.01}$	&$2.11\pm0.01$			& $4.49^{+0.02}_{-1.04}$		& $4.42^{+0.07}_{-1.03}$	   \\
 $ n_{\rm e} [{\rm cm}^{-3}]  $
					&   $16^{+11}_{-9} $		&$64^*$ 				& $4.1\times 10^3 *$				& $1.1\times 10^5 *$	\\
 $ p $				&  $2.87^{+0.12}_{-0.08}$	&$2.1^{+2.9}_{-0.05}$	& $4.74^{+0.26}_{-2.73}$		& $2.02^{+2.98}_{-0.01}$	    \\
 $ \Gamma_r $		&   $1.32 \pm 0.05$		&$1.0\pm 0.03 $ 		& $2.24^{+0.16}_{-0.13}$		& $1.00^{+0.07}_{-0.0}$ 	    \\

 $\chi^2$/dof		& $ 113.6 /  98$			&$129.7 / 98 $			& $107.2 /  97$ 				& $ 107.1 /  97$		    \\
\hline
\end{tabular}
\caption{Model parameters for the time-integrated fits with a magnetically and a
baryonically dominated model for two bursts. Stars mark parameters for which we could not determine the errors. }
\end{center}
\end{table}

\noindent
{\it Parameters constrained:}
In a situation when the spectrum is described by either an extremal magnetic
($\mu=1/3$) or baryonic ($\mu=1$) model,  the varying parameters we set out to
constrain are: the break energy ($\varepsilon_b$, see Eq. \ref{eq:ebreak}) of the prompt
emission, the low- and high energy indices ($\alpha$ and $\beta$) of a Band
function ($N_{\varepsilon}(\varepsilon) \propto \varepsilon^{-\alpha} e^{-\varepsilon(2-\alpha)/\varepsilon_{b}}$ for
$\varepsilon < \varepsilon_b$ and $N_{\varepsilon}(\varepsilon) \propto \varepsilon^{-\beta}$ if $\varepsilon > \varepsilon_b$), the
relative Lorentz factor of the reconnecting regions in the photosphere
$\Gamma_r(\gtrsim 1)$, the total luminosity $L_t$, interstellar density
($n_e$), index of the electron power-law distribution ($p$, where $N_e(\gamma)
\propto \gamma^{-p}$) and the fraction of the energy in the magnetic field and
electrons at the deceleration radius ($\epsilon_{\rm B,FS}$ and $\epsilon_e$).
We do not to include variations in e.g. magnetic field in FS and RS because it
would not be possible to constrain them and we assume that $\epsilon_{\rm
B,FS}=\epsilon_{\rm B,RS}$. The components included in this work are: prompt
Band component (PR), blackbody component (BB), forward- and reverse shock
synchrotron (FS, RS), forward and reverse shock synchrotron self Compton
(FS-SSC, RS-SSC), external inverse Compton of the prompt photons on the forward
and the reverse shock (FSEIC, RSEIC).

The prompt emission is represented as a Band function, the forward and reverse
shock  synchrotron and inverse Compton components as multiple broken power
laws.  We have applied a phenomenological smoothing at the breaks of the type
$( (x/x_{\rm break})^{i_1 b}+(x/x_{\rm break})^{i_2 b})^{1/b}$.  Here $b$ is an
arbitrary constant describing the smoothness of the break (we set $b=-4$),
$i_{1,2}$ are the power law indices and the break occurs at $x_{\rm break}$.
The prompt emission may also have a high energy spectral cut-off, at an energy
depending on the acceleration mechanism at the photosphere, and whether or not
pair creation is present.  Here we have taken a phenomenological cut-off value
at approximately $\Gamma_{\rm ph} m_e c^2 /(1+z)$ (as also in
\citet{thompson94}).  Both without pair formation and for pair formation
estimated from the photons above threshold estimated from an extrapolation of
the observed spectral slopes in the bursts considered here, a cutoff close to
the above value is obtained. The energy coverage is generally sparse at a
putative prompt cutoff energy in the range of  tens of MeV, making it hard to
discern whether there is a pair cutoff or a cutoff due to the particle
acceleration mechanism.  For simplicity, here we use the model without pair
creation for the purposes of the numerical calculation.  This cutoff can
constrain the photospheric Lorentz factor ($\Gamma_{\rm ph}$).	Concerning the
GeV component, we identify the high energy cutoff (and peak) with the
Klein-Nishina cutoff.

The strength of the GeV component is significantly influenced  by the amplitude
of the prompt component (responsible for the intensity of seed photons for the
EIC scattering), $\epsilon_e$ and $n_{\rm ext}$, while	$\epsilon_{\rm B,FS}$
plays a role in some of the cases.  The shape of the GeV emission is defined by
a combination of $\alpha$ and $p$, which makes it possible to constrain these
parameters. The peak energy of the GeV emission depends on $\eta$, $\epsilon_e$
and $p$. The coasting Lorentz factor, $\eta$ can be constrained by the prompt
emission cutoff at tens of MeVs.

The model in general has difficulty in constraining $\epsilon_{\rm B,FS}$, as
this parameter affects only the  low-energy ($\simeq 10 \hbox{~keV}$) behavior of the
spectrum  through the FS synchrotron (and to a lesser extent by RS-SSC
emission). As a consequence this energy range can be used to constrain the
parameter, but in many cases, however, there are no signs of an extra component
emerging at low energies.  We note that, considering the small number of
high-energy photons, the fact that many of the combined GBM-LAT spectra are
consistent with a single Band function does not necessarily exclude additional
components.  The two components (photospheric and external shock-related) can
align to form a continuous component without  excessive fine tuning.  This
requirement is relaxed when we need to fit bursts  with additional components.

\section{Individual GRBs}
\label{sec:individual}

\subsection{GRB 090510 fit for time integrated spectrum}
\label{sec:090510}
This is a short burst \citep{Abdo+09-090510} with  GeV emission.  According to
our modelling the time integrated spectrum can be described by a Band function
and the high energy photons are dominated by the RSEIC emission. The low-energy
excess is best modelled by a FS component in the magnetic case and by a
blackbody component in the baryonic case. See Figs. \ref{fig:090510m} and
\ref{fig:090510b}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.7\columnwidth]{fig1a.eps}\\
\includegraphics[width=.7\columnwidth]{fig1b.eps}\\
\caption{GRB 090510 fit with a magnetic model ($\chi^2/$dof=107.6/97). In the
upper figure lines represent the model. The points with the error bars can be
thought of as measured values, though they are themselves model dependent (i.e.
for different models they are not the same). They are calculated from the model
convolved with the detectors' responses (step curve on lower figure) and their
relation to the observed counts (crosses on lower figure). For brevity we do
not show the plot with the convolved spectra for other bursts.	}
\label{fig:090510m} \end{center} \end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.7\columnwidth]{fig2a.eps}\\
\includegraphics[width=.7\columnwidth]{fig2b.eps}
\caption{GRB 090510 fit with a baryonic model ($\chi^2/$dof=107.1/97). Other
notations are the same as Fig. \ref{fig:090510m}  }
\label{fig:090510b} \end{center} \end{figure}

\subsection{GRB 080916C fit for time integrated spectrum}
\label{sec:080916c}
GRB 080916C is one of the most intensive bursts ever detected
\citep{Abdo+09-080916c}. Its spectrum is consistent with the Band function
throughout the Fermi energy range and in all the analyzed intervals (though
there was weak evidence of an extra component \citep{Abdo+09-080916c}).  The
{extremal} magnetic and baryonic fits are shown in Fig. \ref{fig:080916cm}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.7\columnwidth]{fig3a.eps}\\
\includegraphics[width=.7\columnwidth]{fig3b.eps}
\caption{
Fit for GRB 080916C, for the entire interval, with a magnetic model (above,
$\chi^2=113.6$, d.o.f.=$98$) and baryonic model (below, $\chi^2=129.7$,
d.o.f.=$98$). The figures show the models and the data points obtained by
scaling the measured count rates with respect to the  convolved models.
}\label{fig:080916cm} \end{center} \end{figure}

\subsection{GRB 090926A fit for the time resolved spectrum}
\label{sec:090926}

\citet{Ackermann+11-090926} fitted the	spectrum of this burst	with a Band
function and an extra power-law with a high energy cutoff.  Using data from
\citet{Zhang+11-latgrb} the whole interval of the burst cannot be well fit with
a simple analytical model, suggesting strong spectral evolution. We thus use
the time resolved data for this GRB.
We divided the burst in 4 intervals (with limits at 0, 4, 7, 11.5 and 19 s with
respect to the trigger time). The first interval has no GeV range emission, and
we drop it from our analysis.

We found that both the magnetic and the baryonic photosphere model can give an
accurate description of the resolved data (see Fig. \ref{fig:090926Amultim} for
a superposition of magnetic and baryonic models and the models and data for the
first interval on Fig. \ref{fig:090926A1} ).

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.70\columnwidth]{fig4a.eps}\\
\includegraphics[width=.70\columnwidth]{fig4b.eps}
\caption{GRB 090926A with three intervals fitted using the magnetic (top) and
baryonic (bottom) photosphere model. We have displayed only a small number of
components. The numbers indicate the intervals: 1: 4-7 s, 2: 7-11 and 3:
11.5-19 s. The 0-4 s period could not be fitted, because no LAT photons were
detected.  } \label{fig:090926Amultim} \end{center} \end{figure}

Reduced $\chi^2_r$ values are in the range $1.4-1.8$ for the magnetic model,
and $1.2-2.5$ in the baryonic case. Generally, $\epsilon_e$, $\epsilon_B$ and
terminal Lorentz factors ($\eta$) are systematically lower in the baryonic case,
the external density is systematically higher, while the other parameters are
roughly in the same range for both cases. We note as a general trend, the
magnetic model yielded better fits for the data.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.70\columnwidth]{fig5a.eps}\\
\includegraphics[width=.70\columnwidth]{fig5b.eps}
\caption{Examples of magnetic (above, $\chi^2=182.9$, d.o.f.=$115$) and
baryonic (below, $\chi^2=184.5$, d.o.f.=$115$) models with the data overplotted
for GRB 090926A's  4-7 s interval (see also Fig. \ref{fig:090926Amultim}). }
\label{fig:090926A1} \end{center} \end{figure}

\subsection{GRB 090902B fit for the time resolved spectra}
\label{sec:090902}

GRB 090902B \citep{Abdo+09-090902B} has a relatively sharp Band spectrum with
an underlying power law component. While a simple Planck function is too narrow
to fit the spectrum, a multicolor blackbody and an underlying power law can fit
the data equally well \citep{Ryde+10phot090902}. The fit for the entire
duration of the burst did not give an acceptable fit for neither of the models,
possibly due to the strong spectral evolution.	It has been shown that this
burst involves a strong spectral evolution \citep[e.g.][]{Zhang+11-latgrb}

In order to account for spectral evolution, we have defined $10$ intervals
bracketed by the $0, 5, 7.3, 9.7, 12, 14.3, 16.7, 19, 21.3, 23.7$ and $26$ s
marks relative to trigger (see Fig. \ref{fig:090902B}).

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\columnwidth]{fig6a.eps}\\
\includegraphics[width=0.7\columnwidth]{fig6b.eps}
\caption{GRB 090902B magnetic models (upper) and baryonic (lower). The time
intervals depicted: 0-5 s: gray-dotted, 7.3-9.7 s: black-long dashed, 9.7-12 s:
gray-dash-dotted, 12-14.3 s: gray-dash-dot-dot-dot, 14.3-16.7 s:
black-continuous, 16.7-19 s: gray-continuous, 19-21.3 s: black-dashed,
21.3-23.7 s: light gray-dash-dotted.  }\label{fig:090902B} \end{center}
\end{figure}

$\epsilon_e$ varies in the range $\approx 0.001-0.7$, and it is similar for
both models.  The terminal Lorentz factor is systematically higher in the
magnetic case and it is above $400$ in all intervals, while in the baryonic
case, it is below $250$. The external density, is systematically lower in the
magnetic case, while the magnetic parameter at the deceleration radius and the
semirelativistic shocks' Lorentz factor in the photosphere are systematically
higher in the magnetic case.	These parameters are all within the normal
ranges found previously from fitting afterglow and prompt emission data on a
wider sample of bursts.  Here, we do not find any obvious trend for an
evolution of the parameters throughout this burst. For an example fit for this
burst's first interval, see Fig. \ref{fig:090902B1}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.70\columnwidth]{fig7a.eps}\\
\includegraphics[width=.70\columnwidth]{fig7b.eps}
\caption{Examples of magnetic (above, $\chi^2=102.3$, d.o.f.=$92$) and
baryonic (below, $\chi^2=96.7$, d.o.f.=$92$) models with the data overplotted
for GRB 090909B's  0-5 s interval (see also Fig. \ref{fig:090902B} -
gray-dotted curve). } \label{fig:090902B1} \end{center} \end{figure}

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

{\it Possible relation between $\mu$ and observed trends.-}
In our fits we have purposefully avoided introducing any of the observed
phenomenological correlations claimed for GRBs, since these have errors and
biases which require extensive consideration and would complicate the analysis.
However, it is interesting to speculate on the possible values of $\mu$ in the
context of some of these phenomenological correlations.

{In the $\eta>\eta_T$ case it seems simple considerations do not yield a
consistent value for $\mu$, and this needs to be examined further.  In the
other case ($\eta<\eta_T$), there is no $\mu$ dependence, but we only need to
assume one of the relations to obtain approximately the other one; e.g. if we
assume $\eta\propto L^{0.3}$ \citep{Lu+12corr}, we get $\varepsilon_b\propto L^{0.4}$,
fairly close to the Yonetoku relation.}

Although we did not find any cases where the main MeV episode was provided by
the thermal component, {such cases have been discussed based based on more
elaborate radiation transfer models, e.g. by \citet{Peer+11multicolor,
Beloborodov12peak, Vurm+12therm} and others, typically for baryonic dynamics.
In our simplified treatment the thermal component was taken as a simple
blackbody function, instead of a broader multicolor blackbody as in, e.g.
\citep{Ryde+10phot090902}. Such multicolor fits including our other variables
would be much more cumbersome, but in principle, if we make a quick calculation
similar to that in the previous paragraph, using Eq.  \ref{eq:temp} instead of
Eq. \ref{eq:ebreak}, this would imply $\mu\approx 0.5$.  }

It is interesting to note that for the magnetically dominated model
($\eta>\eta_T$ and $\mu=1/3$), the peak energy does not depend on either the
luminosity and the terminal Lorentz factor, and has a value of $\sim 300 \hbox{~keV}$,
corresponding to the clustering, observed by e.g. BATSE
\citep{Kaneko+06BATSEcat}.  There is a weak dependence on $r_0$ and redshift.

{\it Comparison of the baryonic photosphere model with the magnetic model.-}
The main difference between the magnetic and the baryonic model in terms of the
physical parameters involved in the GRBs is the initial dependence of the
Lorentz factor on the radius (see equ. \ref{eq:accel}).

The main manifestation of this difference in terms of the observed spectrum is
the strength of the thermal component from the photosphere. Because of the
different dynamics, the magnetic model results in a sub-dominant thermal
component, while in the case of the baryonic dynamics the black body component
can have a more prominent contribution.

The main goal of this article has been to present a generalized formalism for
models with an arbitrary acceleration of the GRB ejecta, valid for magnetic and baryonic
dynamics. By means of an analysis of the brightest bursts, and properly accounting
for the instrumental effects, we have shown that both an {extremal}  magnetic and a
baryonic model can provide adequate fits to the spectrum. Although in some of the cases
the magnetic model provides a marginally better fit to the data, the difference is
not significant enough to decide in favor of either extremal model.

{\it Constraints on the parameters.-}
As we can see from Table 1. we can constrain some of the physical
parameters of the models, except for the magnetic parameter at the deceleration
radius, and in some cases the interstellar density is also difficult to constrain.

We found that the burst spectra are consistent with both the magnetic and the
baryonic models presented here.  A possible discrimination would require further
theoretical considerations, or else other observational considerations; e.g. the
length of the deceleration time should be the same order as the observed delay of the GeV
emission. Another discriminant would be how well a model can reproduce the observed
correlations between the rest-frame peak energy and e.g. luminosity would.  As we
indicated above, the extreme magnetic photosphere model appears to be  better in
accommodating these constraints.

One of the main differences between the two models is the systematically higher
value of the coasting Lorentz factor (LF) in the magnetic model. This can be understood
as follows: by applying two models to the same data the photospheric LF is the
best constrained by the observations. While in the baryonic case this is just
the terminal LF, in the magnetic case the ejecta is still accelerating at the
photosphere, having a larger coasting LF.

In GRB 090510 \citet{Abdo+09-090510} noted a correlation between the keV range
radiation and the GeV radiation (but no correlation with the MeV radiation).
Based on their spectrum, the additional power-law would be the responsible for
explaining this correlation.
In the framework of this generalized photospheric model, such a correlation can be
explained if these two components have a common source: e.g.  the keV emission
is ascribed to the forward shock (FS) synchrotron and the GeV to the forward shock
synchrotron-self Compton (FS-SSC); or if a reverse shock is present (as proposed in
this work), the keV radiation can be due to the reverse shock synchrotron-self
Compton (RS-SSC) and the GeV is produced by the reverse shock external inverse
Compton (RS-EIC).

{\it Conclusions.-}
In summary, we have shown that a simple dissipative photosphere model of the
prompt emission of GRB based on synchrotron radiation can fit well the joint
Fermi GBM and LAT data for the four brightest bursts with sufficiently detailed
information. We have introduced a general formalism for treating the jet
acceleration dynamics ranging between an extreme magnetically dominated to a
baryon dominated dynamics regime. Both types of dynamics give acceptable fits
to the data on these four bursts. Additional considerations involving statistical
trends over a much larger sample of bursts might favor slightly the magnetic
models, but firm conclusions would require detailed fits over much larger samples
than undertaken here.

We acknowledge partial support from NASA NNX09AL40G and OTKA K077795 (PV, PM)
and NASA SAO SV4-74018 and NASA SAO GO1-12102X (BBZ). We thank Kazumi Kashiyama, Yi-Zhong Fan and the anonymous referee for comments and suggestions.