\section{Introduction}
The Green Bank North Celestial Cap (GBNCC; \citealt{slr+14}) pulsar survey began in 2009 and, when complete, will cover the entire sky accessible to the 100\,m Robert C. Byrd Green Bank Telescope (GBT; $\delta\geq-40\degr$, or 85\% of the celestial sphere) at 350\,MHz. As of mid-2019, the survey is 85\% complete and 161 pulsars have been discovered, including 25 millisecond pulsars (MSPs) and 16 rotating radio transients \citep[RRATs;][]{mll+06}. Timing solutions for these discoveries have been published in \cite{slr+14}, \cite{kkl+15}, \cite{kmlk+18},  \cite{lsk+18}, and \citet{acd+19}, and more are forthcoming.  As such, this constitutes one of the largest and most uniform pulsar surveys to date.

In addition to the newly discovered pulsars, the uniform coverage of GBNCC allows a robust re-assessment of the known pulsar population with reliable flux density measurements.  Here we present a detailed search for all known pulsars in the GBNCC footprint.
We find that 572 previously published pulsars and 98 unpublished pulsars have been re-detected by the survey pipeline and visually confirmed, comprising 670 detections in total, the largest low-frequency, single-survey sample.  Similar to previous efforts based on results from the Parkes Multibeam Pulsar Survey (PMPS) and the Pulsar Arecibo L-band Feed Array (PALFA) survey \citep[e.g., see][]{lfl+06,slm+14,lbh+15}, we conduct a detailed  analysis of the GBNCC pulsar survey and compare its sensitivity with that of other surveys in overlapping regions of sky. Flux densities at 350\,MHz ($S_{350}$) are presented for all detections, as well as
pulse widths and profiles.

In \S\ref{sec:methods}, we outline the process used to generate a comprehensive list of pulsars as well as predicting and measuring signal-to-noise ratios (S/N) of detections in the survey. In \S\ref{sec:results}, we present the recovered S/N and flux density measurements for all detected pulsars as well as measurements of pulse width, dispersion measure, and spectral index. We also present the profiles for all of these pulsars. In \S\ref{sec:sensitivity}, we discuss how the GBNCC survey is performing compared to expectations and RFI characteristics of the survey, and remark on interesting detections and notable non-detections. We also discuss the implications of our results for the Galactic pulsar population. Finally, in \S\ref{sec:conc}, we summarize the main conclusions of this analysis.

\section{Sample Assembly and Data Reduction}

\label{sec:methods}
The GBNCC data set as of late fall 2018 included $\sim$108,000 120\,s pointings, each tagged with a unique beam number. Each dual-polarization observation was taken with the GBT over the past $\simeq$10 years. The survey utilizes the GUPPI backend, with a sampling time of 82\,$\mu$s and 100\,MHz of bandwidth centered at 350\,MHz \citep[for more information on the observing setup for the GBNCC survey, see][]{slr+14}. We began by organizing a comprehensive list of all known pulsars with parameters that were available for use, whether they were published or not. By utilizing the Australia Telescope National Facility (ATNF) pulsar catalog\footnote{\url{http://www.atnf.csiro.au/research/pulsar/psrcat}} \citep[v1.59,][]{mhth05}, we amassed the bulk of the sources from the list of all published pulsars and their positions on the sky as well as their spin parameters and other relevant quantities (dispersion measure, etc.). Discovery parameters are also available for additional pulsars that have not been published but were detected in a number of other recent or ongoing surveys. Many of these surveys, including AODrift \citep{dsm+13}, SUPERB \citep[][Spiewak et al., 2019, in prep.]{kbj+18}, GBT 350\,MHz Drift \citep{blr+13}, PALFA \citep{cfl+06,lbh+15}, LOTAAS \citep{scb+19}, and HTRU-South \citep{kjs+10} include pulsars that are in GBNCC survey area, and so were included in the list. More information on these surveys is included in Table \ref{tbl:surv}. Furthermore, we included the list of pulsars that had been discovered in the search pipeline for the GBNCC survey. We then limited this list to pulsars within the range of the survey, i.e., pulsars with $\delta>-40\degr$.  In total, this list contained 2299 pulsars. We determined which pulsars were within 30\arcmin\ (FWHM of GBT at 350\,MHZ) of completed GBNCC pointings, adjusting when necessary to compensate for large ($>$30\arcmin) uncertainties in pulsar position. This reduced the total number of pulsar candidates to 1413. We could then match each pulsar with the GBNCC beams closest to its position before beginning to process the data.

Radio frequency interference (RFI) excision is the first step of GBNCC data analysis, and is done primarily with the \texttt{rfifind} tool from the \texttt{PRESTO}\footnote{\url{http://www.cv.nrao.edu/~sransom/presto/}} pulsar data analysis software package \citep{smr+01} as described in \S3.1 of \citet{slr+14}. We also performed an analysis of the \texttt{rfifind} output files spanning the lifetime of the GBNCC survey up to late 2018 (roughly 83\% of the total survey) to characterize the effects of RFI over the course of the survey. These files contain information about which frequency channels were masked due to RFI for every 120 second scan in the survey. For a particular scan, the effective bandwidth $\Delta\nu$ is the total 100\,MHz bandwidth of the GBT 350\,MHz receiver multiplied by the ratio of unmasked to total channels for that scan, minus an additional 20\,MHz for rolloff.

In some cases, the \texttt{rfifind} masks were insufficient to remove additional RFI that was either narrow in frequency space or brief in time. The latter often appears as a very bright burst at $\sim$0 DM for portions of the observation. To mitigate this, we employed some additional narrowband flagging in the \texttt{PRESTO} \texttt{prepfold} command as well as removing corrupted portions of the scan in the time-domain. Note that these changes also alter the values for $\tau_{\rm obs}$ and $\Delta\nu$ which consequently change the measured S/N for a given observation. For this reason, we calculate the fraction of data points from the observations that were not omitted in processing and multiply the total bandwidth by this fraction.

After removal of RFI, we dedispersed and folded the observations at each pulsar's rotational period and integrated the profiles to obtain a single average profile for each observation. For the vast majority of sources included in this analysis, a precise ephemeris from the ATNF catalog was used to perform the folding. In all other cases, only the discovery parameters (period, DM and, if known, period derivative) were used. We also repeated this process while allowing dispersion measure to vary and, in some cases, also allowing variations in period and period derivative. This second iteration allows for fine-tuning previously published parameters at the cost of potentially finding bright RFI, which will often occur when attempting to detect low-DM pulsars as sources of RFI have DM = 0\,pc\,cm$^{-3}$. The 120 second observation times utilized in the GBNCC survey limit sensitivity to period refinement, so fitting for period was only used to increase the S/N of detections of pulsars for which only discovery parameters were used, and no further timing analysis was done as a part of this study. All folded data were visually inspected to determine likelihood of an actual detection. In cases where RFI still existed in the data, we removed high order ($>$5) polynomials from the off-pulse regions of the profile. With folded profiles, we calculated a measured signal-to-noise ratio (S/N) \citep{lk04},
\begin{equation}
{\rm S/N}_{\rm meas} = \sum_{i=0}^{N_{\rm bin}} \frac{p_i-\bar{p}_{\rm off}}{\sigma_{\rm off} \sqrt{W N_{\rm bin}/P}}~\gamma,
\label{eq:snm}
\end{equation}
where $N_{\rm bin}$ is the number of bins across the pulse profile,
$p_i$ is the value of bin $i$,
$\bar{p}_{\rm off}$ is the mean of the off-pulse bins,
$\sigma_{\rm off}$ is the standard deviation of the off-pulse bins, $W$ is the on-pulse width in seconds, $P$ is the pulsar spin period in seconds, and $\gamma$ is a correction factor. When continuous signals are assigned to a finite number of bins in the profile during the folding process in \texttt{PRESTO}, their intensity is ``smeared" over the neighboring bins, resulting in correlations in the bins' intensities. This correction, dubbed $\gamma$, depends on the sampling time and the number of bins in the profile, which (for this study) is dependent on the pulsar spin period. Typical values are close to 0.95. The number of bins $N_{\rm bin}$ was determined by the pulsar period as follows: profiles for pulsars with periods shorter than 1.7\,ms had 28 bins, periods shorter than 10\,ms had 50 bins, periods shorter than 50\,ms had 128 bins, and all others had 200 bins. This prescription retains sensitivity to long-period pulsars but avoids bin widths corresponding to time intervals smaller than the sampling time of 82\,$\mu$s. Pulse widths were determined with a standard process. First, sigma-clipping was used to find the off-pulse region. Then, the peak value above the noise floor was identified, and bins on either side of the peak were added to the on-pulse width. This process was repeated, adding bins on the sides of the peak until we reached bins within 2$\sigma$ of the mean of the noise. The edges of the pulse were found by fitting lines to the two bins on either side of the pulse and finding the fraction of the outermost bins that were above the noise floor. At this point, we consider the full on-pulse width to be determined. Each profile was then checked by eye, and corrections to the on-pulse region were made. Any components of the pulse width that were distinct from the main pulse were determined using the same algorithm. To determine the sensitivity of uncertainties in S/N from the choice of the number of on-pulse bins, noisy Gaussian pulses were simulated and various width choices were used to measure the fractional error on S/N. From this test, it was found that on-pulse widths that exceed at least one $\sigma$ beyond the Gaussian mean were sufficient to greatly reduce the fractional uncertainty on S/N. Beyond this, adding bins had little effect on this fractional uncertainty - so, pulse widths were chosen to encompass all of the pulse visible above the noise.  In some cases, additional RFI features were removed prior to the determination of $W$ to minimize errors in $W$ and S/N (see \S\ref{sec:rfi}).

Characteristic measurements of pulse width include measurements at both 50\% and 10\% of the pulse profile's maximum amplitude (hereafter $W_{50}$ and $W_{10}$, respectively). These widths are dependent on both pulse period and observing frequency, so measurements at 350\,MHz help to fill out the low-frequency regime for a wide range of pulse periods. However, the noise floor in some pulsars limits the ability to determine $W_{10}$ robustly. Note also that $W_{50}$ and $W_{10}$ are distinct from $W$, which includes all bins that contain the pulse signal, and so $W$ is generally slightly larger than $W_{10}$.

The expected S/N of a pulsar  can be estimated as \citep{dtw+85,lk04}
\begin{equation}
{\rm S/N}_{\rm exp} = \frac{ S_{350}G\sqrt{N_{\rm pol}\tau_{\rm obs}\Delta\nu}} {T_{\rm sys}\beta} \sqrt{\frac{P-W}{W}} f(\theta),\label{eq:sne}
\end{equation}
where $S_{350}$ is the flux density at 350\,MHz,
$G=2$\,K/Jy is the gain of the Green Bank Telescope \citep{slr+14},
$N_{\rm pol}=2$ is the number of polarizations recorded,
$\tau_{\rm obs}=120$\,s is the length of the observation,
$\Delta\nu$ is the bandwidth in MHz after removing RFI (see \S\ref{sec:rfi}),
$T_{\rm sys}$ is the system temperature (including the sky temperature at the source position, receiver temperature $\simeq$ 20\,K, and CMB temperature $\simeq$ 3\,K),
$\beta\simeq1.1$ is an instrument-dependent correction factor due to downsampling the data to 2 bits \citep{lk04},
and $f(\theta)$ is a radial Gaussian factor accounting for sensitivity degradation as a function of angular offset from the center of the circular beam $\theta$. The sky temperature in the direction of each pulsar was determined by using the measurements made by \cite{hks+81} for the beam positions, scaled to 350\,MHz using with the spectral index therein, $-$2.6.

Where possible we use flux densities at other frequencies and previous measurements of spectral index ($\alpha$, with $S_\nu \propto \nu^{\alpha}$) from the ATNF catalog to determine an expected flux density at 350\,MHz and the expected S/N \citep{mhth05}. In cases where there was no published value for $\alpha$ but flux densities at both 400\,MHz and 1400\,MHz were published, we determine a spectral index using a simple power law. In all other cases, we assume a spectral index of $-1.4$ \citep{blr+14} to estimate the flux density at 350\,MHz. We also calculate the measured flux density of each pulsar by inverting Equation \ref{eq:sne} and using measured values for S/N (determined from Equation \ref{eq:snm}) and pulse width. Comparing the expected flux density to our measurements can both roughly confirm our current models for pulsar emission as well as aid in explaining non-detections.

\section{Pulsar Flux Density Census at 350\,MHz}
\label{sec:results}

\begin{figure*}[ht]
\centering
\plotone{psr_map_Gal_012820.pdf}
\caption{Sky map with pulsars from overlapping surveys, plotted in Galactic coordinates as a Mollweide projection. The shaded regions indicate completed GBNCC observations. Detected pulsars from the ATNF catalog and pulsars that were detected using discovery parameters from overlapping surveys are differentiated by marker type, with green plus symbols indicating pulsars from the catalog and red triangles indicating pulsars from the surveys listed in Table \ref{tbl:surv}. Pulsars that were not detected are plotted as blue ``x" symbols.}\label{fig:skymap}
\end{figure*}

We detected 670 pulsars out of a total of 1413 in the survey area, and these detections are listed in Table \ref{tbl:det} in Appendix 1. For all following analysis, the beams corresponding to the brightest detections (highest S/N) were used, as these are most likely to represent the pulsars' flux density. Along with pulsar names, we provide several relevant quantities: dispersion measure from searching with \texttt{PRESTO} \citep{smr+01}, MJD of the brightest detection, angular offset from the center of the beam, $W_{50}$, $W_{10}$ (when S/N was large enough), detection S/N, 350\,MHz flux density measured from the GBNCC data, and measured spectral index $\alpha$ (see \S\ref{sec:spindx}). Uncertainties on the S/N and flux densities were calculated using standard error propagation from equations \ref{eq:snm} and \ref{eq:sne} and uncertainties on bandwidth, temperature, and $\theta$ of 5\,MHz, 10\,K, and 0.5 degrees, respectively. Among these are 66 millisecond pulsars (MSP), defined here as pulsars with spin periods shorter than 30\,ms. The integrated pulse profiles for all of the brightest detections are shown in Figures \ref{fig:prof0} $-$ \ref{fig:prof11} along with pulsar names, dispersion measure, and flux density. Figure \ref{fig:skymap} shows all detected pulsars plotted by their Galactic positions, and different markers indicate whether or not the pulsars were from the ATNF catalog or were a part of one of the other survey lists mentioned above.

\subsection{Comparison Between the GBNCC and Overlapping Pulsar Surveys}

\begin{deluxetable*}{lcccc}
      \tablecaption{Pulsar Survey Comparison \label{tbl:surv}}
      \tablehead{
\colhead{Survey} & \colhead{Central Frequency} & \colhead{Limiting Flux Density\tablenotemark{a}} & \colhead{Detections\tablenotemark{b}} & \colhead{Reference}  \\
& \colhead{(MHz)} & \colhead{(mJy)} & & }
\startdata
AODrift \dotfill & \phn327 & 0.59  & 7/13 & \cite{dsm+13}\\
HTRU$-$S (low-lat) \dotfill & 1352 & 0.40   & 0/9\phn & \cite{kjs+10}\\
HTRU$-$S (med-lat) \dotfill & 1352 &  0.95   & 3/27 & \cite{kjs+10}\\
HTRU$-$S (high-lat) \dotfill & 1352 & 1.2\phn   & 1/8\phn & \cite{kjs+10}\\
SUPERB \dotfill & 1352 & 0.4\phn  & 2/15  & \cite{kbj+18}, Spiewak et al. (in prep)\\
LOTAAS  \dotfill& \phn134 & 0.63 & 10/39\phn & \cite{scb+19} \\
PALFA  \dotfill& 1400 & 0.23 & 0/29 & \cite{lbh+15}  \\
GBT350 \dotfill & \phn350 & 0.59  & 3/6\phn & \cite{blr+13} \\
GBNCC \dotfill & \phn350 & 0.70  & 72/72\phn & \cite{slr+14} \\
\enddata
\tablenotetext{a}{Averaged over the survey area and scaled to 350\,MHz.}
\tablenotetext{b}{Number of detections of pulsars from this survey by GBNCC/number of pulsars from this survey within the GBNCC survey area.}
\tablecomments{Information about individual detections is reported in Table \ref{tbl:det}.}
\end{deluxetable*}

Out of the 210 pulsars with discovery parameters that are not currently listed in the ATNF catalog, 98 were detected. Names, central frequencies, scaled limiting flux densities, and the ratio of detected to processed pulsars are given for each survey in Table \ref{tbl:surv}. It should be noted that there are many pulsars from these surveys (excluding GBNCC) in regions of the sky where the GBNCC survey has yet to observe, and so they may be detected in the future; these pulsars are not included in the counts listed in Table \ref{tbl:surv}. Three of these surveys (SUPERB, HTRU-S, PALFA) were conducted at higher frequencies,, where average sky temperature (especially near the Galactic plane) is much lower. This reason and the increased sensitivity to high DM pulsars at high frequency is useful for diagnosing missed detections. Because these pulsars have neither published flux densities nor spectral indices, reasons for missed detections cannot be determined more robustly than those due to sky temperature, position relative to the survey, extreme nulling/intensity variation, and high DM/short periods. It is also possible that for some of these pulsars, the discovery parameters may not be precise enough to be found in this analysis.

The most surprising missed detections come from the GBT350, AODrift, and LOTAAS surveys, which all have comparable sensitivities and frequencies. In an effort to explain why these pulsars were missed, all of the discovery plots were checked against our results, and acceleration searches were run. Three pulsars (J0100+69 and J0121+14 from LOTAAS, and J1854+36 from AODrift) that were originally missed were found on the second trial, as the DM used in the first run was not close enough to the DM at which the pulsar was discovered. For the majority of pulsars that were not detected after re-running the pipeline, the discoveries were quite dim. The LOTAAS survey also has much longer integration times (60 minutes), which significantly improves the chances of the survey detecting pulsars which may be faint and/or nulling. When checking the discovery plots, it became clear that both of these effects were common to many of the missed pulsars. Some pulsars even appeared to exhibit nulling with `off' times as large as 100 seconds. Nulling behavior was also seen in many cases for the AODrift survey. For the GBT350 missed pulsars, all three of those that were missed were faint, and several GBNCC beams in which the pulsars were most likely to be found had RFI that spanned the entire 100\,MHz band.

Eight binary pulsars that were originally discovered in the GBNCC survey were not detected in the first pass of this pipeline. These pulsars required acceleration searches, which are automatically performed as a part of the search pipeline, but not here. As a part of the missed pulsar analysis, we ran an additional acceleration search using \texttt{ACCELSEARCH} from within the \texttt{PRESTO} package, and they were all detected. We also reprocessed data for 15 binary pulsars from the ATNF catalog with short ($\leq$ 0.5 day) orbital periods that were not detected in the first pass using acceleration searches; none of these were detected.

Pulsars with long periods (greater than 2.5\,s) were also followed up with a search for single pulses. Because these pulsars would only be observed for at most 48 pulses, non-detections are more common. To address this, we implemented \texttt{single\_pulse\_search.py} from the \texttt{PRESTO} package, which searches a range of dispersion measures to find bright single pulses in the data and characterize them by their S/N. In this way, a pulsar that is not detected via a periodicity search may be found by individual pulses. However, we were still unable to find these pulsars using this method.

\subsection{Spectral Indices}
\label{sec:spindx}
\begin{deluxetable}{cccc}
\tablecaption{Broken Powerlaw Spectral Indices \label{tbl:spindx}}
\tablehead{\colhead{PSR} & \colhead{$\alpha_{l}$\tablenotemark{a}} & \colhead{$\alpha_{h}$\tablenotemark{a}} & \colhead{Break Frequency} \\
 & & & (MHz)}
\startdata
J0034-0534 & 0.6(3) & $-$3.1(2) & 181 \\
J0218+4232 & 1.15(7) & $-$2.7(4) & 149 \\
J1900-2600 & 0.2(4) & $-$1.89(15) & 204 \\
J2002+4050 & 0.2(16) & $-$1.51(18) & 378 \\
\enddata
\tablenotetext{a}{Spectral indices below ($\alpha_l$) and above ($\alpha_h$) the break.}
\tablecomments{Quantities in parentheses are uncertainties in the last digit.  See Figure~\ref{fig:spindx} for the corresponding plots.}
\end{deluxetable}

Many previously published spectral indices were determined from flux measurements from high-frequency surveys \citep[e.g., see][]{jsk+18}. Therefore, low frequency surveys like the GBNCC survey provide more stringent constraints on these calculations.  Results from this analysis are listed in Table \ref{tbl:det}.
The majority of the pulsars in this data set follow a single power law, or do not have enough ($>$2) flux density measurements to fit multiple power law functions.  However, there are a small number of cases where the emission is better fit by a broken power law, defined instead as a piecewise function composed of two power laws.
All 339 pulsars for which we measured spectral index had three or more flux measurements (including our 350\,MHz measurements) and were checked by eye to determine whether or not a broken line fit was appropriate. Four pulsars fit these criteria. For these pulsars, we fit two lines, one for high frequency flux density measurements and one for low frequency. The breaking point for the power law was determined by finding the maximum change in the derivative of flux density with respect to frequency. A similar analysis was done in \cite{mkb+17}. Plots of these cases are provided in Figure \ref{fig:spindx} with both indices included. These plots also display the best-fit line to all measured flux densities. The measured values of $\alpha_{l}$ and $\alpha_{h}$ are reported in Table \ref{tbl:spindx}.

\begin{figure*}[ht]
\epsscale{1.}
\plotone{broken_spindx.pdf}
    \caption{Pulsars with broken power-law spectral indices. We plot all available measurements of flux density in the ATNF catalog as well as the 350\,MHz measurements made in this study against observing frequency. We  fit two disjoint lines to the low- and high-frequency measurements (orange solid lines). The red dashed line indicates the frequency of the turnover in the spectrum, determined by finding the point at which the two lines match up. Information for these measurements is presented in Table \ref{tbl:spindx}.}
    \label{fig:spindx}
\end{figure*}

\subsection{Comparison of Dispersion Measure with Catalog Values}

The relatively low frequency of the GBNCC survey allows much higher precision DM measurements than typical 1400-MHz surveys, as dispersion across the band scales as $\nu^{-2}$. As pulses propagate through the interstellar medium, this dispersion results in a frequency dependent delay that smears out the arrival time of the pulse. Tools within the \texttt{PRESTO} package adjust for this, shifting the low frequency portion of the signal back in time to line up the pulse across the band. Using the \texttt{dmsearch} flag contained within the \texttt{PRESTO} command \texttt{prepfold}, we processed each of the pulsars and recovered more accurate values of DM. The program adjusts for dispersion and then folds the data at the pulsar's period to line up the pulses in both time and frequency. When \texttt{dmsearch} is off, the program does not tune the DM to maximize S/N; otherwise, the DM which aligns the pulses in frequency is returned as a new DM. In some cases, RFI caused the DM searching algorithm to return erroneous values for DM, and so we were unable to refine the dispersion measure. For these pulsars, we include the previously published DM in Table \ref{tbl:det} and mark them with a double dagger. More often, we were able to improve upon the previously published values of DM. Most of the discrepancies were small, but in some cases, our more precise DM measurement differed from the previous value significantly. For the pulsars with significant changes to their previously catalogued DM, we followed up with \texttt{TEMPO}\footnote{\url{http://tempo.sourceforge.net}} (maintained and distributed by Princeton University and the Australia Telescope National Facility). We split each detection into four subbands and created precise pulse times-of-arrival (TOAs) which can then be utilized to fit for DM. This method provides marginally more precise measurements, and so was only performed on pulsars with significant changes to previous DM measurements ($\geq$3$\sigma$). All newly measured DMs are presented in Table \ref{tbl:det}, and Table \ref{tbl:dm} highlights the pulsars which were followed up with \texttt{TEMPO} timing.

\section{Survey Sensitivity}
\label{sec:sensitivity}

\begin{figure*}
\plotone{GBNCC_senscurve_012820.pdf}
\caption{Flux density sensitivity in the GBNCC as a function of pulse period. Assuming a duty cycle of 6\% and an average unmasked bandwidth of 67\,MHz (which incorporates a 20\,MHz rolloff in the bandpass), we plot the predicted lower limit on the flux density of detectable pulsars for dispersion measures of 20, 50, 100, 150, 200, and 300\,pc\,cm$^{-3}$. To determine the sky temperature for the curves, we found the average sky temperature as a function of DM using the sky temperatures at the positions of all detected pulsars. We then drew from this function the temperatures at each DM for which a curve is plotted. For the above DMs, the function returns 95, 126, 171, 208, 237, and 273\,K. We glean the minimum detectable S/N for the survey by matching the curves to the faintest detection. This was found to be $\sim$3.8. Higher DM pulsars are more susceptible to smearing, and so the likelihood of detection is decreased for high DM, short period pulsars. We also plot both the detections (plus symbols) and non-detections (triangles), which are colored by their DM.}\label{fig:senscurve}
\end{figure*}

\begin{figure}
\epsscale{1.3}
\plotone{snm_vs_sne_120919.pdf}
\caption{Measured S/N vs.\ expected S/N for detections in the GBNCC survey. Extrinsic contributions to expected S/N include system temperature, telescope gain, scintillation, and offset from the beam center (newer pulsars without full timing solutions may have significant uncertainties in position). Errors in these quantities, previous flux measurements, and spectral indices increase the spread about unity, as does variable pulsar emission, i.e. nulling. }\label{fig:sn_mve}
\end{figure}

\subsection{Efficiency of GBNCC Survey}\label{subsec:efficiency}
In total, there were 5633 unique beams analyzed, yielding 1328 unique detections of the 670 pulsars. Given that there were 102948 beams that had been observed at the beginning of this project, this corresponds to an average number of detections per beam of $\sim$\,0.013 (0.063 detections per square degree), and $\sim$\,0.38 detections per hour of observing.
The ability to detect pulsars at 350\,MHz is limited most stringently by sky temperature and scattering in the interstellar medium (which correlates with dispersion). The expected S/N for detections is inversely proportional to system temperature, which is dominated by sky temperature near the Galactic plane. At 350\,MHz, this effect is quite significant, with temperatures approaching 1000\,K in this region. Scattering is especially detrimental in the detection of pulsars with short periods, as even a few milliseconds of smearing can eliminate the pulse entirely. Given a particular spin period and the estimated DM smearing, we can estimate the minimum flux density that will be detected by the survey. This relationship comes from solving Equation \ref{eq:sne} for flux density and assuming both an average sky temperature and duty cycle for the pulsars in the survey. Plotted in Figure \ref{fig:senscurve} are curves corresponding to a number of trial values of DM, showing the  sensitivity floor at those values. Because DM and sky temperature are correlated, we determined the average sky temperature for each curve that is plotted, resulting in an increase in minimum detectable signals for higher DM pulsars. Also plotted are flux density measurements for detections made by this survey and expected flux density measurements for the pulsars which were not successfully detected. The colors in the plot correspond to the dispersion measure of each pulsar, showing how pulsars that may be intrinsically bright enough to be detected can still be missed because of dispersive smearing and/or scattering. The minimum flux density expected to be measured in the survey (regardless of spin period) can be determined to be the asymptotic value of the DM curve corresponding to the faintest detection. This value is directly proportional to the minimum S/N which results in a detection, hereafter ${\rm S/N}_{\rm cut}$, which was found to be $\sim$\,3.8. For all detections, we plot both the expected S/N at 350\, MHz as well as the measured S/N of the detection. These are plotted in Figure \ref{fig:sn_mve} along with a line marking unity. There is a large spread about this line, due mostly to stochastic noise sources in the data (telescope noise, temperature fluctuations, scintillation, and variable pulsar emission). When examining these results, several of the more significant outliers were analyzed in closer detail. One of the three significant outliers in the lower right portion of the plot was found to be a new nulling candidate, and the other two were initially labeled as possible nullers that could not be verified without higher resolution observations.

Low frequency observations can result in significant deterioration of the pulse due to scattering and scintillation effects, as residual dispersive time delay within a frequency channel with finite width increases as $\nu^{-3}$ and scattering roughly as $\nu^{-4}$ \citep{lk04}. Both of these phenomena result in a broadened pulse and subsequently a reduction in S/N.
To shed light on the causes for some of the missed pulsars, we calculate the expected S/N using information from both the catalog and information about the beams in which we expect to detect them. We predict flux density at 350\,MHz calculated as described in \S 2, determine the masked fraction of the closest beam to the pulsar's position (when measured), and determine $T_{sys}$ for the corresponding sky position. To determine $W$, we fit a line to our measurements of $W_{10}$ as a function of spin period and draw from this function. This allows for a measurement of the spin period-pulse width relation at 350\,MHz, supplementing previous measurements at other frequencies. This best fit line was measured to be $W_{10}$ = $18.5\degr(4)P^{0.270(10)}$, which is consistent with the relation determined in \cite{jk19} modulo a frequency-dependent scaling factor \citep[for a more in-depth analysis, see][]{cw14}.  This fit is shown in Figure \ref{fig:pvw}.

\begin{figure}[ht]
    \centering
    \epsscale{1.2}
    \plotone{pvw.pdf}
    \caption{Pulse width at 10\% of the pulse maximum as a function of spin period. The solid line shows the line of best fit through the data, described by $W_{10}$ = $18.5\degr(4)P^{0.270(10)}$. The dashed line shows the minimum bin width as a function of period, as described in \S \ref{sec:methods}. }
    \label{fig:pvw}
\end{figure}

After drawing widths from either the catalog or the above function (based on availability of previous measurements of $W_{10}$ near 350\,MHz), we determined the expected S/N for all non-detections. These are plotted along with the measured S/N for all of the detections in Figure \ref{fig:sn_hist}. The detections have been divided between those found from the catalog and those discovered by the GBNCC survey, and non-detections are divided based on Galactic latitude. These divisions allow for direct comparison between the survey's ability to detect pulsars blindly as well as the limits placed on the survey by high temperatures and scattering near the Galactic plane. Included in the plot are three different S/N cutoffs placed during different stages of the survey. The least stringent cutoff of S/N = 12 comes from \cite{slr+14}, where it was used as an estimated cutoff for detection to predict the survey's sensitivity. At this S/N, $\simeq$75\% of non-detections are not expected to be detected. Pulsars close to the plane generally have lower S/N as the temperature is so high, while pulsars outside of the plane generally have smaller DM and temperature but more scintillation. The two detection curves show that the GBNCC is sensitive to intrinsically fainter pulsars, as the histogram is skewed toward lower measured S/N than those from the catalog. Note that there was one pulsar discovered by the GBNCC search pipeline with S/N = 5.98, which is the bin to the left of the search S/N cutoff.

\begin{figure*}[ht]
    \epsscale{1.2}
    \centering
    \plotone{sn_hist.pdf}
    \caption{Histograms of measured S/N for detections and expected S/N for non-detections. Detections are differentiated by GBNCC discovery/catalog pulsars (green/cyan lines), and non-detections by distance from the Galactic plane (the red line indicates pulsars that are within 10$\degr$ from the plane, and the orange line indicates pulsars outside of this region). The dashed lines indicate three different S/N cutoffs: the first line, in black, show the minimum S/N detected in the survey; the second, in grey, indicates the significance down to which candidates are folded in the GBNCC search pipeline; and the third, in blue, shows the predicted S/N limit used in \cite{slr+14} to predict sensitivity of the survey.}
    \label{fig:sn_hist}
\end{figure*}

In Figure \ref{fig:p0vsdm}, we plot all pulsars' periods against their dispersion measure. Each point's color and shape describe whether or not the pulsar was detected, and if not, whether we expect to have detected it. Missed detections that were unexpected are plotted with point sizes reflecting the expected flux density (calculated as described in \S\ref{sec:methods}) normalized by the value of the effective sensitivity curve for that pulsar, so larger points indicate pulsars with expected flux density much higher than the minimum detectable flux density at the pulsar's position.

\begin{figure*}[ht]
\epsscale{1}
\plotone{p0vsDM_120919.pdf}
    \caption{Period vs.\ DM for all included pulsars. Blue symbols indicate detections made by the survey, and red symbols indicate non-detections. Red triangles indicate missed pulsars that were not expected to be detected, in that they lie below the expected sensitivity of the survey. Red circles indicate missed pulsars that lie above their expected sensitivity, and so were unexpected non-detections (see \S\ref{subsec:efficiency} for details). Blue circles indicate detections that were expected, and blue x symbols indicate detection of pulsars with expected flux densities that were below our sensitivity limit. The area of these points is given by the ratio of expected flux density to the limiting flux density at the pulsar's position.}
    \label{fig:p0vsdm}
\end{figure*}

In total, there are 116 undetected pulsars plotted in Figure \ref{fig:p0vsdm} that have been classified as ``unexpected" by the logic above. Many of these pulsars are quite close to the sensitivity line, and so small errors in other flux density measurements and spectral indices may change them to ``expected." Because the effective sensitivity curve includes temperature and bandwidth (RFI, by proxy) information, reasons for missed detections are limited to effects that are harder to characterize. The most likely contributors include scintillation, abnormal pulsar behavior (i.e. nulling), and imprecise previous measurements of pulsar parameters resulting in inflated expected flux densities. Scintillation depends on DM \citep{cl91}, with increased timescales for smaller DM. Many of the non-detected pulsars that are outside of the Galactic plane are in this low-DM high-scintillation regime, and are likely to have been obscured (the expected number of scintles in the observation are on the order of $\sim$10). Many of the other missed detections, especially those from surveys with comparable limiting fluxes, were inspected individually. Many of these were obscured by significant RFI across the band. For example, PSR J0108$-$1431 (spin period of $\simeq0.81$\,s and DM of 2.38\,pc\, cm$^{-3}$, to the right of the bottom center of Figure \ref{fig:p0vsdm}) should be easily detected but was obscured by RFI. When examining a number of the other sources, it was found that many of the published spectral indices came from a 1400\,MHz study conducted by \cite{hwx+17}, and were unusually steep. This steepness results in high expected values of flux at 350\,MHz, which are not reflected in our results.

\subsection{RFI Analysis}\label{sec:rfi}
To visualize how RFI affects the efficiency of the survey, we determined the limiting flux density for each beam based on a S/N cutoff of 3.8, the temperature at the sky position of the beam, and the bandwidth available after RFI excision. Figure \ref{fig:flux_lim_hist} displays a histogram of the beams by their limiting flux, and Figure \ref{fig:flux_lim_map} shows these same data projected onto their sky positions. The sky map depicts a few important characteristics of the survey: the most obvious is the decreased sensitivity near the Galactic plane, but also visible are many individual pointings within the completed regions where significant RFI masking has reduced sensitivity. To mitigate this, these beams will be scheduled for reobserving. There is a small discrepancy between the number of observed beams displayed in Figures \ref{fig:skymap} and \ref{fig:flux_lim_map} due to a backlog of data which has yet to processed, and so mask fractions have not been determined for these beams.

\begin{figure}
\epsscale{1.2}
\plotone{GBNCC_fluxdensity_limit_hist.pdf}
    \caption{Cumulative histogram of limiting flux density for GBNCC. The mean and median limiting flux densities in the histogram 0.74\,mJy and 0.62\,mJy, and the values range from 0.42\,mJy to 47.\,mJy. All flux density values are given in mJy.}
    \label{fig:flux_lim_hist}
\end{figure}

\begin{figure*}
\plotone{limiting_flux_skymap.pdf}
    \caption{Sky map of GBNCC beams, colored by limiting flux density. The map is plotted in Galactic coordinates on a Mollweide projection, and the flux density is given in mJy.}
    \label{fig:flux_lim_map}
\end{figure*}

\subsection{Nulling/Mode-Changing Candidates}
The large set of data analyzed in this study as well as the ``by-eye" verification of all detections allowed for easy identification of potential nulling/mode-changing candidates in the results. This way, we are sensitive to nulling timescales between that of the pulsar spin period and the observation time (120 seconds). These cases were first identified by the appearance of missing pulses in the time-phase plots from processing using the \texttt{PRESTO} package. When a pulsar was noted as a candidate, we followed up using the \texttt{dspsr}\footnote{\url{http://dspsr.sourceforge.net/index.shtml}} package. We folded the time series data in 10 second integrations, zapped remaining RFI by hand, and integrated across frequency using the \texttt{pav} and \texttt{pam} commands within \texttt{PSRCHIVE}\footnote{\url{http://psrchive.sourceforge.net/index.shtml}}. When it was possible to discern on- and off-pulse regions by eye (i.e. significant changes in intensity for some rotations), the candidates were considered likely to be nulling. Some pulsars exhibited behavior similar to mode-changing, where multiple components of the averaged profile were found to be on during different portions of the observation. These pulsars were not treated differently than other nulling candidates $-$ we folded for single pulses to determine the likelihood that different components were visible. All of these sources will be followed up in later works regarding these data. In total, 223 pulsars were found to exhibit intensity variations similar to nulling or mode-changing during their observations, 62 of which have not previously been found to to do so. These candidates' names are marked in Table \ref{tbl:det} with an asterisk.

\subsection{The Galactic Pulsar Population}
Given its overall sky coverage and the large number of pulsar detections reported here (670), the GBNCC survey will play an important role in future understanding of the Galactic pulsar population. To date, the GBNCC survey has detected 571 non-recycled (long-period) pulsars in the Galactic field and 70 Galactic MSPs, which have undergone recycling and have spin periods, $P<30$\,ms. Remaining detections are either associated with globular clusters (3) or are recycled pulsars with spin periods, $P>30$\,ms (26), and have been intentionally ignored for the following analysis, since our current models do not adequately describe the features of this sub-population.

To estimate expected numbers of non-recycled/millisecond pulsar detections in the GBNCC survey, Galactic populations were simulated using {\sc PsrPopPy2}\footnote{https://github.com/devanshkv/PsrPopPy2}, a more recent and currently-maintained version of {\sc PsrPopPy} \citep[][and citations within]{blr+14}. Pulsar populations were generated using {\sc PsrPopPy2}'s {\tt populate} function, which simulates pulsars by drawing parameters from predefined distributions until some condition is met. Due to its large sample size, population estimates from the Parkes Multibeam Pulsar Survey (PMPS) provide the best-known sample parameters. For this reason, these results were used to set a limit on the number of pulsars simulated by {\tt populate}. For the non-recycled pulsar population, pulsars were generated until a synthetic PMPS ``detected" 1038 sources; for MSPs, the desired population size was set to 30,000 sources. Specific parameters defining pulsars' Galactic radial distribution, as well as scale height, spin period, luminosity, and duty cycle can be found in \citet{slm+14}. However, an updated model for the MSP $P$-distribution \citep{lem+15} was implemented in simulations here.

Synthetic surveys were conducted with 100 realizations each of the Galactic non-recycled/millisecond pulsar populations using {\tt survey} and a GBNCC model file, including survey parameters identical to those presented in \S\ref{sec:methods} and lists of completed/remaining GBNCC pointing positions. In the first round of simulations, we fixed the S/N limit for detections to ${\rm S/N}_{\rm cut}=3.8$ (as determined in \S\ref{subsec:efficiency}). This simulation predicted 1442/126 simulated detections for non-recycled/millisecond pulsar populations, respectively (on average; compared to 571/70 actual detections). We then fixed the number of simulated non-recycled/millisecond pulsar detections to their actual values (571/70) and found nominal S/N thresholds for each sub-population, ${\rm S/N}_{\rm cut}=15.3/9.1$. The discrepancies between simulated and actual yields suggest uncertainties in population parameters informed primarily by the PMPS survey, which targeted the Galactic plane and was conducted at 1.4\,GHz. Population parameters determined by these previous surveys produce over-estimates for GBNCC pulsar yields. As an all-sky, low-frequency search, the GBNCC survey (when complete) will be a valuable counterpoint to further refine non-recycled/millisecond pulsar population parameters. As we will show below, positional and rotational parameters of the simulated populations do not match the detected population when these thresholds are set.

\begin{deluxetable}{lcccc}
\tablecaption{K-S test statistics and $p$-values resulting from comparisons between actual/simulated parameter distributions for non-recycled/millisecond pulsars. In cases where the $p$-value is $<1\%$, the null hypothesis (that the two distributions are the same) is rejected.\label{tbl:ks}}
\tablehead{\colhead{Parameter} & \multicolumn{2}{c}{Normal\tablenotemark{a}}& \multicolumn{2}{c}{MSP\tablenotemark{b}}\\
\cline{2-3} \cline{4-5}
& \colhead{K-S} & \colhead{$p(\%)$} & \colhead{K-S} & \colhead{$p(\%)$}
}
\startdata
Spin Period ($P$) & 0.20 & $\ll1$ & 0.14 & \phn$10$ \\
Dispersion Measure (DM) & 0.21 & $\ll1$ & 0.26 & $<1$ \\
Flux Density ($S_{350}$) & 0.13 & $\ll1$ & 0.21 & $<1$ \\
Galactic Latitude ($b$) & 0.04 & \phn\phn$41$ & 0.17 & \phn\phn$3$ \\
\enddata
\tablenotetext{a}{For simulated non-recycled pulsars, ${\rm S/N}_{\rm cut}=15.3$.}
\tablenotetext{b}{For simulated MSP population, ${\rm S/N}_{\rm cut}=9.1$.}
\end{deluxetable}
To test the validity of underlying non-recycled/millisecond pulsar populations, we compared cumulative distribution functions (CDFs) of simulated ({\tt sim}) pulsar parameters ($P$, DM, $S_{350}$, and $b$) with those of the actual ({\tt act}) detections using a {\tt scipy} implementation of the 2-sample Kolmogorov–Smirnov (K-S) test. For each parameter, the K-S test statistic and $p$-value were computed over a range of ${\rm S/N}_{\rm cut}$. When $p<1\%$, the null hypothesis (that {\tt act}/{\tt sim} parameters are drawn from the same underlying distribution) is rejected. Figures \ref{fig:norm_cdfs} and \ref{fig:msp_cdfs} illustrate these comparisons for non-recycled and millisecond pulsar population parameter distributions, and Table \ref{tbl:ks} summarizes K-S test results when the nominal ${\rm S/N}_{\rm cut}$ values for non-recycled/millisecond pulsar sub-populations (15.3/9.1) are implemented, though we measured these p-values for a range of imposed ${\rm S/N}_{\rm cut}$ values.

\begin{figure*}
\centering
\includegraphics[width=0.75\textwidth]{NORM-Hists+CDFs.pdf}
    \caption{Normalized histograms showing comparisons between (a) spin period, $P$, (b) Galactic latitude, $b$, (c) flux density, $S_{350}$, and (d) dispersion measure, DM, distributions for simulated non-recycled pulsars (blue) and actual detections (orange). The rightmost panel in each row compares actual/simulated CDFs for each parameter. K-S tests comparing these CDFs (see Table \ref{tbl:ks} for details) show disagreement between {\tt act}/{\tt sim} $P$, $S_{350}$, and DM distributions, but $p=41\%$ for $b$ distributions.}
    \label{fig:norm_cdfs}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=0.75\textwidth]{MSP-Hists+CDFs.pdf}
    \caption{Normalized histograms showing comparisons between (a) spin period, $P$, (b) Galactic latitude, $b$, (c) flux density, $S_{350}$, and (d) dispersion measure, DM, distributions for simulated millisecond pulsars (blue) and actual detections (orange). The rightmost panel in each row compares actual/simulated CDFs for each parameter. K-S tests comparing these CDFs (see Table \ref{tbl:ks} for details) show disagreement between {\tt act}/{\tt sim} $S_{350}$ and DM distributions, but distributions for $b$ and $P$ have $p=3\%$ and $10\%$, respectively.}
    \label{fig:msp_cdfs}
\end{figure*}

Comparing {\tt act}/{\tt sim} parameters for the non-recycled pulsar population, we find broad agreement between $b$ distributions, regardless of ${\rm S/N}_{\rm cut}$. Results for other non-recycled pulsar parameters in Table \ref{tbl:ks} show significant inconsistencies between {\tt act}/{\tt sim} samples. DM distributions are clearly different for ${\rm S/N}_{\rm cut}>4$, likely due to an over-abundance of low-DM simulated detections. For ${\rm S/N}_{\rm cut}=15.3$, we find twice as many {\tt sim} detections with ${\rm DM}<35$\,pc\,cm$^{-3}$. Presumably due to the prevalence of nearby {\tt sim} sources, this sample also has a larger fraction of high-flux density sources, so $S_{350}$ distributions are statistically different for ${\rm S/N}_{\rm cut}=15.3$. However, there is a small window ($10.25<{\rm S/N}_{\rm cut}<12.25$) where {\tt act}/{\tt sim} $S_{350}$ distributions become statistically similar, with $p>1\%$. The null hypothesis is rejected for $P$ due to {\tt act}/{\tt sim} log-normal distributions having different mean values: $\langle\log P_{\rm act}\rangle = 2.88$ versus $\langle\log P_{\rm sim}\rangle = 2.72$ (see Figure \ref{fig:norm_cdfs}). This discrepancy persists, regardless of chosen ${\rm S/N}_{\rm cut}$.

Because the simulated versions of the non-recycled pulsar population were primarily informed by PMPS \citep[e.g.][]{lfl+06}, which was conducted at 1.4\,GHz and exclusively covered regions of sky near the Galactic plane ($|b| < 5\,\degr$), we expect there to be bias toward highly dispersed pulsars near the plane. Due to more uniform sky coverage and \--- near the Galactic plane \--- higher sky temperatures and more significant scattering at 350\,MHz, the majority of GBNCC detections ($67\%$) are away from the plane ($|b| > 5\,\degr$). Young pulsars are typically born in the plane and tend to be found nearby, therefore GBNCC's reduced sensitivity to low-latitude sources means that relatively few detections are young pulsars. The $P$-$\dot{P}$ diagram in Figure \ref{fig:ppdot} nicely illustrates this shortage of pulsars detected with characteristic ages, $\tau \leq 1$\,Myr. By imposing an age cutoff on non-recycled pulsars in the ATNF catalog, $\tau > 1$\,Myr, the resulting simulated spin period distribution is statistically similar to that of GBNCC detections (K-S $p>1\%$). This selection effect accounts for the apparent differences between {\tt act}/{\tt sim} $P$-distributions, but can not explain discrepancies in $S_{350}$ and DM distributions for non-recycled pulsars.

K-S tests comparing {\tt act}/{\tt sim} parameter distributions for the MSP population show better agreement (see Table \ref{tbl:ks} and Figure \ref{fig:msp_cdfs}). For MSPs, selection effects based on Galactic latitude and spin period do not come into play since MSPs are more isotropically distributed and model parameters for this sub-population are based on results from multiple Parkes Telescope surveys (see \cite{lem+15}, and references therein). For these reasons, the simulated population's spin periods are statistically similar to the sample detected by GBNCC. This conclusion does not change, regardless of the chosen ${\rm S/N}_{\rm cut}$ value. For $b$, the null hypothesis is still not rejected by our criteria ($p<1\%$). Based on the $b$ histograms themselves, there appears to be an absence of detections in the {\tt act} sample in/near the Galactic plane, which is not the case for {\tt sim} sources. The null hypothesis is rejected for $S_{350}$ due to the over-abundance of high-flux-density sources in the {\tt sim} sample compared to those present in the {\tt act} sample. Median flux densities for {\tt act}/{\tt sim} detections are 4.9/6.8\,mJy respectively. Comparing {\tt act}/{\tt sim} DM distributions, the {\tt sim} sample consists of a higher fraction of high-DM MSPs and 12\% of simulated detections have DMs in excess of the {\tt act} maximum value, $104.5$\,pc\,cm$^{-3}$. This is likely related to the bias toward high-flux-density detections noted in $S_{350}$ for {\tt sim} MSPs mentioned above.

Based on discrepancies between predicted yields from simulations and actual numbers of detections by the GBNCC survey, it appears that model parameters need to be further refined in order to generate more realistic Galactic pulsar populations in the future. For now, we proceed with nominal ${\rm S/N}_{\rm cut}$ values in order to estimate the GBNCC survey's future yields. In the remaining $\approx$21,000 pointings, we expect an additional 160/16 non-recycled/millisecond pulsar detections, or $\approx60/5$ discoveries, accounting for detectable known pulsars in regions of sky remaining \citep{mhth05}.

\begin{figure}
\epsscale{1.2}
\plotone{ppdot.pdf}
    \caption{Period vs.\ period derivative for pulsars in GBNCC survey area. Shown in grey are pulsars that were not detected, and blue plus symbols show detections.}
    \label{fig:ppdot}
\end{figure}

\section{Conclusions}
\label{sec:conc}
We have provided all detections of currently known pulsars that exist within the area of the 350\,MHz GBNCC pulsar survey and performed some preliminary analysis of the resulting data set. Specifically, we have provided new flux density and pulse width measurements as well as pulse profiles for the 670 detections. When possible, we used our flux density measurements with previous measurements at different frequencies to refine spectral index. We also made a measurement of the spin period-pulse width relation, observing a powerlaw correlation of the form $W_{10} \propto P^{-0.27}$. The low frequency of the survey provides increased sensitivity to dispersion, allowing for more precise measurements of DM for many pulsars that have only been measured in high frequency surveys. Using all of this information, we have made quantitative measurements of the survey's efficacy and the RFI environment at 350\,MHz, with a minimum detectable S/N of $\sim$3.8 and a mean limiting flux density of 0.74\,mJy. These measurements have allowed us to make realistic predictions about the survey's yield when complete based on the detectability of known pulsars in the dataset, and we expect to detect on the order of 160 non-recycled pulsars and 15 MSPs. The simulations from which these expectations come uncovered discrepancies in DM, spin period, and spatial distribution in the Galaxy for the simulated populations which will be addressed in a future study. Combing through the data following processing has brought many interesting characteristics of pulsars in the survey to light, including 223 pulsars exhibiting evidence of variable intensities suggestive of nulling/mode-changing and 4 showing evidence for broken power-law spectral energy distributions. These kinds of qualitative observations pave the way for follow-up quantitative analyses of these data and the remaining beams that will be observed in the next few years.

\section*{Acknowledgements}
We thank our anonymous referee for their suggestions and guidance. This work was supported by the NANOGrav Physics Frontiers Center, which is supported by the National Science Foundation award 1430284. The Green Bank Observatory is a facility of the National Science Foundation (NSF) operated under cooperative agreement by Associated Universities, Inc.
R.S.\ acknowledges support through the Australian Research Council grant FL150100148.
WF acknowledges the WVU STEM Mountains of Excellence Graduate Fellowship. MM and MS acknowledge the National Science Foundation OIA award number 1458952. JvL acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 617199 (`ALERT'), and from Vici research programme `ARGO' with project number 639.043.815, financed by the Netherlands Organisation for Scientific Research (NWO). VMK acknowledges the NSERC Discovery Grant, the Herzberg Award, FRQNT and CRAQ, Canada Research Chairs, CIFAR and the Webster Foundation Fellowship, and the Trottier Chair in Astrophysics and Cosmology. Computations were made on the supercomputer Guillimin at McGill University\footnote{\url{www.hpc.mcgill.ca}}, managed by Calcul Quebec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), NanoQuebec, RMGA and the Fonds de recherche du Quebec - Nature et technologies (FRQ-NT). The CyberSKA project was funded by a CANARIE NEP-2 grant. PC acknowledges the FRQNT Doctoral Research Award. SMR is a CIFAR Senior Fellow. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canada Foundation for Innovation.

\software{\texttt{Astropy} \citep{astropy:2018}, \texttt{PRESTO} \citep{smr+01},
\texttt{PsrPopPy2} \citep{blr+14},
\texttt{SciPy} \citep{scipy},
\texttt{NumPy} \citep{numpy},
\texttt{dspsr} \citep{sb11},
\texttt{PSRCHIVE} \citep{hsm04},
\texttt{TEMPO} (\url{http://tempo.sourceforge.net/})}

\facilities{Robert C. Byrd Green Bank Telescope (GBT)}

\bibliographystyle{apj}


\section{Background and Motivation}

In the field of ElectroMagnetics (EM), Boundary Value Problems (BVPs) are the problems for which the EM field in a given region of space is determined from a knowledge of the field over the boundary of the region \citep{BVP}. Traditional numerical methods such as Finite Difference Time Domain (FDTD) and Finite Element Methods (FEM) are typically employed to solve BVPs. However, these methods involve discretization of the domain to reduce it to higher-order system of linear algebraic equations and solving for them. As such, these methods are not local i.e. they do not give the value of the solution directly at an arbitrary point, where the field needs to be determined, but its value should be extracted from the complete field solution and hence are not amenable to parallel processing.

    Neural Network (NN) based numerical method provides an alternate approach to solving BVPs \citep{Lagaris},\citep{Mcfall}. The principal advantages of the NN based numerical method are the discrete data points where field is computed, can be unstructured and therefore, the issues of meshing (uniform/non-uniform) are not a factor; the solutions are in a differentiable, closed analytic form which avoids the need to interpolate between data points where solutions are obtained using other methods. Also, recent advances in the field of Machine Learning (ML) and Artificial Intelligence (AI) has jump started the design and implementation of computer architectures that are optimized to implement training and inference tasks more efficiently. Since, this NN based numerical method is inherently parallel and hence can be efficiently implemented on parallel architectures,
 it stands to gain from advances in AI based computer architectures.
	In this study, NN based field computation is formulated and presented for BVPs with Dirichlet boundary condition and BVPs with mixed boundary conditions on a uniform rectangular grid. Multiple non-uniform training grid distributions are also explored in this study  to showcase that this method is not limited by the domain discretization like traditional methods.

\section{Problem Statement}
NN based method is developed, implemented and investigated for solving Laplace equation with Dirichlet boundary condition and Poisson's equation with mixed boundary condition. Numerical investigations are carried out to
understand efficacy
of this method and to provide qualitative understanding
of various parameters involved.

\section{Formulation}
BVPs, where the EM field $\psi(x)$ is given by the linear second order partial differential equation (PDE), are considered in this study \cite{Lagaris}
\begin{eqnarray}
G(x,\psi(x),\nabla\psi(x),\nabla^2\psi(x))&=& 0,\quad \forall x\in D
\label{l2pde}
\end{eqnarray}
subject to boundary condition (B.Cs). Here $x =(x_1,x_2,\cdots,x_n) \in\mathbb{R}^n$ and $D\subset{\mathbb{R}^n}$. To solve this BVP using NN based method, a trial form of the solution is constructed, which is written as sum of two parts: the first part satisfies the boundary conditions and contains no adjustable parameters and the second part which involves a feed forward neural network with adjustable parameters is constructed so as not to contribute to boundary conditions.
\begin{eqnarray}
\psi_t(x,W,b)&=&\hat{\psi}(x)+F(x)N(x,W,b)
\end{eqnarray}
$N(x,W,b)$ is a feed forward NN with weights $W$ and biases $b$. $F(x)$ is chosen such that second part does not contribute to boundary conditions.

In such a case, the task of field computation reduces to learning the NN parameters, $W$ and $b$, which is done by first transforming the equation (\ref{l2pde}) to a discretized version and applying it at discretized domain, $\hat{D}=\{x^{(i)}\in D; i=1,\cdots,m\}$.
\begin{eqnarray}
G(x^{(i)},\psi(x^{(i)}),\nabla\psi(x^{(i)}),\nabla^2\psi(x^{(i)}))&=& 0,\quad \forall x^{(i)}\in D \nonumber
\end{eqnarray}
and then training the NN, where the PDE error or cost corresponding to $x^{(i)}
 has to become zero.
\small
\begin{eqnarray}\small
W,b=\argminA_{W,b}{} G(x^{(i)},\psi_t(x^{(i)},W,b),\nabla\psi_t(x^{(i)},W,b),\nabla^2\psi_t(x^{(i)},W,b))^2 \nonumber
\end{eqnarray}
\normalsize
Note that the computation of this cost function involves not only the network output, but also the derivatives of the output with respect to any of its inputs.

\subsection{Neural Network Architecture}
The neural network that will be implemented for this solution is a $3$-layer network with an input layer ($n+1$ nodes), a hidden layer ($H$ nodes) and an output layer ($1$ node) as shown in Fig.\ref{fig:NN}. Sigmoid function $(\sigma(.))$ will be used as activation function for the hidden layer.
\begin{eqnarray}
h&=&W^{[1]}x+b^{[1]} \nonumber \\
N&=&W^{[2]}\sigma(h) \nonumber
\end{eqnarray}
where $W^{[1]} \in \mathbb{R}^{H\times n}, W^{[2]} \in \mathbb{R}^{1\times H}$ and $h,b^{[1]} \in \mathbb{R}^{H\times 1}$. $W^{[1]}_j$ is used to denote the $j^{th}$ column of $W^{[1]}$.

\begin{figure}[!hbt]

		\centering

		\includegraphics[width=.7\linewidth]{NN_Architecture.png}

		\caption{Neural Network Architecture with $n$ input nodes, $H$ hidden nodes and $1$ output node.}

		\label{fig:NN}

	\end{figure}

\subsection{Network Output Derivative Computation}
For this neural network, it can be shown that
\begin{eqnarray}
\frac{\partial^{\lambda_1}}{\partial x_1^{\lambda_1}}\frac{\partial^{\lambda_2}}{\partial x_2^{\lambda_2}}\cdots \frac{\partial^{\lambda_n} }{\partial x_n^{\lambda_n}} N&=&\sum_{i=1}^H W^{[2]}_i\left(\prod_{j=1}^n(W^{[1]}_{ij})^{\lambda_j}\right)\sigma^{(\lambda)}(h_i) \nonumber
\end{eqnarray}
where $\lambda=\sum_{i=1}^n\lambda_i$ and $\sigma^{(\lambda)}(h_i)$ denotes the $\lambda$th order derivative of the sigmoid.

Therefore, first order derivate of $N$, with respect to any input parameter $x_j$ is given by
\begin{eqnarray}
\frac{\partial N }{\partial x_j} &=& \sum_{i=1}^H W^{[2]}_i(W^{[1]}_{ij})\sigma^{(1)}(h_i) = \left(W^{[2]}\circ W^{[1]}_j\right)) \sigma^{(1)}(h_i) \nonumber
\end{eqnarray}
which can be interpreted as the output of the feedforward neural network of same schematic, where the activation function for the hidden layer is given by the first order derivative of sigmoid, instead of sigmoid and the $W^{[2]}$ replaced by $W^{[2]}\circ W^{[1]}_j$.

Similarly, the second order derivative of $N$ with respect to $x_j$ can be interpreted as the output of a feedforward neural network with the same architecture, where the activation function for the hidden layer is given by the second order derivative of sigmoid and the $W^{[2]}$ replaced by $W^{[2]}\circ (W^{[1]}_j \circ W^{[1]}_j)$.
\begin{eqnarray}
\frac{\partial^2 N }{\partial x_j^2} &=& \sum_{i=1}^H W^{[2]}_i(W^{[1]}_{ij})^2\sigma^{(2)}(h_i) \nonumber \\&=& \left(W^{[2]}\circ W^{[1]}_j\circ W^{[1]}_j\right)) \sigma^{(2)}(h_i) \nonumber
\end{eqnarray}
For the cost function, which includes network output as well as derivatives of network output, these interpretations become extremely useful during training, when cost function is being optimized.

\section{Implementation}
The formulated method is implemented to find solution of two BVPs - Laplace equation with Dirichlet boundary condition and Poisson's equation with mixed boundary condition. In both examples, the domain is chosen to be a square $D=[0,1]\times [0,1]$.
Neural Network has been implemented using Tensor Flow framework and optimized using Stochastic Gradient Descent (SGD) with annealing learning rate and regularization.

\subsection{Laplace Equation with Dirichlet boundary condition}
Electrostatic potential distribution inside a rectangular region where the potential on the boundary is specified, is given by the Laplace equation with Dirichlet boundary condition in a $2D$ rectangular region \cite{janaswamy}. NN based method is used to compute the solution and compared with the analytical solution \cite{stanford}.
\begin{eqnarray}
\nabla^{2}\psi(x)=0, \quad \forall x \in D
\end{eqnarray}
The boundary conditions are
\begin{eqnarray}
\psi(x)&=&0, \quad \forall x \in \{(x_1,x_2)\} \in \begin{cases}
 \partial D| x_1=0\\
 \partial D| x_1=1\\
 \partial D| x_2=0
\end{cases} \nonumber \\
\psi(x)&=&sin(\pi x_1),\quad \forall x \in \{(x_1,x_2)\in\partial D|x_2=1 \} \nonumber
\end{eqnarray}

The analytical solution is
\begin{eqnarray}
\psi_a(x)&=&\frac{1}{e^\pi-e^{-\pi}}sin(\pi x_1)(e^{\pi x_2}-e^{-\pi x_2})
\end{eqnarray}
The trial solution constructed for NN based method is
\begin{eqnarray}
\psi_t(x)&=&x_2sin(\pi x_1)+x_1(1-x_1)x_2(1-x_2)N(x,W,b) \nonumber
\end{eqnarray}
In this case, the cost function is given by
\begin{eqnarray}
&&{-\pi^2 x_2 sin(\pi x_1)+} \nonumber \\
&&{x_2(1-x_2)\left(x_1(1-x_1)\frac{\partial^2 N}{\partial x_1^2} +(2-4x_1)\frac{\partial N}{\partial x_1}-2N\right)+} \nonumber\\
&&{x_1(1-x_1)\left(x_2(1-x_2)\frac{\partial^2 N}{\partial x_2^2} +(2-4x_2)\frac{\partial N}{\partial x_2}-2N\right)}\nonumber
\end{eqnarray}
The analytical solution (Fig. \ref{fig:laplace_a}) and the NN based solution (Fig. \ref{fig:laplace_nn}) computed by minimizing the cost function when $K=16$ and $H=15$ are shown below
The absolute value of the delta between the two solutions, $|\psi_a(x)-\psi_t(x)|$ is also plotted in Fig. \ref{fig:delta_1}.

\begin{figure}[h]
\centering
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{laplace_analytical.png}
  \captionof{figure}{Analytical solution}
  \label{fig:laplace_a}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1.1\linewidth]{NN_H_15_K_16.png}
  \captionof{figure}{NN based solution}
  \label{fig:laplace_nn}
\end{minipage}
\end{figure}

\begin{figure}[!hbt]

		\begin{center}

		\includegraphics[width=.56\linewidth]{Delta_H_15_K_16.png}

		\caption{$E_{abs}$, Laplace equation}

		\label{fig:delta_1}
		\end{center}
	\end{figure}
\subsection{Poisson's equation with mixed boundary condition}
Electrostatic potential in the presence of charge distribution inside a rectangular region where the potential is specified on a section of the boundary and the gradient of the potential is specified on the rest is considered here.
\begin{eqnarray}
\nabla^{2}\psi(x)=(2-\pi^2x_2^2)sin(\pi x_1), \quad \forall x \in D
\end{eqnarray}
The mixed boundary conditions are
\begin{eqnarray}
\psi(x)&=&0, \quad \forall x \in \{(x_1,x_2)\in
\begin{cases}
 \partial D| x_1=0\\
 \partial D| x_1=1\\
 \partial D| x_2=0
\end{cases} \nonumber \\
\frac{\partial \psi(x)}{\partial x_2}&=&2sin(\pi x_1),\quad \forall x \in \{(x_1,x_2)\in\partial D_n|x_2=1 \} \nonumber
\end{eqnarray}
The analytical solution is
\begin{eqnarray}
\psi_a(x)&=&x_2^2 sin(\pi x_1)
\end{eqnarray}
The trial solution constructed for NN based method is
\begin{eqnarray}
\psi_t(x)&=&2x_2sin(\pi x_1)\nonumber \\
&&+x_1(1-x_1)x_2\left[N(x_1,x_2,W,b)\right] \nonumber \\
&&-x_1(1-x_1)x_2\left[ N(x_1,1,W,b)+\frac{\partial N(x_1,1,W,b)}{\partial x_2}\right] \nonumber
\end{eqnarray}
In this case, the cost function is given by
\begin{eqnarray}
&&-(2-\pi^2 x_2)sin(\pi x_1)-2\pi^2 x_2 sin(\pi x_1)+ \nonumber \\
&&x_2x_1(1-x_1)\left[\frac{\partial^2 N(x_1,x_2)}{\partial x_1^2}-\frac{\partial^2 N(x_1,1)}{\partial x_1^2} -\frac{\partial^3 N(x_1,1)}{\partial x_1^2 \partial x_2}\right]+\nonumber\\
&&2x_2(1-2x_1)\left[\frac{\partial N(x_1,x_2)}{\partial x_1}-\frac{\partial N(x_1,1)}{\partial x_1} -\frac{\partial^2 N(x_1,1)}{\partial x_1 \partial x_2}\right]-\nonumber\\
&&2x_2\left[N(x_1,x_2)-N(x_1,1)-\frac{\partial N(x_1,1)}{\partial x_2}\right]+\nonumber \\
&&x_2x_1(1-x_1)\left[\frac{\partial^2 N(x_1,x_2)}{\partial x_2^2}\right]
+2x_1(1-x_1)\left[\frac{\partial N(x_1,x_2)}{\partial x_1}\right
\end{eqnarray}
The analytical solution and the NN based solution computed by minimizing the cost function when $K=16$ and $H=15$ are shown below
The absolute value of the delta between the two solutions, $|\psi_a(x)-\psi_t(x)|$, with a maximum value of $0.0056$, plotted in Fig. \ref{fig:delta} showcases good agreement between the two solutions.
\begin{figure}[!hbt]

		\begin{center}

		\includegraphics[width=.56\linewidth]{Poisson_Delta_H_15_K_16.png}

		\caption{$E_{abs}$, Poisson's equation}

		\label{fig:delta}
		\end{center}
\end{figure}
\section{Numerical Investigation on Error Properties}
Let $\hat{D}_{train}=\{x^{(i)}_{train}\in D; i=1,\cdots,m_{train}\}$ and $\hat{D}_{test}=\{x^{(i)}_{test}\in D; i=1,\cdots,m_{test}\}$ represent the set of training and testing dataset respectively. For the case of uniform discretization of $2D$ domain with a training grid resolution, $K$, $m_{train}=K^2$ and for the non-uniform discretization as well, the term training grid resolution, $K$, is meant to indicate that $m_{train}=K^2$
The testing grid resolution for all the experiments is kept fixed at $21$, therefore the testing dataset is of the size $m_{test}=441$, where the $2D$ domain is discretized uniformly.

The performance metrics that will be used in this study are the absolute value of the delta between analytical solution and NN based solution
\begin{eqnarray}
E_{abs}=|(\psi_a(x)-\psi_t(x))|
\end{eqnarray}
and relative error norm as given by
\begin{eqnarray}
E_{norm}=\frac{\sqrt{\sum_x(\psi_a(x)-\psi_t(x))^2}}{\sqrt{\sum_x(\psi_a(x))^2}}
\end{eqnarray}

$E_{abs}$ gives us information about the spatial distribution of the divergence between analytical and NN based solution, while the $E_{norm}$ paints an aggregate picture over the entire domain.

\begin{figure}[!hbt]

		\begin{center}

		\includegraphics[width=.56\linewidth]{Convergence_Laplace_500Epoch.png}

		\caption{Progression of PDE error and $E_{norm}$ during training}

		\label{fig:convergence}
		\end{center}
	\end{figure}

There is a subtle, but important point that needs to be understood when solving BVPs using NN, $i.e$ training is driven by the cost function (or PDE error), while performance is gauged by the error in the computed field \cite{Mcfall}.

During the optimization of NN, the PDE is more closely satisfied at the training points for each successful training epoch. But the reduction of PDE error at training points has no obvious guarantee that the NN based solution is converging to the analytical solution, even at the training points. In addition to that, NN with poor generalization because of over fitting may not satisfy the PDE in the regions of the domain away from the training points, and hence NN based solution diverges from analytical solution at such points. Convergence theorem outlines that the difference between the analytical and NN based solutions will decrease everywhere in the domain, when the PDE error is also reduced everywhere in the domain \cite{Mcfall2}. In practice, however, PDE error cannot be guaranteed to decrease everywhere in the domain during training. However, convergence plots such as one shown in Fig. \ref{fig:convergence} for the Laplace equation when $H=15,m_{train}=256,m_{test}=441$, should be used to monitor the progression of PDE error and $E_{norm}$ during training to make sure both training and test set error decrease along with PDE error during training. Training grid resolution, Number of hidden nodes, Training dataset distribution and size can be tweaked accordingly, if needed to ensure $E_{norm}$ as well as PDE error decreases during training.

\begin{figure}[!h]
\centering
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=\linewidth]{Error_K.png}
  \captionof{figure}{Error Dependence on Training Grid \\ Resolution, K}
  \label{fig:Error_K}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=\linewidth]{Error_H.png}
  \captionof{figure}{Error Dependence on  Number of \\ Hidden Nodes, H}
  \label{fig:Error_H}
\end{minipage}
\end{figure}

\begin{figure*} [!hbt]
\centering
\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{Grid_1.png}
\caption{Grid 1: Uniform}
\label{fig:grid1}
\end{subfigure}%
\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{Grid_2.png}
\caption{Grid 2: Non-Uniform}
\label{fig:grid2}
\end{subfigure}

\bigskip

\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{Grid_3.png}
\caption{Grid 3: Non-Uniform}
\label{fig:grid3}
\end{subfigure}%
\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{Grid_4.png}
\caption{Grid 4: Random}
\label{fig:grid4}
\end{subfigure}

\caption{Four different spatial distribution of the training data points investigated}
\label{fig:Grid Experiment}
\end{figure*}

 \begin{figure}[!hbt]

		\begin{center}

		\includegraphics[width=.56\linewidth]{Error_Grid.png}

		\caption{Error dependence on training dataset distribution}

		\label{fig:Error_Grid}
		\end{center}
	\end{figure}
\subsection{Error Dependence on Training Grid Resolution, K}
Neural network was trained for various combinations of the number of hidden nodes, $H$ and training grid resolution, $K$ for the BVP with Laplace equation with Dirichlet boundary condition. For a fixed $H=15$, the variation of test set error when training grid resolution $K={8,16,24,30 \text{ and } 40}$ is shown in Fig.\ref{fig:Error_K}. As $K$ increases from $8$ to $16$, the $E_{abs}$ becomes relatively narrower and $E_{norm}$ decreases from $0.015$ to $0.002$, but plateaus beyond that resolution. Though it is intuitively expected that the test set error improves with increase in the training set size, as the learning model generalizes better, in this particular case, it could also be inferred in the context of the convergence theorem as well. As the PDE error is minimized at more training points, the convergence between analytical and NN based solution gets better.
\subsection{Error Dependence on Number of Hidden Nodes, H}
To understand the dependence of the test set error on number of hidden nodes, $H$, NN based solution is computed for $H=5,10,15,25,35 \text{ and } 45$, for a fixed $K=16$. Fig.\ref{fig:Error_H} shows how the error improves significantly as $H$ is increased from $5$ to $15$, but stays about the same beyond that. As the computational requirements grow with increase in the number of hidden nodes and/or training grid resolution, and only marginal gains beyond $H=15,K=16$; Hence, this choice of number of hidden nodes and training grid resolution is deemed appropriate for the problem at hand.
\subsection{Error Dependence on Non-Uniform Training Grid}
As mentioned previously, one of the main advantages of NN based method is the training and testing dataset can be totally unstructured and hence they do not suffer from the issue of meshing like other numerical methods. To showcase this, four different grid structures, shown in the Fig.\ref{fig:Grid Experiment} where each dot shows a training point in the training dataset, $m_{train}=256$, are investigated in this study. Fig \ref{fig:grid1} shows ``Grid-1'', the uniform discretization that will be used as baseline in this experiment; Two structured, but different spatial distributions are considered in ``Grid-2" shown Fig \ref{fig:grid2}, which has training points densely distributed closer to the boundary and coarsely away from it and the spatial distribution showing opposite trend in ``Grid-3", shown in Fig \ref{fig:grid3}. Fig \ref{fig:grid4} shows ``Grid-4", the training dataset generated by uniform random distribution in $[0,1]\times[0,1]$.

Fig.\ref{fig:Error_Grid} shows the computed test set error for the four grid structures. ``Grid-4" has the best error performance, followed by ``Grid-1", while ``Grid-2" and "Grid-3" show similar performance. This may be due to fact that the training dataset distribution ``Grid-4" closely resembles the test dataset distribution, compared to the other distributions, as the testing grid is a uniform discretization with a grid resolution of $21$. This demonstrates the NN based method's merits when dealing with unstructured data and its potential to be applied to BVPs with curved or irregular boundaries.

\section{Conclusion and Future work}
In this study, a NN based numerical method has been formulated, implemented and validated for two examples with Dirichlet and mixed boundary conditions. Numerical experiments were carried out to assess the dependence of the error performance on training grid resolution and number of hidden nodes. To showcase the merit of this method, numerical results for structured and unstructured training dataset are presented as well.

Future efforts could explore application of this method to BVPs with curved or irregular boundaries. Some of the commonly found BVPs in EM are given by PDEs with complex coefficients, which could be the focus of future work as well.

\small