We present a systematic search for optical counterparts to 13 gravitational wave (GW) triggers involving at least one neutron star during LIGO/Virgo's third observing run. We searched binary neutron star (BNS) and neutron star black hole (NSBH) merger localizations with the Zwicky Transient Facility (ZTF) and undertook follow-up with the Global Relay of Observatories Watching Transients Happen (GROWTH) collaboration. The GW triggers had a median localization of 4480 deg^2, median distance of 267 Mpc and false alarm rates ranging from 1.5 to 1e-25 per yr. The ZTF coverage had a median enclosed probability of 39%, median depth of 20.8mag, and median response time of 1.5 hr. The O3 follow-up by the GROWTH team comprised 340 UVOIR photometric points, 64 OIR spectra, and 3 radio. We find no promising kilonova (radioactivity-powered counterpart) and we convert the upper limits to constrain the underlying kilonova luminosity function. Assuming that all kilonovae are at least as luminous as GW170817 at discovery (-16.1mag), we calculate our joint probability of detecting zero kilonovae is only 4.2%. If we assume that all kilonovae are brighter than -16.6mag (extrapolated peak magnitude of GW170817) and fade at 1 mag/day (similar to GW170817), the joint probability of zero detections is 7%. If we separate the NSBH and BNS populations, the joint probability of zero detections, assuming all kilonovae are brighter than -16.6mag, is 9.7% for NSBH and 7.9% for BNS mergers. Moreover, <57% (<89%) of putative kilonovae could be brighter than -16.6mag assuming flat (fading) evolution, at 90% confidence. If we further account for the online terrestrial probability for each GW trigger, we find that <68% of putative kilonovae could be brighter than -16.6mag. Comparing to model grids, we find that some kilonovae must have Mej < 0.03 Msun or Xlan>1e-4 or phi>30deg to be consistent with our limits. (Abridged)

Cosmic neutrinos provide a unique window into the otherwise-hidden mechanism of particle acceleration in astrophysical objects. A flux of high-energy neutrinos was discovered in 2013, and the IceCube Collaboration recently associated one high-energy neutrino with a flare from the relativistic jet of an active galaxy pointed towards the Earth. However a combined analysis of many similar active galaxies revealed no excess from the broader population, leaving the vast majority of the cosmic neutrino flux unexplained. Here we present the association of a radio-emitting tidal disruption event (AT2019dsg) with another high-energy neutrino, identified as part of our systematic search for optical counterparts to high-energy neutrinos with the Zwicky Transient Facility (ZTF). The probability of finding any radio-emitting tidal disruption event by chance is 0.5%, while the probability of finding one as bright in bolometric energy flux as AT2019dsg is 0.2%. Our electromagnetic observations can be explained through a multi-zone model, with radio analysis revealing a central engine, embedded in a UV photosphere, that powers an extended synchrotron-emitting outflow. This provides an ideal site for PeV neutrino production. The association suggests that tidal disruption events contribute to the cosmic neutrino flux. Unlike previous work which considered the rare subset of tidal disruption events with relativistic jets, our observations of AT2019dsg suggest an empirical model with a mildly-relativistic outflow.

(Abridged) Using the Zwicky Transient Facility alert stream, we are conducting a large campaign to spectroscopically classify all transients occurring in galaxies in the Census of the Local Universe (CLU) catalog. The aim of the experiment is to construct a spectroscopically complete, volume-limited sample of transients coincident within 100" of CLU galaxies out to 200 Mpc, and to a depth of 20 mag. We describe the survey design and spectroscopic completeness from the first 16 months of operations. We present results from a systematic search for Calcium rich gap transients in the sample of 22 low luminosity (peak absolute magnitude $M > -17$), hydrogen poor events found in the experiment (out of 754 spectroscopically classified SNe). We report the detection of eight Calcium rich gap transients, and constrain their volumetric rate to be at least $\approx 15\pm5$% of the SN Ia rate. Combining this sample with ten events from the literature, we find a likely continuum of spectroscopic properties ranging from events with SN Ia-like features (Ca-Ia objects) to SN Ib/c-like features (Ca-Ib/c objects) at peak light. Within the Ca-Ib/c events, we find two populations of events distinguished by their red ($g - r \approx 1.5$ mag) or green ($g - r \approx 0.5$ mag) spectral colors at $r$-band peak, wherein redder events show strong line blanketing signatures, slower light curves, weaker He lines and lower [Ca II]/[O I] in the nebular phase. Together, we find that the spectroscopic continuum, volumetric rates and striking old environments are consistent with the explosive burning of He shells on low mass white dwarfs. We posit that Ca-Ia and red Ca-Ib/c objects are consistent with the double detonation of He shells with high He burning efficiency, while green Ca-Ib/c objects could arise from less efficient He burning scenarios such as detonations in low density He shells or He shell deflagrations.

Possessing the strongest magnetic fields in the Universe, magnetars mark an extremum of physical phenomena. The strength of their magnetic fields is sufficient to deform the shape of the stellar body, and when the rotational and magnetic axes are not aligned, these deformations lead to the production of gravitational waves (GWs) via a time-varying quadrupole moment. Such gravitational radiation differs from signals presently detectable by the Laser Interferometer Gravitational-Wave Observatory. These signals are continuous rather than the momentary 'chirp' waveforms produced by binary systems during the phases of inspiral, merger, and ringdown. We evaluate the sensitivity requirements of future iterations of GW detectors to continuous GW signals resulting from magnetars. Here, we construct a computational model for magnetar stellar structure with strong internal magnetic fields. We implement an n = 1 polytropic equation of state (EOS) and adopt a mixed poloidal and toroidal magnetic field model constrained by the choice of EOS. We utilize fiducial values for magnetar magnetic field strength and various stellar physical attributes. Via computational simulation, we measure the deformation of magnetar stellar structure to determine upper bounds on the strength of continuous GWs formed as a result of these deformations inducing non-axisymmetric rotation. We compute predictions of upper limit GW strain values for sources in the McGill Magnetar Catalog, an index of all detected magnetars.

We present a survey of the early evolution of 12 Type IIn supernovae (SNe IIn) in the Ultra-Violet (UV) and visible light. We use this survey to constrain the geometry of the circumstellar material (CSM) surrounding SN IIn explosions, which may shed light on their progenitor diversity. In order to distinguish between aspherical and spherical circumstellar material (CSM), we estimate the blackbody radius temporal evolution of the SNe IIn of our sample, following the method introduced by Soumagnac et al. We find that higher luminosity objects tend to show evidence for aspherical CSM. Depending on whether this correlation is due to physical reasons or to some selection bias, we derive a lower limit between 35% and 66% on the fraction of SNe IIn showing evidence for aspherical CSM. This result suggests that asphericity of the CSM surrounding SNe IIn is common - consistent with data from resolved images of stars undergoing considerable mass loss. It should be taken into account for more realistic modelling of these events.

While tidal disruption events (TDEs) have long been heralded as laboratories for the study of quiescent black holes, the small number of known TDEs and uncertainties in their emission mechanism have hindered progress towards this promise. Here present 17 new TDEs that have been detected recently by the Zwicky Transient Facility along with Swift UV and X-ray follow-up observations. Our homogeneous analysis of the optical/UV light curves, including 22 previously known TDEs from the literature, reveals a clean separation of light curve properties with spectroscopic class. The TDEs with Bowen fluorescence features in their optical spectra have smaller blackbody radii, as well as longer rise times and higher disruption rates compared to the rest of the sample. The Bowen fluorescence mechanism requires a high density which can be reached at smaller radii, which in turn yields longer diffusion timescales. Thus, the difference in rise times suggests the pre-peak TDE light curves are governed not by the fallback timescale, but instead by the diffusion of photons through the tidal debris. The small subset of TDEs that show only helium emission lines in their spectra have the longest rise times, the highest luminosities and the lowest rates. We also report, for the first time, the detection of soft X-ray flares from a TDE on day timescales. Based on the fact the flares peak at a luminosity similar to the optical/UV blackbody luminosity, we attribute them to brief glimpses through a reprocessing layer that otherwise obscures the inner accretion flow.

The Zwicky Transient Facility (ZTF) is performing a three-day cadence survey of the visible Northern sky (~3$\pi$). The transient candidates found in this survey are announced via public alerts. As a supplementary product ZTF is also conducting a large spectroscopic campaign: the ZTF Bright Transient Survey (BTS). The goal of the BTS is to spectroscopically classify all extragalactic transients brighter than 18.5 mag at peak brightness and immediately announce those classifications to the public. Extragalactic discoveries from ZTF are predominantly Supernovae (SNe). The BTS is the largest flux-limited SN survey to date. Here we present a catalog of the761 SNe that were classified during the first nine months of the survey (2018 Apr. 1 to 2018 Dec. 31). The BTS SN catalog contains redshifts based on SN template matching and spectroscopic host galaxy redshifts when available. Based on this data we perform an analysis of the redshift completeness of local galaxy catalogs, dubbed as the Redshift Completeness Fraction (RCF; the number of SN host galaxies with known spectroscopic redshift prior to SN discovery divided by the total number of SN hosts). In total, we identify the host galaxies of 512 Type Ia supernovae, 227 of which have known spectroscopic redshifts, yielding an RCF estimate of $44\% \pm1\%$. We find a steady decrease in the RCF with increasing distance in the local universe. For z<0.05, or ~200 Mpc, we find RCF=0.6, which has important ramifications when searching for multimessenger astronomical events. Prospects for dramatically increasing the RCF are limited to new multi-fiber spectroscopic instruments, or wide-field narrowband surveys. We find that existing galaxy redshift catalogs are only $50\%$ complete at $r\approx16.9$ mag. Pushing this limit several magnitudes deeper will pay huge dividends when searching for electromagnetic counterparts to gravitational wave events.

We present DeepStreaks, a convolutional-neural-network, deep-learning system designed to efficiently identify streaking fast-moving near-Earth objects that are detected in the data of the Zwicky Transient Facility (ZTF), a wide-field, time-domain survey using a dedicated 47 sq. deg camera attached to the Samuel Oschin 48-inch Telescope at the Palomar Observatory in California, United States. The system demonstrates a 96-98% true positive rate, depending on the night, while keeping the false positive rate below 1%. The sensitivity of DeepStreaks is quantified by the performance on the test data sets as well as using known near-Earth objects observed by ZTF. The system is deployed and adapted for usage within the ZTF Solar-System framework and has significantly reduced human involvement in the streak identification process, from several hours to typically under 10 minutes per day.

We present photometry and spectroscopy of four hydrogen-poor luminous supernovae discovered during the two-month science commissioning and early operations of the Zwicky Transient Facility (ZTF) survey. Three of these objects, SN2018bym (ZTF18aapgrxo), SN2018avk (ZTF18aaisyyp) and SN2018bgv (ZTF18aavrmcg) resemble typical SLSN-I spectroscopically, while SN2018don (ZTF18aajqcue) may be an object similar to SN2007bi experiencing considerable host galaxy reddening, or an intrinsically long-lived, luminous and red SN Ic. We analyze the light curves, spectra, and host galaxy properties of these four objects and put them in context of the population of SLSN-I. SN2018bgv stands out as the fastest-rising SLSN-I observed to date, with a rest-frame g-band rise time of just 10 days from explosion to peak -- if it is powered by magnetar spin-down, the implied ejecta mass is only ~1 M$_{\odot}$. SN2018don also displays unusual properties -- in addition to its red colors and comparatively massive host galaxy, the light curve undergoes some of the strongest light curve undulations post-peak seen in a SLSN-I, which we speculate may be due to interaction with circumstellar material. We discuss the promises and challenges of finding SLSNe in large-scale surveys like ZTF given the observed diversity in the population.

The third observing run by LVC has brought the discovery of many compact binary coalescences. Following the detection of the first binary neutron star merger in this run (LIGO/Virgo S190425z), we performed a dedicated follow-up campaign with the Zwicky Transient Facility (ZTF) and Palomar Gattini-IR telescopes. The initial skymap of this single-detector gravitational wave (GW) trigger spanned most of the sky observable from Palomar Observatory. Covering 8000 deg$^2$ of the initial skymap over the next two nights, corresponding to 46\% integrated probability, ZTF system achieved a depth of $\approx$\,21 $m_\textrm{AB}$ in $g$- and $r$-bands. Palomar Gattini-IR covered 2200 square degrees in $J$-band to a depth of 15.5\,mag, including 32\% integrated probability based on the initial sky map. The revised skymap issued the following day reduced these numbers to 21\% for the Zwicky Transient Facility and 19\% for Palomar Gattini-IR. We narrowed 338,646 ZTF transient "alerts" over the first two nights of observations to 15 candidate counterparts. Two candidates, ZTF19aarykkb and ZTF19aarzaod, were particularly compelling given that their location, distance, and age were consistent with the GW event, and their early optical lightcurves were photometrically consistent with that of kilonovae. These two candidates were spectroscopically classified as young core-collapse supernovae. The remaining candidates were ruled-out as supernovae. Palomar Gattini-IR did not identify any viable candidates with multiple detections only after merger time. We demonstrate that even with single-detector GW events localized to thousands of square degrees, systematic kilonova discovery is feasible.

We report the discovery of six active galactic nuclei (AGN) caught "turning on" during the first nine months of the Zwicky Transient Facility (ZTF) survey. The host galaxies were classified as LINERs by weak narrow forbidden line emission in their archival SDSS spectra, and detected by ZTF as nuclear transients. In five of the cases, we found via follow-up spectroscopy that they had transformed into broad-line AGN, reminiscent of the changing-look LINER iPTF 16bco. In one case, ZTF18aajupnt/AT2018dyk, follow-up HST UV and ground-based optical spectra revealed the transformation into a narrow-line Seyfert 1 (NLS1) with strong [Fe VII, X, XIV] and He II 4686 coronal lines. Swift monitoring observations of this source reveal bright UV emission that tracks the optical flare, accompanied by a luminous soft X-ray flare that peaks ~60 days later. Spitzer follow-up observations also detect a luminous mid-infrared flare implying a large covering fraction of dust. Archival light curves of the entire sample from CRTS, ATLAS, and ASAS-SN constrain the onset of the optical nuclear flaring from a prolonged quiescent state. Here we present the systematic selection and follow-up of this new class of changing-look LINERs, compare their properties to previously reported changing-look Seyfert galaxies, and conclude that they are a unique class of transients well-suited to test the uncertain physical processes associated with the LINER accretion state.

We present the discovery of an optical transient (OT) in Messier 51, designated M51 OT2019-1 (also ZTF19aadyppr, AT 2019abn, ATLAS19bzl), by the Zwicky Transient Facility (ZTF). The OT rose over 15 days to an observed luminosity of $M_r=-13$ (${\nu}L_{\nu}=9\times10^6~L_{\odot}$), in the luminosity gap between novae and typical supernovae (SNe). Spectra during the outburst show a red continuum, Balmer emission with a velocity width of $\approx400$ km s$^{-1}$, Ca II and [Ca II] emission, and absorption features characteristic of an F-type supergiant. The spectra and multiband light curves are similar to the so-called "SN impostors" and intermediate-luminosity red transients (ILRTs). We directly identify the likely progenitor in archival Spitzer Space Telescope imaging with a $4.5~\mu$m luminosity of $M_{[4.5]}\approx-12.2$ and a $[3.6]-[4.5]$ color redder than 0.74 mag, similar to those of the prototype ILRTs SN 2008S and NGC 300 OT2008-1. Intensive monitoring of M51 with Spitzer further reveals evidence for variability of the progenitor candidate at [4.5] in the years before the OT. The progenitor is not detected in pre-outburst Hubble Space Telescope optical and near-IR images. The optical colors during outburst combined with spectroscopic temperature constraints imply a higher reddening of $E(B-V)\approx0.7$ mag and higher intrinsic luminosity of $M_r\approx-14.9$ (${\nu}L_{\nu}=5.3\times10^7~L_{\odot}$) near peak than seen in previous ILRT candidates. Moreover, the extinction estimate is higher on the rise than on the plateau, suggestive of an extended phase of circumstellar dust destruction. These results, enabled by the early discovery of M51 OT2019-1 and extensive pre-outburst archival coverage, offer new clues about the debated origins of ILRTs and may challenge the hypothesis that they arise from the electron-capture induced collapse of extreme asymptotic giant branch stars.

Efficient automated detection of flux-transient, reoccurring flux-variable, and moving objects is increasingly important for large-scale astronomical surveys. We present braai, a convolutional-neural-network, deep-learning real/bogus classifier designed to separate genuine astrophysical events and objects from false positive, or bogus, detections in the data of the Zwicky Transient Facility (ZTF), a new robotic time-domain survey currently in operation at the Palomar Observatory in California, USA. Braai demonstrates a state-of-the-art performance as quantified by its low false negative and false positive rates. We describe the open-source software tools used internally at Caltech to archive and access ZTF's alerts and light curves (Kowalski), and to label the data (Zwickyverse). We also report the initial results of the classifier deployment on the Edge Tensor Processing Units (TPUs) that show comparable performance in terms of accuracy, but in a much more (cost-) efficient manner, which has significant implications for current and future surveys.

We investigate ZTF18aalrxas, a double-peaked Type IIb core-collapse supernova (SN) discovered during science validation of the Zwicky Transient Facility (ZTF). ZTF18aalrxas was discovered while the optical emission was still rising towards the initial cooling peak (0.7 mag over 2 days). Our observations consist of multi-band (UV, optical) light-curves, and optical spectra spanning from $\approx0.7$ d to $\approx180$ d past the explosion. We use a Monte-Carlo based non-local thermodynamic equilibrium (NLTE) model, that simultanously reproduces both the $\rm ^{56}Ni$ powered bolometric light curve and our nebular spectrum. This model is used to constrain the synthesized radioactive nickel mass (0.17 $\mathrm{M}_{\odot}$) and the total ejecta mass (1.7 $\mathrm{M}_{\odot}$) of the SN. The cooling emission is modeled using semi-analytical extended envelope models to constrain the progenitor radius ($790-1050$ $\mathrm{R}_{\odot}$) at the time of explosion. Our nebular spectrum shows signs of interaction with a dense circumstellar medium (CSM), and this spetrum is modeled and analysed to constrain the amount of ejected oxygen ($0.3-0.5$ $\mathrm{M}_{\odot}$) and the total hydrogen mass ($\approx0.15$ $\mathrm{M}_{\odot}$) in the envelope of the progenitor. The oxygen mass of ZTF18aalrxas is consistent with a low ($12-13$ $\mathrm{M}_{\odot}$) Zero Age Main Sequence mass progenitor. The light curves and spectra of ZTF18aalrxas are not consistent with massive single star SN Type IIb progenitor models. The presence of an extended hydrogen envelope of low mass, the presence of a dense CSM, the derived ejecta mass, and the late-time oxygen emission can all be explained in a binary model scenario.

We report the discovery of non-stellar hydrogen Balmer and metastable helium absorption lines accompanying a transient, high-velocity (0.05$c$) broad absorption line (BAL) system in the optical spectra of the tidal disruption event (TDE) AT2018zr ($z=0.071$). In the HST UV spectra, absorption of high- and low-ionization lines are also present at this velocity, making AT2018zr resemble a low-ionization broad absorption line (LoBAL) QSO. We conclude that these transient absorption features are more likely to arise in fast outflows produced by the TDE than absorbed by the unbound debris. In accordance with the outflow picture, we are able to reproduce the flat-topped H$\alpha$ emission in a spherically expanding medium, without invoking the typical prescription of an elliptical disk. We also report the appearance of narrow ($\sim$1000~km~s$^{-1}$) NIII$\lambda$4640, HeII$\lambda4686$, H$\alpha$, and H$\beta$, emission in the late-time optical spectra of AT2018zr, which may be a result of UV continuum hardening at late time as observed by Swift. Including AT2018zr, we find a high association rate (3 out of 4) of BALs in the UV spectra of TDEs. This suggests that outflows may be ubiquitous among TDEs and may be less sensitive to viewing angle effects compared to QSO outflows.

We present Zwicky Transient Facility (ZTF) observations of the tidal disruption flare AT2018zr/PS18kh reported by Holoien et al. and detected during ZTF commissioning. The ZTF light curve of the tidal disruption event (TDE) samples the rise-to-peak exceptionally well, with 50 days of g- and r-band detections before the time of maximum light. We also present our multi-wavelength follow-up observations, including the detection of a thermal (kT~100 eV) X-ray source that is two orders of magnitude fainter than the contemporaneous optical/UV blackbody luminosity, and a stringent upper limit to the radio emission. We use observations of 128 known active galactic nuclei (AGN) to assess the quality of the ZTF astrometry, finding a median host-flare distance of 0.2" for genuine nuclear flares. Using ZTF observations of variability from known AGN and supernovae we show how these sources can be separated from TDEs. A combination of light-curve shape, color, and location in the host galaxy can be used to select a clean TDE sample from multi-band optical surveys such as ZTF or LSST.

The Zwicky Transient Facility (ZTF) is a new optical time-domain survey that uses the Palomar 48-inch Schmidt telescope. A custom-built wide-field camera provides a 47 deg$^2$ field of view and 8 second readout time, yielding more than an order of magnitude improvement in survey speed relative to its predecessor survey, the Palomar Transient Factory (PTF). We describe the design and implementation of the camera and observing system. The ZTF data system at the Infrared Processing and Analysis Center provides near-real-time reduction to identify moving and varying objects. We outline the analysis pipelines, data products, and associated archive. Finally, we present on-sky performance analysis and first scientific results from commissioning and the early survey. ZTF's public alert stream will serve as a useful precursor for that of the Large Synoptic Survey Telescope.

The Zwicky Transient Facility is a large optical survey in multiple filters producing hundreds of thousands of transient alerts per night. We describe here various machine learning (ML) implementations and plans to make the maximal use of the large data set by taking advantage of the temporal nature of the data, and further combining it with other data sets. We start with the initial steps of separating bogus candidates from real ones, separating stars and galaxies, and go on to the classification of real objects into various classes. Besides the usual methods (e.g., based on features extracted from light curves) we also describe early plans for alternate methods including the use of domain adaptation, and deep learning. In a similar fashion we describe efforts to detect fast moving asteroids. We also describe the use of the Zooniverse platform for helping with classifications through the creation of training samples, and active learning. Finally we mention the synergistic aspects of ZTF and LSST from the ML perspective.

Wide-field optical surveys have begun to uncover large samples of fast (t_rise < 5d), luminous (M_peak < -18), blue transients. While commonly attributed to the breakout of a supernova shock into a dense wind, the great distances to the transients of this class found so far have hampered detailed investigation of their properties. We present photometry and spectroscopy from a comprehensive worldwide campaign to observe AT2018cow (ATLAS18qqn), the first fast-luminous optical transient to be found in real time at low redshift. Our first spectra (<2 days after discovery) are entirely featureless. A very broad absorption feature suggestive of near-relativistic velocities develops between 3-8 days, then disappears. Broad emission features of H and He develop after >10 days. The spectrum remains extremely hot throughout its evolution, and the photospheric radius contracts with time (receding below R<10^14 cm after 1 month). This behaviour does not match that of any known supernova, although a relativistic jet within a fallback supernova could explain some of the observed features. Alternatively, the transient could originate from the disruption of a star by an intermediate-mass black hole, although this would require long-lasting emission of highly super-Eddington thermal radiation. In either case, AT2018cow suggests that the population of fast luminous transients represents a new class of astrophysical event. Intensive follow-up of this event in its late phases, and of any future events found at comparable distance, will be essential to better constrain their origins.

A fraction of active galactic nuclei (AGN) exhibit dramatic variability, which is observed on timescales down to minutes in the X-ray band. We introduce the case study of 1H 1934-063 (z = 0.0102), a Narrow-line Seyfert 1 (NLS1) among the brightest and most variable AGN ever observed with XMM-Newton. This work includes spectral and temporal analyses of a concurrent XMM-Newton and NuSTAR 2015 observation lasting 130 kiloseconds, during which the X-ray source exhibited a steep (factor of 6) plummet and subsequent full recovery of flux level, accompanied by deviation from a single log-normal flux distribution. We rule out Compton-thin obscuration as the cause for this dramatic variability observed even at NuSTAR energies. In order to constrain coronal geometry, dynamics, and emission/absorption processes, we compare detailed spectral fitting with Fourier-based timing analysis. Similar to other well-studied, highly variable Seyfert 1s, this AGN is X-ray bright and displays strong reflection features. We find a narrower broad iron line component compared to most Seyfert 1s, and constrain black hole spin to be < 0.1, one of the lowest yet discovered for such systems. Combined spectral and timing results are consistent with a dramatic change in the continuum on timescales as short as a few kiloseconds dictating the nature of this variability. We also discover a Fe-K time lag, measuring a delay of 20 seconds between relativistically-blurred reflection off the inner accretion flow and the hard X-ray continuum emission.

We present the results of a population synthesis of radio and gamma-ray millisecond pulsars (MSPs) from the Galactic Disk (GD). Using 92 radio millisecond pulsars detected in 13 radio surveys and 54 Fermi MSPs detected as point sources in the first point source catalog, we establish six free parameters corresponding to the overall factor and the exponents of the period and period derivative dependence for each of the radio and gamma-ray empirical luminosity models. We test three high-energy emission models described by the two pole caustic Slot Gap, Outer Gap, and Pair Starved Polar Cap geometries. The simulated distributions of pulsar properties adequately describe the distributions of detected MSPs from the GD. We explore the gamma-ray emission from groups of MSPs in globular clusters and in the Galactic Bulge. The simulation predicts reasonable numbers of Fermi MSPs detected in the other point source catalogs and anticipates a bright future for Fermi observations of MSPs, expecting a total of ~ 170 MSP detections from the GD within ten years. Our numbers of simulated MSPs in globular clusters are in agreement with those derived from Fermi detections. The simulation predicts about 11,000 MSPs in the Galactic Bulge are required to explain the gamma-ray Galactic Center Excess.