X-ray detection of the most extreme star-forming galaxies at the cosmic noon via strong lensing

Hyper-luminous infrared galaxies (HyLIRGs) are the most extreme star-forming systems observed in the early Universe, and their properties still elude comprehensive understanding. We have undertaken a large XMM-Newton observing program to probe the total accreting black hole population in three HyLIRGs at z = 2.12, 3.25, and 3.55, gravitationally lensed by foreground galaxies. Selected from the Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts (PASSAGES), these HyLIRGs have apparent infrared luminosities>E14 Lsun. Our observations revealed X-ray emission in each of them. PJ1336+49 appears to be dominated by high-mass X-ray binaries (HMXBs). Remarkably, the luminosity of this non-AGN X-ray emission exceeds by a factor of about three the value obtained by calibration with local galaxies with much lower star formation rates. This enhanced X-ray emission most likely highlights the efficacy of dynamical HMXB production within compact clusters, which is an important mode of star formation in HyLIRGs. The remaining two (PJ0116-24 and PJ1053+60) morphologically and spectrally exhibit a compact X-ray component in addition to the extended non-AGN X-ray emission, indicating the presence of Active Galactic Nuclei (AGNs). The AGN appears to be centrally located in the reconstructed source plane images of PJ0116-24, which manifests its star-forming activity predominantly within an extended galactic disk. In contrast, the AGN in the field of PJ1053+60 is projected 60 kpc away from the extreme star-forming galaxy and could be ejected from it. These results underline the synergistic potential of deep X-ray observations with strong lensing for the study of high-energy astrophysical phenomena in HyLIRGs.


INTRODUCTION
Hyperluminous infrared galaxies (HyLIRGs) are the brightest galaxies in the Universe and represent the most extreme version of dusty star-forming galaxies (DSFGs) (e.g.Irwin et al. 1998;Rowan-Robinson 2000;Blain et al. 2002;Casey et al. 2014).Identified by their large intrinsic infrared (IR) luminosity > 10 13 L ⊙ , HyLIRGs are typically found at z > ∼ 2 and are rare (a few per deg 2 according to Gruppioni et al. (2013)).If the IR luminosity is due to star forma-⋆ Contact e-mail:wqd@umass.edution (SF), then the corresponding SF rate (SFR) of such a HyLIRG is > ∼ 10 3 M ⊙ yr −1 .However, the IR luminosity could be significantly contributed by an Active Galactic Nucleus (AGN), as observed especially at z < 1 (e.g., Farrah et al. 2002;Ruiz et al. 2007).In this case, the galaxy may be detected as a so-called Hot Dust-Obscured Galaxy (Hot DOG; e.g., Wu et al. 2018).In any case, HyLIRGs are natural laboratories for studying the physical processes that govern the rapid growth and quenching of massive young galaxies and/or the coevolution between SF and supermassive black hole (SMBH) accretion in the early Universe (e.g.Blain et al. 2002;Casey et al. 2014).Since there are no local analogs to their large inferred molecular gas mass (up to 10 11 M ⊙ ; Harrington et al. (2018Harrington et al. ( , 2021))) and gas mass fraction (40-80%), as well as their extreme luminosity (e.g., Carilli & Walter 2013), a detailed multi-wavelength study of HyLIRGs is the only way to probe the underlying physical processes.

Studying HyLIRGs in the high-energy astrophysics context
High-energy astrophysical processes are thought to play a central role in regulating the formation and evolution of galaxies and even their large-scale environment.Such processes are mostly associated with SF and SMBH accretion, which are probably coordinated processes, especially during the early rapid co-evolution stage (e.g., Hopkins et al. 2008;Mancuso et al. 2016).However, due to the limited sensitivity and resolving power of existing X-ray observatories -a primary tool for studying high-energy phenomena, it remains unclear how this regulation works and whether relevant key scaling relations derived locally can be applied, for example, to high-z extreme star-forming galaxies evolving under very different conditions.Here we focus on two key questions: 1) How might the ratio of the non-AGN X-ray luminosity (dominated by high-mass X-ray binaries or HMXBs) to SFR, L X /SFR (hereafter X HMXB ), in HyLIRGs differ from that in more normal star-forming galaxies?2) How do AGNs co-evolve with SF in HyLIRGs?

Uncertainty in the X HMXB calibration
The calibration of the X HMXB factor is important for multiple reasons.First, if X HMXB is known, then one may use it to infer the SFR of a galaxy from its L X , which may be measured more reliably than other SF tracers, especially when dust obscuration is severe (e.g., Ghosh & White 2001).However, when the AGN contribution is potentially significant, L X gives only an upper limit to the SFR (e.g., Mineo et al. 2012Mineo et al. , 2014)).Second, one may infer the presence of an AGN in a galaxy if its L X is significantly greater than that estimated from the SFR (if independently measured), together with a known X HMXB .Third and probably more importantly, the determination of a systematic X HMXB variation with certain galaxy properties could provide deep insights into the formation process of non-AGN X-ray sources.In active star-forming galaxies, for example, L X is mainly determined by HMXBs, especially so-called ultra-luminous X-ray sources (ULXs) with L 0.5−10 keV .> ∼ 1 × 10 39 ergs s −1 (e.g., Kaaret 2014;Lehmer et al. 2021).Here we consider HMXBs to include intermediate-mass X-ray binaries with donor stars whose lifetimes are up to a few ×10 8 yrs (e.g., Poutanen et al. 2013;Hunt et al. 2023) -the characteristic starburst age of HyLIRGs.Each HMXB consists of a compact object (neutron star or black hole) and a donor star whose current mass is greater than 3 M ⊙ 1 .The black hole (BH) may have stellar mass (< 10 2 M ⊙ ) or intermediate mass (IMBH with the BH mass > ∼ 10 2 M ⊙ but < 10 5 M ⊙ ).
The probably most commonly used calibration is obtained by Mineo et al. (2012) primarily through a careful spatially resolved Xray study of star-forming galaxies within a distance of < 40 Mpc.In these galaxies, the HMXB contributions are statistically determined.The best-fit L X −SFR relation (Mineo et al. 2012) is L X (erg s −1 ) = 2.61 × 10 39 SFR(M ⊙ yr −1 ), (1) which is equivalent to a X HMXB calibration (Fig. 1).The scatter around this relation (or X HMXB calibration) is large (RMS = 0.43 dex), at least partly due to both the limited sampling of the star formation history and the variation in the number and luminosity of the brightest sources (Mineo et al. 2012).The SFRs of the local sample galaxies are all < 20 M ⊙ yr −1 and have a median value of only 5.3 M ⊙ yr −1 , leading to large stochastic uncertainties.For comparison, as in Mineo et al. (2012), Fig. 1 also includes the X HMXB data obtained for more distant galaxies for which only the integrated total X-ray luminosities are measured, extending the SFR range to higher values by about two orders of magnitude.
The inclusion of such X-ray-unresolved integrated measurements in the calibration of X HMXB would introduce substantial systematic effects (e.g., Mineo et al. 2014;Riccio et al. 2023).In fact, using only the unresolved sample galaxies yields a coefficient value of ∼ 1.5 higher than in Eq. 1.This statistical X-ray emission excess is also apparent for the data points with SFRs > 20 M ⊙ yr −1 in Fig. 1 and is expected because of the inclusion of unresolved X-ray contributions (e.g., from diffuse hot gas, LMXBs, and AGNs; Mineo et al. (2012)).The contamination from diffuse hot plasma is typically important at < ∼ 1.5 keV (e.g., Lehmer et al. 2022;Wang et al. 2021) and is thus sensitive to both the X-ray absorption and the actual energy band used to measure the X-ray emission.AGNs are present in ∼ 25 − 50% of galaxies with high SFRs (Kirkpatrick et al. 2017) and contribute more at higher energies.Such varying contamination remains difficult to quantify, but could contribute significantly to the RMS of the galaxies, as well as the bias toward an over-estimated X HMXB in such a calibration.
The use of a different IMF must also be considered.We here use the IR luminosity to SFR conversion that assumes the Kroupa IMF (Kennicutt 1998;Calzetti 2013).This IMF is more realistic than the Salpeter IMF assumed to get Eq. 1 when low-mass stars (down to 0.1 M ⊙ ) are also considered, which contribute the bulk of the stellar mass.To account for this difference, we apply a correction factor of 1.5 (Madau & Dickinson 2014) to the aforementioned calibration, resulting in the following Kroupa IMF-based calibration, X HMXB,r ≈ 3.9 × 10 39 erg s −1 /(M ⊙ yr −1 ), (2) which we use as the reference or unit for our measurement of X HMXB .The use of this calibration, instead of the one based on the total X-ray luminosity of a galaxy in the 0.5-8 keV range (Mineo et al. 2014), is appropriate here because we measure X HMXB or the apparent non-AGN luminosity µL X by fitting the spectra of our sample HyLIRGs with z > 2 in the rest-frame energy range above 1.5 keV ( § 4), where the diffuse hot plasma contamination can be neglected.The LMXB contribution should also be negligible in such extreme star-forming galaxies.
There are other systematics that affect existing X HMXB calibrations, most of which are based on the total X-ray luminosities of galaxies.Both theoretical argument and observational evidence suggest that X HMXB increases with decreasing metallicity and/or decreasing total age of stellar populations (e.g., Fragos et al. 2013;Kaaret 2014;Lehmer et al. 2021).In addition, the dynamical forma-   2012), normalized to their best-fit mean or to the coefficient of the L X −SFR relation (as given in Eq. 1) for those galaxies with resolved X-ray emission.Our X HMXB measurements for the three HyLIRGs (Table 4) with the SFR corrected for the lensing magnifications are included for comparison.
The required N H would be on the order of > ∼ 10 24 cm −2 .The compactness of starbursts in (U)LIRGs is also related to their extreme youth (with starburst ages less than a few Myr;  2023)), in which massive stars contribute strongly to the IR emission, but the number of HMXBs is small.Alternatively, the compactness could lead to the rapid destruction of clusters (e.g., Linden et al. 2021Linden et al. , 2023)).These effects are still difficult to quantify for (U)LIRGs.Furthermore, recent X-ray stacking analyses of "normal" galaxies up to z ∼ 2.4 (e.g. using data from the Chandra Deep Field and Chandra COSMOS surveys) suggest that X HMXB increases moderately (e.g, by a factor of ∼ 2) with increasing redshift and decreasing metallicity, although the confusion with low luminosity AGNs can be serious (e.g., Kaaret 2014;Lehmer et al. 2019;Fornasini et al. 2019Fornasini et al. , 2020;;Lehmer et al. 2021Lehmer et al. , 2022;;Gilbertson et al. 2022).These statistical and systematic effects on the existing calibrations illustrate the importance of a spatially resolved measurement of X HMXB , preferably in a harder energy band and with a sufficiently large total SFR.Incidentally, this measurement could also help us to determine how BH binaries of masses > ∼ 30−60 M ⊙ may be formed, which have been detected via their gravitational wave emission from coalescence (e.g., Abbott et al. 2020;Belczynski 2020).

Uncertainty in the AGN and HyLIRG co-evolution
To date, information on the presence of AGNs in HyLIRGs at z > ∼ 2 is very limited.While some HyLIRGs have been identified as Hot DOGs via their enhanced mid-IR emission observed in Wide-field Infrared Survey Explorer (WISE) bands (e.g., Wu et al. 2018), they probably represent only a special AGN accretion stage and do not tell a complete story about the SF and SMBH co-evolution in such galaxies.An intriguing example is the XMM-Newton detection of the unlensed HyLIRG, HATLAS J084933.4+021443 at z = 2.41 (Ivison et al. 2019).Its XMM-Newton spectrum suggests the presence of an AGN, with its high X-ray luminosity indicating a potential energy release that could rival the starburst.Interestingly, however, the apparent X-ray absorption (N H < ∼ 5 × 10 21 cm −2 ) suggests that this AGN is not deeply embedded.Because the XMM-Newton data provide little spatial information about the X-ray emission from this unlensed HyLIRG, its actual physical relationship to the AGN is not clear.Several scenarios remain to be tested.One possibility is that the AGN is observed through a cavity of the galaxy it excavates.Alternatively, if the AGN is not co-spatial with the starburst, it may have been ejected from the galaxy or has evolved into the naked quasar phase in a separate faint companion galaxy.

X-ray HyLIRGs through strong gravitational lensing foregrounds
We have been exploring the use of strong gravitational lensing -nature's magnifying glass -to overcome the observational limitations in observing dusty star-forming galaxies (DSFGs), especially HyLIRGs (e.g., Harrington et al. 2021;Berman et al. 2022;Kamieneski et al. 2023).Strong lensing provides a unique opportunity to amplify the flux of a galaxy, making it easier to detect.
In addition, the increased angular size potentially allows us to resolve its structure and separate different components, while an AGN that remains point-like in multiple images can then be relatively easily identified.Existing studies have already shown that HyLIRGs can be complex systems of multiple components: e.g., a dominant DSFG accompanied by separate galaxies, which may or may not be physically associated.Aided by the strong lensing, we can have a better chance to pinpoint the responsible component for an observed signature.
We report here initial results from a large XMM-Newton program observing three strongly lensed HyLIRGs at the redshifts 2.12, 3.25, and 3.55 (Table 1) with a total exposure of 526 ks (Table 2).The excellent sensitivity and broad energy coverage of XMM-Newton, combined with the intrinsic brightness and lensing magnification of HyLIRGs, provide us with a unique opportunity to probe the highenergy activity in these extraordinary galaxies.The observations are sensitive to HMXBs and AGNs in the ∼ 1.5 − 36 keV rest frame.In this energy range, both the diffuse emission contamination and the X-ray absorption of the interstellar medium (ISM) are substantially smaller than in the 0.5-1.5 keV band (e.g., Lehmer et al. 2022;Wang et al. 2021).We may reasonably assume that the HMXB emission statistically follows the rest-frame far-IR emission since both typically trace recent massive SF on similar spatial and temporal scales (down to ∼ 1 kpc and 10 7 yr; Persic & Rephaeli (2007); Mineo et al. (2012); Cochrane et al. (2019)).Under this assumption, both the observed X-ray and dust continuum emissions are expected to have similar mean magnification factors.We can then measure X HMXB directly in the image plane without the need for the magnification correction.We discuss how the assumption can be violated, especially in the extreme case of HyLIRGs, in § 5.2.We empirically measure the X HMXB factor in such individual galaxies, as well as find evidence for AGNs, to provide an important check on how the locally calibrated X HMXB,r may be applied to galaxies with the most extreme SFRs observed exclusively at high-z and to offer insights into high-energy processes in HyLIRGs.
The organization of the rest of the paper is as follows.We describe our selection of targets (including their key multi-wavelength properties) in § 2, the XMM-Newton observations and data analysis in § 3, and present the results in § 4. In § 5 we compare the results with those obtained for local star-forming galaxies, including (U)LIRGs, and discuss the origin of the enhanced X HMXB and the implications for the detection of the AGNs observed in our targets.§ 6 gives a summary of our results, conclusions, and near-future prospects.Additional multi-wavelength data reduction and modeling are detailed in Appendices A, B, C and D. All images presented in this paper have north up and east to the left.We assume the ΛCDM cosmology 2 with H o = 69.6 km s −1 Mpc −1 , Ω M = 0.286 and Ω Λ = 0.714.

TARGET DESCRIPTION
Our targets for the XMM-Newton observations (Table 1; Figs.2-5) are selected from the Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts (PASSAGES).This survey provides a sample of 30 gravitationally-lensed DSFGs representing the brightest IR galaxies observed in the Universe (Fig. 2 Cañameras et al. 2015;Harrington et al. 2016;Frye et al. 2019;Trombetti et al. 2021;Berman et al. 2022;Kamieneski et al. 2023).A wealth of (sub)arcsecond multi-wavelength data have been taken for these DSFGs, using facilities including the Hubble Space Telescope (HST), the Submillimeter Array (SMA), the Atacama Large Millimeter/submillimeter Array (ALMA), the Northern Extended Millimeter Array (NOEMA), and the Karl G. Jansky Very Large Array (VLA).These data have helped us to probe various stellar and interstellar components of the underlying DSFGs.
Although our intended measurements as presented here are largely based on the observed (apparent) emission from the lensed DSFGs, the ultimate understanding of the underlying astrophysics will still be related to their intrinsic properties.Detailed lens modeling has been performed for many of the PASSAGES DSFGs (e.g.Frye et al. 2019Frye et al. , 2023b;;Kamieneski et al. 2023), showing that most of them are moderately to highly magnified by µ = 2 − 28 (see also Bussmann et al. 2013;Geach et al. 2018;Cañameras et al. 2018).The lensing magnification depends on the location and hence the  Harrington et al. (2016)), but with the three XMM-Newton HyLIRG targets highlighted here with filled red circles.This plot mainly shows the PASSAGES DSFGs (blue dots Harrington et al. 2016;Berman et al. 2022), compared with other sources identified by the WISE (Tsai et al. 2015), South Pole Telescope (SPT) (Vieira et al. 2013;Weiß et al. 2013;Reuter et al. 2020), and Herschel (Bussmann et al. 2013;Wardlow et al. 2013;Bakx et al. 2018Bakx et al. , 2020;;Urquhart et al. 2022;Cox et al. 2023).
intensity distribution across a galaxy at a wavelength.The quoted mean magnification factors for our targets in Table 1 are best estimated with the ALMA CO(J=3-2) emission data for PJ0116-24 (comparable to that with the dust continuum emission and close to µ ≈ 16 with the HST 1.6-µm data; Liu et al. (2023)), with the VLA 6 GHz continuum data for PJ1336+49 (Kamieneski et al. 2023), and with the SMA 850 µm data for PJ1053+60 (Appendix C).Our targets are mostly intrinsically luminous [L IR ≈ 1 × 10 13 L ⊙ ], hence HyLIRGs, rivaling the luminosities of even the brightest known unlensed objects at their respective redshifts.De-lensed images show that SF in the DSFGs tends to be widely distributed with the intrinsic sizes of far-infrared continuum regions (R e = 1.7 − 4.3 kpc; median 3.0 kpc; Kamieneski et al. (2023, and references therein)).In contrast, low-z ULIRGs, which tend to form stars in compact nuclear regions, are often unresolved and have their core sizes < ∼ 1.5 kpc (e.g., Díaz-Santos et al. 2010).
Our main selection criteria for the XMM-Newton targets are as follows: 1) high apparent 1.1-mm flux > ∼ 60 mJy and SED-inferred far-IR luminosity > ∼ 10 14 L ⊙ (e.g., Fig. 2) to maximize detection of the X-ray emission associated with extreme SF; 2) no evidence for AGN in the radio-infrared spectral energy distribution (SED) analysis (Berman et al. 2022) to minimize its potential contamination, and 3) lensing mainly by a single galaxy or a small group of galaxies without any evidence for AGN (e.g., radio jets) to minimize X-ray contamination from them.Nevertheless, our selected three targets can still be complex systems with multiple galaxy components.We describe the main characteristics of the targets in Table 1 and below: PJ1336+49 (or G104.43+66.26):This HyLIRG has the smallest Einstein radius (2 ′′ diameter) and lensing amplification among our targets (Table 1).The lensed image of the galaxy is most vividly seen as a nearly complete ring in the VLA 6 GHz continuum and consists of several emission clumps larger than the synthesized beam (Fig. 3).Three doubly imaged families of these radio clumps are identified and labeled in the left panel of Fig. 3.These clumps are used to construct the lens model of the source (Kamieneski et al. 2023).The two relatively bright and unlabeled clumps between '1b' and '2b' are not used in the lens modeling; their counterparts on the east side of the ring are not clear and may be confused with the unresolved peaks '1a' and '2a'.The ring of the DSFG is also seen in the SMA 1.2-mm continuum image at ∼0.8 ′′ resolution (Fig. 3 right panel), further confirming the lensed star-forming regions (see Appendix A).The lensing is dominated by a foreground compact galaxy labeled as 'L1' in the HST/WFC3 1.6 µm image (Fig. 3 middle panel), which shows a faint counterpart of the lensed DSFG ring, more apparent after the foreground light is subtracted with GALFIT (Kamieneski et al. 2023;Lowenthal 2023).However, this counterpart is quantitatively very uncertain.As a result, the re-constructed source-plane image in the 1.6 µm band may be problematic and is presented here just for completeness.): This HyLIRG has the largest apparent SFR among our targets (Table 1; Berman et al. 2022;Liu et al. 2023).Its HST/WFC3 1.6 µm image shows a nearly complete double Einstein ring with an average diameter of 6 ′′ , while the VLA 6 GHz and the ALMA 1.1-mm continuum and CO images show two major resolved arcs, which lie mostly in the middle of the ring (Fig. 4; Kamieneski et al. 2023).The results from detailed modeling of the VLT/ERIS near-IR IFU data, as well as the ALMA CO data, have been reported by Liu et al. (2023).The detection of the Hα, Hβ , [NII], and [SII] lines shows a slightly supersolar gas-phase metallicity and a high mean Balmer decrement (observed Hα/Hβ flux ratio of 9.2 ± 1.2, corresponding to A V ≈ 3).However, the stellar continuum spectrum suggests a lower value (A V ≈ 1.5), indicating that the attenuation is highly structured.The galaxy is well resolved in high resolution (∼ 0. ′′ 2) ALMA continuum and CO observations, after correction for lensing, revealing a secu-larly evolving, massive galaxy ecosystem with, unlike most known HyLIRGs, no major merger signatures.The systematic rotation in both cold and ionized gas tracers is consistent with a massive baryonic disk (M baryon > ∼ 10 11.1 M ⊙ ) living in a massive dark matter halo (M DM ∼ 10 13.2 M ⊙ ).The massive disk of HyLIRG forms stars at an intrinsic rate of ∼ 1.5 × 103 M ⊙ yr −1 and has a Sérsic profile halflight radius of ∼ 2.9 kpc.The disk also shows indications of two grand spiral arms.The compact bluer peak seen within the critical curve in the inset forms the four images labeled 2a, 2b, 2c, and 2d in the middle panel of Fig. 4. By measuring the half-light radius of these images as ∼ 4 ′′ in the image plane and dividing by the linear magnification factor of µ 0.5 , we estimate the effective radius of the peak as ∼ 0.9 − 1.3 kpc in the source plane.Accounting for about 1/7 of the rest-frame g-band flux of the HyLIRG, the peak likely represents a giant star-forming region -a typical signature expected for high Toomre instability in massive high-z disks.In contrast, the stellar bulge of the galaxy is largely devoid of cold gas, suggesting inside-out quenching.PJ1053+60 (or G145.25+50.84):This HyLIRG (Fig. 5) has the highest redshift and has the largest gravitational amplification among our XMM-Newton targets (Table 1).It is one of the PAS-SAGES sources also studied with Herschel (Harrington et al. 2016;Berman et al. 2022), independently identified by Cañameras et al. (2015), and followed-up with HST observations (Frye et al. 2019).Though not particularly bright, intrinsically, the source shows highly excited CO emission and very broad line profiles (e.g., Harrington et al. 2021).The SMA ∼200-230 GHz continuum image of the source (see Appendix A) shows an almost complete elongated Einstein ring with an average diameter of 11 ′′ (see also Cañameras et al. 2015).The ring is also well detected in the Spitzer/IRAC 3.6 and 4.5 µm imaging data (Appendix B).In contrast, the HST 1.6-µm image shows little evidence for the ring, indicating strong dust attenuation.The HST 1.6-µm image, however, reveals additional lensed features from galaxies in the background other than the HyLIRG.These features are included in the lens modeling detailed in Appendix C. The result is illustrated in Fig. 5.

XMM-Newton observations
The relevant data sets of the XMM-Newton observations are obtained from a set of three X-ray CCD cameras comprising the European Photon Imaging Camera (EPIC) 3 .Two of the cameras are MOS (Metal Oxide Semiconductor) CCD arrays, while the other uses pn CCDs and is referred to as the pn camera.The angular resolution of the EPIC pn is about 6.6 ′′ (FWHM) or 15 ′′ (80% energyencircled radius), while the MOS detectors have slightly better resolutions (e.g. 6 ′′ FWHM).Table 2 summarizes the key parameters of the XMM-Newton observations.Limited by the orbital period, each observation of a target was split into multiple exposures.For this study, we use data from the European Photon Imaging Camera (EPIC) instruments, MOS1/2 plus pn cameras.

Imaging Data
We use the pipeline product (PPS) directly to visualize the data.To obtain maximum counting statistics of the imaging data, we combine the PPS MOS1/2+pn count intensity images produced for each exposure in the 0.3-0.7 keV, 0.7-1.2keV, and 1.2-7.0keV bands.According to the XMM-Newton specifications (XMM-SOC-GEN-ICDD-0024, Issue 4.6), these PPS images have been corrected for vignetting and exposure, while events from MOS1/2 cameras are given weights 5.41, 4.64, and 2.96 in the three bands for the pn ones, from which minimum instrument backgrounds and out-of-time events have also been removed.We further merge these combined images of individual exposures, excluding pixels labeled NaN, accounting for the slightly different pointing directions, and using the PPS (band 8000) exposure images as weights.The merged images in the three bands are finally added together to form a broadband (0.3-7 keV) image of the target (Fig. 6).
We also perform a preliminary spatial analysis of the XMM-Newton data mainly to check how point-like the X-ray emission is for each target (Table 3).Fig. 7 shows the close-up pn count images of the three targets in the 0.5-4.5 keV band.The signal to noise ratios (S/N) in these images tend to be optimal.The imaging data of PJ1336+49 look quite noisy, although there is obviously a net excess of X-ray emission.This is more apparent in the 1-D radial intensity profile of the emission shown in Fig. 8A.The other two targets are much brighter.We fit the profile of each target to obtain the normalization of the point spread function (PSF) and the local constant background.The PSF is approximately constructed at an effective photon energy of 3 keV and is insensitive to a small change in this energy.We check this procedure with SDSS J105315.15+605145.7 -a QSO at z = 1.1908, 58 ′′ west of PJ1053+60 (Fig. 6).This source is quite bright and is expected to be point-like, as confirmed by the comparison shown in Fig. 9 and Table 3, suggesting that our adopted PSF is good.Fig. 8 presents the comparison of the profiles with the best-fit normalizations and local backgrounds of the targets.The profile of PJ1053+60 is clearly inconsistent with being point-like at a high confidence level (Table 3).Based on the present analysis of the X-ray data, we cannot reject the point-like nature of the emission for both PJ0116-24 and PJ1336+49 with > ∼ 80% confidence.However, a close examination of the image of PJ0116-24, as well as a preliminary 2-D fit, suggests that this target may not be point-like in the 2-D emission distribution (e.g., Fig. 10A).While more careful spatial analysis is required to make a firm conclusion, we focus here on the spectral analysis and results.

Spectral Data
For spectral analysis, we independently process the XMM-Newton data with the Science Analysis System (SAS), primarily for our X-ray spectral analysis of each target.This process includes creating event files separately for the MOS and pn data using the emproc and epproc routines; filtering out time intervals with high background activity4 ; and applying the edetect chain to detect discrete sources.We also use evselect to extract the MOS1/2 and pn spectra and their corresponding local background spectra, and rmfgen and arfgen to obtain the corresponding instrument response and effective area files.Spectra from the same instrument but different exposures of the same target are combined with epicspeccombine.The two MOS spectra are then combined due to their similar properties.The centers of the source extraction regions are determined from the edetect chain source detection results.The X-ray centroid coordinates (R.A./Dec.J2000) of PJ0116-24, PJ1336+49, and PJ1053+60 are 1:16:46.63/-24:37:03.0,13:36:35.13/49:13:12.9,and 10:53:22.08/60:51:43.2,respectively.The on-source spectral extraction aperture of each source is 20 ′′ radius with a corresponding annulus background extraction region with an inner radius of 35 ′′ and an outer radius of 175 ′′ .The aperture size is large enough to obtain the majority of the spectral data, while regions are excluded for other detected X-ray sources in each background annulus.

Spectral Modeling
The modeling of the background subtracted spectral data must account for potential X-ray contributions of several types: (discrete and diffuse) X-ray sources in the foreground lenses, and the AGN/non-AGN components of the lensed HyLIRGs.In Appendix D we estimate the foreground lens contributions and find them to be negligible (count rates < ∼ 1% in the 0.5-2 keV band and less in broader bands).Therefore, we focus here on modeling contributions from the HyLIRGs.
Our spectral modeling uses a combination of the X-ray spectral fitting package Xspec (Arnaud 1996) and the Bayesian X-ray Analysis (BXA) software (Buchner et al. 2014) 5 .We further use the python package, chainconsumer (Hinton 2016) to improve the presentation of the corner plots.The BXA takes the initial bestfit parameters from Xspec to perform the MCMC analysis.We use the High Energy Astrophysics Science Archive Research Center (HEASARC) routine ftgrouppha to group spectral channels of the resulting combined MOS and pn spectra to obtain certain counting statistics.For the use of the C-statistic (called cstat in Xspec) or more precisely the W-statistic (wstat), appropriate for dealing with the Poisson data (https://heasarc.gsfc.nasa.gov/xanadu/Xspec/manual/node319.html),we group the spectral data to obtain a minimum of seven counts per bin, as recommended in Xspec (https://asd.gsfc.nasa.gov/Xspecwiki/low_count_spectra).For ease of the visualization of the spectra, we group them to obtain background-subtracted S/N > 1 per bin in Xspec.
The sophistication of our spectral modeling is limited by the counting statistics of the data.We start with simple models with the minimum number of fitting parameters and add complications only when justified, either statistically or physically.All spectral components are subject to the Galactic foreground absorption, modeled with tbabs with column density fixed at N H,G (Table 1).We fit an effective intrinsic (rest-frame) X-ray absorbing gas column density (N H,HMXB ) with the Xspec model ztbabs.The Xspec default solar metal abundances are assumed for all absorbing gases.We model the intrinsic HMXB contribution with a redshifted power-law zpowerlw; the corresponding Xspec model syntax is const*clumin(zpowerlw), where the convolution model clumin, allows us to conveniently pre-set the intrinsic rest-frame luminosity of the contribution to a reference luminosity (which we set to be X HMXB,r × the SFR of a target), while the model const or X HMXB serves as the normalization of the component to be fitted.Taken together, we have a basic set of spectral models: tbabs * (ztbabs * const*clumin(zpowerlw)).We find that this model is sufficient for the characterization of PJ1336+49, where the three free parameters to be constrained by the spectral fit are: X HMXB , the photon index of the power law (Γ HMXB or Γ nonAGN ), and the intrinsic X-ray absorbing column density (N H,HMXB or N H,PJ1336 used later in the two-source joint fit).However, PJ0116-24 and PJ1053+60 show X-ray evidence for the presence of AGN (see § 4).Therefore, we have included another redshifted power-law component, ztbabs * clumin(zpowerlw), with the corresponding X-ray absorbing column density (N H,AGN ), AGN power-law photon index (Γ AGN ), and luminosity (L X,AGN ) to be fitted.We adopt a prior for the AGN photon index as a normal distribution with a mean of 1.95 and a standard deviation of 0.15 (Buchner et al. 2014, see also Migliori et al. (2023); Zappacosta et al. ( 2023)).From the fitted spectral models we derive the observed absorbed 0.5-8 keV flux for each spectral component ( f X,HMXB or f X,AGN ).

RESULTS
Fig. 6 presents the merged MOS+pn images of our targets in the stacked broadband, as well as in individual energy bands.These images give an overview of the XMM-Newton data quality of individual targets, with respect to their local environment.They are all clearly detected in the broadband images.Fig. 7 shows the unsmoothed XMM pn data close-ups of the targets in the 0.5-4.5 keV band.It is clear that the XMM-Newton counting statistics of PJ1336+49 is too limited to give useful 2-D spatial information.Fig. 8 compares the radial intensity profiles of the targets with that of the PSF (Table 3).We further present the close-ups of the smoothed broadband images of PJ0116-24 and PJ1053+60 in Fig. 10.The Xray emission of PJ1053+60 is clearly not point-like, even in the radial profile analysis (Table 3).In the smooth 2-D image, PJ0116-24 seems to be elongated, as indicated by the higher intensity contours, although this is yet to be quantified via detailed spatial analysis.The results from our spectral analysis are summarized in Table 4.In particular, our measured X HMXB values are compared with the reference calibration data in Fig. 1.In the following, we focus on describing the X-ray spatial and spectral properties of the targets.
We start with PJ1336+49 -probably the simplest case here.The detected net number of XMM-Newton counts of the target is very limited (Fig. 11A; Table 4), which is expected if the X-ray emission from the HyLIRG is dominated by HMXBs.Their spectrum can be well characterized by a power law with an intrinsic absorption column density of ∼ 3×10 21 cm −2 (Figs.11A and 12; Table 4).The best-fit model gives X HMXB ∼ 4 or an apparent X-ray luminos- ity L X ≈ 3 × 10 44 erg s −1 .The relatively steep spectrum (e.g., power law index Γ HMXB ∼ 2.8) is consistent with the hypothesis that the Xray emission represents the collective contribution from the HMXBs of the galaxy.Although both the power law index and the absorption are not tightly constrained, we see no evidence for a significant AGN component, which should typically have a harder spectrum, consistent with the lack of a point-like source component in radio (Fig. 3A).
The spectra of PJ0116-24 and PJ1053+60 (Fig. 11B and C) do not fit as well with the simple power-law model with wstat/   .. 3.9(0.1,7.1) 35.9(26.6,40.9) ../1.00(0.05,6.06) The goodness value of a fit represents the percentage of 1000 data simulations that have a smaller wstat/d.o.f., while the 90% confidence interval of each parameter is included in the parenthesis next to the best-fit value.The unit of X HMXB is X HMXB,r in Eq. 2. The AGN luminosity is calculated in the rest frame 0.5-8.0keV band of each source, while the absorbed fluxes of the HMXB and AGN components ( f X,HMXB and f X,AGN ) are inferred (hence the marked '-' ) from the best-fit model in the observed 0.5-8.0keV band.
d.o.f=189.65/152(where d.o.f is the degree of freedom) and 817.94/699, respectively.We find that the double power-law model gives better fits to the spectra of these two sources (Fig. 11B-C, 13, and 14; Table 4).Their estimated luminosities are far too large to be consistent with the values expected from the X HMXB estimated from the above fit to PJ1336+49 or from the reference calibration X HMXB,r (Eq.2), indicating a substantial AGN contribution.Indeed, a follow-up analysis of the data available in other wavelength bands supports this hypothesis (see § 5.3).Here, the fits with the double power-law model allow us to decompose these two contributions in PJ0116-24 and PJ1053+60.Studies have shown that ULXs on average tend to have spectra that steepen substantially at rest frame energies of a few keV (e.g., Bachetti et al. 2013;Brightman et al. 2018;Walton et al. 2020), whereas an AGN spectrum intrinsically has a flatter power-law-like spectrum and tends to be flattened due to the Compton hump at rest frame energies > ∼ 10 keV (e.g., Reynolds et al. 1995;Brightman et al. 2018).This makes a crude spectral decomposition of the two components possible.
We further jointly fit the non-AGN spectral components of the PJ0116-24 and PJ1336+49 to better constrain X HMXB and isolate the AGN contribution in PJ0116-24.The reason for doing this twosource joint fit (instead of a joint fit to the three sources) is that the AGN contribution in PJ0116-24 is relatively small, while the PJ1053+60 spectrum is more AGN-dominated.Furthermore, the lopsided highest X-ray intensity contour seen in Fig. 10A is consistent with the AGN being at the nucleus of PJ0116-24, forming the two lensed peaks 1a and 1b in Fig. 4A (Kamieneski et al. 2023).The results of the joint fit are included in Table 4, as well as shown in Figs. 15 and 16.The fitted X HMXB is consistent with that from the fit to the PJ1336+49 spectrum alone but the uncertainty interval is significantly tightened.The obtained X HMXB is 3.4 [with the 90% confidence uncertain interval of (2.2,7.1)].
The presence of an AGN in the PJ1053+60 field is now evident.Its X-ray emission in projection extends much further away from the lensed DSFG towards the southwest, where an apparent mid-IR counterpart of this AGN is found (Appendix B).We decompose the XMM-Newton spectra of PJ1053+60 using two power law models for the AGN and HMXB contributions and assuming that the redshift of the AGN is the same as that of the DSFG (Fig. 11C and 14; Table 4).The fitted X HMXB of 14.6 (8.7,23.8) is surprisingly high, although the reliability of this measurement cannot be fully judged with the limited quality of the XMM-Newton spectral data.4).
The AGN redshift assumption should have little direct effect on the power-law index of the AGN component, but would change the bestfit values of its intrinsic absorption and luminosity if its actual redshift is different, which in turn could affect the constraint on X HMXB .Future anticipated Gemini-N observations will help to measure the redshift of the AGN directly, and further spatial/spectral analysis of PJ1053+60 will be reserved for Diaz et al. (2023).What we can conclude here is that the spectral data of HyLIRGs are all consistent with an enhanced X HMXB .

DISCUSSIONS
The goal here is to better understand our results and their implications in a broader context.We first compare our X HMXB estimate with existing ones based primarily on local star-forming galaxies and   4).
then examine the origin of the enhanced X HMXB in the HyLIRGs.We also discuss the implications of the AGN detections in PJ0116-24 and PJ1053+60 for the AGN and extreme SF co-evolution.

Uncertainties in the X HMXB Estimate
Our best constraint on the non-AGN component comes from the joint fit of the PJ0116-24 and PJ1336+49 XMM-Newton spectra (Figs. 15 and 16; Table 4) 4).The X-ray spectral model is a power law for PJ1336+49, representing the collective contribution of the HMXBs in the galaxy, while for PJ0116-24 a separate power law is included to account for the contribution from the AGN contribution.These two components are also separately plotted.component, Γ HMXB = 2.46 (2.09, 3.32), is consistent with what is expected from the collective X-ray spectrum of HMXBs, which can be characterized by a power law with Γ HMXB = 2.1±0.1 in the 0.25-8 keV range, after accounting for their intrinsic diversity of emission and absorption properties (Sazonov & Khabibullin 2017).The overall contributions from hard, soft, and supersoft sources are comparable in the energy range.While the soft part (below 2 keV) of the spectrum is mostly due to soft and supersoft sources, classical harder ULXs dominate at energies above a few keV, which is most relevant here since we are observing the rest-frame > ∼ 2 keV emission of HyLIRGs at z > 2. The spectrum is expected to steepen with increasing energy, as typically seen for individual ULXs (e.g., Bachetti et al. 2013;Brightman et al. 2018;Walton et al. 2020).
Our best estimate of the X HMXB factor (Table 4) is a factor of ∼ 3.4 larger than expected from X HMXB,r (Fig. 1).A possible systematic uncertainty is in the X-ray absorption correction.The X HMXB estimate is correlated with the absorption in the X-ray spectral fits.Our best-fit intrinsic N H,HMXB (Table 4) is quite moderate; the absorption correction for X HMXB or µL X is only 3.1% and  4).0.8% for PJ0116-24 and PJ1336+49.Our X-ray spectral modeling approximates the absorption as being the foreground of the X-ray emission.Any non-uniformity of the absorption towards the starforming regions in a galaxy should lead to an underestimate of the actual N H,HMXB and hence X HMXB in our spectral modeling.Accounting for this underestimate, which is typically hard to quantify, would make X HMXB even higher.For PJ0116-24, we can compare our X-ray absorption estimate with the optical reddening measurements via the VLT/ERIS spectroscopy: E(B − V ) * = 0.49 ± 0.10 and E(B − V ) neb = 0.99 ± 0.21 towards stars and HII nebulae, respectively (Liu et al. 2023).These measurements, together with the conversion of N H = 4.8 × 10 21 E(B − V ) cm −2 mag −1 for the ISM in our Galaxy (Bohlin et al. 1978), give the corresponding N H = (2.3 ± 0.5) × 10 21 cm −2 and = (4.5 ± 1.0) × 10 21 cm −2 , which are slightly larger than N H,HMXB ∼ 1.6 × 10 21 cm −2 obtained from our X-ray spectral modeling (Table 4).Such uncertainty at the N H,HMXB ∼ 10 21 cm −2 level has a negligible effect on our measurement of X HMXB (especially in the rest-frame energy range above 1.5 keV), compared to other systematical and statistical uncertainties in our spectral analysis.Therefore, we conclude that the large X HMXB from our modeling of PJ0116-24 and PJ1336+49 cannot be due to an overestimation of the X-ray absorption.
In fact, we could have underestimated the absorption, which may be caused by the presence of soft X-ray emission in addition to the power-law components used in our spectral modeling.Such a soft excess can arise in an AGN.The origin of the soft excess is one of the major open questions in AGN research (e.g., Turner & Miller 2009;Nandi et al. 2023).Common in Type 1 AGNs (e.g., Piconcelli et al. 2005), the soft excess can be comparable to the luminosity of the AGN power-law component in the rest-frame 0.5-2 keV band and can typically be fitted with a thermal model with a characteristic temperature ∼ 0.1 − 0.3 keV.Its expected peak emission is thus clearly outside the rest-frame energy range used in our spectral analysis.However, the presence of the high-energy tail of the excess could potentially lead to an underestimation of the X-ray absorption.
A similar effect can be expected from diffuse soft X-ray emission around an AGN and in the general ISM.Soft X-ray emission from the general ISM can be expected in star-forming galaxies, primarily because of stellar feedback from massive stars via fast stellar winds and core-collapsing SNe.These diffuse soft X-ray emission components are again typically important at photon energies < 1.5 keV, which are not included in our spectral analysis.Similarly, the interaction between the energetic feedback of an AGN (via jets, winds, and radiation) and its surrounding medium can produce an extended emission that tends to be soft.The physical scale of the emission is typically significant on sub-kpc scales (e.g., Travascio et al. 2021;Fabbiano & Elvis 2022).
Extended hard X-ray emission can also be produced by AGNs, mostly due to scattering of AGN radiation (including both reflection and fluorescence of photons) by the ISM, although heating by fast shocks can occasionally be considerable.Such emission typically accounts for ∼ 1% of the intrinsic X-ray luminosities of local AGNs (e.g., Ricci et al. 2017).For such an emission to explain the enhanced X HMXB , the AGNs in PJ1336+49, PJ0116-24, and PJ1053+60 need to have apparent luminosities ∼ 3, 4, and 10 ×10 46 erg s −1 , much greater than those inferred from our spectral analysis (Table 4).However, the emission can be important for AGNs that are strongly obscured.Since the dense clouds responsible for the scattering are typically more concentrated at smaller galactocentric radii, the intensity of this emission is expected to decrease rapidly with off-AGN distance.Indeed, such emission has been detected up to several ×10 2 pc radii around nearby Compton-thick AGNs and with total 3-7 keV luminosities up to several ×10 40 erg s −1 (Fabbiano & Elvis 2022;Trindade Falcao et al. 2023).The angular to physical scale conversion is 7. 659, 8.46, or 7.433 kpc per arcsecond at redshifts 3.25, 2.12, or 3.55 of PJ1336+49, PJ0116-24, or PJ1053+60.Therefore, even taking into account the lensing amplification (e.g., 6.5 and 6.0 for the AGN in PJ0116-24), a kpc-scale feature in the galactic core region would likely be observed on the sub-arcsecond scale in the image plane, smaller than the resolved X-ray emission extent.
It may be possible that the bulk of the enhanced X-ray emission arises from the hot intergalactic medium (IGM) surrounding individual HyLIRGs.The situation may be similar to the Spiderweb protocluster (J1140-2629) at z = 2.156), for which a Chandra ACIS-S observation (700 ks) reveals the presence of significant extended Xray emission out to a radius of ∼ 12 ′′ (100 kpc), as well as an AGN and radio jets (e.g., Tozzi et al. 2022;Di Mascolo et al. 2023).The luminosity of the thermal emission is comparable to those of the non-AGN component of the HyLIRGs reported here.However, in this scenario one should not expect any significant X-ray absorption towards the diffuse emitting regions, which seems to be marginally inconsistent with the intrinsic N H measurements (or N H,HMXB values in Table 4).Furthermore, protoclusters like the Spiderweb are predicted to be very rare at z > 2 (e.g.Saro et al. 2009).It is therefore unlikely that our target HyLIRGs, especially those PJ1336+49 and PJ1053+60 with z > 3, are all inside such protoclusters.Nevertheless, the protocluster scenario of enhanced X-ray emission in the HyLIRGs needs to be explored in future work.
Finally, the large X HMXB in HyLIRGs is in great contrast to the Xray deficiency observed in nearby (U)LIRGs (e.g., Persic & Rephaeli 2007;Iwasawa et al. 2009;Lehmer et al. 2010;Torres-Albà et al. 2018).This is presumably due to the fact that the SF in HyLIRGs is typically more distributed than in (U)LIRGs, which minimizes both the severe absorption and youth/rapid cluster destruction effects of compact starbursts on the X-ray emission ( § 1.1.1).In particular, both soft X-ray contamination and absorption effects should be much less important in our X-ray luminosity or X HMXB estimate made in the rest-frame energy range above 1.5 keV.As discussed above, the moderate X-ray absorption toward star-forming regions in the HyLIRGs seems to be reasonably well accounted for in our spectral modeling.Taking all this into account, we conclude that X HMXB is probably indeed ∼ 3 or could be even higher in the PJ1053+60 case.

Origin of the X HMXB Enhancement
It has been shown that a significant non-AGN L X or X HMXB enhancement can be caused by the very low metallicity ( < ∼ 10% solar, due mostly to low stellar mass loss) and/or the youth of a stellar population (e.g., Fragos et al. 2013;Gilbertson et al. 2022).However, this scenario does not apply to HyLIRGs.They tend to have extreme but distributed SF (Swinbank et al. 2015;Rujopakarn et al. 2016;Iono et al. 2016;Kamieneski et al. 2023), instead of a very localized starburst as typically seen in nearby (U)LIRGs (e.g., Díaz-Santos et al. 2010).Such distributed SF should typically have lasted longer than the dynamic time scale of a galaxy (∼ 10 8 yr).Although the SF could have started in a low metallicity environment, the chemical enrichment is expected to have a short time scale in HyLIRGs; so their metallicity is expected to be close to solar, or at least not particularly low, consistent with their dusty nature.Indeed, the metallicity of PJ0116-24 is shown spectroscopically to be slightly supersolar (Liu et al. 2023).In comparison, HMXBs, which should dominate the non-AGN component in such galaxies and trace the very recent star formation history over a period probably mostly shorter than the dynamic time scale of our target HyLIRGs.Therefore, we don't expect the large X HMXB to be due to the very low metallicity or youth of their stellar population.
An X HMXB enhancement could also be caused by an underestimation of µSFR.For each of our target HyLIRGs, µSFR (Table 1) is converted from the IR (rest-frame 8-1000 µm) luminosity derived from the best fit of the spectral energy distribution (SED) (Berman et al. 2022;Liu et al. 2023).The accuracy of the SED depends mostly on the photometry, the uncertainty of which is largely systematic and difficult to quantify (Berman et al. 2022), but could be up to 30-40% or even more.For PJ0116-24, Liu et al. (2023) have recently obtained an alternative estimate of µSFR = (2.11± 0.45) × 10 4 M ⊙ yr −1 from the H α emission corrected for extinction by the Balmer decrement H α /H β ratio.This estimate is only slightly less than µSFR= (2.51 ± 0.65) × 10 4 M ⊙ yr −1 from the SED fit (Table 1) and within their quoted 1σ error bars.An under-estimate of the µSFR is expected from the H α emission, since the decrement correction most likely underestimates the extinction because of its mixing with the emission on the galaxy scale.Therefore, the comparison gives us some confidence in concluding that our SED-inferred µSFRs of the HyLIRGs are not systematically underestimated to explain the large X HMXB (see also Mineo et al. 2012).
One may wonder if the differential lensing of the HyLIRGs might play a role in the measurement of X HMXB .We have so far assumed that L X and L IR (or other SFR tracers) of a galaxy are magnified equally, with their ratio essentially unchanged in the image plane of the lensed galaxies § 1.2).However, the assumption may not hold if the dynamical formation of HMXBs is important in HyLIRGs, which may be expected in dense clusters and probably in galactic nuclear bulges (e.g., Garofali et al. 2012;Rangelov et al. 2012;Kremer et al. 2020;Rizzuto et al. 2022).If this is the case, then the X HMXB depends on local clustering properties and is thus intrinsically nonuniform across a galaxy.This non-uniformity can, in principle, be enhanced by differential lensing magnification.X-ray emission from a dense stellar cluster, if located at the lensing caustic and much more compact than its surrounding far-IR emission region, could be disproportionally enhanced.However, because the source plane area covered by very high magnification (µ > 20) is very small, we don't expect that such a differential magnification would systematically enhance X HMXB .
We do suspect that the dynamically formed population of HMXBs itself is responsible for the enhancement of the non-AGN X-ray emission in HyLIRGs.This scenario is similar to that used to explain the population of LMXBs observed in nearby globular-clusterrich elliptical galaxies (e.g., Gilfanov et al. 2022).Dynamical interactions can cause mass segregation, tidal capture and disruption, and multibody exchange and binary mergers (e.g., via Lidov-Kozai cycles or by stellar collisions; Kozai (1962); Eggleton & Kiseleva-Eggleton (2001); Baumgardt & Klessen (2011)) in dense stellar clusters (e.g., Fabian et al. 1975;Hills 1976;Stone et al. 2017) and possibly in galactic nuclear bulges (e.g., Voss & Gilfanov 2007;Zhang et al. 2011).More massive stars, more compact stellar remnants (neutron stars and BHs), and more X-ray binaries are thus expected to form (e.g., Stone et al. 2017;Rizzuto et al. 2023).In extremely dense stellar environments (e.g., when the central number density of stars > ∼ 10 7 pc −3 ), initial stellar mass BHs may even grow into IMBHs via tidal capture and disruption events (e.g., Rizzuto et al. 2023;Stone et al. 2017;Di Matteo et al. 2023), which could explain BHs with masses of a few 10s (even > ∼ 100 M ⊙ ), detected by the LIGO-Virgo-KAGRA interferometers (e.g., Mahapatra et al. 2021).Such formed IMBHs have also been considered to be potential seeds of SMBHs.Furthermore, the dynamical interactions can also kick out a good fraction of the binaries from the clusters, eventually forming X-ray binaries that are observed outside individual clusters.In short, the dynamical effect in young and intermediateage dense clusters (e.g., Kremer et al. 2020;Rizzuto et al. 2022;Hunt et al. 2023) is expected to produce disproportionately more compact stars (neutral stars and BHs) and more massive ones and more in accreting binaries, which can be responsible for generating more non-AGN X-ray emission and hence the large X HMXB .Indeed, an indication for the increasing X HMXB with the surface SFR, traced by the observed sub-mm emission, is seen in PJ1053+60.The XMM-Newton image of this HyLIRG is sufficiently resolved to allow for a spatial decomposition of the non-AGN and AGN components (Fig. 10B).A preliminary analysis of the image shows an excess of the X-ray emission associated with the northeast CO intensity peak of the lensed ring of the DSFG, which can naturally be explained by a surface-SFR-dependent enhancement of the X HMXB .The quantification of this dependence will be attempted in an upcoming paper by Diaz et al. (2023).
We expect that HyLIRGs are the sites where the dynamical effects manifest greatly in generating a large population of luminous HMXBs such as ULXs.It has long been proposed that gas-rich galaxies such as HyLIRGs tend to have higher molecular gas densities, resulting in higher SFRs and clustering efficiencies, and can form cluster populations that extend to higher initial cluster masses (e.g., Adamo et al. 2020).In particular, galaxy mergers are considered to be the most efficient producer of dense clusters.A recent high-resolution simulation of galaxy mergers (Li et al. 2022) shows that the final coalescence of two large gas-rich galaxies is very turbulent and clumpy, which can lead to the formation of a large number of dense clusters whose masses dominate the total SFR.This can be seen in a PASSAGES DSFG at z = 2.24 in the PLCK G165.7+67.0galaxy cluster, for which spatially-resolved JWST NIRSpec spectroscopy shows the galaxy to separate into several interacting high mass clumps (Frye et al. 2023a).There is also a strong positive correlation between the efficiency of cluster formation and the surface SFR, as well as the turbulent gas pressure.Thus it is natural to predict that the BH merger rate in stellar clusters is highest at cosmic noon (Di Matteo et al. 2023), where the population of HyLIRGs also peaks.Furthermore, starbursts in HyLIRGs seem to be similar to other clumpy star-forming galaxies observed at high-z (e.g., Elmegreen et al. 2009;Swinbank et al. 2015;Rujopakarn et al. 2016;Iono et al. 2016;Kamieneski et al. 2023) and substantially more spatially distributed than in local (U)LIRGs (Kamieneski et al. 2023), which minimizes the effects of the extreme extinction/X-ray absorption and stellar cluster destruction (see § 1.1.1).
The stellar dynamical effect has also been proposed to explain other properties of high-z galaxies.Perhaps most notable are the anomalously high [N/O] ratios observed by JWST, similar to those known in globular clusters of the Milky Way (Belokurov & Kravtsov 2023).Stellar collisions or other dynamical processes in such clusters can naturally increase the formation of close binaries, in addition to the formation of very massive stars, which may be responsible for high UV continuum emission observed (e.g., Cameron et al. 2023).Thus, newly synthesized nitrogen can be exposed and ejected into the ISM (e.g. via Wolf-Rayet phases; Harrington et al. 2019).Massive close binaries could further evolve into HMXBs.We conclude that stellar dynamics may play a very important role in determining the amount and spectrum of the radiation from stars in extreme star-forming galaxies, which are common at high-z.

AGN and SF co-evolution in HyLIRGs
Our detection of AGN in PJ0116-24 and in the immediate vicinity of PJ1053+60 was not expected based on our selection criteria stated in § 2. The PASSAGES galaxies were selected to be apparently starburst-dominated.This selection is based on the comparison in the mid-and far-IR color-color diagrams for starburst/AGN templates (Harrington et al. 2016(Harrington et al. , 2018) ) and the degree of radio and/or mid-IR excess in the SED fits (e.g., Yun et al. 2001;Berman et al. 2022).Two of our three XMM-Newton targets, PJ0116-24 and PJ1336+49, are included in the study by Kamieneski et al. (2023).They show that both have radio/FIR correlation parameters q FIR consistent with the median for their selected PASSAGES DSFGs (q FIR ∼ 2.4) and apparently above the q FIR < 1.8 threshold for an AGN-powered object (see their Table 4 and Fig. 9).Our X-ray detection of AGN in PJ0116-24 and PJ1053+60 thus indicates that AGN may be quite common in HyLIRGs, although this is not evident in wavelength bands other than X-ray.Now that we know of the presence of the AGNs, let us examine their multi-wavelength properties, which helps to determine their nature and potential role in the AGN and SF co-evolution.
The association of an AGN with the DSFG in PJ0116-24 is com-pelling, especially now with spatially resolved multi-wavelength spectroscopic data.Although the dust continuum and CO emissions do not peak at the center of the DSFG, the radio continuum emission does (Fig. 4A), and most likely represents this AGN.This scenario is consistent with a high [NII]λ 6584/Hα ratio of 1.35 ± 0.20 at the center (compared to the galaxy-wide [NII]λ 6584/Hα ratio of 0.36 ± 0.04) measured with the near-IR ERIS IFU data, as well as with the X-ray spatial and spectral properties of the galaxy (e.g., § 4).The apparent lack of dust emission in the central region of the DSFG as seen in its source plane sub-mm image suggests that much of the gas may have been ejected from the central region, probably by the AGN (e.g., Stacey et al. 2022).The northeast and southwest images of the AGN have lensing magnifications of 6.5 and 6.0 in the HST image (Liu et al. 2023).Therefore, the total lensing magnification factor is about ∼ 13.With this factor, we infer that the intrinsic luminosity of the AGN is ∼ 10 43 erg s −1 .However, the uncertainty in the luminosity is large, mainly because of the uncertainty in the foreground absorption, which could be enhanced by the presence of the AGN torus (e.g., Yamada et al. 2023).
The presence of an AGN in the vicinity of PJ1053+60 was first indicated by its X-ray morphology and spectrum ( § 4).A followup multi-wavelength examination described in Appendix B led to the identification of this AGN based on its near-to mid-IR colors.
It is located at R.A./Dec.(J2000) = 10 : 53 : 21.75/60 : 51 : 41.4, or about 6.6 ′′ offset from the lensed DSFG (Fig. A1).The AGN appears point-like in the HST 1.6-mm image.To firmly establish the relationship between the AGN and the DSFG, we will first need to measure the redshift of the AGN.
If the AGN is indeed at the same redshift as PJ1053+60, then we will have to seriously consider several intriguing questions.One of them is whether the AGN was ejected from the DSFG (e.g., van Dokkum et al. 2023;Smole & Micic 2023).Such an ejection can occur during the merger of a binary SMBH.Gravitational radiation is emitted anisotropically due to asymmetries in the merger configuration.This anisotropic radiation leads to a gravitational wave kick, or recoil velocity, as large as a few 10 3 km s −1 (e.g., Campanelli et al. 2007;Smole & Micic 2023;Di Matteo et al. 2023).At such a speed, the separation between the AGN and the DSFG can probably be achieved well within 10 8 yr.This ejection scenario is consistent with the apparent lack of evidence for a counterpart to the AGN in the existing sub-mm continuum and CO data, indicating the dearth of the ISM.Alternatively, the AGN and DSFG may represent a pair of two rare galaxies, separated by < ∼ 50 kpc in projection, in the same forming galaxy cluster.Interestingly, the DSFG has a similar color, most likely due to the strong reddening of stellar light.
While a detailed analysis of PJ1053+60 will be presented in a later paper (Diaz et al. 2023), we note that such associations of strongly lensed extreme DSFGs and AGNs provide us with excellent laboratories for understanding similar systems that are mostly unlensed and thus difficult to study in a spatially resolved fashion.One example is the unlensed HATLAS J084933.4+021443(Ivison et al. 2019, see also the discussion in § 1.1.2).Existing studies also show that extreme SFRs ( > ∼ 300 M ⊙ yr −1 ) tend to be seen in the hosts of AGNs with intermediate X-ray luminosities (10 42.5 erg s −1 < ∼ L(8 − 28 keV) < ∼ 10 44 erg s −1 ), and that many dusty and fainter AGNs remain to be detected in distant galaxies (z > ∼ 2; e.g., Barger et al. (2019)).Furthermore, Kirkpatrick et al. (2017) point out that there is a rare population of weakly obscured AGNs that appear to coexist with galaxies having large amounts of cold gas or extreme SFR, which may represent a short-lived phase in the co-evolution of SMBHs and galaxies (see also Sokol et al. (2023)).The strongly lensed HyLIRGs allow for the spatially resolved study of such extreme DSFGs associated with AGNs, especially with high angular resolution X-ray observations which could be provided by Chandra and potentially AXIS (Advanced X-ray Imaging Satellitea NASA concept probe mission) (Mushotzky & AXIS Team 2019).

SUMMARY, CONCLUSIONS AND FUTURE PROSPECTS
We have presented the XMM-Newton observations of three HyLIRGs (Table 1) in the redshift range z = 2.12 − 3.55, selected from the PASSAGES sample of strongly lensed DSFGs, the ratio to sub-mm SEDs of which show no sign for AGNs.The primary goal of our XMM-Newton observations is to measure the collective X-ray emission of HMXBs and to constrain the X HMXB factor in our target HyLIRGs.The combination of the lensing magnification and the deep XMM-Newton observations, together with the extreme intrinsic infrared brightness of the HyLIRGs, has allowed us to detect the X-ray emission from each of them.This detection is complemented by the analysis and modeling of multi-wavelength data, some of which are presented for the first time.We have obtained the following main results and conclusions: • The 0.5-8 keV luminosity of PJ1336+49 is significantly greater than that expected from the reference calibration X HMXB,r of 3.9 × 10 39 erg s −1 /(M ⊙ yr −1 ).Its X-ray spectrum is consistent with that of ULXs, which are expected to dominate the HMXB emission.This, together with the absence of multi-wavelength evidence for an AGN in the lensed DSFG, suggests that the observed X-ray emission is primarily non-AGN in origin.
• The X-ray emission from PJ0116-24 is partly contributed by an AGN.The X-ray morphology is consistent with this AGN being at the center of the DSFG, as concluded from a multiwavelength analysis of the HyLIRG.The DSFG consists of a central stellar spheroid with little dusty gas and a disk of extreme SF with a compact blue clump.A joint spectral analysis of the non-AGN X-ray emission from PJ0116-24 and PJ1336+49 allows us to better constrain its properties.We find that X HMXB is about 3.4 or 2.2-7.1 (90% confidence uncertain interval) in units of X HMXB,r .
• The X-ray emission of PJ1053+60 morphologically shows the presence of an AGN.It has a near-to mid-IR counterpart, ∼ 6.6 ′′ (or ∼ 50 kpc in projection) from the lensed image centroid of the DSFG component, and is largely naked, lacking dust continuum and CO emission in existing observations.A spectral decomposition of the X-ray emission from PJ1053+60 suggests that the AGN has a restframe 0.5-8 keV luminosity of ∼ 10 46 erg s −1 and is obscured with log[N H (cm −2 )] ∼ 23.7 if it is at the same redshift as the DSFG.In this case, the AGN could be kicked out of the DSFG in the recent past (e.g.< ∼ 10 8 years ago), probably during a major galaxy merger.
• Our results suggest that AGNs may be quite common in HyLIRGs, although a study of a large sample of them is highly desirable.The strong lensing of such galaxies, as provided by the PASSAGES sample, provides a powerful tool to examine the AGN co-evolution with extreme SF.
• Our measured X HMXB factor provides a tight constraint on the HMXB population in the HyLIRGs.The large X HMXB (∼ 3; Fig. 1) can naturally be explained by an enhanced population of HMXBs (due to increased numbers and BH masses) that could dynamically form in dense clusters, some of which may eventually evolve into today's globular clusters.HMXBs formed in this way also produce BH binaries, the eventual merger of which can be detected by gravitational waves.The inclusion of the dynamical formation could significantly change the HMXB population synthesis modeling (e.g., Fragos et al. 2013), which has so far neglected such an effect.
While this XMM-Newton study demonstrates the scientific potential of X-ray observations of strongly lensed HyLIRGs, more work is clearly needed.For example, to further test our postulated scenario for the enhanced X HMXB , we need high angular resolution X-ray data, which can be provided by deep Chandra observations.Such data will allow a clear separation of the AGN and non-AGN components.We can also study the spatially resolved X HMXB distribution.If the dynamical formation of HMXBs is indeed important, X HMXB should statistically increase with the surface SFR or correlate with the presence of dense clusters, which can be probed by JWST observations.Such studies will allow us to learn about important highenergy processes, in particular the co-evolution of stellar to massive BHs with extreme SF in HyLIRGs.

ACKNOWLEDGMENT
We thank the referee, Giuseppina Fabbiano, for prompt and constructive comments, which helped to improve the paper.This work is largely supported by NASA under the grant 80NSSC22K1923, based on data obtained from XMM-Newton.This research is also based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.These observations are associated with programs GO-14653 and GO-14223.Support for program GO-14653 was provided by NASA through a grant from the Space Telescope Science Institute.Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute.This paper further makes use of the following ALMA data: ADS/JAO.ALMA # 2017.1.01214.S, as well as data from VLA and SMA.ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile.The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.E.F.-J.A. acknowledges support from UNAM-PAPIIT project IA102023, and from CONAHCyT Ciencia de Frontera (project ID: CF-2023-I-506).

APPENDIX A: SUBMILLIMETER ARRAY OBSERVATIONS (SMA) OF PJ1336+49 AND PJ1053+60:
We present data from two separate SMA observations for PJ1336+49 and PJ1053+60 to complete the mm-wavelength perspective of these objects in the context of the XMM-Newton observations.PJ1336+49 was observed in the extended configuration 2022 June 12 (6 antennas; PWV<2.5 mm) during the 2021B semester (Project ID: 2021B-S011; P.I.Harrington, K.) using the SWARM correlator and dual receiver tuning at ∼243 GHz and 343 GHz.The data were calibrated using the latest SMA pipeline6 implementation in CASA and will be further detailed in a forthcoming manuscript (Harrington et al. 2023) analyzing the CO(9-8) and [NII]205 µm emission lines.PJ1053+60 was observed in the extended configuration 2016-Nov 22 (8 antennas; PWV<3.5 mm) during the 2016B semester (Project ID: 2016B-S062; P.I.Yun, M.) using the SWARM correlator during its initial upgrade to wider bandwidth capabilities, here tuned to ∼200 GHz and 220 GHz.The data were calibrated using the MIR package and further converted to a CASA measurement set before imaging.More information will be presented in a detailed multi-wavelength analysis for PJ1053+60 including this SMA data set (covering CO(8-7) and CO(9-8) emission lines) (Diaz et al. 2023).APPENDIX B: Spitzer / IRAC DATA ANALYSIS OF PJ1053+60 In the preparation of the observing program for PJ1053+60, we find that Spitzer/IRAC data of this target 7 are useful and are not published.We here present a preliminary analysis of the data in the multi-wavelength context of this HyLIRG.Our primary interest is in the mid-IR properties of the AGN candidate source southwest of the DSFG, as indicated in the XMM-Newton and HST images (Figs.10B and A1B).We clearly see the counterparts of the source in the IRAC 4.5 µm and 3.6 µm images (Fig. A1D-F).We estimate its specific energy flux to be 14.8 ± 0.2 and 23 ± 2 µJy in the IRAC 3.6 and 4.5 µm bands.The color [3.6]-[4.5](AB)= 0.48 ± 0.09, corresponding to a power-law SED f ν ∝ ν −2 , is consistent with the AGN nature (e.g., Fig. A2; Stern et al. 2005;Assef et al. 2018).
The IRAC data also show an apparent counterpart to the DSFG with a red color (similar to that of the AGN), indicating strong dust reddening of the underlying stellar light.The DSFG and AGN may represent two members of a proto-cluster of galaxies, which may further include a linear feature off from both DSFG and AGN, as 7 https://sha.ipac.caltech.edu/applications/Spitzer/SHA/seen in the HST 1.1-µm image (Fig. A1B).Because no radio or dust continuum counterpart has been found for this feature, it likely represents an evolved galaxy strongly lensed by a foreground galaxy group complex (Appendix C).

APPENDIX C: LENS MODELING FOR PJ1053+60
The main purpose of this modeling is to provide a reasonable estimate of the lensing mass of PJ1053+60.Fig. C1 shows a complex foreground lensing system of galaxies together with the identified multiple-imaged systems.The green contour displays the SMA 850 µm emission at z = 3.549 which is lensed into a partial Einstein ring.The centroids of three emission peaks, as marked in Fig. C1A, are used for the lens modeling.We also find a long arc consisting of two sets of image systems: 2ab and 3ab (Fig. C1A), which appear to be caustic crossings and thus help constrain the lens model (Fig. C1B).However, these systems have no measurements of redshifts, which are left as free parameters in our lens modeling.These sets of lensed image systems are used to constrain the lens model.
Also marked in Fig. C1 are the foreground lensing galaxies that are included in our lens-mass modeling.They are identified by searching for the red cluster sequence in a color-magnitude diagram constructed with HST f140W and F160W (Gladders & Yee 2000).Unfortunately, the HST images are heavily contaminated by a foreground star, potentially affecting photometric and color measurements of galaxies in the field.We therefore make the foreground membership particularly strict, setting a tight color cut and including only those galaxies with the most influential lensing potential.
We characterize the lensing mass model of the PJ1053+60 foreground galaxy group, following the semi-parametric Light Traces Mass (LTM) approach.The details regarding the general implementation of the LTM model, as well as the initial lens modeling of PJ1053+60, are presented in Frye et al. (2019).The modeling assumes that the gravitational halo mass of the lens is traced by the observed foreground galaxies (Zitrin et al. 2009(Zitrin et al. , 2015)).The surface mass density profile of each galaxy is characterized by a 2-D power law with a given index or a free parameter left for the model to optimize.The 2-D mass distributions of the lensing galaxies are combined and smoothed with a Gaussian kernel of a given size, which is another fitting parameter.We also implement an external shear to allow some flexibility in modeling the overall ellipticity of the lens mass distribution, which introduces two additional fitting parameters, a shear amplitude and a position angle.More freedom is given to the brightest group galaxy by fitting its core radius, ellipticity, and position angle.The redshifts of the identified galaxies are also left as free parameters in absence of reliable measurements.
The lens model is constrained via MCMC minimization by the positions of the identified multi-lensed images, which provide strong lensing information such as parity.We find that the best-fit model reproduces the image positions of the lensed sources well with an RMS∼0.3 ′′ with a χ 2 /d.o.f = 24.02/12.The redshifts of image systems 2 and 3 cannot be constrained well, but are consistent with that of the DSFG.Fig. C1C shows the best-fit mass distribution of the lens, which is consistent with a galaxy group morphology with distinct multiple mass peaks.The masses enclosed by the circles around the mass peaks A, B, and C are (7.7,7.8, and 2.8) ×10 12 M ⊙ , respectively, while the total mass enclosed by the critical curve in Fig. C1B is ∼ 1.4 × 10 13 M ⊙ .The mass of peak A is most relevant here and is listed in Table 1 for PJ1053+60.The mean magnification of the DSFG as detected by the SMA 850 µm emission is µ ∼ 24 which is also listed in the table.In addition, we present a magnification gradient map in Fig. C1D to show the differential magnification ranging from µ ∼ 3-100 over the range of its images.The foreground lens (or eastern peak of the mass distribution) of PJ1053+60 has an effective Einstein radius of 5.9 ′′ , enclosing a mass of ∼ 1.4 × 10 13 M ⊙ (included in Table 1).With the limited constraints from the three systems providing only local constraints, the global foreground structure remains poorly understood.Nevertheless, the modeling should be sufficient for our lens mass estimate (hence its X-ray contribution) and a simple reconstruction of the source plane image in the SMA band of PJ1053+60 (Fig. 5 middle panel), as described in the main text ( § 3).

Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .
Figure 3. Multi-wavelength montage of the PJ1336+49 field, showing VLA 6 GHz continuum, HST/WFC3 F160W, and SMA 1.2-mm continuum.Zoomed-in insets show the reconstructed source-plane structure.Light from the foreground lensing galaxy is subtracted before reconstructing the HST image.White dotted curves show the model-derived caustic and critical curves from Kamieneski et al. (2023).The synthesized beams for interferometric observations are shown in cyan in the lower left.Identified image systems used in the lens modeling are labeled on the image where they are most prominent.VLA 2" HST

Figure 10 .
Figure 10.XMM-Newton 0.3-7 keV intensity contours overlaid on the (A) ALMA 1.1-mm image of PJ0116-24 and (B) VLA 6-GHz image of PJ1053+60.These X-ray intensity contours (all in units of 10 −2 counts s −1 arcmin −2 ) are at 0.4, 0. 6, 1, 1.6, and 2.4 above a local background of ∼ 1.4 for PJ0116-24 and at 1, 2, 4, and 8 above a local background of ∼ 4.5 for PJ1053+60, while the X-ray images are constructed in the same way as for Fig 6, but smoothed with a Gaussian with FWHM= 3 ′′ .The position of the AGN, identified in the PJ1053+60 field, is marked with a plus in panel B.

Figure 11 .
Figure 11.Spectral data and fitting models of individual targets: simple power law fit for PJ1336+49 (A); double power law fit for PJ0116-24 (B) and PJ1053+60 (C).The two components are also separately plotted in (B) and (C).

Figure 12 .
Figure 12.Corner plot of the model parameters from the fit to the spectrum of PJ1336+49 (Table4).

Figure 13 .
Figure 13.Corner plot of the model parameters from the fit to the spectrum of PJ0116-24 (Table4).

Figure 14 .
Figure 14.Corner plot of the model parameters from the fit to the spectrum of PJ1053+60 (Table4).

Figure 15 .
Figure 15.Joint fit to the PJ0116-24 and PJ1336+49 MOS spectra (A) and pn spectra (B) (Table4).The X-ray spectral model is a power law for PJ1336+49, representing the collective contribution of the HMXBs in the galaxy, while for PJ0116-24 a separate power law is included to account for the contribution from the AGN contribution.These two components are also separately plotted.

Figure 16 .
Figure 16.Corner plot of the best-fit model parameters from the joint-fit to the spectra of PJ0116-24 and PJ1336+49 (Table4).

Figure A2 .
Figure A2.NIR-MIR spectral shape of the AGN candidate in the vicinity of PJ1053+60.This shape is estimated from the flux densities observed in the HST WFC3/IR F160W and the Spitzer/IRAC 4.5 µm, 3.6 µm bands.

Figure
Figure C1.(A)The HST f140W/F160W color image of the PJ1053+60 field, together with the SMA 850 µm intensity (pink) contours.The lensed image systems 2 and 3 are labeled, while the cluster member galaxies used in the lens modeling are marked with magenta stars.(B) The same color image as in panel A, but together with the (magenta) critical curve of the lensing.(C) The same color image as in (A), but overlaid with the marked surface mass density contours in units of the critical mass surface density.The white dashed circles enclose the regions used to estimate the gravitational masses of A, B, and C peaks as (7.7, 7.8, and 2.8) ×10 12 M ⊙ yr −1 .(D) The WFC3/IR F160W image, together with the same SMA 850 µm intensity contours as in (A), but further overlaid with our best-fit magnification contours at the marked levels.

Table 1 .
Parameters of the targets

Table 2 .
XMM-Newton Observation Log Individual exposure numbers are enclosed in the parentheses.
a b Background flare-cleaned exposures of the individual instruments.

Table 3 .
PSF comparison results a The probability for a source to be point-like.