TRAPUM search for pulsars in supernova remnants and pulsar wind nebulae − I. Survey description and initial discoveries

We present the description and initial results of the TRAPUM (TRAnsients And PUlsars with MeerKAT) search for pulsars associated with supernova remnants (SNRs), pulsar wind nebulae and unidentified TeV emission. The list of sources to be targeted includes a large number of well-known candidate pulsar locations but also new candidate SNRs identified using a range of criteria. Using the 64-dish MeerKAT radio telescope, we use an interferometric beamforming technique to tile the potential pulsar locations with coherent beams which we search for radio pulsations, above a signal-to-noise of 9, down to an average flux density upper limit of 30 μ Jy. This limit is target-dependent due to the contribution of the sky and nebula to the system temperature. Coherent beams are arranged to overlap at their 50 per cent power radius, so the sensitivity to pulsars is not degraded by more than this amount, though realistically averages around 65 per cent if every location in the beam is considered. We report the discovery of two new pulsars; PSR J1831 − 0941 is an adolescent pulsar likely to be the plerionic engine of the candidate PWN G20.0+0.0, and PSR J1818 − 1502 appears to be an old and faint pulsar that we serendipitously discovered near the centre of a SNR already hosting a compact central object. The survey holds importance for better understanding of neutron star birth rates and the energetics of young pulsars.


INTRODUCTION
Pulsars are rapidly spinning, highly magnetised neutron stars (NSs), detected if the radiation that is emitted from their magnetic poles sweeps across the line of sight of an observer.NSs are predominantly thought to be born from the core-collapse supernova explosions of fuel-exhausted early spectral type main sequence stars Baade & Zwicky (1934).Pulsars emit spin-powered radiation for far longer than their supernova remnants (SNRs) are visible (e.g.Braun et al. 1989), therefore identifiable supernova sites offer the best prospects for uncovering young pulsars.
SNRs are the expanding front of ejecta and swept up interstellar medium (ISM) resulting from supernovae (SNe).It is now understood that about 75 per cent of SNe are Type Ib/c or Type II (Cappellaro ★ E-mail: james.turner-13@postgrad.manchester.ac.uk et al. 1999), all of which are triggered by the core-collapse of stars of initial mass ≳ 9 M ⊙ (Heger et al. 2003) and form a NS.The remainder are Type Ia supernovae, which are the result of either a thermonuclear deflagration (Nomoto et al. 1984) of an accreting white dwarf star at the Chandrasekhar mass limit of 1.4 M ⊙ (Chandrasekhar 1931) or a collision between two white dwarfs (Raskin et al. 2009;Rosswog et al. 2009).Type Ia SNe do not produce a neutron star, therefore an effective pulsar search strategy should exclude SNRs that show evidence of a thermonuclear origin, which are primarily distinguished from core-collapsed remnants by a strong presence of iron in the optical spectra (see e.g.Vink 2020, for a review).It is possible that the most massive 1 per cent of massive stars undergo so-called pair-instability SNe which does not produce a neutron star but does result in a SNR (Heger et al. 2003).
There are 303 confirmed SNRs listed in the 2022 version of the Galactic SNR Catalogue 1 (Green 2022, henceforth referred to as G22).Radio continuum imaging surveys such as those with the Molonglo Observatory Synthesis Telescope (MOST, Clark et al. 1975;Gray 1994;Whiteoak 1996) or the Multi-Array Galactic Plane Imaging Survey (MAGPIS, Helfand et al. 2006;Brogan et al. 2006) with the Very Large Array (VLA) have reported regular batches of new sources.New instruments and surveys with a high sensitivity and resolution, even at lower frequencies, have identified a plethora of new candidates, particularly small diameter shells of low surface brightness on the Galactic plane.Nearly 200 of these are from three surveys; the 1-2 GHz HI, OH, Recombination line survey of the Milky Way (THOR, Anderson et al. 2017), the GaLactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM, Hurley-Walker et al. 2019), and the GLOSTAR (GLObal view of STAR formation in the Milky Way) survey (Dokara et al. 2021).The success of these surveys is demonstrated in Figure 1, which shows an accelerating discovery rate for new SNRs.Based on an all-type SN rate of one event per 40 yr (Tammann et al. 1994) and a dissipation timescale of approximately 60 kyr (Frail et al. 1994), there should be around 1500 visible radio SNRs in the Galaxy.Figure 1 includes 200 new candidates identified by MeerKAT and reported by Anderson (2023), though we do not include these unpublished targets in the survey reported here.
Most SNRs are shell-type at radio and X-ray wavelengths (Green 2019), where only emission from the shock front is visible.A smaller number of shells encompass non-plerionic radio and non-thermal Xray features of increasing brightness towards the remnant centre, and are wholly referred to as filled centre SNRs.These are not to be confused with mixed-morphology SNRs, which exhibit central soft X-ray emission in an apparent delineation from the Sedov-Taylor description of expansion (Sedov 1946;Taylor 1950) as applied to SNR evolution, due to a dense but uniform ISM (Rho & Petre 1998).If a shell encloses a plerion then they are conjointly described as a composite-type SNR.A plerion, or pulsar wind nebula (PWN), is formed when the magnetised pulsar wind is confined behind a termination shock, where the wind and pressures equalise (e.g.Torres et al. 2013).The pressure confining the wind can be due to the over-dense SNR interior medium, or if the pulsar is moving supersonically.The post-shock particle radiation is dominated by synchrotron emission at radio wavelengths and X-rays, and includes γ-rays from 1 http://www.mrao.cam.ac.uk/surveys/snrs/ inverse Compton scattering of low energy photons (see e.g.Slane 2017;Mitchell & Gelfand 2022, and references therein).Supersonic pulsar motion can be supplied by a kick from an asymmetric core collapse (Podsiadlowski et al. 2005).The kick velocity, which averages somewhere between 250 and 440 km s −1 (Hansen & Phinney 1997;Hobbs et al. 2005), is occasionally sufficient for a NS to depart a still-visible SNR, e.g. the 'Frying Pan' (Camilo et al. 2009b), the 'Mouse' (Camilo et al. 2002) and the 'Mini Mouse' (Motta et al. 2023).The bow-shock directs the energetic particles behind the pulsar, forming a tail (Kargaltsev et al. 2017).
PWNe are sources of very high energy (VHE) γ-rays.This emission is well described by models of energetic leptons injected by the pulsar wind (see e.g.Slane 2017, and references therein), and can extend into the PeV regime γ-rays (Cao et al. 2021).Due to the much slower cooling timescale than for pure synchrotron X-rays, the TeV component is more representative of the cumulative energy supplied by the spin-down luminosity of the pulsar over its lifetime, , i.e.
∫  0  ()  (Kargaltsev et al. 2013), whereas the X-ray power is some fraction of the instantaneous output, .The High Energy Stereoscopic System (H.E.S.S.), located in Namibia, is an array of Čerenkov telescopes detecting γ-rays from 0.1 to 100 TeV (Bolmont et al. 2014).A full survey of the Galactic plane has identified a total of 95 VHE sources (Abdalla et al. 2018a), the majority of which are firmly or potentially associated with known PWNe (Abdalla et al. 2018b).Sources of non-plerionic origin include hadronic VHE emission due to SN shocks colliding with molecular clouds, or from cataclysmic binary systems.According to the catalogue of TeV sources tevcat2 , 25 TeV sources seen by H.E.S.S. remain unidentified.Aside from PWNe, H.E.S.S. has also helped uncover an interesting new population of TeV halos, regions of VHE γ-rays extending over a degree away from a pulsar.It is possible this emission is from a relic wind nebula where the energetic leptons would, in the case of a pulsar that has left the SNR shell (López-Coto et al. 2022), not be bound within the magnetised wind bubble as is expected for young PWNe (Gaensler & Slane 2006).Thus TeV halos could be the final stage in the life cycle of PWNe.As many as 100 halos could be detectable in the near future (Sudoh et al. 2019;Martin et al. 2022).TeV sources, therefore, are useful indicators for the presence of undetected young or adolescent pulsars.
All the sources discussed here constitute a zoo of targets that could contain many undiscovered young pulsars.The discovery rate for PSR-SNR systems has slowed in recent years; 63 out of the 303 G22 SNRs have evidence for a neutron star counterpart.It should be noted that already 20 years ago there were 51 associations (Kaspi & Helfand 2002).Additionally, 63/303 is consistent with the fraction of pulsars estimated to be visible from Earth assuming a beaming fraction between 1/6 to 1/5 (Taylor & Manchester 1977;Biggs & Lyne 1996).Both points suggest that a ceiling on finding new pulsars in known SNRs is approaching.Nevertheless, the jump in new SNR candidates (Figure 1), new VHE sources and the high sensitivity of current generation telescopes like MeerKAT motivate a new, large survey to discover many more young pulsars.The establishment of a more statistically significant sample is necessary for deepening our understanding of the pulsar population, the energy budgets of PWN systems, the coupling of pulsars to the evolution of their remnants, and constraining the magnitude of effects that prohibit detections of some young pulsars.
Many SNRs and PWNe are located close to, or on the Galactic plane, and therefore are covered by the footprint of many non-targeted pulsar surveys.Targeted searches are usually able to go much deeper by centring the beam pattern directly on the remnants and by integrating longer than the non-targeted survey dwell times.However, of the 40 pulsars associated with SNRs or PWNe for which radio pulsations were detected first, a minority of 12 were discovered in targeted searches while 28 were found by non-targeted surveys (Manchester et al. 2005).On several occasions, searches with sensitivities of order 1 mJy of many SNRs have produced no new discoveries; Kaspi et al. (1996) searched 39 3 Southern SNRs using Murriyang, the 64-m radio telescope at Parkes, Gorham et al. (1996) searched 17 remnants with the 305 m Arecibo telescope, Sett et al. (2021) focused on 5 remnants with the Green Bank telescope and searched deeper than any previous attempts, but did not see anything new, and Straal & van Leeuwen (2019) were fruitless in their different approach of using LOFAR to detect steep-spectrum pulsars.In a search of 33 Northern remnants with the Lovell telescope, Lorimer et al. (1998) found only one pulsar that is likely to be associated with the target (Tian & Leahy 2006).
It might be argued that it is preferable to focus on targeting pulsar tracers such as compact VHE and radio emission from PWNe, as these searches have often succeeded in detecting new radio pulsars (e.g.Wolszczan et al. 1991;Camilo et al. 2009a).The difficulty in explaining the non-detection of a radio pulsar in any given SNR is closely linked to an incomplete understanding of selection effects.For example, a radio pulsar may escape detection if the pulses are substantially affected by the circumstellar medium and the ISM.Frequency-dependent dispersion due to the cold plasma, which delays radio waves proportionally to  −2 , risks smearing pulsar signals within frequency channels beyond recovery.Scattering, which is caused by a turbulent medium, can broaden a pulse and spread the signal too thinly, masking it in a similar fashion to dispersion but with a weaker frequency dependence of  −4 .Alternatively, a pulsar may also not be intrinsically luminous enough to surpass the signalto-noise threshold.These explanations for a non-detections are not mutually exclusive, nor is the size of their effect on pulsar discovery truly constrained.
New telescopes like MeerKAT provide increased sensitivity to dimmer radio pulsars.Taking the (pseudo-) luminosity,   =    2 , where   is the pulsar flux density at a frequency, , the distribution  (  ) is a power law of negative index (Bagchi 2013).The unavoidable bias of pulsar searches against distant and intrinsically dim sources, despite being more populous, limits modelling the behaviour of  (  ) below 0.1 mJy kpc 2 (Faucher-Giguere & Kaspi 2006;Lorimer et al. 2006). (  ) must turn over and fall to zero at some unknown value of .MeerKAT provides the opportunity to probe the low-luminosity population, improve models of  (  ) and maximise pulsar discovery.
Increasing the sample of young radio pulsars would also provide constraints on the Galactic NS birth rate.The lowest estimate by Keane & Kramer (2008), amongst all the model and methoddependent formation rates they considered, is 5.7 +4.1 −2.7 century −1 for all NS manifestations combined.This value exceeds the independently measured CCSN rate of 1.9±1.1 century −1 from the Galactic content of 26 Al (Diehl et al. 2006), though the uncertainties just about overlap.Keane & Kramer (2008) suggest two possible solutions to this inconsistency.Either there is an evolutionary link between some classes of NSs which would shrink the total birth rate, or this is possibly a case of underestimated uncertainty or overestimated rates.
3 40 are listed in the survey but one was later found to be the north-western segment of the catalogued G350.0-2.0 (Green 2019) The most constrained birth rate between the NS types is that of the radio pulsar, from which the others are extrapolated.Therefore more informed estimates of this birth rate will reduce the uncertainty propagating through to the total.
Another manifestation of young NSs are the central compact objects (CCOs).There are ten known CCOs and a further four objects are candidates.These objects have only been seen in Xrays, which are of predominantly thermal origin and of steady flux.The rotation periods of three CCOs has been measured; PSR J1210−5226/1E 1207.4−5209 which rotates every 424 ms (Zavlin et al. 2000), PSR J1852+0040 in Kes 79 with a period of 105 ms (Gotthelf et al. 2005) and PSR J0821−4300 with a period of 113 ms (Gotthelf & Halpern 2009).This cements CCOs as neutron stars, but the characteristic age derived from assuming the rotational evolution of a magnetic dipole would be >100 Myr, seemingly at odds with the youth of their remnants.Transient X-ray emission from many CCOs can be explained as surface hot spots (see De Luca 2017 for a review).Deep radio searches that detect radio pulses or otherwise set deep upper limits would help to understand whether CCOs are faint, quiescent or non-emitting at radio wavelengths.
This survey is a part of the TRAPUM (TRAnsients and PUlsars with MeerKAT) Large Survey Project (Stappers & Kramer 2016), which has collectively discovered over 200 pulsars 4 using the MeerKAT telescope.Having presented the motivation for deep searches of many SNRs and PWNe with MeerKAT, we continue in Section 2 by describing the front and backend infrastructure and how the design of the survey seeks to minimise some of the aforementioned selection effects.The search pipeline is described in Section 3, the first discoveries of the survey are presented in Section 4 and a discussion of these pulsars follows in Section 5.

Beamforming with MeerKAT
MeerKAT is the most sensitive radio telescope in the Southern hemisphere (Camilo et al. 2018;Jonas 2018).It is comprised of 64 dishes of Gregorian-offset design, each 13.5 m in diameter and installed with three dual polarisation receivers that collectively offer a wide frequency coverage; the L-band receiver (856-1712 MHz) (Lehmensiek & Theron 2012), the Ultra High Frequency (UHF)-band receiver (544-1088 MHz) (Lehmensiek & Theron 2014) and the S-band receiver (1750-3500 MHz) (Barr 2018).Table 1 lists the specifications for our time-domain observations which so far have exclusively used the L-band receiver.The presence of narrowband radio frequency interference (RFI) sometimes causes losses of up to 20 per cent of the band.
Voltage-based beamforming functions by coherently summing the signals arriving at each receiver, using a set of geometric delays that correct to a reference point in the array.This forms a tied-array beam, or coherent beam (CB), with a size dependent on the projected baselines of the telescopes.Up to 780 CBs5 can be packed hexagonally into a tiling pattern bounded by any -sided polygon(s) using the Filterbanking BeamFormer User Supplied Equipment (FBFUSE), an on-site, dedicated, high-performance cluster (Barr 2018).To do this, FBFUSE receives the channelised voltages from the correlator/beamformer's (CBF) F-engines (van der Byl et al. 2022) and exe- Using interferometric beamforming for targeted pulsar searches has a number of advantages over previous survey methods.The 'small diameter, long baseline' design of MeerKAT provides a wide field of view (FoV) on the angular scale of SNR diameters, while also providing accurate localisation of pulsars without the need for a timing solution.By compartmentalising searches into many CBs, the trade-off between field of view, telescope time and sensitivity are much more navigable.Additionally, multiple targets in the field can be searched simultaneously.At 1284 MHz, the sensitivity drops to 50 per cent at around 30 arcminutes from the FoV centre (Asad et al. 2021), setting a limit on the size of a source, or distance between sources, that can be observed in one instant.Beyond providing a flexible geometry, CBs are less susceptible to RFI (Chen et al. 2021).Furthermore, they do not 'see' the full signal sky temperature contribution from the source, which is particularly useful for the ring shape of a SNR shell.This provides a sensitivity boost compared singledish searches as the noise is reduced by a factor of approximately √  d when  d telescopes are used, compared to the incoherent beam (e.g.Thompson et al. 2017).

Sensitivity
The sensitivity of our searches is equivalent to the minimum flux density,  min we can detect for pulses of period,  and width, .It is described in terms of the system gain, , temperature,  sys , the observing time,  obs and bandwidth, Δ by the modified radiometer equation (Dewey et al. 1985) as (1) We have appended an efficiency factor,  of 0.7 as determined by Morello et al. (2020) to account for additional sensitivity loss from searching with a Fast Fourier Transform (FFT)-based search pipeline (see Section 3). is the factor accounting for the raw amplitude information lost during digitisation, which for 8-bit digitisation is <1 per cent (Kouwenhoven & Voûte 2001).For our calculations we assume this value to be 1.At L-band, the system temperature,  sys is a summation of the sky temperature,  sky , the 18 K receiver temperature and the 4.5 K spillover contribution at 45 degrees elevation (Bailes et al. 2020).Assuming a minimum signal-to-noise, S/N of 9, 40 minute integration time and an upper bound on the duty cycle,  = / of 10 per cent (Posselt et al. 2023), our observations achieve a flux density upper limit of 21 μJy or 31 μJy for the 64 or 44 dishes respectively.These calculations refer specifically to the centre of the coherent beam.The average dilution in sensitivity considering the pulsar could be located at any point in a coherent beam is around 35 per cent, though this varies depending on the baselines used and the position of the source in the sky.
Figure 2 displays the flux density sensitivity curve from Equation 1for this survey, and also for some non-targeted pulsar search surveys of the Galactic plane.For this survey, a  sky of 5.3 K is estimated by scaling the cold-sky temperature from Haslam et al. (1982) using a spectral index of −2.6 for diffuse radio synchrotron emission (Reich et al. 1988;Zheng et al. 2017).We scale the quoted  sky of the other survey publications to 1284 MHz in the same way.We set  equal to 1 for all of the flux density curves in Figure 2 to allow for a better comparison, as search pipeline efficiencies are sometimes not assigned for the non-MeerKAT surveys.The sensitivity degradation due to intra-channel dispersion smearing affects shorter pulse periods more strongly.The intrinsic width of the pulse,  i is broadened by sampling, scattering and DM smearing such that the measured width, (Dewey et al. 1985), where  accounts for the finite pulse profile resolution due to the sampling time,  samp , though for our searches we take it to be 1 (Men et al. 2023).For  scatt , we use the expression derived by (Lewandowski et al. 2015), where log  scatt (ms) = −6.344+ 1.467 log DM − 0.509 (log DM) 2 .Our choice of 4096 channels ensures that at a DM below 1000 pc cm −3 , the dominant factor limiting our sensitivity for  > 10 ms is  sky .Also shown are the flux densities of pulsars associated with SNRs and/or PWNe for which this has been measured (Manchester et al. 2005), also scaled from 1400 MHz using the mean pulsar spectral index of −1.60±0.03(Jankowski et al. 2018).The Galactic Plane Pulsar Snapshot (GPPS) survey (Han et al. 2021) with the Five-hundredmeter Aperture Spherical Telescope (FAST, Li & Pan 2016) achieves flux density limits of 7 μJy, consistently more sensitive than our observations across all periods at which young pulsars are known to rotate.The sensitivity achieved by the PALFA survey with the Arecibo telescope (Cordes et al. 2006) is fairly competitive with our searches for  >100 ms, but quickly deteriorates at higher DMs.
Further to this, we expect many young pulsars in the Galactic plane to have fast periods in the range of 30-100 ms.We would expect to detect all 45 of the SNR-associated pulsars, with the exception of PSR J1907+0602, the pulsar with the lowest flux density in Figure 2.This pulsar was discovered following a 1.8 hr search with the Arecibo telescope (Abdo et al. 2010).Consideration of such cases, where exceptionally deep observations uncovered a very faint pulsar, informed our decision to efficiently balance telescope time between the depth of each search and the number of targets.

Sources targeted to date
So far, we have observed 55 targets, of which 32 are G22 SNRs, 17 are candidate SNRs, 2 are isolated PWNe and 4 are unidentified TeV sources.Due to the extremely competitive sensitivity of the FAST GPPS, we avoid the footprint of that survey, therefore the targets of this survey are located at Galactic longitudes <30 • .Information on each source and its observation are displayed in Table 2.The sky temperature for each source was estimated from the GSM2016 model (Zheng et al. 2017) using the pygdsm7 package (Price 2016) to obtain more accurate estimates of system temperatures for each source (Price 2021).The sizes quoted are the extent of each source, not the explicit area searched with CBs.
There is a trade-off in sensitivity for the more extended sources due to the CB limit.If a source is too large to be tiled with 780 beams, there are two options; i) change the CB overlap level or ii) increase the beam size by using only the innermost 44 dishes, which have a maximum baseline of 1 km.To calculate the number of beams required to fully tile a target, we use the multibeam simulation package mosaic8 , described in Chen et al. ( 2021) to simulate the point spread function (PSF) of the CBs and their tiling positions.The PSF is time-dependent due to the movement of the source across the sky (Chen et al. 2021), so the PSF is approximated at the midpoint of the planned observation.We fix the radius at which the CBs overlap to be the 50 per cent power radius of the PSF.FBFUSE can only beamform using multiples of four dishes.
Among the first sources to be observed was SNR G28.7−0.4,which we wished to observe with the full 64-dish gain.Due to the reduced beam size of the 8 km baseline, we mitigated the prohibitively high data rate of 684 CBs by reducing the time resolution to 306 μs.In all subsequent observations of more extended sources, however, we have optioned to use the core 44-dish baseline to provide a eight-fold increase in CB coverage, rather than halve the sampling rate which only increases the coverage by a factor of two.
All the sources in Table 2 can be covered by the half-power beam width of the incoherent beam, with the exception of G28.8+1.5.A handful of SNRs could not be fully tiled with coherent beams by just using the core baselines of MeerKAT.In these cases, which are indicated in Table 1, we have targeted regions of interest.For G27.8+0.6 and G28.8+1.5, we placed coherent beams on two and seven X-ray point sources, respectively, that were identified by Misanovic et al. (2010) as potential PWN or NS candidates.The candidate remnant G350.8+5.0 (Hurley-Walker et al. 2019) contains the steep spectrum source GLEAM J170154−333857 which could be a pulsar.A circular region 21 ′ in diameter and centred on the GLEAM source has been searched.For the remnants G296.5+10.0,G347.3−0.5 and G353.6−0.7,we placed 480 CBs centred on the associated CCO, tiling a circular area 9.7 ′ , 6.1 ′ and 9.1 ′ across, respectively, instead of the entire shell.G355.4+0.7 contains a possible associated γray source 4FGL J1731.2−3235(Acero et al. 2016), around which we tiled a region 5 ′ across with CBs.In future work we intend to produce maps of all beam tiling regions.
Table 2. Sources searched so far in the survey and ordered by observation date.Some sources have alternative names which are provided at the bottom of the table.The type is defined as follows; 'S' for shell-type G22 SNR, 'C' for composite-type G22 SNR, 'F' for filled centre G22 SNR, 'P' for isolated pulsar wind nebula, 'cand' for candidate SNR, 'U' for unidentified TeV source.'?' denotes ambiguity in the literature.Source sizes are taken from the following; [1] Green (2019) or Green (2022) in the case of SNRs observed before or after 2023-09-01 respectively, [2] GLOSTAR (Dokara et al. 2021), [3] (Brogan et al. 2006), [4] MAGPIS (Helfand et al. 2006), [5] GLEAM (Hurley-Walker et al. 2019), [6] the catalogue of TeV sources tevcat for H.E.S.S. sources, [7] (Filipović et al. 2022). min is calculated using Equation 1, where  = 0.7, and specifically refers to the centre of a coherent beam.The upper limit on the radio luminosity of the pulses,  min =  min  2 , where the distance,  is the most recent measurement.Previous flux density limits from targeted searches have been scaled to 1284 MHz using a spectral index of −1.6, and their original values are taken from the following; † Target was observed simultaneously with one or more other sources, identifiable by their equivalent Observation Start times.
The associated object is a CCO or candidate CCO.
Due to the angular size of the source, only regions of interest were targeted with coherent beams.These are described in the text. The associated object has a known period from X-ray timing.

Timing observations
We continue to observe newly discovered pulsars in order to accurately measure their changing rotation rate.If the pulsar has been detected in more that one beam in the discovery observation, we can begin timing immediately by using SeeKAT9 (Bezuidenhout et al. 2023), a multibeam localisation tool that uses the S/N and a PSF simulated by mosaic, to obtain a position with arcsecond accuracy.If the pulsar is seen only seen in one beam, we attempt a multibeam detection using the first timing observation.The timing campaign is pseudo-logarithmic, consisting of nine epochs where the first two observations are spaced about half a day apart.This cadence is essential for keeping track of the rotation of the pulsar, as young pulsars tend to have a high rate of spin-down.A minimum of 58 of the 64 telescopes are used with FBFUSE and the APSUSE backend, to which the data are written for offline analysis.
We also make use of the UHF-band receivers for timing.The negative spectral index of pulsar radio emission may result in a higher signal-to-noise than at L-band for an equivalent integration time.The extended frequency range allows better opportunity for investigating frequency-dependent profile changes, which provide insights into the emission mechanism (e.g.Xu et al. 2021).However, the effects of pulse broadening due to scattering and dispersive smearing are more pronounced in this band.To mitigate these effects, UHF observations use the Pulsar Timing User Supplied Equipment (PTUSE, Bailes et al. 2020) backend in search mode.PTUSE is a dedicated pulsar timing backend, providing much finer sampling in time and recording full Stokes polarisation information.Data are written in the PSRFITS format (Hotan et al. 2004), sampled every 7.5 μs, recorded across 4096 channels and coherently dedispersed at the DM of the source.We did not use PTUSE for the L-band observations, because this could introduce a systematic offset in L-band pulse arrival times between the two backends.
The MeerKAT timing observations are supplemented with observations from Murriyang.These data are recorded with the Ultra Wide L-band (UWL) receiver covering 704-4032 MHz (Hobbs et al. 2020), and are stored in PSRFITS format.We use the search mode, channelising the data into 3328 channels across 26 subbands, and sampling every 64 μs.

Searching
To search for periodic signals, we use peasoup10 (Barr 2020; Morello et al. 2019), a GPU-based C++/CUDA library for time-domain linear acceleration searching for pulsations with frequencies between 0.1-1100 Hz.While we do not expect to find pulsar binaries in SNRs and wind nebulae, we elected to search accelerations between ±5 m s −2 , assuming constant acceleration.A non-constant acceleration could result in a more accelerated pulsar being smeared below detection threshold if the entire integration time is searched at once, especially for tight-orbit binaries.We do not attempt to segment our searches as our integration times of 20-40 minutes would not risk much loss in sensitivity to mildly accelerated systems.Dedispersion is performed incoherently per DM trial using the dedisp library (Levin 2012;Barsdell et al. 2012).As many sources in the survey lie on the Galactic plane, it is necessary to search up to a high DM.The contributions to pulse smearing due to DM step size and sampling time become dominated by the intra-channel smearing at higher DMs, so one can make use of a dedispersion plan that downsamples in DM and time to reduce the number of trials.We therefore generated a plan using presto/DDPlan.py(Ransom et al. 2002(Ransom et al. , 2011)), consisting of 8784 trials up to 3000 pc cm −3 .Before being searched with peasoup, the data are cleaned using the filtool RFI cleaning algorithm from pulsarX11 (Men et al. 2023) which eliminates narrowband RFI using statistical thresholds and broadband RFI using a zero-DM matched filter (Men et al. 2019).It also replaces intra-channel samples of interference with noise.A separate channel mask is also applied, which masks 64 MHz of the band in order to mitigate common sources of narrowband RFI.After cleaning, the data are delivered to peasoup where they are dedispersed.For each DM trial, the cleaned and dedispersed time-series is de-reddened to remove low-frequency RFI before being re-sampled to account for the acceleration trial value.After this, the corrected series is Fourier transformed and harmonics are summed up to the 8 th harmonic.To facilitate efficient use of the Fast Fourier Transform operation, the data are zero-padded from the final sample up to the next instance of 2  .If the amplitude at a sampled frequency exceeds the S/N threshold, it is checked against a list of previously identified RFI frequencies before being written as a candidate.
The candidates, which normally number 10 3 per beam, are stored in XML format files and require sifting.We use candidate_filter.py12, which performs further RFI mitigation and multibeam spatial clustering (a detailed description is found in Voraganti Padmanabh 2021).This reduces the number of candidates by a factor of approximately 100.The filtered candidates are folded with psrfold_fil, another application of the pulsarX suite.psrfold_fil optimises the period, period derivative and DM.The time-series to be folded is dedispersed using an algorithm that avoids redundant DMs, and is thus more efficient than brute force methods (Men et al. 2023).Optimised candidates are classified using two Pulsar Image-based Classification System (PICS) Machine Learning classifiers (Zhu et al. 2014) that are trained on PALFA or TRAPUM data sets respectively.Candidates with a pulsar likelihood below a lenient score of 10 per cent are rejected.The final set of candidate diagnostic information and metadata are inspected by eye in the specialised candidate viewer CandyJar13 .Candidates are designated either Tier 1 (high probability of novelty), Tier 2 (potentially novel but uncertain), RFI, noise or as a known pulsar.Upon identification of a Tier 1 candidate, we check that the discovery period is not a harmonic of the true period by refolding the filterbank using psrfold_fil at integer multiples of the period; 2 and 3, and unit fractions; /2 and /3.
We investigate all the Tier 2 candidates by directly refolding the coherent or incoherent beam data.To do this, the data are cleaned with filtool then folded at the candidate's topocentric period with DSPSR14 .A second round of cleaning is done using clfd15 .Finally, we use pdmp to optimise the S/N of the signal by searching over period, DM and acceleration.Both DSPSR and pdmp are tools from the pulsar data analysis package psrchive16 (Hotan et al. 2004;van Straten & Bailes 2011).The priors used by pdmp are informed by the uncertainty on the period, period derivative and acceleration of the candidate.The diagnostics are then inspected to see if the candidate is real.A higher S/N can potentially be extracted from the refolded data compared to the search due to the different optimisation methods; pdmp selects the highest S/N across every element in the parameter grid, whereas psrfold_fil instead maximises the  2 of the profile against the noise before calculating the S/N of the best profile.Any latent RFI that was not cleaned by filtool or clfd is flagged interactively using the psrchive/psrzap tool to flag, and psrchive/psrsh to mask the affected channels or subintegrations, which can further boost the S/N of the refolded data.
Two sources in the survey are associated with neutron stars of a known rotational period but remain without a radio detection.The data of the coherent beam at the position of the neutron star is folded directly using an ephemeris from the known rotation parameters, and use both DM values predicted by the ne2001 (Cordes 2004) and ymw16 (Yao et al. 2016) electron density models.To obtain the DM values, the pygedm17 (Price et al. 2021) package was used.The data are cleaned and folded as before, except this time pdmp searches over a DM range of ±10 pc cm −3 .

Timing
Both the CB data from APSUSE and search mode PTUSE data are folded and cleaned by way of the same method as just described, but without searching in period or DM space, as these are already well constrained.Instead, the psrchive/pam command is used to dedisperse and then combine the folded data in frequency and time to obtain the summed pulse profile.To generate the topocentric timeof-arrival (TOA) of the average pulse, the profile is cross-correlated with a noise-less profile template using the psrchive/pat command.Noise-less templates are generated from the discovery observation using psrchive/paas, in which we interactively fit one or more scaled von Mises functions to return a high S/N mean profile.Separate templates are generated for each receiver band.For the template, we use either the profile of the discovery observation, or where possible, a high S/N summation of multiple profiles that have been aligned using a phase-connected timing solution.With TEMPO218 (Hobbs et al. 2006), the TOAs are corrected by referring them to the barycentre of the Solar System and fitted with a timing model.For UWL data from Murriyang, the same process for MeerKAT data is followed for cleaning folded data and generating TOAs.Any systematic offset in TOAs due to differing backends is corrected by fitting for a temporal jump.

RESULTS
Of the candidates inspected from the searches of all the sources in Table 2, two Tier 1 candidates were identified and subsequently confirmed as new pulsars.This difference is less than one bin width for both pulsars but is nonetheless accounted for.The standard deviation of the noise is shown as an error bar at the bottom of each pulse.We see evidence for a precursor component for PSR J1818−1502 at a phase of approximately 0.45.
candidate PWN G22.0+0.0 (Ueno et al. 2006;Yamauchi et al. 2016, hereafter YSB16).The pulsar appears to be isolated, and was detected at a topocentric rotation period of 300.85 ms, a dispersion measure of 377 pc cm −3 , and with a S/N of 16.8.The pulse profile generated from folding the data from this observation is shown in the upper left panel of Figure 3. PSR J1831−0941 was detected in an additional 3 beams adjacent to the discovery beam.The S/N of the four profiles, each with 128 flux bins, was measured using the pdmp option of psrchive/psrstat and using SeeKAT we obtained a 1-sigma positional uncertainty of about 10 ′′ , or 1/6 th the width of the CB.In total we have TOAs from 26 observations spanning 408 days.The residuals of these arrival times against the model of timing parameters in Table 3, are shown in the upper panel of Figure 4.The residuals have a root mean square (rms) of 2.2 ms.We measure a period derivative of 8.56 × 10 −14 s s −1 , from which we derive a characteristic age of 56 kyr and a spin-down luminosity of 1.2 × 10 35 erg s −1 .This places PSR J1831−0941 within the population of young pulsars in SNRs.

PSR J1818−1502
PSR J1818−1502 was discovered in a coherent beam with a S/N of 14.6 in a search of the shell-type remnant G15.9+0.2.It has a period of 572.5 ms and a DM of 435 pc cm −3 .The profile of PSR J1818−1502 from this observation was produced in the same way as for PSR J1831−0941, and is shown in the upper right panel of Figure 3. Using the same procedure as in Section 4.1, SeeKAT re-Table 3. The properties of both new pulsars, including the measured and derived quantities from their coherent timing solutions after fitting TOAs in TEMPO2.
The numbers in parentheses are the 1-sigma uncertainty on the final digit.The inferred distances are derived from the ne2001 and ymw16 electron density models.The upper limit on the period derivative for PSR J1818−1502 is the 1-sigma uncertainty on the value returned when fitting for  in TEMPO2.The position for PSR J1818−1502 has not been fitted, instead it is held at the SeeKAT position.turned a 4-sigma position of 18 h 18 m 54.s 3 +0.2−0.5 , −15 • 02 ′ 04.′′ 5 +4.7 −5.9 .The probability distribution of the localisation is shown in the lefthand side of Figure 5.The pulsar's position within G15.9+0.2 is shown in the map on the right-hand side.The quoted uncertainty from SeeKAT is more representative of the total position error than 3-sigma due to the low signal-to-noise of the detections (Bezuidenhout et al. 2023).
As of now we have 8 observations spanning 118 days, for which TOAs have been generated using a template made from the discovery observation.The timing residuals from a fit of the SeeKAT position, and the period and period derivative in Table 3, are shown in the lower panel of Figure 4.They have a rms of 7.6 ms.We are not yet able to make a significant measurement of , and so present an upper limit equivalent to the 1-sigma uncertainty of 1.2 × 10 −16 s s −1 on the fitted values.We are confident that the actual magnitude of the spin-down rate is much smaller than this value, as it would otherwise have been resolved with the cadence and baseline of these TOAs.The upper limit on  from the fit leads us to conclude that PSR J1818−1502 is a very old canonical pulsar with a characteristic age greater than 78 Myr.We will continue timing PSR J1818−1502 in order to determine the true  and constrain the age, spin-down luminosity and magnetic field strength.It is possible that in fitting for , we are absorbing an apparent change in period due to a position error.The maximum position error at 4-sigma is 7 ′′ .If the true position differs by this amount then the apparent change in rotation period, if the period were measured at the near and far side of the Earth's orbit, induces a period derivative of 8 × 10 −17 s s −1 , which is the same order of magnitude the uncertainty of the fit.Therefore a baseline of at least one year will be required to ensure the position and period change of the source are disentangled.

Pulse profiles
The widths of the profiles at 50 per cent of the peak intensity, W50 at 1284 MHz have been calculated using a high S/N profile of each pulsar by combining data from all MeerKAT observations where a coherent beam was placed directly on the best position of the pulsar.For PSR J1831−0941 and PSR J1818−1502 there are 6 and 7 such observations, respectively.The pulse profiles of these high signalto-noise data are shown in Figure 3 where the lower left panel the profile of PSR J1831−0941 and the lower right panel shows that of PSR J1818−1502.We obtained these profiles by folding the data from each timing observation with DSPSR, using ephemerides composed of the measured parameters in Table 3.They were cleaned with clfd before being summed in frequency and time.psrchive/psradd was used to combine the profiles into an average profile.To calculate the width of each profile, a noise-less template was generated by modelling the pulse component(s) of the profile with von Mises functions in psrchive/paas.To determine W50, the template had noise added to it using the noise statistics of the average profile's off-pulse region, before being fitted using the fitvonMises function in PSRSALSA (Weltevrede 2016) which considers the concentration parameter of the best combination of von Mises fits.This process is repeated 1000 times.Using the concentration parameter of the fitvonMises function from all 1000 fits, we get a distribution of widths and obtain W50 and its uncertainty from the mean and standard deviation.The mean widths from each distribution are W50 = (25.7 ± 0.5) ms for PSR J1831−0941, and W50 = (24.1 ± 2.7) ms for PSR J1818−1502.These are converted to a duty cycle,  and presented in Table 3.The pulse width for PSR J1831−0941 is larger than most pulsars of a similar period (Posselt et al. 2021).
The pulse shape of PSR J1831−0941 appears single-peaked and approximately symmetric, whereas PSR J1818−1502 appears to have a main pulse preceded by a weaker and slightly narrower precursor.The evidence for the reality of the precursor is two-fold.Two von Mises functions are required to adequately model the profile of PSR J1818−1502 in psrchive/paas.Secondly, the distribution of W50 from the iterative profile fits is not Gaussian; aside from a peak at approximately 25 ms, there is another peak at about 21 ms, which likely corresponds to narrower von Mises functions not capturing the precursor.This asymmetry inflates the uncertainty on .
We investigated the effect of interstellar scattering on the profile of PSR J1831−0941 at L-band.The same high signal-to-noise profile with which we measured the pulse width was used.The scattering analysis tool scatfit19 (Jankowski et al. 2023) was used to split the 856 MHz band into four subbands, each of 256 channels, and fit an exponentially modified Gaussian distribution to the flux.This function is a robust description of a pulse scattered by a thin scattering screen (Cordes & Lazio 2001;Oswald et al. 2021).We find little evidence for pulse broadening due to scattering across the band.A scattering timescale of (12 ± 4) ms is measured between 856-1070 MHz.The dispersion smearing timescale dominates broadening at higher frequencies.
The flux densities of the new pulsars are calculated using Equation 1 where  = , omitting  and using Δ eff from Table 1.The sky temperatures at both positions were determined using the GSM2016 model (Zheng et al. 2017).We find the flux density at 1284 MHz,  1284 to be (119 ± 12) μJy for PSR J1831−0941, and (33 ± 3) μJy for PSR J1818−1502.The uncertainties on these values of 10 per cent are determined from the error on both the sky temperature and the signal-to-noise.Both are fainter than could reasonably be detected by the vast majority of previously targeted pulsar searches of SNRs.Using the DM-derived distances in Table 3, PSR J1831−0941 has a radio luminosity,  1284 of 2.1-2.9 mJy kpc 2 , whereas PSR J1818−1502 has a luminosity of just 0.6-0.8mJy kpc 2 .The range is the difference between the ne2001 and ymw16 models.PSR J1818−1502 has a pseudo-luminosity comparable to some of the more intrinsically faint pulsars that have been found by the FAST GPPS survey (Han et al. 2021).

Non-detections
Many of the observations produce a handful of Tier 2 candidates.So far, none of these have been found to be convincing signals from new pulsars after follow up analysis described in Section 3. We regularly detect known pulsars present in the FoV, mainly in the incoherent beam, but occasionally in clusters of coherent beams at the edge of the tiling arrangements.Indeed, analysis by Padmanabh et al. (2023) find that known pulsars are detected as expected in the MPIfR-MeerKAT Galactic Plane Survey at L-band, which also uses the peasoup pipeline.The upper limits on the targets we have searched are 21-52 μJy at 1284 MHz, and are a significant improvement on the depth of previous targeted searches or of non-targeted searches of the plane.

DISCUSSION
5.1 Association between PSR J1831−0941, G22.0+0.0 and G22.045−0.028G22.0+0.0 was first identified by Ueno et al. (2006) during a reanalysis of ASCA Galactic Plane Survey (Yamauchi et al. 2002) data to search for synchrotron emission from low latitude SNRs with a low radio surface brightness (Sugizaki et al. 2001).They identified an extended region of diffuse, hard X-ray emission about 10 ′ across.YSB16 re-detected the emission and measured the distance, for which they used Galactic HI column density and CO velocity measurements, to be 5.2±2.2kpc.The uncertainty range overlaps with the DM-derived distances to PSR J1831−0941 in Table 3.The discovery of G22.045−0.028(Dokara et al. 2021), a 15.4 ′ across shell encompassing the X-ray position, provided more evidence for a plerionic scenario.The position of PSR J1831−0941 relative to the radio shell and X-ray emission is shown in Figure 6.We used the reprocessed 20 Suzaku data from the observation by YSB16 (Obs ID: 505025010 data from the three detectors were added and then smoothed using a 1-sigma Gaussian kernel.YSB16 estimated a characteristic age of 10 4 -10 5 yr for the putative pulsar, based on its clumpy morphology and also the relation between age and size in the X-ray band (Bamba et al. 2010), to which our characteristic age of 55.7 kyr agrees.The position, inferred distance and adolescence of PSR J1831−0941 provides strong evidence that G22.0+0.0 and G22.045−0.028form one composite-type SNR.
At 5.2±2.2kpc and 10 ′ across, the nebula is 15±6 pc in diameter, and the shell is 23±9 pc across.If we assume the age of the remnant is equivalent to  c from Table 3, then we find G22.0+0.0 follows the observed relation between X-ray PWN size and age (Bamba et al. 2010), which predicts a size of 10-21 pc.Given an X-ray flux from YSB16 of 8.36 × 10 −9 erg s −1 m −2 , we derive an X-ray conversion efficiency of 2.2% at 5.2 kpc, which is higher than other pulsars of a similar  (Possenti et al. 2002), though contained within the spread of values for any given .A high photon index of the X-ray spectra is expected for more efficient winds, as the spectrum of the radiation depends on the injected energy spectrum and wind magnetisation, which themselves are driven by the spin-down power of the pulsar (e.g.Torres et al. 2013).The photon index of 1.7±0.3measured by YSB16 is typical for PWNe and agrees with this trend (Kargaltsev & Pavlov 2008).
As Figure 6 shows, PSR J1831−0941 is located centrally within the shell whilst the nebula is concentrated in the South-Western quadrant.Given the estimated age of the system, it is probable that the nebula has already been disrupted by the reverse shock (Gelfand et al. 2009), which would have temporarily counteracted the expansion of the nebula and broken up the wind as it mixed with the reheated interior (Blondin et al. 2001).It is possible that this process is ongoing, which would perhaps explain the higher than average conversion efficiency as a consequence of a surge in pressure and magnetisation of the pulsar wind.The displacement of the PWN could be explained by a westwards gradient in the ISM density, making G22.0+0.0 a similar system to SNR G327.1−1.1 (Temim et al. 2015).The cluster of Hii regions to the West of the shell (Dokara et al. 2021) could be related to the density asymmetry.
c is not a fully reliable proxy for the true age, a point which is discussed this further ahead in this section.If, for now, we use  c , we calculate the transverse component of the kick velocity,  ⊥ .The geometric centre of the G22.045−0.028from Dokara et al. (2021) is assumed as the explosion site.The PWN is not an indicator of birth location, due to the effects of the reverse shock interaction and the movement of the NS.Assuming a single system, we use the distance to the PWN from YSB16, we find  ⊥ = 32 km s −1 , with a range of 18 km s −1 to 45 km s −1 given the ±2.2 kpc uncertainty.Depending on the magnitude of the radial velocity, the  kick vector is likely small compared to the young pulsar population (Hobbs et al. 2005).The observed change in period induced by  ⊥ across the celestial sphere, known as the Shklovskii effect (Shklovskii 1970), would increase the observed period derivative by 7 × 10 −21 , thus is not significant compared to the uncertainty on the  we have obtained.We have investigated whether a bow-shock scenario is realistic.No radio emission, nor any GeV or TeV excess has been identified at position of the PWN.However, it is not unusual for bowshock nebulae to only show up in one frequency regime due to the different radiative cooling timescales.Pulsars within shells are not expected to cause bow-shocks due to the low density of the medium after being swept up by the SN shock, nevertheless it is possible the pulsar may have escaped the far side of the remnant.A pulsar with  kick ⩾ 650 km s −1 would be able to escape its shell before the remnant dissipates (van der Swaluw et al. 2003), and in doing so would form a bow-shock that emits copious radio synchrotron photons (Shull et al. 1989).We note both the lack of observational evidence for a bow-shock scenario, and also the small likelihood of a dominant radial component of  kick that would be required.Thus, we do not consider G22.0+0.0 to be a bow-shock nebula.Dokara et al. (2021) were not able to identify any radio emission from the wind nebula itself.We can estimate the expected angular size of this emission to determine whether it is not luminous enough, or has been resolved out.Stappers et al. (1999) derived the angular size of the region shocked by the wind to be where  kpc is the distance in kpc,  1000 = v kick × 1000 km s −1 ,  is the energy conversion fraction from the spin-down energy to the wind and  0 is the local number density of particles in cm −3 .If we assume a spherical wind in the face of reverse shock disruption, and also assuming  kick ≈  ⊥ ,  0 = 0.1 and  = 1%, then  = 2 ′′ , much smaller than the 18 ′′ resolution of the GLOSTAR observations.There is no mention of a point source, thus we can estimate an upper limit on the radio luminosity, using the methodology of Frail & Scharringhausen (1997) and scaling the frequency using a Crab-like spectral index of −0.3.Given the upper value for the rms noise is 0.15 mJy at an effective central frequency of 5.8 GHz (Dokara et al. 2021), the radio luminosity,  R given by  R = 9.1 × 10 28  2 kpc  5.8 (3) is 3.7 × 10 29 mJy kpc 2 , corresponding to a conversion efficiency,  R /  of about 3 × 10 −6 .This is a higher efficiency upper limit than the peak of the distribution from non-detections in the study of pulsars with 31.88 ⩽ log( ) ⩽ 36.53 by Frail & Scharringhausen (1997), which is not unexpected given that the radio nebula luminosity appears to fall with age.It is possible that the nebula may be a TeV source, but falls below current sensitivity limits given the distance and age of the pulsar, though may be detectable by next generation TeV observatories such as the Čerenkov Telescope Array (e.g.Mitchell & Gelfand 2022).As mentioned, these calculations come with the caveat that the model-dependent distances are subject to large and unquantified uncertainties, and that the characteristic age is not a robust evaluation of the true age. c is contingent on the assumption that the pulsar spindown is purely driven by the magnetic dipole radiation, and thus the braking index,  is equal to 3, where Ω ∝ Ω  (Manchester & Taylor 1977), where Ω = 1/.It also assumes the birth period is much less than the present day period.A realistic birth period of 30 ms would only cause  c to change by one per cent for PSR J1831−0941.A much larger uncertainty is imposed by the reality that  < 3 (Noutsos et al. 2013) and may not be constant.For  = 2.5 (Lyne et al. 2015), the true age would be 30 per cent greater that  c .Apart from seeking an age from which to derive  kick , an accurate velocity could instead be obtained from measuring the proper motion.The contribution due to the natal kick would be only about 13 mas year −1 , inducing an apparent  of 8 × 10 −21 s s −1 , corresponding to a TOA drift of less than 10 −4 turns per year.
It is reasonable to expect that high  pulsars might be detectable at γ-rays, with one reason being that the wider beams in the VHE band gives a greater sky coverage, perhaps a factor of 2 higher than radio beams, even for pulsars of moderately high  (Ravi et al. 2010).
PSR J1831−0941 lies within the localisation region of a Fermi Large Area Telescope (LAT) γ-ray source, 4FGL J1830.8-0947(Abdollahi et al. 2020;Ballet et al. 2023), however this source has already been identified as being powered by the energetic PSR J1831−0952 (Laffon et al. 2015).It is possible that γ-ray emission from PSR J1831−0941 is masked by this source, but unfortunately the validity period of the timing solution and the precision on the spin parameters are insufficient for applying the methods of Bruel (2019) and Smith et al. (2019) to search for γ-ray pulsations in the absence of a possibly associated point source.γ-ray pulsations from PSR J1831−0941 may yet be detected in the Fermi-LAT data with a future, more precise timing ephemeris, and with a detailed study of the γ-ray emission from this busy region of the Galactic plane.As a young pulsar, PSR J1831−0941 is a member of a population of pulsars known to glitch, i.e. suddenly increase their spin frequency counter to their spin-down.Pulsars are more likely to glitch when young, possibly because the higher spin-down rate is conjunctive to the trigger mechanisms that might be responsible for disrupting the angular momentum distribution of the NS interior (see e.g.Antonopoulou et al. 2022).We see no evidence for a glitch of any magnitude in the timing residuals we have obtained so far.Using the power-law relation between  c and the glitch rate from Basu et al. (2022), PSR J1831−0941 could be expected to glitch at a rate between 0.2-0.4yr −1 using the 1-sigma uncertainty range.We would therefore expect to wait 4.5 yrs before we could expect a glitch to occur with a confidence of 68 per cent.Furthermore, one may also expect a handful of glitches within the 15-yr interval spanned by the Fermi-LAT data, which would further complicate efforts to detect γ-ray pulsations from.

PSR J1818−1502
G15.9+0.2 is a shell-type SNR of small diameter first identified in the radio band (Clark et al. 1973).The detection of a radio pulsar at this position is interesting given that the SNR already has an associated neutron star.The central compact object CXOU J181852.0−150213, also shown in Figure 5 is best candidate for the collapsed core of the progenitor star the assumed explosion site (Mayer & Becker 2021).We see no evidence of radio pulsations at the position of the CCO to an upper limit of 29 μJy.The distance to PSR J1818−1502 from a DM of 435 pc cm −3 is 5.4 kpc (ne2001) or 4.2 kpc (ymw16) though these are subject to large uncertainties.The distance to the shell is 7-16 kpc from HI absorption (Tian et al. 2019), suggesting that PSR J1818−1502 is an unrelated pulsar in the foreground.
The CCO is offset by about 0.6 ′ from the geometric centre of the SNR at 18 h 18 m 54.s 2, −15 • 01 ′ 55 ′′ .However, the pulsar is much closer to this position at only 0.2 ′ away, and only 0.5 ′ from the CCO.The height of PSR J1818−1502 above the Galactic plane would be 23 pc at the lowest estimate of the distance to G15.9+0.2 of 7 kpc.If we take the scale height for pulsars with a characteristic age of >78 Myr to be 400-600 pc (Mdzinarishvili & Melikidze 2004), approximately 95 per cent of pulsars as old as PSR J1818−1502 would lie >23 pc from the Galactic plane.However, the Galactic plane is densely populated with known pulsars, and becoming increasingly more so, thus it is not completely improbable that we should discover an old pulsar in this part of the sky.We therefore note that despite the offset between the CCO and PSR J1818−1502 being unusually small, this does not by itself suggest a evolutionary connection between the two.

Upper limits on sources with known periods
Two of our targets are young neutron stars with known rotation periods measured from X-ray timing.One is the 38.5 ms radio-quiet pulsar PSR J1849−0001 (Gotthelf et al. 2011).This pulsar emits non-thermal X-rays and powers the wind nebula IGR J18490−0000/HESS J1849−000.The nebula has been detected at both X-ray (Calas et al. 2018) and VHE γ-ray (Terrier et al. 2008;Abdalla et al. 2018a) wavelengths.This is an archetypal PWN system but not easily studied due to its faintness at a distance of about 7 kpc (Gotthelf et al. 2011).We speculated that the distance and proximity to the Galactic plane may have inhibited the detection of radio pulsations of previous surveys, for which our search would be more resilient due to the sensitivity and wide bandwidth extending above 1700 MHz where scattering effects are weaker.However, to a flux limit of  min = 34 μJy, we do not see any pulsations from PSR J1849−0001.This could be due to a differing geometry between the X-ray and radio beams, such that the radio beam is directed too far away from the line of sight to Earth to be detectable.Alternatively, the radio beam may be incident but the emission too faint due to a insufficient radio conversion efficiency,  R .We can estimate an upper limit on  R using our flux limit and the relation between luminosity and flux density from Lorimer & Kramer (2005), where the radio luminosity,  R = 2.7 × 10 27 (  1400 mJy ) 2 (  kpc ) 2 erg s −1 .Using the estimated distance of 7 kpc, it follows that the upper limit for  R =  R  = (3±1) × 10 −11 , where  min has been scaled from 1284 MHz to 1400 MHz using a spectral index of −1.6.A distance uncertainty of 1 kpc was assumed to account for the depth of the Scutum-Centaurus Arm in which the system is thought to located.At 7 kpc, the X-ray conversion efficiency,  X is ≈0.003 for flux integrated between 2-10 keV (Gotthelf et al. 2011), so the ratio between  R and  X is ∼10 −8 .Both  R and  R / X are approximately 2 orders of magnitude lower than for pulsars of high  for which both radio and X-ray emission is seen (Szary et al. 2014), which hints at the possibility that the radio emission could be detected, but is instead beaming away.
We also targeted the CCO associated with the supernova remnant G296.5+10.0,PSR J1210−5226/1E 1207.4−5209, which has been found to emit at X-ray wavelengths with a period of 424 ms (Zavlin et al. 2000).PSR J1210−5226 is one of only three CCOs to have a measured rotation rate.We see no evidence of radio pulsations above a flux density of 21 μJy.The depth of this limit suggests that the CCO is radio-quiet at the time of the observation.The pulsed X-rays likely come from thermal emission from hotspots in the polar cap region, incidentally the same zone where radio emission would be generated for canonical pulsar radio emission (Ruderman & Sutherland 1975).The opening angle of the X-ray emission is likely larger than the radio beam, so it would not be extraordinary to only detect the former.However, the presence of X-rays does not itself compel there to be radio emission.The very low surface magnetic field of around 10 11 G (Halpern & Gotthelf 2015) that could be suppressing spin-powered radio emission may also be insufficient for heating to be provided from the impact of rotation-powered particles.Nevertheless, heating is still possible via internal cooling focused on the polar car by a strong crustal magnetic field (Shabaltas & Lai 2012).Despite our deep upper limit we are still unable to rule out a geometrical explanation for a lack of radio emission.A pulsed radio detection or otherwise deeper upper limits are essential for understanding the physics governing multiwavelength emission from CCOs, and how this differs from other classes of neutron stars (De Luca 2017).

SUMMARY
We have introduced and described a new large survey of supernova remnants, pulsar wind nebulae and unidentified TeV sources from the H.E.S.S. catalogue.The high sensitivity of MeerKAT has allowed upper limits on radio pulsations between 856-1712 MHz of around 30 μJy.The 55 sources that have been targeted so far are composed of 32 supernova remnants, 17 candidate supernova remnants, two isolated pulsar wind nebulae and four unidentified TeV sources.As a result of our searches, two new pulsars have been discovered.The positional alignment of the new pulsar PSR J1831−0941, the candidate radio shell G022.045−0.028, the PWN G22.0+0.0, a characteristic age of 56 kyr and the pulsar being energetic enough to power the plerionic X-ray emission confirms within reasonable certainty that these are associated.PSR J1818−1502 is a faint, slow pulsar that we have discovered close to the centre of the shell-type SNR G15.9+0.2.Given the distance from DM, the slow rotation period and the presence of CCO J181852.0−150213, this pulsar is likely not associated with the remnant.We do not anticipate an evolutionary connection between the pulsar and the progenitor of the CCO, despite their close proximity on the sky.We do not detect radio pulsations from the CCO.We were unable to detect radio pulsations from five other CCOs and three candidate CCOs.Both pulsar discoveries continue to be timed with MeerKAT and the Parkes radio telescope, and the search for new pulsars is ongoing.Once completed, we intend to use the statistics of the survey to perform a robust analysis by way of a population synthesis of young pulsars born from supernovae.With this, we aim to explain the factors that dictate the detectability of this population, such as beaming fractions, the influence of the interstellar medium, the Galactic neutron star birth rate and the luminosity distribution of radio pulsars.

Figure 1 .
Figure 1.Histogram of new SNRs identified in the radio regime binned by their date reported in the literature.Bins are 5 years in width.The blue columns constitute the confirmed SNRs listed in G22.

Figure 2 .
Figure 2. Lower limit for detectable pulse flux densities,  min as a function of pulse period for two representative DMs. is set to be unity and full bandwidth availability is assumed for all flux density lines.Two lines are shown for this work, each corresponding to the core 44-dish or full 64-dish array configuration.Also shown are the same sensitivity curves for the HTRU low-lat (Keith et al. 2010), PALFA survey with Arecibo (Cordes et al. 2006), MPIfR-MeerKAT Galactic Plane Survey at L-band (Padmanabh et al. 2023) and the FAST Galactic Plane Pulsar Survey (Han et al. 2021).Flux densities of pulsars in SNRs and PWNe are plotted, though most have fluxes beyond the maximum shown and are therefore displayed as up-arrows.All curves and pulsar fluxes have been scaled to the central frequency of the MeerKAT L-band reciever, using a spectral index of −1.60±0.03.The error bars account for the uncertainty on the spectral index.

Figure 3 .
Figure3.The integrated pulse profiles at L-band for the PSR J1831−0941 discovery observation (upper left) and for combined timing observations (lower left) and equivalently for PSR J1818−1502 (upper right and lower right, respectively).The data are folded onto 64 flux bins using the best ephemerides from timing, and are centred using the psrchive/pdv centring function.The effective time resolution is greater than the sampling time due to DM smearing within frequency channels.This difference is less than one bin width for both pulsars but is nonetheless accounted for.The standard deviation of the noise is shown as an error bar at the bottom of each pulse.We see evidence for a precursor component for PSR J1818−1502 at a phase of approximately 0.45.

Figure 4 .
Figure 4. Residuals of the arrival times from the timing models for PSR J1831−0941 (upper panel) and PSR J1818−1502 (lower panel).The TOAs were fitted in TEMPO2 resulting in the timing solutions in Table 3.The folded data is split into multiple subintegrations in time.There are 58 TOAs shown for PSR J1831−0941 and 49 for PSR J1818−1502.

Figure 5 .
Figure 5. Left: Close up of the localisation from SeeKAT.The red contours are the two (bold), three and four-sigma levels of the position likelihood.The likelihood function is plotted on the opposite sides of the position axes.The 50 per cent power level contour of the PSF simulated by mosaic is plotted for each beam as solid black lines.The PSF is displayed in the upper right with an elliptical approximation of the 50 per cent power contour.Right: Image of G15.9+0.2 from MAGPIS at 20cm with the VLA(White et al. 2005).The seven beams in which PSR J1818−1502 was detected are shown in yellow, and the best position from SeeKAT is marked with a green cross.The white dashed ellipse defines the shape of the remnant (G22,Dubner et al. 1996).The white cross is the position of CXOU J181852.0−150213fromMayer & Becker (2021).The white arrow points North.

Figure 6 .
Figure 6.Image of the field towards G22.045−0.028from MAGPIS at 20cm with the VLA(White et al. 2005).The blue contours are the X-ray emission seen by YSB16 (reprocessd Suzaku data).The yellow ellipses and green cross are the four discovery beams and best position of PSR J1831−0941, respectively.The dashed circle is the shape of G22.045−0.028(Dokara et al. 2021).The white arrow points North.

Table 1 .
A summary of specifications that are utilised for searching with the L-band receivers of MeerKAT.