A Search for Pulsars around Sgr A* in the First Event Horizon Telescope Data Set

In 2017 the Event Horizon Telescope (EHT) observed the supermassive black hole at the center of the Milky Way, Sagittarius A* (Sgr A*), at a frequency of 228.1 GHz (λ = 1.3 mm). The fundamental physics tests that even a single pulsar orbiting Sgr A* would enable motivate searching for pulsars in EHT data sets. The high observing frequency means that pulsars—which typically exhibit steep emission spectra—are expected to be very faint. However, it also negates pulse scattering, an effect that could hinder pulsar detections in the Galactic center. Additionally, magnetars or a secondary inverse Compton emission could be stronger at millimeter wavelengths than at lower frequencies. We present a search for pulsars close to Sgr A* using the data from the three most sensitive stations in the EHT 2017 campaign: the Atacama Large Millimeter/submillimeter Array, the Large Millimeter Telescope, and the IRAM 30 m Telescope. We apply three detection methods based on Fourier-domain analysis, the fast folding algorithm, and single-pulse searches targeting both pulsars and burst-like transient emission. We use the simultaneity of the observations to confirm potential candidates. No new pulsars or significant bursts were found. Being the first pulsar search ever carried out at such high radio frequencies, we detail our analysis methods and give a detailed estimation of the sensitivity of the search. We conclude that the EHT 2017 observations are only sensitive to a small fraction (≲2.2%) of the pulsars that may exist close to Sgr A*, motivating further searches for fainter pulsars in the region.


INTRODUCTION
The first test of strong-field gravity came from measurements of the relativistic orbital decay in a binary pulsar system (Taylor & Weisberg 1982, 1989), with sub-sequent tests of increasing precision using other binary (or even triple) pulsar systems (see e.g.Wex & Kramer 2020, for a review).So far, all gravity experiments using pulsars conform with the predictions of general relativity (GR), but it is expected that the most constraining tests will come from a pulsar black hole binary system (Wex & Kopeikin 1999;Liu et al. 2014), and in particular from pulsars orbiting Sagittarius A* (Sgr A*) itself (Kramer et al. 2004;Liu et al. 2012).Discovery and timing observations of pulsars near Sgr A* could provide measurements of the spin and quadrupole moment of the supermassive black hole (SMBH), in addition to unique information on Galactic Centre (GC) stellar populations, dark matter, and the γ−ray excess (Ajello et al. 2016;Daylan et al. 2016;Bartels et al. 2016), along with measurements of the magnetoionic plasma (Eatough et al. 2013;Desvignes et al. 2018).
Together with the groundbreaking measurements of binary black hole mergers with gravitational wave detectors (Abbott et al. 2016) and high-precision astrometry of stars orbiting Sgr A* (Gravity Collaboration et al. 2018, 2019, 2020;Do et al. 2019), millimeter interferometric imaging of SMBHs by the Event Horizon Telescope (EHT) Collaboration represents a transformation in the way black holes can be observed and studied (see Event Horizon Telescope Collaboration et al. 2019a, 2022a, and references therein for an overview of key results).Comparison of EHT image properties with synthetic images derived from general relativistic magnetohydrodynamic (GRMHD) simulations, and general relativistic ray tracing, provide a new framework in which to measure the fundamental properties of black holes and test theories of gravity in the strong-field regime (e.g., Mizuno et al. 2018;Event Horizon Telescope Collaboration et al. 2019b,c, 2022b;Özel et al. 2022;Younsi et al. 2023).
EHT images of the SMBHs Messier 87* (M87*) and Sgr A* are consistent with predictions for the black hole shadow of a spinning Kerr black hole in GR, whilst certain alternative theories of gravity, which deviate from the Kerr metric in the strong-field regime, have been ruled out (Psaltis et al. 2020).For Sgr A*, the EHT images can be synergized with the measurements from an orbiting pulsar.For example Psaltis et al. (2016) showed that the spin and quadrupole moment of Sgr A* from the motion of orbiting stars (e.g.Weinberg et al. 2005) and pulsars have correlated uncertainties that are almost orthogonal to those from black hole shadow images, thereby increasing overall measurement accuracy.
Multipath propagation by scattering from electrondensity inhomogeneities typically broadens radio pulse emission, reducing the signal-to-noise ratio of periodic or single-burst radio signals (Cordes & Lazio 1997).The severity of this effect depends on the observing frequency ν, (∼ ν n ; where n ≈ −4), spin period (P ), pulse or burst width, and also on the scattering environment toward the GC.The NE2001 model of the Galactic distribution of free electrons (Cordes & Lazio 2002) predicts a large scattering time scale that is inconsistent with the combined measurements of PSR J1745−2900 and Sgr A* (Spitler et al. 2014;Bower et al. 2014Bower et al. , 2015) ) but may apply to other lines of sight toward (≲ 0.5 • ) the GC (Dexter et al. 2017).Therefore, one rationale is to search the GC at mm and sub-mm wavelengths where multipath scattering can be avoided altogether.However, given the average steep spectrum of pulsars' radio emission where flux density ∼ ν α , with α typically between −1 and −2.5 (see e.g., Jankowski et al. 2018), many pulsars are expected to be faint at frequencies where scattering is negligible.Nevertheless, a subset of pulsars and their high magnetic field counterparts, magnetars, show flat or inverted spectra 1 with detectable mm and sub-mm emission (Camilo et al. 2007;Torne et al. 2020;Chu et al. 2021;Torne et al. 2022), including from PSR J1745−2900 itself (Torne et al. 2015(Torne et al. , 2017)).
Additional physical processes may contribute to emission in these bands and could enhance the probability of detection.For example, induced inverse Compton scattering of low-frequency radio photons at frequency ν can produce higher-frequency photons at ν ′ ∼ γ 2 ν in some pulsar models (e.g.Blandford & Scharlemann 1976;Philippov & Kramer 2022), where γ is the Lorentz factor.For example, with ν = 100 MHz and γ = 50, scattered photons will boost emission at ν ∼ 250 GHz.Since a wide range of Lorentz factors is possible, it is necessary to search in a corresponding wide range of observing bands with sensitive mm and sub-mm telescopes.
The first attempts at searching the GC for pulsars at mm wavelengths were undertaken with the IRAM 30 m 1 Around 13% of known pulsars with a measured spectral index have α ≥ −1.0 (Manchester et al. 2005, PSRCAT in 2023).
telescope (λ ≃ 2, 3 mm) and the Atacama Large Millimeter/submillimeter Array (λ ≃ 3.5 mm) (Torne et al. 2021;Liu et al. 2021).In this work we report the first 1.3 mm searches for GC pulsars using data from the EHT 2017 campaign on Sgr A*.As well as utilizing the most sensitive mm/sub-mm telescopes, the simultaneous EHT observations from multiple telescope sites enable efficient rejection of false positive detections.
In Section 2 we describe the observations and data processing; in Section 3 the results of our pulsar search of the EHT data are detailed; Section 4 is a discussion of the findings; Section 5 presents potential improvements for the future and Section 6 summarizes our conclusions.

Observations
The EHT 2017 observing campaign was scheduled on five nights during 2017 April 5-14, where three nights (April 6,7 and 11) included exposures to Sgr A* (α J2000 = 17 h 45 m 40 s .0361,δ J2000 = −29 • 00 ′ 28 ′′ .168).At each epoch, the track on Sgr A* was divided into individual scans of approximately 7−10 min each, switching between Sgr A* and an active galactic nucleus calibrator source, J1924−2914 and/or NRAO 530.Baseband voltage data were recorded with 2-bit sampling at a total rate of 32 Gbps in two 2-GHz IF-bands centered at 229.1 and 227.1 GHz, labeled "high" and "low", respectively.More details of the EHT array, the observations and data recording can be found in Event Horizon Telescope Collaboration et al. (2019dCollaboration et al. ( ,e, 2022c)).
Since system sensitivity is one of the primary considerations for a pulsar search of the GC, we chose to analyze data only from the three most sensitive telescopes in the EHT 2017 observations: the phased Atacama Large Millimeter/submillimeter Array (ALMA), the Large Millimeter Telescope Alfonso Serrano (LMT) and the IRAM 30 m telescope (PV).While the sensitivity of phased ALMA is significantly higher than the other two stations, its field of view, given by the full width at half maximum (FWHM) of its synthetic beam during the EHT 2017 observations on Sgr A*, is only 1-2 ′′ (Goddi et al. 2021); offering comparatively limited sky coverage.LMT and IRAM 30 m, both of which are used here as single-dish telescopes with a beam size of approximately 10 ′′ at 1.3 mm, are useful in supplementing the sky coverage of ALMA for the pulsar search (see Section 4.3 for more details) and covering the position of PSR J1745−2900.The total length of the tracks on Sgr A* of these telescopes varied from 5 to 10 hours on different nights, and are summarized in Table 1.

Data conversion and preparation
The baseband voltage data were captured by Mark 6 recorders2 at the telescopes, and the disk packs shipped to the correlator at the Max Planck Institute for Radio Astronomy in Bonn, Germany, for post-processing.The voltage data selected for the pulsar search were reduced into multi-channel intensity time series (commonly refered to as a filterbank) written in search-mode PSR-FITS format, a standard FITS specification for pulsar data (Hotan et al. 2004) 3 .This used the software tool MPIvdif2psrfits4 , an upgraded version of the original toolkit developed under the ALMA Pulsar Mode Project5 (Liu et al. 2019) that incorporates parallel processing capability using OpenMPI6 .The resulting properties of the PSRFITS products are presented in Table 1.
After conversion to PSRFITS the data consisted of a number of short scans for each of the two frequency sub-bands (high and low), corresponding to observations of Sgr A*, a calibrator, or another science target.From the scans on Sgr A*, those showing potential issues were flagged to be excluded from the analysis.The main reasons to flag scans were large variations or jumps of the mean power level, and, specifically for the case of ALMA, those scans not showing the expected array phasing noise features (see Section 2.2.1).
Once the selection of Sgr A* scans to be analyzed were identified for ALMA, LMT and IRAM 30 m, the data were further prepared for pulsar searching.The steps for each telescope are presented in the following subsections.

Atacama Large Millimeter/submillimeter Array
When ALMA is used as a phased array the scans are partitioned into subscans with updates of the phase corrections every 18.192 s to enable the coherent summation of signals from individual antennas (Goddi et al. 2019).This phasing cycle introduced a periodic feature in the time series, consisting in a decrease of the power level with a duration slightly below two seconds, followed by a large negative peak (we hereafter refer to this feature as the "broad" feature).Additionally, a forest of negative narrow spikes occur a few seconds after the broad feature (see Appendix in Liu et al. 2021, for details).The first step in the preparation of the ALMA data was therefore to try to remove those noise features, which could introduce undesired power in the Fourier domain and decrease the sensitivity to pulsars in the search.To locate Table 1.Details of observations and data analyzed in this paper.The first seven columns indicate the EHT station, the real (or equivalent for the phased ALMA array) single dish diameter (Ds), the effective beam FWHM, the central observing frequency (ν obs ), the total effective bandwidth (∆ν), the number of frequency channels (n ch ), and initial sample time (tsamp), in the resulting PSRFITS files ‡ .The remaining columns show, for each epoch, the total time span of the dataset constructed for the pulsar search (i.e., the time between the start of the first scan and the end of the last scan) and the net integration time on Sgr A*.The difference between the total time span and the net integration time is due to the interleaved observations of calibrators, slew time of the telescope, time for pointing and focus adjustments and the flagged data.The data with the total time span are the ones searched, but only the net integration time on Sgr A* is considered for the search sensitivity analysis.‡ During data preparation the number of frequency channels and sample time is reduced to varying degrees (see Section 2.2). * Due to slightly-overlapping channels in the digital filterbank, the effective bandwidth is 2×1.875 GHz (see Goddi et al. 2019).

DS
and clean the features, we first produced a time series collapsing all the frequency channels of the PSRFITS file with the routine prepdata of the software package presto (Ransom 2011) 7 , downsampling the time resolution to 32 µs, and visually inspected each scan, marking the center of the first appearance of the broad feature.Then, two-second wide windows around the first broad feature position, and every 18.192 s thereafter, were replaced with the mean level of the time series calculated outside the marked windows.Next, we use a moving window of width 10 s that clips negative narrow spikes below a threshold of 5 times the standard deviation of the samples within the window.The flagged data due to the phasing features correspond to ≈11% of the total.
Before combining the two frequency subbands (low and high) to obtain a full-bandwidth dataset, we remove slow variations of the mean power level in the time series with a "running-fit" filter that fits and subtracts a first order polynomial in a running window of 10 s.This step reduces the excess of power in the low-frequency end of the Fourier-domain data (also known as 'red noise'), which tend to be large in mm-observations due to atmospheric opacity variations during observations.We also normalize the time series by its standard deviation to combine the low and high subbands with a similar root mean square (RMS) level.At this stage, the cleaned time series of each individual scan and sub-band is saved as a sigproc8 -format filterbank file (Lorimer 2011) with one single frequency channel of 2 GHz and a time sampling of 32 µs, except for April 7 where we downsample to 64 µs to reduce the data size and speed up the process-ing of this particularly-long epoch.The production of the sigproc filterbank files is made with custom Python code supported by tools from presto and sigpyproc 9 .
Then, for each sub-band, we concatenate all the scans for each observing track into a single, continuous sigproc filterbank file, padding the gaps in between scans with the mean value.This concatenation is necessary to maintain the coherence of the data and maximize the sensitivity to periodic signals like those expected from pulsars.Before the next step, we ensure that the start time and length are equal between the concatenated files in the low-and high-band datasets, adjusting the length if necessary.Finally, we splice together the two frequency sub-bands into one file with the splice routine of sigproc.The result is a two-channel filterbank file containing all the scans of the observing track with the full bandwidth.This process is repeated for each night, and the final filterbank files are the ones passed to the searching pipelines (see Section 2.3).
In a separate step, we Fourier transformed and analyzed with presto's accelsearch routine three scans on calibrators (one for each observing night) to obtain a list of locally-generated periodic signals through an analysis of the corresponding Fourier spectrum.A list of periodicities -most likely produced in the receiving chain -is created, and later used to flag them when applying the searching pipeline to the observations on Sgr A* by zero-weighting the corresponding Fourier bins.

IRAM 30m Telescope
The IRAM 30m telescope data were mainly affected by two features: a strong ripple in the time series with a frequency of about 210.7 Hz, and regular power drop-offs, part of them synchronized with the 210.7 Hz cycle (see Figures 6 and 7 in the Appendix).The origin of these noise features is under investigation, suspected to be related to oscillations in the bias circuits of the first mixer of the receiver.To reduce the impact in the subsequent analysis, we first zero-weight a few Fourier bins centered on the 210.7 Hz signal, including 3 harmonics.We also zero-weight a few bins around 1 Hz, to remove another strong periodic signal related to the cryogenic pump cycle.To reduce the number of power drop-offs, we visually inspect the resulting time series, and manually mark a threshold below which all samples are substituted with the median of the remaining data.The percentage of data flagged related to the power dropoffs is on average 2.8%.The following steps are analogous as for the ALMA data: a running-fit filter with a window of 10 s to subtract slow mean level variations, a normalization of the data by their standard deviation, and a concatenation of scans plus splicing into a single filterbank file.Five scans showing larger instabilities in the form of ripples or considerably more signal dropoffs than on average, amounting to 0.72 hr in total, were flagged and excluded from the analysis.

Large Millimeter Telescope
The data from the LMT show strong periodicities associated with the local receiver.The most prominent signals have a frequency of 1.2 Hz, probably related to the cryogenics refrigeration cycle, and around 185 Hz, which has an unknown origin.The data also show some power drop-offs, but less frequently than the ALMA and IRAM 30 m data.A variation of the mean power level is usually present, related to opacity variations during the observations (see Figures 6 and 7 in the Appendix).To reduce these undesired signals, we zero-weight the strongest peaks related to the 1.2 and 185 Hz signals in the Fourier domain, and after an inverse Fourier transform, use a moving window of 10 s that clips the powerdrops substituting them with the median value inside the window.The mean power level variations are minimized using the running-fit filter with a window of 10 s as done for the ALMA and IRAM 30 m data.The remaining steps -normalization, concatenation and splicing -are analogous to the ALMA and IRAM 30 m datasets.For the LMT, 16 scans in total (amounting to 2.32 hr) were flagged because they still showed clear artifacts after the cleaning, like jumps in the mean power level and regular trains of strong, wide-pulse-like features.
We remark that although the undesired signals are highly reduced in the data from the three stations by our cleaning algorithms, they are not fully removed.We show examples of the time series before and after the cleaning procedures in Appendix A.

Pulsar searching
Searching for pulsars via their inherent periodicity is done with two independent pipelines: a Fourier-domain search (Ransom et al. 2002) and a search using the Fast-Folding-Algorithm (FFA, Staelin 1969).Fourier-domain methods are widely used in pulsar searches and have proven successful in discovering pulsars with a variety of spin parameters (Lyne 2003).The FFA works particularly well at detecting long-period pulsars that often show narrow pulse profiles (Cameron et al. 2017;Morello et al. 2020).In addition, we searched the data for singlepulse burst-like transient emission10 as seen in pulsar giant pulses, Rotating Radio Transients (RRATs) and Fast Radio Bursts (FRBs) (see e.g.Keane et al. 2011).In the following subsections, we detail the search algorithms used.The parameters utilized for each pipeline are summarized in Tables 2 and 3.
Pulsar surveys at low radio frequencies typically search over a range of dispersion measures (DM) -the integrated column density of free electrons along the line-of-sight.However, there is no dedispersion in our pipelines because the effect of even large DMs is negligible at the EHT observing frequencies and for the target time resolution of 32 µs.Interstellar dispersion delays the pulses' arrival times scaling with DM ν −2 (see e.g., Lorimer & Kramer 2004), but a DM=10000 pc cm −3 would produce a pulse smearing of only ∆t DM ≃ 28 µs across our ∼4-GHz band.The highest DM known to date is for the magnetar PSR J1745−2900, located close to Sgr A*, at DM=1778 pc cm −3 (Eatough et al. 2013).The negligible dispersion has some disadvantages (see Sections 2. 3.3, 5.3), but reduces considerably the computing costs, which can be concentrated on the acceleration search.

Periodicity search in the Fourier domain
The Fourier-domain pipeline predominantly searches for pulsars by detecting and analysing peaks in the Fourier spectrum of the time series.We use the software presto, which includes the capability to search for accelerated pulsars, like those orbiting Sgr A* or other companions.The acceleration parameter space is searched using a template-matching algorithm by recovering the power spread over contiguous Fourier bins (Ransom et al. 2002).This power spread is the result of the Doppler effect when a pulsar changes its radial velocity during an observation.
A limitation of the template-matching algorithm is that it loses efficiency for observations covering a substantial fraction of the binary period (P b ).This limiting fraction is about 0.1P b (given a companion mass of 1.4 solar mass), i.e., the algorithm works well when the total observing time does not exceed about 10% of the binary period (e.g., Ransom 2001;Ng et al. 2015).To improve the sensitivity for longer observing spans, one can add a search in the line-of-sight (LoS) acceleration derivative (commonly referred to as "jerk") parameter space (Andersen & Ransom 2018), although this increases greatly the computational costs11 .The LoS acceleration (a l ) and its derivative (j l ) of a pulsar in orbit can be calculated as (e.g., Bagchi et al. 2013): and where a p , i, e, ω and A T are the semi-major axis, inclination angle, eccentricity, longitude of periastron and true anomaly of the orbit, respectively.Our presto-based pipeline is based on the one used by Torne et al. (2021).It first produces a time series from the filterbank file and then uses a Fast Fourier Transform (FFT) to transform the data to the frequency domain.Next, it searches for periodic signals as discussed above.To increase the sensitivity to long-period pulsars, the rednoise filtering routine of presto is applied to further decrease the red noise.Up to N h = 32 harmonics of each detected periodicity are summed.
To be sensitive to a wide range of different pulsar systems, two passes on each observation are made.In the first pass, highly accelerated pulsars are searched, and a second pass searches for accelerated pulsars including a search in the jerk parameter space.In both cases searches for isolated pulsars (i.e., with no acceleration) are included.The control of the ranges for surveying acceleration and jerk are made through the parameters zmax=max(|z|) and wmax=max(|w|) of the routine accelsearch, which in turn represent the maximum shift in the Fourier frequency bin (z) and its derivative (w) within an entire observation of length T obs that the pipeline will search.Following Andersen & Ransom (2018), these two terms can be written as: and where f , c, h represent the pulsar rotational frequency, speed of light, and the index of the Fourier harmonic.
Here, we use zmax=1200 and wmax=0 for the pass aiming for highly accelerated pulsars, and zmax=300 and wmax=900 for the pass including a jerk search.For isolated pulsars both parameters are set to zero.
The exact relationship between sensitivity, pulsar spin period, acceleration (or orbital period), jerk, pipeline parameters, and observing time is complex (see e.g., Liu et al. 2021;Eatough et al. 2021).Nonetheless, with the parameters zmax and wmax chosen here, we can correct for the LoS acceleration of virtually any pulsar orbiting Sgr A* down to orbital periods P b ≥ 2.5 yr (assuming circular orbits).Pulsars in orbits with shorter orbital periods can still be detected, depending on their strength and spin period.A more detailed discussion on the coverage of pulsar orbits by our search is presented in Section 4.3.
After all the pulsar candidates found by accelsearch are saved, a sifting step removes duplicated and harmonically related ones.This step admits pulsar candidates with σ sift ≥ 2.0 as calculated by the sifting.pycode 12 .Such a σ sift threshold allows for the detection of weak pulsars while maintaining a manageable number of candidates (see Torne et al. 2021, and Section 2.4).A final step uses prepfold to fold13 the data with the information of each candidate, producing plots that can be inspected to decide if any of the candidates corresponds to a real pulsar.
The folding step produces four plots from each candidate.Two of these plots correspond to the raw data with and without an optimization of candidate parameters from prepfold.The other two plots arise from a filtered version of the raw data using the rednoise filter from presto.The rednoise filter is very effective in reducing some local, interfering periodic signals in some datasets, and can in some cases be key to enabling a detection (for an example, see Figure 2 in Torne et al. 2021).This multiple candidate plot production, in exchange of increasing the number of candidate plots to inspect, decreases the risk of missing weak pulsars by an insufficient cleaning or by interference from the locally-generated periodic signals.

Periodicity search with Fast Folding Algorithm
Unlike the Fourier-domain pipeline, the FFA searches for pulsars by folding the time series at different trial periods, forming a sequence of profiles, and testing the significance of each profile using boxcar matched filters.The pulse width and arbitrary total pulse power produced by the filters are used to calculate the signalto-noise ratio (S/N FFA ) of the pulse profiles (see e.g., Cameron et al. 2017, for details).We used the software riptide 14 (Morello et al. 2020) as the basis for the FFA search pipeline.We include an acceleration search by resampling the time series using sigpyproc at a series of acceleration trials using a custom pipeline (Wongphechauxsorn et al., in prep.).This implementation of the FFA pipeline assumes a constant acceleration during an observation.
The riptide software uses a matched filter to detect the pulse in the folded data; making the number of folding profile bins (N bin ) an important parameter.If the duty cycle is less than one bin, the sensitivity will be reduced.Morello et al. (2020) demonstrated that the number of trials is proportional to N 3 bin .Folding with many profile bins consequently results in more trials and longer processing times.Furthermore, the acceleration step between trials is proportional to N bin , thus the total number of trials is proportional to N 4 bin .Our FFA pipeline uses the cleaned time series generated following the description in Section 2.2.We applied the FFA acceleration search to the time series for periods from P min = 1 s to P max = 1025 s, using N bin = 128.As a result, the pipeline provides full sensitivity to duty cycles down to about 1% 15 ).We used an acceleration range of a l = ±50 m s −2 , i.e., approximately equivalent to zmax=1400 for detecting a 1-s pulsar in the Fourier 14 https://github.com/v-morello/riptide 15The duty cycle was calculated using the period pulse width from PSRCAT (Manchester et al. 2005).N.b., less than 7% of known pulsars have a lower duty cycle.
domain acceleration search pipeline with N h = 32 and T obs of 4.5 hr.The FFA can therefore probe the same type of binary orbit as the Fourier-domain acceleration search for long-period pulsars.Furthermore, the FFA pipeline is sensitive to an acceleration derivative up to approximately 2.6×10 −5 m s −3 without needing any jerk search.The range of acceleration derivative covered is equivalent approximately to a wmax=10 for a 1-s pulsar in the Fourier domain acceleration search pipeline with the same parameters (N h = 32 and T obs = 4.5 hr).
In addition, to search for very narrow duty cycles, the FFA was repeated using N bin = 512 with no acceleration, meaning that the pipeline is sensitive to pulsars with a duty cycle larger than 0.19% (only ≈ 0.2% of currently-known pulsars have a lower duty cycle).All candidates detected with S/N ≥ 6 generated by either the accelerated or non-accelerated pipeline were folded to be inspected visually.

Single pulse search
While a pulsar signal is often characterized as a series of stable pulsations, individual single pulses -which can be orders of magnitude brighter than the averaged pulse emission -are observed in some cases (McLaughlin & Cordes 2003) and are detectable with alternative search techniques (Cordes & McLaughlin 2003).For this task we employed single pulse search.py in presto.This routine identifies pulses by matchedfiltering the time series with boxcar filters of different pulse widths up to a maximum of 300 samples; corresponding to 9.6 ms and 19.2 ms for our datasets resampled at 32 µs and 64 µs, respectively.To avoid excessive automatic time series flagging, we adjusted the single pulse search.pyinternal parameters lo std and hi std to 0.88 and 3.68 respectively.For the EHT data at 1.3 mm pulse dispersion is negligible, so it cannot be used as the marker of a genuine astrophysical pulse.Instead the simultaneous nature of observations conducted at multiple sites is utilized to match coincident pulses from each telescope in the common rest frame of the solar system barycentre.
In practice single pulse search.pygives information on the single pulse significance, σ SP , its arrival time along with the corresponding number of samples relative to the beginning of the observation and the pulse width (determined by the optimum boxcar filter size).By transforming all the single pulse arrival times from each site to their equivalent barycentric Modified Julian Date (MJD) − determined from presto using the NASA Jet Propulsion Laboratory planetary ephemeris DE40516 in tempo (Nice et al. 2015) 17 − we search for pulses that are coincident to an accuracy of one time sample (32 or 64 µs depending on dataset).
In addition, any pulse with σ SP > 12 is visually inspected (regardless of coincidence) to check for other astrophysical markers such as scatter broadening of the profile or 'multi-component profiles' as seen in many pulsars.This also accounts for different Sgr A* visibility windows at each station, their various individual singlepulse sensitivities, and the different spatial coverages due to different beam sizes (see Tables 1 and 3).

Verification of periodicity search pipelines with synthetic signal injections
Similarly to the method in Torne et al. (2021), before searching the data a number of pre-processing tests were carried out to verify and tune the pipelines.This was done by injecting synthetic pulsars into the real data and checking their correct detection.The signals were produced with a custom-version of sigproc's fake and injected to the filterbank files.The synthetic pulsars included slow-spinning pulsars (P ≥ 2 s), canonical pulsars (P ∼ 500 ms), and fast-spinning millisecond pulsars (MSPs; 1 ms < P ≤ 30 ms).Within each category, different intensities and duty cycles were used, and both isolated and binary pulsars in a range of configurations were injected.The simulated companions included neutron stars, stellar-mass black holes and a supermassive black hole with the mass of Sgr A*.
The tests served two main purposes: first, to allow fine adjustments of the pipeline parameters so that injected signals were recovered and secondly, to verify that even extreme pulsar systems, like MSPs in tight orbits around massive companions, were detected when the characteristics are within the theoretical limitations of the searching algorithms (see Section 4.3).Here we confirmed the need of a low value of presto's sifting.pysigma threshold, which is a significance threshold to accept or reject candidates after a sifting step (set to σ sift > 2.0), together with a weak requirement on the minimum number of harmonics required per candidate, that is set to one (Torne et al. 2021).

Fourier-domain acceleration search
The analysis with the Fourier-domain pipeline of the three epochs of ALMA resulted in 3146 pulsar candi-dates.However, no new pulsars are discovered.We identify a large number of candidates with a round number in the spin period or frequency (e.g.ν s = 437.500000, 111.111111, 2875.000000 Hz, etc.).Another remarkable characteristic of the detected signals is a very short period, with 56.3% of the candidates having P < 0.5 ms.Those characteristics, in particular the round numbers, suggest a human-generated origin.Furthermore, we observe similar signals in the analysis of off-source scans (see Section 2.2.1) and we therefore relate these found periodicities to the local receiving chain and the properties of the Superconductor-Insulator-Superconductor (SIS) receivers used in the observations (see also Torne et al. 2021).
The single epoch with the IRAM 30 m telescope resulted in 99 pulsar candidates.No signal resembling a real pulsar was found.The candidates are similar in properties to those detected for ALMA, with a significant number showing round values in frequency, and the majority with short periods, with P ≤ 4 ms.The IRAM 30 m telescope used also a SIS receiver, and similarly to the case of ALMA, we relate the candidates to local oscillations inside the receiver or data transport chains.
Searching the three epochs from the LMT yielded 984 candidates.One candidate stands out after the analysis, with a spin period very close to 32 ms.After a careful examination we conclude that this candidate is related to locally-generated periodic signals or a digitization artifact, because the spin period is almost exactly a thousand times the sampling time of 32 µs.Other reasons to classify this candidate as local are its high power, the fact that the signal suffers at least one significant jump in rotational phase during two of the three observed epochs, and because the same signal is not detected in the ALMA nor the IRAM 30 m datasets.Other strong periodicities at 100, about 180, and 200 Hz are detected in the LMT data, likely related to the power supply.Similarly to the other stations, the data show a significant number of candidates found at round periodicities (both in period or frequency).We conclude that no new pulsars are discovered in the LMT data.
Table 2 provides details of the number of candidates from each epoch, station, and searching pass, and we represent visually the spin period of the candidates in Figure 1.

Fast-Folding-Algorithm with acceleration search
The FFA pipeline applied to the ALMA data resulted in 3044 candidates.The strongest candidate had a period around 18.2 s, which is close to the feature found in the raw ALMA data due to the phasing cycle.The rest  2).The blue circles, green crosses, and orange squares represent the candidates from ALMA, LMT and the IRAM 30 m telescope, respectively.The size of the marker is proportional to the σ sift , and S/N ratio values, for the Fourier-domain and FFA cases, respectively.None of the candidates is a convincing pulsar signal after the inspection of their candidate evaluation plots.The black arrows in the top panel indicate the P = 0.5 ms candidate, the only periodicity present in all datasets and related to a locally-generated signal at ν = 2000 Hz (see Section 3.1.3).
of the candidates have a period close to the harmonics of this signal.This candidate and its harmonics likely originates from residual contamination from the known 18.2-s period phasing features.
The FFA pipeline on the LMT data resulted in 300 candidates.Most of the candidates with S/N > 7 showed repetitions of the pulse in the folded profile, indicating that the fundamental periods of these signals are below our period search range.The candidates also have periods with highly-round numbers, e.g., 1.0000 s, 2.5000 s, and 1.6000 s.We conclude that these candidates have a local origin.
For the IRAM 30 m the FFA search with acceleration yielded over 16000 candidates.This number is much higher than the number of candidates detected for ALMA and the LMT, indicating that our S/N threshold is likely too low resulting in too many false candidates.To limit the number of accepted candidates, the minimum S/N ratio was increased to 7, reducing the candi-dates to 4185 while still maintaining a low significance limit.
A strong candidate at 3.0 s was found without acceleration, with detections of harmonics (i.e., period of 6.000 s and even 150.0s).This candidate's period is close to an integer to the third decimal, suggesting that it is not astrophysical.We confirmed that this candidate likely originates locally by a repeat of the search in topocentric frame, which yields the same candidate with an exact period of 3.000 s.
We conclude that no pulsars were found by the FFA pipeline.Table 2 includes a summary of the candidates, and Figure 1 a visual presentation of the candidates' spin periods.

Comparison of candidates among stations
Apart from an individual review of all candidates, we cross-checked repetitions of periodicities between epochs and stations.For the Fourier-domain pipeline, we use a To match candidates with periodicities that may not coincide exactly (this can occur for example by a slightly different result of optimization of parameters from prepfold), we cross-check the periodicities with precision down to 100 µs.From the 442 reviewed coincidences, none resembles a real pulsar signal; the results are again dominated by local round-number periodicities.
In the few cases where repetitions of periodicities with no round values are found, the characteristics of the candidate (e.g., significance or profile shape) made us classify them as false, produced again likely by local sig-nals or the noise.The most repeated candidate is this step was the P = 0.5 ms, identifiable in Figure 1.This P=0.5 ms candidate is clearly detected in the three stations at all epochs.It is however related to a local signal at a frequency of ν s =2000 Hz.The fact that it exists in all of the datasets reinforces a relationship with the EHT-specific backend (e.g., its clock signal generator, digitizers and/or Mark 6 recorders), that were common hardware at the three telescopes.For the FFA pipeline the periodicities were compared to within 1% of precision, yielding no coincidences between the candidates.

Single pulse search
The single-pulse pipeline resulted in a total of 6893 candidate pulses; the majority of which we determine to be erroneous events.Table 3 shows the number of pulses detected per station epoch, σ SP thresholds and representative flux density limits.The large discrepancy in pulse numbers corresponds to a combination of the non-Gaussian noise properties of each dataset and the different observing lengths.The ALMA data on April 6 and 11 yield the lowest number of single-pulse candidates, however, on April 7 where there are more instrumental instabilities and residual artifacts in the time series (see Section 2.2 and Appendix A), combined with the longer integration length, a significantly higher number of candidates is observed.In general the LMT and the IRAM 30 m datasets also produced larger numbers of candidates for the same reasons and despite the higher σ SP thresholds used.The analysis of simultaneous events among stations yielded 42 single-pulse candidates coincident within time windows of 86.4 ms or less.Further inspection of the time series around the detected pulses showed that these candidates were produced by baseline level jumps or residual time series artifacts.From the visual inspection of all single-pulse candidates with σ SP > 12, no convincing real pulses were found.
Using Equation 1in Karako-Argaman et al. ( 2015), and the system sensitivity parameters outlined in Table 5 we estimate σ SP = 6 or σ SP = 10 on source sensitivity limits to representative 1, 10 and 50 ms duration single pulses.These limits are outlined in Table 3.Several potential factors contribute to the nondetection of pulsars in this search, including the typical steep spectrum of radio pulsars and the large distance to the center of our Galaxy.The sensitivity of a pulsar search is usually quantified by the minimum detectable mean flux density of a pulsar, S min , given a certain set of observational and pulsar properties.Based on the radiometer equation, the theoretical sensitivity of a Fourier-domain pulsar search can be expressed following (Cordes & Chernoff 1997) where β is a correction factor due to imperfections in digitization (β = 1.136 in our case for 2-bit sampling, Cooper 1970), η ph is the phasing efficiency of a phased array (η ph = 1 for single-dish telescopes, and η ph ≃ 0.95 for phased ALMA observations during the EHT 2017 observations, Goddi et al. 2019), σ is the requested detection significance of the pulsar signal, T * sys stands for the system temperature (including the contribution due to Earth's atmopshere), G is the telescope gain, n p is the number of summed polarizations, ∆ν is the effective observing bandwidth, t int is the integration time, α ≡ 1 − π/4 is a coefficient related to the RMS noise in the Fourier transform of the intensity, N h is the optimum number of harmonics summed 18 and R l is the amplitude ratio of the lth harmonic of the signal to the zero frequency in the Fourier transform.Assuming a Gaussian profile shape of the pulsar signal with a width w (the pulse width at half height of the pulse), R l can then be written as where ϵ ≡ w/P is the duty cycle of the pulse.Based on these, the theoretical sensitivities of our searches can be computed using the telescope and data properties summarised in Tables 2, 4 and 5.These results are shown in Figure 2. In practice, Equation 5 is known to overestimate the realistic sensitivity of a pulsar search.This is particularly true for slow pulsars with large duty cycles, when red noise is prominent in the data (Lazarus et al. 2015;Liu et al. 2021;Eatough et al. 2021).Therefore, a more practical accurate estimate is obtained empirically by injecting synthetic pulsar signals each time of a different flux density into the real data, conduct the search and obtain the minimum detectable flux density.Such 18 The N h parameter refers to the number of summed harmonics that minimizes S min .It depends on pulse duty cycle and shape (see Appendix A in Cordes & Chernoff 1997, for details).N h may not coincide with the maximum number of harmonics that the pipeline can sum.The presto pipeline searches for the optimal number of summed harmonics for each detected periodicity up to the limit, which in our case is 32 harmonics.a procedure has been carried out using the longest and most sensitive (see Table 4) dataset from each of the three telescopes (all on Apr.7) and the results are presented in Figure 2. It can be seen that for fast spinning pulsars, the estimate from injection into real data gives a minimum detectable flux density close to the theoretical estimate, compared with those for slow pulsars.As the pulsar period increases, e.g., from 1 ms to 1 s, the empirical sensitivity of the search drops drastically, by up to more than an order of magnitude.
The most sensitive search comes as expected from phased ALMA, where S min ∼ 0.01 mJy.For the singledish telescopes, where the field of view is significantly larger, the best sensitivity is yielded by the IRAM 30 m telescope, where S min ∼ 0.4 mJy.For the LMT, S min ∼ 1 mJy.The GC magnetar PSR J1745−2900 had a flux density of 0.39 mJy at 86 GHz on Apr. 3, 2017; the closest published observation, in both epoch and frequency, to those of the EHT (Liu et al. 2021).Under the assumption of a flat spectrum (e.g., Torne et al. 2015Torne et al. , 2017)), PSR J1745−2900 would fall below the empirically-derived sensitivity limits of IRAM 30 m and LMT and would explain our non-detection in the data from these two telescopes (see Figure 2).Although phased ALMA may have had enough sensitivity to detect PSR J1745−2900 during the EHT 2017 observations, its field of view of 1-2 ′′ was not wide enough to cover this source.
The sensitivity estimates given above are primarily for the Fourier-domain search methods.The sensitivity of the FFA is in theory on a similar level.In Morello et al. (2020), the theoretical sensitivity difference between these two methods is characterized by the search efficiency factor (E), a function of the number of harmonics summed and duty cycle of the pulsar.With a summation of up to 32 harmonics, this factor is approximately 0.65, 0.76, 0.83 for signals with duty cycle of 2%, 5%, 10% (used in our estimates above), respectively.Assuming E FFA = 0.93 as in Morello et al. (2020), the theoretical sensitivity of the FFA search is better than the FFT search by a factor of 1.4, 1.2, 1.1, respectively.
Finally, to illustrate the advantage of our search at very high frequency to overcome the scattering effect, we consider a pessimistic (but still possible) scenario of temporal scattering toward the Galactic Center as predicted by the NE2001 model, i.e., τ s = 2000 ν −4 s (ν in unit of GHz), where τ s is the exponential characteristic time on the pulses due to the scattering, and ν is the observing frequency.Even in this case, at 228 GHz the scattering time would be τ s ≈ 740 ns, roughly 40 times smaller than our smallest sampling interval, three orders of magnitude less than the spin period of any pulsar known, and ∼ 50 times less than the narrowest pulse width known to date (Manchester et al. 2005).We note that the scattering measured for the currently-known closest pulsar to Sgr A*, PSR J1745−2900 (τ s ≃ 1.3ν −3.8 ; Spitler et al.  2014), is much smaller than the predicted value from the NE2001 model.However, we cannot fully discard a patchy environment in the innermost region of the Milky Way, with different scattering for different lineof-sights due to multiple screens (e.g., Schnitzeler et al. 2016;Dexter et al. 2017).In any case, due to the veryhigh observing frequency of the EHT, our search sensitivity is unaffected by interstellar scattering.

Potential of the search to detect a Galactic Center pulsar population
Once the flux density limits of the observations are known, we can estimate the potential to detect pulsars around Sgr A* with this search.We do this by calculating how many pulsars, from a simulated population in the GC, would be detectable with our actual sensitivity.
In the following analysis we focus on the most-sensitive observation from each station as representative of the best limit for GC pulsar detectability, i.e., the data from April 7 (see Table 4).
The flux density limit from Equation 5or through the injection of signals (Section 4.1) can be converted to a pseudo-luminosity limit at the GC by multiplying S min by the square of the distance to the GC (see e. The simulated population in the GC is created from the known pulsar population in the Milky Way, i.e., we assume that the putative pulsars at the center of Table 5. Parameters from each station, corresponding to the most sensitive datasets, used for the sensitivity analysis.Columns indicate the station, observing frequency (ν), average system temperature (T * sys ), telescope gain (G), net integration time onsource (tint), instantaneous bandwidth (∆ν), assumed pulsar duty cycle for the calculations (ϵ), resulting theoretical minimum detectable flux density for σ = 5 (Smin), equivalent pseudo-luminosity at the distance of the Galactic Center (L GC min ), and the percentage of the simulated pulsar population in the Galactic Center that the observations could detect for the theoretically-(ft) and empirically-derived (fe) luminosity limits (see Section 4.2).In all cases the number of summed polarization np = 2.The two last rows show the results from similar recent pulsar searches around Sgr A* at wavelengths of 3.5 mm with ALMA (Liu et al. 2021) and 3 to 2 mm with the IRAM 30 m telescope (Torne et al. 2021) the Galaxy resemble in properties those of the Milky Way.Although the properties of the pulsars that may populate the innermost region of the Galaxy are unknown, we consider this is a valid assumption given that both old and new populations of stars exist in the GC (Pfuhl et al. 2011;Nogueras-Lara et al. 2022).We use the online pulsar catalog PSRCAT 19 (Manchester  et al. 2005) in its version 1.67, extracting the pseudoluminosity of those pulsars for which this value is available (in this catalog version, from 3320 pulsars, 2457, i.e. 74%, have a pseudo-luminosity entry).These pseudoluminosity values are in most cases provided for frequencies of 400 and 1400 MHz.We therefore extrapolate to 228 GHz by using the spectral index of each pulsar that is available in the catalog.When the spectral index is not available, we draw a sample from a normal distribution with mean value ⟨α 1 ⟩ = −1.60 ± 0.54 (Jankowski et al. 2018).Both when the spectral index is known and when we draw from the distribution, a single power law for the extrapolation is used.Finally, we compute the percentage of pulsars with an extrapolated pseudo-luminosity at 228 GHz above the pseudoluminosity limit obtained with Equation 7, both with the theoretically-and empirically-derived S min .Those percentages of the simulated pulsars would be in theory detectable in this search.
Because the exact number of pulsars above the limit will depend on the samples drawn from the normal distribution of spectral indices, we run a Monte Carlo simulation with 10000 executions to obtain the average percentage of detectable population for each station.The summary of the utilized parameters and the results are presented in Table 5. Figure 3 shows a snapshot of the Monte Carlo simulation (representative of the results), together with the empirical luminosity limits for each station.
In addition, to investigate how the spectral index distribution could affect the results, we repeated the analysis in other scenarios for the simulated pulsar population: one with ⟨α 2 ⟩ = −1.4±1.0 (Bates et al. 2013), and with ⟨α 3 ⟩ = −1.8± 0.2 (Maron et al. 2000).Similarly, the impact of the duty cycle parameter was evaluated by repeating the simulation, only for the spectral index distribution ⟨α 1 ⟩, for ϵ = 0.01, and 0.2.
Focusing on the results from the empirical sensitivity, the spectral index distribution ⟨α 2 ⟩ results in a detectable percentage of the simulated population of 8.8%, 1.3%, and 3.4% for ALMA, LMT and the IRAM 30 m telescope, respectively.In the case of using ⟨α 3 ⟩, the percentage decreases to 1.4%, 0.1%, and 0.3%, respectively.These results show that the percentage of detectable population can vary a factor ∼2 to 10 depending on the spectral index distribution utilized; mainly a result of the required extrapolation to 228 GHz and the difference in the standard deviation of the assumed distribution.The simulations show that the mean value itself has a smaller impact on the result than the standard deviation.Even with a steeper mean value, if the standard deviation is sufficiently large, pulsars with flat spectral indices and thus more likely emitting above our sensitivity threshold at 228 GHz, will be drawn from the distribution.In the future, understanding better the distribution of pulsar spectral indices at high radio frequencies (ν ≳ 20 − 30 GHz) would help to obtain more accurate estimations for the potentially-detectable fraction from pulsar searches at very high radio frequencies (see Löhmer et al. 2008).In contrast, the simulations 10 3 10 2 10 1 10 0 10 1 Period (s) In the empirical estimate, each point represents the median from 5 injection iterations.The star in purple stands for the GC magnetar PSR J1745−2900 using the flux density reported by Liu et al. (2021) and assuming a flat spectrum.
Here we used Equation 5(with σ = 5) to calculate the theoretical estimate of minimum detectable flux density.
with additional values of the duty cycle ϵ = 0.01 and 0.2 in the ⟨α 1 ⟩ scenario shows negligible variations on the detectable population.
A first conclusion from the results summarised in Table 5 is that our potential to detect pulsars in the Galactic Center with this search is in general low.This was somewhat expected, primarily due to the effect of the average steep spectral index of pulsars.However, with the observations at 228 GHz we were mainly targeting flat-spectrum pulsars, like radio magnetars, and probing a frequency space where certain emission models include possible emission boosts, making certain pulsars poten-tially visible.Furthermore, the single pulse analysis may detect pulsars through bright single pulses.
A second finding shows that the percentage of detected pulsars can decrease by up to ∼ 50 − 60% when using the empirically-derived sensitivity limits versus the theoretical ones.This sensitivity loss is substantial, and is caused by the noise characteristics, which are far from showing the Gaussian statistics that the radiometer equation assumes.From the three stations, the data from the IRAM 30m seems to be the less affected by this loss of sensitivity.This is apparent in Figure 2, where the empirical limit for this station is flatter, and starts rising at longer spin periods as compared to the ALMA and LMT cases.
The highest fraction of detectable population in the EHT 2017 data comes from ALMA that, despite its superb sensitivity, could detect only around 2.2% of the brightest pulsars in the simulated population with ⟨α 1 ⟩ -that we consider as the most accurate.This fraction is lower by a factor ∼2 compared to the one yielded by ALMA in 2017 Global Millimeter VLBI observations (GMVA), of about 4%, despite the use of only 2 GHz of bandwidth (Liu et al. 2021).In the case of the IRAM 30 m telescope and the LMT, less than 0.4% and 0.1% of the simulated pulsars, respectively, would be detectable with the EHT 2017 data.For the IRAM 30 m, previous searches at 3 and 2 mm yielded theoretical population coverage fractions of 1 − 2% (Torne et al. 2021).The differences arise from a combination of different search parameters (e.g., T sys * , G, t int , or ∆ν), together with the higher average luminosity of pulsars at the longer wavelengths.
Focusing on the EHT 2017 data, the simulations show therefore that only a very small number of very bright pulsars at 228 GHz (1.3 mm), i.e., those that not only are bright but also show a flat (or inverted) spectrum, had chances to be detected in this search.In fact, even bright pulsars at cm wavelengths like Vela, and bright, relatively flat-spectrum pulsars like PSR B0355+54, would not be detected in our survey if located at the Galactic Center.Consequently, it is not surprising that no detections of new pulsars were achieved by our search.Related to our low sensitivity, it is therefore plausible that pulsars emitting below our current luminosity limits exist in the region.This is particular relevant for MSPs, that may be less luminous than normal pulsars (e.g., Bailes et al. 1997;Kramer et al. 1998).From our simulations, our sensitivity to MSPs in this search is close to zero (see Figure 3).Given that the GC may be dominated by an MSP population (Gordon & Macias 2014;Macquart & Kanekar 2015;Schödel et al. 2020), this could also explain the lack of discoveries.We remark that even pulsars emitting above our detection limits and with orbital parameters within the searching pipelines capabilities (see Section 4.3) may be missed by other factors.For example, a pulsar could be highly variable in flux density (as occurs with radioemitting magnetars, see e.g., Torne et al. 2017), or eclipses due to surrounding gas could hide a pulsar emission for some periods of time in binary systems (e.g., Freire 2005).Relativistic precession can move the beam of emission of a pulsar orbiting a companion in and out our line-of-sight (e.g., Kramer 1998;Perera et al. 2010).These reasons alone justify repeated surveys for pulsars around Sgr A* at different epochs separated in time to overcome these time-dependent effects.
Lastly, another case potentially affecting our sensitivity is when a pulsar's spin and orbital parameters are beyond the capabilities of our pipelines to recover the accelerated signals.This mostly affects fast MSPs (P ≲ 5 ms) in tight orbits around massive companions.Nevertheless, the possibility that the center of the Galaxy is populated to a large extent by MSPs cannot be ruled out (e.g., Rajwade et al. 2017).Since the main signals aimed for here are pulsars orbiting the super massive black hole Sgr A*, we should consider our sensitivity to such MSP-(SM)BH systems potentially lower than the limits shown in Table 5 and Figures 2 and 3.In the next Section, we discuss this relationship between sensitivity and binary systems in more detail.

Search capability for pulsar in various orbits
around Sgr A* The phased ALMA field of view is significantly smaller than LMT and IRAM 30 m, sampling a smaller volume and range of potential orbits near Sgr A*.Using the distance and mass of Sgr A* derived from Gravity Collaboration et al. (2021) and assuming a circular orbit, the phased ALMA field of view is able to cover an orbital period of approximately up to P b ≈ 300 yr.This is already enough to cover a large fraction of the stars known to date in the S-star cluster (Sabha et al. 2012;Gillessen et al. 2017), and certainly also the orbits that would enable Sgr A* measurements and gravity tests with pulsars (Liu et al. 2012).Meanwhile, the field of view of LMT and IRAM 30 m can cover orbital periods of up to P b ≈ 4000 yr.So, even though their observations are less sensitive, the LMT and IRAM 30 m telescope are able to cover a significantly larger volume of any putative pulsar population in the GC.
In practice, as already mentioned in Section 2.3, the range of acceleration and its derivative explored in the search would also place a constraint on the types of orbits around Sgr A* that the search is able to detect.For a given time span of the observation, the constraint can be estimated by calculating the maximum absolute z and w values in a range of orbits and compare with those used in the search, as demonstrated in Liu et al. (2021) (see also Eatough et al. 2021).Here as case studies, we carried out the same practice for two time spans of the observation, T obs = 10.2 and 3.2 hr, which correspond to the April-7 observation of ALMA and IRAM 30 m, respectively, being the longest and shortest time span in our dataset.We focus on a full sensitivity signal recovery based on 32 summed harmonics (see Section 2.3.1).The results are summarised in Figure 4, with a selection of two typical spin periods for ordinary pulsars (P ≃ 500 ms) and MSPs (P ≃ 5 ms), respectively.
With the settings detailed in Section 2.3.1, our search using 10.2-hr long data is capable of detecting ordinary pulsars well within one-year orbits around Sgr A*, likely down to P b ≳ 0.5 yr.For MSPs, we are able to detect those with an orbital period longer than 10 yr.Meanwhile, the search using 3.2-hr long data can detect ordinary pulsars in much closer orbits, i.e., down to 0.2 yr (≈ 2.4 months).For MSPs, it can detect orbits with as short as approximately 2-yr periods.These are also inferred by Equation 3 and 4, where the parameter space in acceleration and jerk to be explored scales down for short observations and slow pulsars.Note that as shown in Table 1, the ALMA observations from April 6 and 11 have similar time span and thus coverage of pulsar orbits as the IRAM 30 m observation.
As discussed in Liu et al. (2021), the above limit estimate on the orbits corresponds to the boundary where the optimal sensitivity can always be retained throughout the entire orbit.If in practice the search conducted does not cover the recorded max(|z|) and max(|w|), it would start to lose sensitivity in the worse scenario but could still detect the pulsar.For instance, for highly eccentric orbits, the pulsar spends most of the time near the apoastron, where the acceleration and its derivative are considerably lower than those near periastron.In this case, searches that do not cover max(|z|) and max(|w|) of the entire orbit, could still be able to find the pulsar at a large fraction of the orbital phase.

Constraints on a putative Galactic Center pulsar population
Theoretical estimations on the number of pulsars that may exist in the inner part of our Galaxy differ considerably, with typical assumed values from around ∼200 (Chennamangalam & Lorimer 2014) to thousands (Wharton et al. 2012).The wide range relates mainly to the complexity and peculiarity of the star forming history and evolution in the Galactic Center, and the dynamics that the high-density environment and the SMBH introduce (see e.g., Figer 2009;Morris 2023).Using our survey results, we can try to constrain the number of pulsars that may exist in the Galactic Center.
We will assume a population of pulsars located at the distance of Sgr A* .For simplicity in the calculations, the population will have a Gaussian radial distribution, with two parameters: n 0 as the number density of pulsars, and σ R the RMS radial scale, for a 3D Gaussian distribution.Following a similar approach as in Cordes & Chernoff (1997), integrating over the angles sampled by each telescope we calculate the mean expected number of pulsar detections from this assumed population (N d ), for a given telescope, as where f s is the fraction of detectable pulsars with a given telescope (we use the empirical detection fraction f e for our calculations, see Table 5), d gc is the distance to the Galactic Center, and σ θ is the RMS angular scale of a Gaussian beam pattern for the telescope.
To constrain the parameters of the population, we assume that our probability for detections follows a Poisson distribution.Given our non-detections, we build our survey likelihood function as L 0 = P 0 = e −N d .We note that the two parameters for the assumed GC pulsar population are not statistically independent, so the upper limits for n 0 will depend on σ R .We calculate an upper limit for N d (and so, also for n 0 ) at a given confidence level, that we choose here as 99.7%.In order to obtain the most stringent constraint on the expected number of pulsars, we select the lowest n 0 for each σ R from the values obtained from the three used telescopes.Finally, we compute an upper limit for the expected mean number of total pulsars in the assumed population as N GC psr = n 0 V pop , where V pop = (2π) 3/2 σ 3 R is the volume enclosing the assumed population.The results are plotted in Figure 5.Although we have not detected any pulsar in the EHT 2017 observations, we know that at least one pulsar exist near Sgr A* : the magnetar PSR J1745−2900, just three arcsec away (i.e., a projected distance of ≈ 0.1 pc, Rea et al. 2013).This pulsar was not detected because it is highly variable and was likely on a weak-emission state during the EHT 2017 campaign.PSR J1745−2900 would have been detected by the IRAM 30m telescope, and even phased ALMA (even offset from the center of its synthetic beam), if it were in a bright-emission episode (PSR J1745−2900 can reach S ∼ 6 mJy at 1.3 mm, Torne et al. 2017).Taking into account an alternative scenario with a detection of PSR J1745−2900, we also present the upper limits considering our survey likelihood function L 1 = P 1 = N d e −N d in Figure 5.The derived upper limits show that, given the assumptions made on pulsar population properties (i.e., similar in properties to the known pulsar population, and with the chosen spectral index distribution ⟨α 1 ⟩ = −1.60 ± 0.54), and the characteristics of the telescopes used, of the order of 1000s of pulsars can populate the GC if distributed in a region of width ∼1 pc or more around Sgr A*.For instance, if the pulsars are assumed to exist within the inner central pc (σ R = 1 pc), N GC psr ≲ 1600 (N GC psr ≲ 5300 in the case of using L 1 ).A scenario of particular interest is when the assumed pulsar population in the GC is very compact around Sgr A*, within a volume containing (circular) pulsar orbital peri-ods P b ≲ 100 yr around Sgr A* .This smaller volume encloses the most promising pulsars for gravity tests with Sgr A* (Wex & Kopeikin 1999;Pfahl & Loeb 2004;Liu 2012), and is given by a σ R ≈ 0.02 pc, yielding N GC psr ≲ 50 (N GC psr ≲ 160 in the case of using L 1 ).Our results therefore cannot set stringent constraints on the number of pulsars in the inner parsec of the galaxy, mainly due to our low sensitivity at large radial scales (constrained by the wider beams of IRAM 30m and LMT).However, on smaller radial scales (σ R ≲ 0.1 pc), phased ALMA offers better sensitivity.In these smaller scales, we are still compatible with a population of tens to hundreds of pulsars close to the central SMBH and beamed towards us.
We note however that N GC psr is subject to big uncertainties.One of the largest ones arise from assumption made on the properties of the pulsars that may exist in the GC, from which the detection fraction, f s , is derived.To illustrate the impact of just a different mean spectral index taken for the GC pulsar population, the upper limit for N GC psr for a population within 1 pc around Sgr A* varies between N GC psr ≲ 190 and N GC psr ≲ 2100 (N GC psr ≲ 620 and N GC psr ≲ 7100 if L 1 ) for the spectral index distributions ⟨α 2 ⟩ = −1.4± 1.0 and ⟨α 3 ⟩ = −1.8±0.2, respectively.In the case of a compact (Bottom panels:) Upper limit on the number of expected pulsars around Sgr A* within a spherical Gaussian volume with radial width σR.The top axis indicates the Gaussian radial width in arcsec scale, for easier comparison with the telescope beam sizes, represented by the vertical dotted lines in blue, green, and red color, at ≈0.06, ≈0.4 and ≈0.43 pc for phased ALMA, LMT, and the IRAM 30m, respectively.The survey likelihood function is L0 (left) and L1 (right).In all the panels the continuous thick, dotted-dashed and dashed lines marks the upper limits for a pulsar population with a mean spectral index ⟨α1⟩, ⟨α2⟩ and ⟨α3⟩, respectively.population enclosed within 0.02 pc around Sgr A* , N GC psr varies between N GC psr ≲ 10 and N GC psr ≲ 75 (N GC psr ≲ 40 and N GC psr ≲ 250 if L 1 ), again for ⟨α 2 ⟩ and ⟨α 3 ⟩, respectively.The upper limits on N GC psr for the scenarios with mean spectral indices ⟨α 2 ⟩ and ⟨α 3 ⟩, for comparison, are also shown in Figure 5.
Finally, we remark that the upper limit N GC psr refers to pulsars active in radio and beamed towards us.The limits for total number of neutron stars in the assumed regions would be considerably greater depending on the assumed beaming fraction and ratio of radio-active neutron stars over all the existing ones.For instance, for a mean beaming fraction of 0.1 (Faucher-Giguère & Kaspi 2006) and a ratio of radio-active over all neutrons stars of 10 −2 -10 −3 (Faucher-Giguère & Kaspi 2006; Keane & Kramer 2008;Sartore et al. 2010), the upper limits for neutron stars would scale up by a factor O(10 3 -10 4 ).Therefore, despite our non-detections, our results are consistent with a large population of neutron stars in the Galactic Center; and compatible with results from theoretical simulations and population synthesis (Alexander 2017; Baumgardt et al. 2018;Gravity Collaboration et al. 2020;Chen et al. 2023).

Flux density sensitivity
The pulsar search in this paper used EHT 2017 campaign data with a total bandwidth of ∆ν ≈ 4 GHz.Later EHT observations, starting from 2018, began to record an instantaneous bandwidth of ∆ν = 8 GHz, offering a factor of √ 2 improvement in system sensitivity.Additionally, for single-dish telescopes, observing in the VLBI mode requires a notably larger overhead compared to those in stand-alone mode, as a result of the time spent on phase calibrators and the associated slewing time (see Table 1).For the same observing length, dedicated pulsar surveys would largely increase the fraction of time spent on Sgr A*, by possibly more than a factor of 2; this would improve the system sensitivity by roughly another factor of 1.5.
Furthermore, in this work we searched the dataset individually which were in fact collected from simultaneous observations at different telescopes.In principle, it is possible to increase the search sensitivity by coherently summing the data collected with the entire EHT array, forming tied array beams towards areas in the vicinity of Sgr A* (see e.g., Bassa et al. 2016).This may also effectively mitigate some of the systematics present in individual telescope data because they should not correlate between stations.A downside is that a large number of beams needs to be formed, which requires substantial computing power.For instance, given the longest baseline of the EHT 2017 observation, an order of 20 million beams would be needed to cover an area of 0.1 arcsecond around Sgr A*.These beams will also have to be contained by the smallest beam in the array, which is the synthetic beam of ALMA for the EHT 2017 campaign.This could be improved if multi-beaming capability of ALMA becomes available.
To improve the feasibility of coherent beam forming of the EHT array, the beams can be directed towards compact objects that may be identified by imaging observations.Alternatively, incoherently adding the data and/or coherently summing only a subset of the array, may also be considered as an option to increase the system sensitivity.Table 6 summarizes the delivered equivalent diameters of the array under different ways of summing the telescope signals.For coherent addition, the summation of the entire array will deliver a sensitivity equivalent to a 111-m dish, while summing only the most sensitive four telescopes in the array will still yield an equivalent diameter of 107 m, only 7% less in sensitivity.
The sky area where a coherent beam may be formed is restricted by the smallest beam size of the array, which is from one of the interferometers.If summing only the single-dish telescopes, the sky coverage can then be enlarged to the smallest beam size of the single dishes.This will give an equivalent diameter of 64 m, approximately 60% more sensitive than LMT, the largest single dish in the array.For incoherent summation, adding the data from the entire array will give an equivalent diameter of 81 m.However, only a few percent difference is expected from summing only the top three sensitive telescopes, i.e., ALMA, LMT and the NOrthern Extended Millimeter Array (NOEMA)20 .Similarly, incoherent summation of all single-dish telescopes in the array, will deliver an equivalent diameter of 52 m, only 8% improvement in sensitivity from the LMT.This is largely due to the significant difference in gain among the telescopes.
Finally, the sensitivity of the LMT during the EHT 2017 campaign is far from its full capacity, mainly due to the under-illumination of its surface.It is expected that the LMT data from later EHT observations will deliver much higher sensitivity.

Orbital parameters and searching algorithms
As discussed in Section 4.3, the acceleration search performed in this work has optimal sensitivity for pulsars in a certain range of orbits.To cover more compact orbits with Sgr A* optimally, or where a fast MSP is involved, it is possible to simply increase the range of acceleration and jerk used in the search (at the expense of higher computational costs).
Alternatively, new searching schemes could also be introduced to directly explore the space of Keplerian parameters (e.g., Balakrishnan et al. 2022), delivering optimal sensitivity to more types of orbits.Nonetheless, the challenge for these search strategies is the computing power required, considering that pulsar orbits in the high stellar density region of the GC are expected to be significantly eccentric.As such, the full set of five parameters is required to describe the orbit rather than only three for circular orbits.
Another possible approach to enable optimal sensitivity for more compact orbits, is to divide up the entire length of data into segments and to search each segment individually (Ng et al. 2015;Liu et al. 2021;Eatough et al. 2021).Admittedly, this approach suffers from a loss in sensitivity due to shorter integration times in each segment, but it is more computationally feasible.

Data cleaning and quality
The impact of red noise is a major source of sensitivity loss in pulsar searches, particularly detrimental in high radio frequencies, where the instrumentation and opacity variations produce an excess of power in the low-frequency end of the Fourier spectrum.This effect combines with the fact that radio magnetars, typically found slowly spinning, are one of the main targets of pulsar surveys in the short millimeter range, due to their typically-flat spectrum.An optimal red noise reduction scheme can improve the sensitivity of our surveys, although the exact filtering and parameters required are non-trivial and depend on the data properties and the specific properties of the pulsars to be detected (Singh et al. 2022).Future surveys may need different passes with a suite of red noise filtering techniques to optimize the results.
Extra improvements in sensitivity to pulsars could potentially come from the reduction of the locallygenerated signals that exist in the data.The SIS receivers themselves and the downconversion systems in the (sub)millimeter telescopes are complex, with many oscillators and mixers, and therefore can introduce external signals into the final data recorded to disk.These local signals are often highly periodic.This is generally not an issue for spectroscopic or continuum observations, where the scans are shorter in integration length and the data are averaged over a comparatively longer time as compared to the data analyzed here for the pul-sar search21 .Nonetheless, in the case of using the data for a search for pulsars, such local periodic signals can have a substantial impact on the sensitivity.
The detrimental effect comes mainly from large excesses of power in many Fourier bins, which then interfere with the searching algorithms by the extra candidates found.Although local interfering signals (both internal and external to the telescope system) are common in radio observations, they can be particularly harmful at very high radio frequencies, where the dispersion effect is very small or negligible22 .
With risk of having still a strong scattering affecting the centimeter-wavelength observations, and with the large impact on sensitivity by the steep spectrum of pulsars at the short millimeter wavelengths, attempting new surveys with high sensitivity in the short centimeter and long millimeter-wavelength regime is compelling.Such frequency range (around ∼ 10 − 50 GHz) would still diminish the scattering and dispersive effects of the ISM, while observing in a spectral window with stronger emission from the steep-spectrum pulsars.This range may contain the frequency "sweet spot" between scattering mitigation and sensitivity to pulsars near Sgr A* (Macquart et al. 2010;Macquart & Kanekar 2015;Bower et al. 2018).

SUMMARY AND CONCLUSIONS
We carried out a search for pulsars and fast transients in the Galactic Center with the Event Horizon Telescope 2017 observations, the first search at an observing wavelength of λ = 1.3 mm (ν = 228 GHz).The search used data from phased ALMA, LMT, and the IRAM 30 m telescope, the most sensitive telescopes in the EHT array of 2017.Periodicity searches both in the Fourier domain and with a Fast-Folding-Algorithm were conducted to each individual dataset collected on three different nights, with acceleration and jerk search incorporated to cope with the potential orbital motion of pulsars around Sgr A* or a binary companion.Singlepulse searches were also performed to the whole dataset.
These searches did not detect any pulsars or transients.We estimated the search sensitivity both theoretically and in practice by injecting artificial signals into the real data.The practical sensitivities for fastspinning pulsars (P = 1 − 10 ms) are approximately 0.02, 0.4, 1 mJy for ALMA, LMT and the IRAM 30 m telescope, respectively, but are roughly an order of magnitude worse for slow pulsars (P = 1 − 10 s).In addition, we showed the possible pulsar orbits that can be detected with our searching scheme.We explored the detectability of the search towards a simulated pulsar population in the Galactic Center and around Sgr A*, concluding that the sensitivity of the observations is still low.For MSPs in particular, which may be a dominant population in the region, the search sensitivity is close to zero.The lack of discoveries is therefore not indicative that pulsars do not exist in the region; there could still be pulsars in the covered areas that simply emit below our detection thresholds.Finally, we discussed future improvements that can be introduced to optimize the usage of these and similar data, to improve the system sensitivity and, overall, the chances for the detection of pulsars that likely still hide in the vicinity of the supermassive black hole Sgr A*.The top two panels corresponds to data from phased-ALMA, the middle two panels to the LMT, and the bottom two panels to IRAM 30 m telescope data.Each pair of panels corresponding to one station show: (Top) the raw total intensity data just after the conversion to PSRFITS and (Bottom) the same total intensity time series after the preparation and cleaning described in Section 2.2.A description of the main features of the data from each station is presented in Section 2.2.The rightmost panels present a histogram of the distribution of the samples in linear scale.After the cleaning, the histograms tend to a Gaussian shape (although not perfect), indicating that the filtering and flagging schemes result in statistically better-behaved data.The top two panels corresponds to data from phased-ALMA, the middle two panels to the LMT, and the bottom two panels to IRAM 30 m telescope data.Each pair of panels corresponding to one station show: (Top) the raw total intensity data just after the conversion to PSRFITS and (Bottom) the same total intensity time series after the preparation and cleaning described in Section 2.2.A description of the main features of the data from each station is presented in Section 2.2.The panels on the right show an even closer zoom in, to facilitate the identification of the main issues of each dataset.Note that for the ALMA and LMT this extra zoom shows one second, while for the IRAM 30 m data, which suffer of much faster oscillations, the time axis of the right-hand panel encompasses only 0.1 s.The black solid line shows a running mean with resolution of 10 ms.

Figure 1 .
Figure 1.Spin period of the candidates from the two periodicity-search pipelines: Fourier-domain analysis (top panel) and FFA (bottom panel).Each epoch shows two lines of markers, corresponding to the two data processing passes with each pipeline (labeled as FFT1, FFT2, FFA1 and FFA2 in Table2).The blue circles, green crosses, and orange squares represent the candidates from ALMA, LMT and the IRAM 30 m telescope, respectively.The size of the marker is proportional to the σ sift , and S/N ratio values, for the Fourier-domain and FFA cases, respectively.None of the candidates is a convincing pulsar signal after the inspection of their candidate evaluation plots.The black arrows in the top panel indicate the P = 0.5 ms candidate, the only periodicity present in all datasets and related to a locally-generated signal at ν = 2000 Hz (seeSection 3.1.3).

Figure 3 .
Figure 3. (Left:)  Luminosity for the simulated pulsar population at 228 GHz (black dots) with the empirically-derived minimum detectable luminosity for each telescope at the distance to the Galactic Center (d = 8.28 kpc) over-plotted as lines.The minimum detectable luminosity corresponds to the empirical limits on flux density as seen in Figure2for a pulse duty cycle of 5%.The blue line, dashed-dotted green line, and dashed orange line represents the limits for ALMA, LMT and the IRAM 30 m telescope, respectively.Each black dot represents a simulated pulsar, and the dots above each line would in theory be detectable for that telescope.For reference, PSR J1745−2900's average luminosity at 228 GHz(Torne et al. 2015(Torne et al. , 2017)), the Vela pulsar (PSR B0833−45) and PSR B0355+54 are marked with a yellow star, brown triangle and pink diamond symbols, respectively.The vertical line over the yellow star marks the range of known variability of PSR J1745−2900.(Middle:) Histogram of the luminosity distribution, and (Right:) histogram of the luminosity distribution in a cumulative layout.

PFigure 4 .
Figure4.Values of zmax (upper row) and wmax (lower row) required for the Fourier-domain search to be conducted with optimal sensitivity in orbits around Sgr A* with different periods and eccentricity (e).Here we use two pulsar spin periods, which are typical values for ordinary pulsars and MSPs, respectively.The time spans of the observations are assumed to be 10.2 and 3.2 hr, which are in turn the longest and shortest in our analyzed datasets.

Figure 5 .
Figure 5. (Top panels:) Upper limits on pulsar number density (n0) versus Gaussian radial width (σR) constrained by the EHT 2017 observations for a survey likelihood L0 = P0 (left) and L1 = P1 (right).(Bottompanels:) Upper limit on the number of expected pulsars around Sgr A* within a spherical Gaussian volume with radial width σR.The top axis indicates the Gaussian radial width in arcsec scale, for easier comparison with the telescope beam sizes, represented by the vertical dotted lines in blue, green, and red color, at ≈0.06, ≈0.4 and ≈0.43 pc for phased ALMA, LMT, and the IRAM 30m, respectively.The survey likelihood function is L0 (left) and L1 (right).In all the panels the continuous thick, dotted-dashed and dashed lines marks the upper limits for a pulsar population with a mean spectral index ⟨α1⟩, ⟨α2⟩ and ⟨α3⟩, respectively.

Figure 6 .
Figure 6.Example time series of the observations of Sgr A* analyzed in this paper.The top two panels corresponds to data from phased-ALMA, the middle two panels to the LMT, and the bottom two panels to IRAM 30 m telescope data.Each pair of panels corresponding to one station show: (Top) the raw total intensity data just after the conversion to PSRFITS and (Bottom) the same total intensity time series after the preparation and cleaning described in Section 2.2.A description of the main features of the data from each station is presented in Section 2.2.The rightmost panels present a histogram of the distribution of the samples in linear scale.After the cleaning, the histograms tend to a Gaussian shape (although not perfect), indicating that the filtering and flagging schemes result in statistically better-behaved data.

Figure 7 .
Figure 7. Zoom to the example time series of the observations of Sgr A* analyzed in this paper.The top two panels corresponds to data from phased-ALMA, the middle two panels to the LMT, and the bottom two panels to IRAM 30 m telescope data.Each pair of panels corresponding to one station show: (Top) the raw total intensity data just after the conversion to PSRFITS and (Bottom) the same total intensity time series after the preparation and cleaning described in Section 2.2.A description of the main features of the data from each station is presented in Section 2.2.The panels on the right show an even closer zoom in, to facilitate the identification of the main issues of each dataset.Note that for the ALMA and LMT this extra zoom shows one second, while for the IRAM 30 m data, which suffer of much faster oscillations, the time axis of the right-hand panel encompasses only 0.1 s.The black solid line shows a running mean with resolution of 10 ms.
Though the full geometric diameter of LMT is 50 m, during the EHT 2017 campaign only 32.5 m were illuminated.

Table 2 .
Summary of the parameters of the two periodicity-search pipelines.Columns show the station, epoch, pipeline (with subscript indicating different passes with different parameters), parameters used in each case (see Sections 2.3.1 and 2.3.2 for details), and number of pulsar candidates produced in each search pass.A "−" symbol indicates that for that pipeline the parameter is not applicable.

Table 3 .
Parameters and results of the search for single pulses and transients.Columns show the dataset analysed (date and station), the signal-to-noise threshold used for single pulse events (σ ′ SP ), the number of single pulses detected (N SP cand ), and flux density limits (Smin) for single pulses of representative widths 1, 10 and 50 ms.

Table 4 .
Calibration information per station and epoch.Remaining columns show the telescope gain (G), the average effective system temperature (that includes the contribution due to Earth's atmosphere, T * sys ), the net integration time on Sgr A* (tint) and a figure-of-merit of the sensitivity of each dataset, Γ = 1000 • T * sys /(η ph G , for comparison.The "n/a" abbreviation indicates that a result is not available.

Table 6 .
The equivalent diameter delivered by a variety of summation forms of the EHT array.The equivalent diameters of each individual telescope are obtained from Event Horizon TelescopeCollaboration et al. (2019d).