Efficient selection of gravitationally lensed OH megamasers with MeerKAT and the Square Kilometre Array

There has been a recent resurgence in hydroxyl (OH) megamaser research driven by Square Kilometre Array (SKA) precursor/pathfinder telescopes. This will continue in the lead-up to the SKA mid-frequency array, which will greatly expand our view of OH megamasers and their cosmic evolution over $\gtrsim80$ per cent of the age of the universe. This is expected to yield large scientific returns as OH megamasers trace galaxy mergers, extreme star formation, high molecular gas densities, and potentially binary/dual supermassive black hole systems. In this paper, we predict the distortion to the OH luminosity function that a magnification bias will inflict, and in turn, predict the distortion on the OH megamaser number counts as a function of redshift. We identify spectral flux density thresholds that will enable efficient lensed OH megamaser selection in large spectral line surveys with MeerKAT and SKA. The surface density of lensed galaxies that could be discovered in this way is a strong function of the redshift evolution of the OH megamaser luminosity function, with predictions as high as $\sim$1 lensed OH source per square degree at high redshifts ($z \gtrsim 1$) for anticipated SKA spectral line survey designs. This could enable efficient selection of some of the most highly-obscured galaxies in the universe. This high-redshift selection efficiency, in combination with the large survey speed of the SKA at $\lesssim$1 GHz frequencies and the high magnifications possible with compact OH emission regions ($\mu_{\rm OH} \gg 10$), will enable a transformational view of OH in the universe.


INTRODUCTION
Hydroxyl (OH) megamasers are luminous, extragalactic maser sources.As in their Galactic counterparts, the emission in OH megamasers (OHMs) is dominated by the masing lines at 1665 and 1667 MHz, with much weaker satellite lines at 1612 and 1720 MHz.However, unlike Galactic OH masers, the emission line at 1667 MHz is stronger than the emission line at 1665 MHz, with a typical ratio of the line strengths of 9:5 in local thermodynamic equilibrium (Lo 2005, and references therein).Additionally, due to Doppler broadening in massive, often merging galactic systems, OHMs have significantly larger line widths (∼100-1000 km s −1 ; Darling 2005) than Galactic OH masers.When compared to galaxy scale emission components, OHMs are compact sources, with sizes of order ∼100 pc revealed by high-resolution radio imaging (e.g.Pihlström et al. 2001;Rovilos et al. 2003;Lo 2005).
Because OHMs require luminous far infrared (IR) radiation to maintain the population inversion that can produce stimulated emission (Lockett & Elitzur 2008), they are typically found in the nuclear regions of luminous and ultra-luminous infrared galaxies (LIRGs and ULIRGs), many of which are also major merger systems.Indeed, the integrated OH line luminosity of OHMs is strongly correlated with the far-IR luminosity of the host galaxy, where the correlation follows ★ E-mail: charissa@imago-web.co.za a super-linear power law, (1) (Baan et al. 1992;Darling & Giovanelli 2002a;Glowacki et al. 2022;Wang et al. 2023).However, not all (U)LIRGs host OHMs.The fraction of (U)LIRGs that host OHMs is a strong function of the IR luminosity and increases to about one in three for ULIRGs (Lo 2005).
Using HCN and CO observations of a sample of OHMs, Darling (2007) show that high molecular gas densities ( H 2 ≳ 10 4 cm −3 ) are required in addition to strong far-IR radiation for the production of OHMs.Since OHMs reside in LIRGs and ULIRGs, the number density of OHMs should evolve strongly with redshift as the number of (U)LIRGs increases with redshift due to the ≳ 1 dex increase in the cosmic star formation rate density (e.g.Madau & Dickinson 2014).Furthermore, the number density of OHMs should also increase with redshift due to the expected increase in the molecular gas fraction and density at high redshift (Obreschkow & Rawlings 2009).
OHMs have been demonstrated to be useful tracers for classes of galaxies that are important for understanding several aspects of galaxy evolution.Because they are associated with strong IR radiation and high molecular gas densities, they trace extreme star formation (Darling 2007;Lockett & Elitzur 2008).Additionally, since they seem to be produced primarily in major galaxy mergers, they ought to provide an independent probe of the galaxy merger rate (Briggs 1998).They are also likely signposts for dual or binary AGN in what are typically obscured environments, especially at higher redshifts.Some of the first searches for OHMs targeted luminous IR sources with strong radio continuum (Baan et al. 1985).The largest systematic search for OHMs, to date, was conducted with the Arecibo OH megamaser survey (Darling & Giovanelli 2002a).This survey targeted galaxies selected from the IRAS Point Source Catalogue Redshift Survey (PSCz; Saunders et al. 2000) that were within the declination range 0 • <  < 37 • and that fell within the redshift range 0.1 ≤  ≤ 0.23.Due to the flux limit of the PSCz and the lower redshift limit of Arecibo, the target galaxies were primarily luminous IR galaxies (LIRGS) with  FIR ≥ 10 11.4 L ⊙ .The Arecibo survey detected 52 new OHMs, almost doubling the number of known OHMs at the time.
Despite this major step forward by the Arecibo OH Survey over two decades ago, the number of known OHMs today is still at a similar level (for more up to date catalogues of OHMs, see Zhang et al. (2014) and Sotnikova et al. (2022)).Upcoming wide area spectral surveys on MeerKAT and the SKA1-Mid will mark a step-change in the population size.Roberts et al. (2021) predicted that the LADUMA survey (Blyth et al. 2016) alone could detect ∼ 80 OHMs in its single pointing, while Glowacki et al. (2022) reported an OHM detection at  = 0.52 in the LADUMA survey, and Jarvis et al. (2023) report the detection of a  = 0.71 OHM in the MIGHTEE survey, making these the highest redshift detections to date by factors of ∼ 2 − 3, respectively.It is worth highlighting that the latter, MIGHTEE  = 0.71 detection, is likely to be strongly lensed, with a magnification factor of  OH ∼ 3.This highlights the opportunities that the SKA1-Mid and its pathfinders, including MeerKAT, will open up for studying both high-redshift and low-luminosity OHMs, which will further our understanding of OHMs and provide useful tracers for understanding aspects of galaxy evolution.
As the upper redshift limit of cosmic OH is increased by more sensitive instruments, so too is the probability of detecting gravitationally lensed OH megamasers, analogous to what SKA and its precursors/pathfinder will do for the neutral hydrogen line (e.g.Deane et al. 2015Deane et al. , 2016;;Blecher et al. 2019).A targeted approach to lensed OHMs is described in Manamela et al. (in prep.), while in this work we describe the untargeted, statistical approach.This relies on the increase in the observed OHM number density due to magnification bias.This distortion to the luminosity function by the sub-population of lensed objects would be most noticeable at high luminosity values because of the exponential decline in the number density at this end, assuming that the luminosity function follows the form of a Schechter function.This increase in the number counts can be exploited to find a peak flux density threshold above which the lensed population dominates over the unlensed population, thereby enabling an efficient approach to a challenging but scientifically rich objective.This selection technique has been used very successfully to discover strong gravitational lenses in far-IR and sub-mm surveys undertaken with Herschel and the South Pole Telescope (e.g.Negrello et al. 2010;Vieira et al. 2013;Wardlow et al. 2013;Negrello et al. 2017).
In this paper, we investigate the extent to which a statistical selection approach could be applied to MeerKAT and SKA1-Mid spectral line surveys to identify lensed OHMs.In order to apply the statistical selection approach to OHMs, the OH luminosity function has to be extrapolated to a larger range of luminosities for which it was measured.Section 2 briefly reviews the relevant surveys that are ongoing with MeerKAT and that are planned for the SKA1-Mid.Section 3 describes the current constraints on the OH luminosity function and investigates suitable models to extrapolate the function to a larger range in OH luminosity and redshift.Section 4 reviews the necessary steps for calculating the integrated source counts for OHMs, while Section 5 presents and discusses the results of the lens selec-tion approach.Finally, Section 6 presents the conclusions and gives some perspectives on the future outlook for OH studies.Unless otherwise stated, we assume a Planck 2018 cosmological model (Planck Collaboration et al. 2020).

OVERVIEW OF MEERKAT AND SKA SURVEYS
Next-generation radio telescopes, such as the Square Kilometre Array (SKA) and its precursors and pathfinders will have dramatically improved instantaneous sensitivity, bandwidth and field of view which will enable the advance of spectral line studies in both H i and OH.This section describes the surveys that are currently being undertaken with MeerKAT and that have been proposed for the SKA1-Mid.In particular, we detail the expected sensitivity of these surveys.
There are two large survey projects currently being undertaken with MeerKAT that are relevant for spectral line searches for OH emission lines.The first is the MIGHTEE survey (Jarvis et al. 2016;Maddox et al. 2021) that will survey a sky area of 20 deg 2 over four fields at L-band and a smaller region in S-band.Each pointing will have an integration time of ∼ 16 hrs and the survey should reach a sensitivity of ≲ 100 Jy in a 209 kHz channel.The second is the LADUMA survey (Blyth et al. 2016) that will spend ∼ 300 hrs covering a single pointing on Chandra Deep Field South in L-band and another ∼ 3000 hrs on the same field in UHF-band.
A large part of the observing time on SKA1-Mid will be devoted to large survey projects, as has been the case with MeerKAT and ASKAP, the two SKA1-Mid precursors.Staveley-Smith & Oosterloo (2015) outline three prospective tiered surveys, each of a 1000 hours with survey areas ranging from 400 deg 2 to 1 deg 2 .Additionally, Staveley-Smith & Oosterloo (2015) discuss commensal surveys that could still be useful for spectral line science but that could have up to 10 000 hours of observing time, covering an area of up to of order Σ ∼  sr.
In order to calculate the anticipated sensitivity of the SKA1-Mid surveys, we used the estimated  eff / sys values for the SKA1-Mid dishes which are provided by Braun et al. (2019).Additionally, we used the measured mean System Equivalent Flux Density (SEFD) values for a single MeerKAT antenna available on the MeerKAT specifications page to account for the system noise of the MeerKAT dishes at the frequency intervals at which they will contribute to the SKA1-Mid array.From these SEFD values, the sensitivity of the SKA1-Mid for a single pointing can be estimated from the radiometer equation.The estimated sensitivities of the SKA1-Mid surveys, proposed by Staveley-Smith & Oosterloo (2015), are summarised in Table 1.While the survey strategies will only be finalized in the future, these surveys serve as useful, indicative, reference points as we consider searching for lensed OH megamasers in the SKA1-Mid and MeerKAT surveys.

THE OH MEGAMASER LUMINOSITY FUNCTION
The OH luminosity function, Θ( OH ), is the comoving number density of OHMs as a function of their OH luminosity.However, it is often more convenient to express the luminosity function as the comoving number density per logarithmic luminosity interval rather than linear luminosity interval.In this case, the luminosity function is denoted as Φ( OH ).The OH luminosity function can be calculated either directly from OHM number counts, or indirectly from the FIR luminosity function, assuming a  OH - FIR correlation and an OHM abundance in FIR-luminous galaxies.Both approaches were used in  (2015).The field of view is calculated at the maximum frequency quoted for the survey.Column 7 gives the effective integration time per pointing calculated from the number of pointings and the total number of hours for the survey.The flux density sensitivity is calculated for a single pointing at the central frequency of the survey and for a rest-frame velocity width of 200 km s −1 which ranges between 950 and 470 kHz for H i between redshifts of 0 and 1.
the early work on the OH luminosity function; Baan (1991) (2) Later, Roberts et al. (2021) used a Markov chain Monte Carlo approach in a re-analysis of the OH luminosity function using a power law model and the same data as Darling & Giovanelli (2002b).They obtained the following parameter estimates (scaled to the same cosmological model as used in Darling & Giovanelli (2002b)), As mentioned, in their analyses both Darling & Giovanelli (2002b) and Roberts et al. (2021) only used data points within the luminosity range 2.2 ≤ log( OH ) ≤ 3.8.The data outside of this range were excluded as they argued that the OHM detections in these luminosity bins did not uniformly sample the available cosmological volumes.To test for uniformity, Darling & Giovanelli (2002b) used the ⟨/ a ⟩ test from Schmidt (1968) where a uniformly distributed sample has a ⟨/ a ⟩ value of 0.5.Here,  is the survey volume at the redshift of each detected OHM and  a (also referred to as  max by other authors) is the maximum available volume for each OHM given the survey sensitivity and the luminosity of the OHM. Figure 1 shows the ⟨/ a ⟩ values and uncertainties estimated by Darling & Giovanelli (2002b).The points marked by crosses indicate the luminosity bins centred on log( OH ) = 1.6, 2.0, and 4.0 that were excluded in the analyses by Darling & Giovanelli (2002b) and Roberts et al. (2021).The ⟨/ a ⟩ values for the luminosity bins centred on log( OH ) = 1.6, 2.0 and 4.0 are consistent with a value of 0.5 within 1.2; however, each of the bins centred on log( OH ) = 1.6, and 4.0 had only one detection which means that the uncertainties on the ⟨/ a ⟩ values for these bins are on order of 1.
With well-sampled number counts, luminosity or mass functions of galaxy properties are typically well described by a Schechter Average  to  a ratios for each luminosity bin.The ⟨/ a ⟩ ratio is used to test for uniform sampling; a value of 0.5 indicates that the volume available within a luminosity bin is well-sampled.The coloured intervals indicate the different data sets that are used in the model selection tests in Section 3.1, where the Data Set A shows the data points that were used in the analyses performed by Darling & Giovanelli (2002b) and Roberts et al. (2021).Figure reproduced from the data points reported in Darling & Giovanelli (2002b).
function (Schechter 1976).While the OH luminosity function and its evolution are poorly constrained at present, upcoming observations with SKA precursors/pathfinders will significantly increase the number of known OH megamasers across a wider range of luminosity and redshift.
In what follows, we perform Bayesian model selection between a power law and Schechter function model.Taking a Bayesian approach, our analysis includes all the detections from the Arecibo OHM survey.However, we also perform this same model selection with subsets of the data to compare with previous work and the effect of excluding data.We also present a typical parametrization of the redshift evolution of the OH luminosity function.

A Bayesian approach to modelling the OH luminosity function
In this section, we model the OH luminosity function with both a power law and a Schechter function using three combinations of the data presented in Darling & Giovanelli (2002b).We used PyMulti-Nest (Buchner et al. 2014) to constrain the model parameters for both a power law model and a Schechter model, as well as to compute the Bayesian evidence for each model required for model selection.The 2-parameter power law model is given by while the 3-parameter Schechter model is given by (5) Since Darling & Giovanelli (2002b) did not include all data points, we carried out three investigations using different subsets of data points (labelled Data Sets A, B, and C, see Fig. 1).First, we modelled the luminosity function using the same data points as used by Darling & Giovanelli (2002b), i.e. the luminosity bins centred on log( OH ) = 2.4, 2.8, 3.2 and 3.8.This is referred to as Data Set A in Table 2 and Figure 1.Second, in Data Set B we extend the data subset to include the high luminosity point at log( OH ) = 4.0 and repeat the same modelling approach.Third, Data Set C includes all the luminosity bins and we again perform the same modelling with this data set.
Figure 1 shows the different subsets of data points used in the three investigations.
The results of the parameter estimation and model selection are presented in Table 2.It shows, for each of the three investigations, the Bayes factor comparing the evidence for the Schechter model to the evidence for a power law model, as well as the median values of the posterior probability distributions for each parameter in the two models.In the first two investigations that exclude some data points, there is only weak or inconclusive support for a power law model over a Schechter function, using the Jeffrey's scale interpretation (see e.g.Trotta 2017).When considering all the data points, the Schechter model is preferred, although the Bayes factor shows inconclusive support for the Schechter model in this case.
The median posterior models for the three investigations are shown in Figure 2. It is interesting to note that for Data Set A, the nested sampling approach produces a similar power law to that found by Darling & Giovanelli (2002b).Similarly, for Data Set C, the nested sampling approach produces a power law that is consistent to within the uncertainties with the power law produced by the MCMC approach used by Roberts et al. (2021).The fact that the power law model from Data Set C agrees with the model found by Roberts et al. (2021) can be viewed as argument against the exclusion of the data points.The differences in the derived power law parameters for Data Set A are expected to be due to differences in the MCMC algorithm used in Roberts et al. (2021) and the nested sampling algorithm used in this work.The latter has been shown to be a more robust, albeit computationally expensive approach (Skilling 2004).Given that nested sampling computes the Bayesian evidence while MCMC does not, we would expect superior results with our approach.
Overall, the results show that the model selection between a power law model and a Schechter model is inconclusive given the current data.There is weak or inconclusive support that a power law model is preferred over the narrower selected luminosity range used in Darling & Giovanelli (2002b).However, it is unclear that this model can be extrapolated to either lower or higher luminosities, with the latter being a key point of interest in this paper.
In the remainder of this work, the OH luminosity function is taken to be modelled by a Schechter function, with parameter estimation based on all the available data in all luminosity bins.There are several reasons that motivate this choice.First, within a Bayesian framework, all data points should be included even where there are large uncertainties.Second, when all the data points are considered, a Schechter model is marginally favoured even though the Bayes factor is inconclusive.Third, a wide range of luminosity functions of all kinds are generally well modelled by a Schechter function and it is reasonable to expect that when the OH luminosity function is measured over a larger range of luminosities, it will be no different.Fourth, if the  OH - FIR correlation holds to higher luminosities, we expect that the OH luminosity function will follow the far-IR luminosity function which is well modelled by a Schechter function.
We note that the lens selection approach we utilise in this paper is sensitive to the exponential steepening of the luminosity function at high luminosities, which is poorly constrained at this point.However, this will become more tightly constrained as more high redshift OHMs are detected with MeerKAT and the SKA1-Mid.In anticipation thereof, we consider what the possibilities are for detecting both lensed and unlensed OHMs en route to a new era in OHM studies.
The Schechter parameters for the OH luminosity function that are used in this work, converted to the Planck 2018 cosmology (Planck Collaboration et al. 2020) are log  * = −7.36 ± 0.31, log  * OH = 3.68 ± 0.52 and  = −1.18± 0.26.

Parametrization of the redshift evolution
As discussed, the number density of OHMs is expected to evolve strongly with redshift for several reasons (e.g.Briggs 1998;Darling & Giovanelli 2002b).First, since OHMs have been observed to be associated with luminous far-IR radiation, the redshift evolution of the OH luminosity function is expected to be strongly influenced by the evolution of the number density of (U)LIRGs.Although (U)LIRGs in the local universe are predominantly major merger systems, at high redshift the (U)LIRG population also contains normal star forming galaxies (e.g.Daddi et al. 2010).Thus, the number density of (U)LIRGs should evolve with redshift due to both the increasing star formation activity (e.g.Madau & Dickinson 2014) and the evolution of the merger rate.However, the merger rate of galaxies is still highly uncertain (see e.g.Mundy et al. 2017).Some studies find that for massive galaxies, the merger rate increases (e.g.Bluck et al. 2012), while others find that the merger rate is constant or even decreases (e.g.Williams et al. 2011;Newman et al. 2012) with selection biases seen to be a major hurdle to reconcile.Second, since OHMs are also associated with dense molecular gas regions, the redshift evolution of the OH luminosity function should also be influenced by the increased mid-plane pressure and dust temperatures of molecular gas in galaxies at higher redshifts (Darling 2007;Lockett & Elitzur 2008).The molecular gas is expected to have a higher density at higher redshifts due to the galaxy size evolution (Gunn & Gott 1972) and due to the predicted increase in the H 2 /H i mass ratio with redshift (Obreschkow & Rawlings 2009).
Although the redshift evolution of the number density of (U)LIRGs and the increasing density of molecular gas should strongly influence the redshift evolution of the OH luminosity function, at this point it is difficult to quantify this evolution of the OH luminosity function.Therefore, in this work we parametrize the redshift evolution by including a factor of (1 + )  OH , as is commonly used to parametrize the redshift evolution of a variety of luminosity or mass functions (e.g. the H i mass function, Pan et al. 2020).The OH luminosity function is then given by where the exponent,  OH , of the (1 + ) term is referred to as the evolution parameter in this work.We assume this evolution parameter incorporates a number of relevant, introduced effects, including the contribution from the merger rate, the increase in the abundance of (U)LIRGs, and the increase in the molecular gas density mid-plane pressure.Because this evolution is largely unknown at this stage, we investigate values of  OH = 0, 2, 4 and 6 in order to cover a wide range of evolution scenarios.This is similar to the approach taken by Darling & Giovanelli (2002b) to model what they isolated as the merger rate.These different redshift evolution scenarios will be constrained as observations from MeerKAT and SKA1-Mid become available and increase the number of known OHMs to sample sizes with sufficient statistical power, however, for this paper they provide the necessary functional form to explore the lensed OH number counts at high redshift.

Lensing probability
The effect of gravitational lensing on the observed luminosity function depends on the probability that a background source at a given redshift is lensed and by a given magnification factor.Here, we follow the method outlined in Perrotta et al. (2002) in order to calculate this probability, which depends on the mass profile of the foreground lens, on the mass and redshift distribution of the foreground lenses, the redshift of the background source, and the cosmological parameters which we fix to the values measured by Planck Collaboration et al. (2020).Following Perrotta et al. (2002) and Wardlow et al. (2013), we assume that the foreground lenses are dark matter haloes.Wardlow et al. (2013) found that their predictions did not have the necessary precision to discriminate between the singular isothermal sphere (SIS) and Navarro-Frenk-White (NFW) density profiles.Additionally, the lensing quantities of an SIS density profile are simpler than the lensing quantities of an NFW density profile.For these reasons, we adopt SIS profiles here.For an individual lens, the mass surface density will be enhanced near its nucleus by a baryonic component, however, following other authors, we focus on the dark matter lensing potential for this statistical, rather than individual lens analysis.
The probability that a background source is lensed depends on the lens cross-section which is the area in the source plane within which the magnification is larger than a given value,  min , and is given by the equation (Lima et al. 2010), where Ω is the lens cross-section,  L is the redshift of the foreground lens,  S is the redshift of the background source,  vir is the virial mass of the foreground lens and  is the source position.In the case of an SIS profile the cross-section becomes Given the lens cross-section, the probability that a background source is lensed with a magnification factor greater than  is given by the fraction of the area of the source sphere where the magnification is greater than  (Perrotta et al. 2002): Here  A ( S ) is the angular diameter distance at the redshift of the source, d/d is the comoving volume element per unit redshift and solid angle, Ω(,  L ,  S ,  vir ) is the lens cross-section, and d/d vir d L is the comoving number density of the lenses, which, following Perrotta et al. (2002) is given by the Sheth and Tormen dark matter halo mass function (Sheth et al. 2001) and is implemented in the Halomod1 package (Murray et al. 2013(Murray et al. , 2021)).This expression for the probability is only valid in the limit of non-overlapping cross-sections where  ≪ 1.However, strong lensing has a very small probability since it requires the source and the lens to be very closely aligned ( ≲ 1 arcsec) and it is rare that a source is lensed by more than one foreground mass.That said, massive galaxies are typically clustered; however, we do not take this into consideration since systematic modelling uncertainties (e.g.OHM source size) far outweigh this second-order effect.
The magnification probability distribution is given by the differential probability,

Integrated counts of OH sources
The effect of lensing on the observed luminosity function can then be calculated from the magnification probability distribution.Since the magnification is equal to the ratio of the observed flux density to the actual flux density, the observed or apparent luminosity is So, if all the sources were lensed by a factor of , the observed number density will be related to the intrinsic number density by (Pei 1995) where Θ( OH ) is the number density calculated for intervals of  OH and the factor of 1/ takes into account the amount by which the magnification changes the width of the luminosity bins.However, it is often more convenient to calculate the number density on intervals of log  OH , as is done in Equation 5.In this case, the observed number density, Φ ′ ( OH ) is related to the intrinsic number density by since the magnification factor does not affect the width of the logarithmic luminosity intervals.
When the sources are lensed by a range of magnification factors with a probability distribution given by Equation 10, then the observed number density is given by Here, the value of  min is restricted to the strong lensing regime ( min > 2) by the assumption in Equation 9that the lens crosssections do not overlap, while the value of  max is predominantly limited by the maximal solid angle of the source.
The integrated source counts of OH megamasers are given by where  min is the integrated luminosity of the OH spectral line corresponding to a given peak flux density and d/d is the differential comoving volume.In the above equation, if the integrated source counts are calculated for the unlensed population, Φ( OH , ) is the OH luminosity function.On the other hand, if the source counts are calculated for the lensed population, Φ ′ ( OH , ) is the OH luminosity function modified by the magnification bias, as given in Equation 14.Thus, in order to calculate the integrated source counts, we need to specify the maximum magnification, the integrated luminosity corresponding to a given peak flux density, and the redshift interval.
As discussed, OHMs are typically ≲100 pc in extent (Rovilos et al. 2003;Lo 2005), which is considerably smaller than most emission components in their host galaxies.This smaller spatial extent should result in much higher magnification factors compared to the magnification factors considered in, for example the H i case, particularly for galaxy-scale Einstein radii (e.g.Deane et al. 2015;Blecher et al. 2019).In order to take into account this expectation that OHMs can have high magnifications ( ≫ 10) and to be consistent with the modelled magnifications seen in similarly sized emission components in the literature (e.g.H ii regions, Kneib & Natarajan 2011), maximum magnification factors in the range  max = 10-100 are considered here.Figure 3 shows the effect of the magnification bias on the OH luminosity function for different values of the redshift evolution parameter and for maximum magnification factors in the range 10-100.These plots indicate that it is only at high luminosities that the OH luminosity function distorted by the magnification bias for a maximum magnification of 100 differs significantly from the OH luminosity function with a maximum magnification of 10.At the point where the distorted OH luminosity function intersects the unlensed OH luminosity function, the results for the two values of the maximum magnification are similar.A similar result is seen in the integrated source counts that are plotted in Figure 4.This could change as the OH luminosity function constraints improve, however, they show that our predictions are relatively insensitive to this assumption with current models.
The integrated luminosity of the OH spectral line is related to the peak flux density through the integrated flux and the definition of the luminosity distance.Assuming that the spectral line has a rectangular profile, the integrated flux is related to the peak flux density as where  OH is the rest-frequency of the OH spectral line,  is the redshift of the OH megamaser and Δ rest is the rest-frame velocity width.Constraints on the velocity width of higher redshift OHMs do not exist, but for consistency with previous predictions by Darling & Giovanelli (2002a) a constant line width of 150 km s −1 is assumed in this paper, however, note that much broader profiles are seen in the low redshift universe.Then, from the definition of the luminosity distance, the integrated flux is related to the integrated luminosity by Combining Equations 16 and 17, the luminosity can be expressed in terms of the peak flux density as The SKA1-Mid Band 1 receivers will extend down to a frequency of 350 MHz, so in principle, OHMs could be detected out to a redshift of  ∼ 3.7.Hence, the integrated number counts in this work are considered for redshifts up to 3.7.This redshift range is divided into smaller redshift intervals and the number counts are calculated on each interval.It is useful to consider smaller redshift intervals in order to investigate how the number counts change with redshift, even though the number counts are smaller, since a smaller redshift interval corresponds to a smaller volume (at a given central redshift).Here, we consider redshift intervals of width Δ = 0.10 and 0.37.This also accounts for practical data processing considerations, as well as the fact that RFI is often clustered in specific windows.

RESULTS AND DISCUSSION
This section presents the results of applying this selection approach to OHMs.As discussed in the previous sections, the integrated counts of the lensed and unlensed populations are calculated for redshift intervals of width Δ = 0.1 and 0.37 within the range 0 ≤  ≤ 3.7 and for values of the evolution parameter of  OH = 0, 2, 4, and 6.This results in 47 redshift intervals for each value of the evolution parameter.We do this in order to provide a practical sense of how these searches can be carried out, while also accounting for our ignorance on  OH and the OH luminosity function.In Figure 4, we show the lensed and unlensed integrated counts for each value of  OH in an example redshift bin of 1.6 <  < 1.7, which is accessible by observed frame frequencies that lie just inside the MeerKAT UHF band ( obs ∼ 630 MHz).Here, the grey shading indicates the region where the lensed number counts are equal to or exceed the unlensed number counts.For the sake of clarity, the point of intersection between the lensed integrated source counts (for  max = 100) and the unlensed integrated source counts (indicated by the intersection of the grey dashed lines) are referred to as the source count equality points.A sample selected at the flux density threshold indicated by the vertical dashed line would contain lensed and unlensed sources in a ratio of 1:1.
The OHM surface density is a strong function of the evolution parameter and is expected to increase by orders of magnitude in even the more conservative estimates.OHM searches with MeerKAT and the SKA1-Mid should help to constrain the observed OHM surface density as a function of redshift, which will, in turn, enable the redshift evolution of the OHM number density to be constrained.
The extent to which this method could be useful for selecting lensed OHMs in the SKA1-Mid surveys can be explored by plotting the source count equality points as a function of the central redshift and the width of the redshift interval for the different values of  OH .These source count equality points are shown in Figure 5 where the vertical lines indicate the approximate sensitivities of the SKA1-Mid medium wide survey, the MIGHTEE survey and the LADUMA survey (5- ≃ 80 Jy beam −1 , 100 Jy beam −1 , 50 Jy beam −1 , respectively, assuming a rest-frame velocity width of 150 km s −1 These source count equality points occur at flux densities that are easily accessible to the SKA1-Mid medium wide survey and the MIGHTEE survey, which has a similar depth to the SKA1-Mid medium wide survey.Since the LADUMA survey only covers an area of ∼ 3 deg 2 at 580 MHz, its area does not make it an optimal survey to discover lensed OHMs, while the wider areas of the SKA1-Mid medium wide survey and the MIGHTEE survey make them potentially more interesting surveys to consider here.The surface densities of the lensed sources at these intersection points depend strongly on the evolution parameter, but are small (< 10 −4 deg −1 ) at redshifts of  ≲ 1 for all the evolution scenarios.In the case where there is no redshift evolution (i.e. OH = 0), the surface density of lensed OHMs at the source count equality points is close to or below the limit of one lensed source in the whole sky for all redshifts.It only starts to exceed this limit at high redshifts ( ≳ 2.5) and when integrating over the wider redshift interval (Δ = 0.37).However, for  OH = 6, the surface density increases to 1 lensed source per 1 deg 2 at high redshifts.This indicates that if the redshift evolution of the OH luminosity function is relatively strong,  OH ≥ 4, this could be a promising method for selecting lensed OHMs, with lensed OHMs being detectable out to redshifts of  ∼ 3.5 at a relatively high surface density.On the other hand, if the OH luminosity function does not evolve strongly with redshift, this selection method becomes limited by the low surface density of the lensed OHMs; however, non-detections of any OH sources, whether lensed or unlensed, will provide joint constraints on the evolution parameter as well as the high end of the OH luminosity function.The above is, of course, for the most pessimistic scenario, where no multi-wavelength information is taken into account.Doing so will likely make a significant enhancement to this selection technique and is the subject of future work.In the current paper, we simply explore indicative levels of contaminant removal.

Lensed OHM surface densities at higher contaminant ratios
Multi-wavelength information, such as that from large optical and IR surveys, will greatly assist in identifying lensed OHM candidates in a given sample.If the multi-wavelength selection is efficient, it should be possible to select lensed OHMs from samples that contain lensed to unlensed OHMs in a ratio much less than 1:1, greatly enhancing yield with this technique.We, therefore, investigate how many lensed sources per square degree are expected at lower, indicative, ratios.Figures 6 show the surface densities of the lensed OHMs at flux densities where the unlensed source count is 10 times the lensed source count.We use this value as an indicative exploration of what is possible.These results suggest that for a scenario with no evolution, which is highly unexpected, even a contaminant removal accuracy of 1 in 10 will result in a very low sky density (less than 1 lensed source per 1000 deg 2 ) and, hence, have low scientific yield.For stronger evolution parameters ( OH = 2, 4 or 6), being able to select 1 lensed source out of 10 unlensed sources would naturally increase the surface density of the lensed OHMs.For example, for  OH = 6 the surface density increases to ∼ 3 lensed sources per 1 deg 2 , at the highest redshifts.By way of comparison, the lens selection in the Herschel ATLAS survey obtained a lens surface density of 0.13 lens candidates per 1 deg 2 (Negrello et al. 2017).This implies that the surface density of lensed OHMs could be significantly larger than the Herschel far-IR selection for galaxies in a very similar redshift window centred on cosmic noon (1 ≲  ≲ 3).Additionally, selecting lensed OHMs has a significant advantage of immediate spectroscopic confirmation, making follow up observations far more efficient.

CONCLUSION AND FUTURE OUTLOOK
Currently, the number of known OHMs is relatively small (on order of 100) and limited to low redshifts.However, as tracers of extreme star formation and major mergers, these objects will provide useful perspectives on galaxy evolution processes, especially at higher redshifts, as well as providing signposts for dual/binary AGN in obscured environments.The SKA1-Mid and its precursors, especially MeerKAT, will increase the number of OHM detections and advance our understanding of these objects through multi-faceted, multi-wavelength studies of the resultant OHM samples out to significantly larger cosmological distances.In anticipation of these developments, this paper investigates the possibility of selecting lensed OHMs in the ongoing/proposed MeerKAT and SKA1-Mid wide area spectral line surveys.
In the first instance, the OH luminosity function is constrained to both higher and lower luminosities than is carried out in Darling & Giovanelli (2002b) and Roberts et al. (2021).Although the number of current OHM detections are few, using a Bayesian framework with a phenomenological expectation that a Schechter function is an appropriate model for sufficient number counts, we model the OH luminosity function using measurements from Darling & Giovanelli (2002b) with a Schechter function.This chosen model is supported by Bayesian model selection; however, this approach and the resultant luminosity function parameter constraints will be tested as more OHM detections become available in the near future.
Following our Bayesian parameter estimation and model selection, we present the results of the lens selection as applied to OHMs.The prospect of detecting lensed OHMs in wide area spectral line surveys is a strong function of the evolution parameter,  OH , which we define to include a wide range of contributing factors, including the major merger rate, the evolution of the IR source counts and the increase in the density of molecular gas with redshift, amongst others.For no or weak evolution, which is highly unlikely, this lens selection method is not effective as the surface density of the lensed OHMs is very small (< 10 −4 deg −2 ), a scenario easily ruled out by a small sample of detections from a targeted lensed OHM survey (Manamela et al., in prep.).For strong  OH evolution, this selection method becomes promising even without any other information as the lensed OHMs should reach a surface density of 1 in 1 deg 2 at the highest redshifts surveyed with SKA1-Mid Band 1.This surface density can be improved up to ∼ 3 lensed OHMs per 1 deg 2 if ancillary multi-wavelength information can be used to remove contaminants.
Clearly, there is great potential scientific yield, which motivates using sophisticated techniques, including machine learning methods, to remove unlensed contaminants, a subject of current research, with a relevant link reported in Roberts et al. (2021), where OHM are seen as the contaminants in separating H i and OHM sources in the LADUMA survey.These results illustrate how the different redshift evolution scenarios of the OH luminosity function can be tested in the near future as more OHMs are discovered in upcoming widearea spectral line surveys.We predict that our view of the OH in the universe is about to be transformed beyond previous predictions by the SKA and its precursors/pathfinders, through the power of strong gravitational lensing.

Figure 2 .
Figure2.Results from the parameter estimation of the OH luminosity function.The top panel shows the results for the parameter estimation using Data Set A, the middle panel shows the results for the parameter estimation using Data Set B, and the bottom panel shows the results for the parameter estimation using Data Set C. In the top panels, the points that are excluded from the parameter estimation are marked with crosses.In this figure, the same cosmological parameters are assumed as inDarling & Giovanelli (2002b) and the parameters from Roberts et al. (2021) are scaled accordingly.

Figure 3 .
Figure 3.Comparison of the OHM luminosity function modified by the magnification bias with the unlensed luminosity function for four values of the evolution parameter.The lensed sources are assumed to be at a redshift of 1.5.

Figure 5 .Figure 6 .
Figure 5.The source count equality points for all the redshift intervals considered, at the different values of the evolution parameter,  OH , as indicated on each panel.The centre of the redshift interval is indicated by the colour and the size of the points indicates the width of the redshift interval.The approximate sensitivities (assuming a rest-frame velocity width of 150 km s −1 ) for the MIGHTEE and LADUMA surveys of MeerKAT, as well as for the SKA1-Mid Medium Wide survey are shown by the vertical lines.The horizontal dashed line indicates the limiting source surface density of 1 source in the whole sky area.

Table 1 .
Flux density sensitivities of the SKA1-Mid surveys proposed in Staveley-Smith & Oosterloo