© J.-F. Lestrade et al. 2022

1. Introduction

Sub-millimetre galaxies (SMGs) found in submillimetre- and millimetre-wavelength surveys are forming stars at prodigious rates, sometimes exceeding 1000 M_⊙ yr⁻¹, and are enshrouded in dust with bright far-infrared luminosities > 10¹² L_⊙ (e.g., Hodge & da Cunha 2020). Thus, SMGs have only faint or undetectable optical counterparts but probe the most active cosmic epoch of galaxy formation at z ≅ 2.

Theory predicts that galaxy formation preferentially occurs along large-scale dark matter filamentary or sheetlike overdense regions that are related to initial fluctuations of the density field in the primordial universe (Bond et al. 1996) and formed the ‘cosmic web’ as depicted by cosmological hydrodynamical simulations. It is thought that galaxy growth and high star formation at an early stage of cosmic evolution is critically tied to massive haloes of dark matter at the knots of this network of filaments (Béthermin et al. 2013). A large fraction of the neutral hydrogen from the inter galactic medium falls into these filaments and constantly replenishes the gas reservoirs of galaxies through cold gas accretion; the cold gas streams along the filaments, survives the shocks at the virial radii of the dark matter haloes and infalls onto the galaxies to fuel their star formation and growth (Kereš et al. 2005; Dekel et al. 2009).

Definitive observational confirmation of this mechanism − cold gas accretion along filaments of the cosmic web − is actively sought. Photoionized Lyman-α nebulae around quasars have been observed with morphologies and kinematics suggestive of intergalactic gas filaments and gas infall (Cantalupo et al. 2014; Hennawi et al. 2015; Martin et al. 2015, 2019) but alternative scenarios based on outflows cannot be ruled out. Recently, extended Ly-α emission plausibly powered by collisional ionization along filamentary structures has been observed. Umehata et al. (2019) have discovered two Ly-α filaments more than 1 Mpc in extent within the protocluster SSA22 at z = 3.1 and, remarkably, also detected embedded SMGs and X-ray sources that are distributed all along these two filaments. Bacon et al. (2021) report the discovery of diffuse extended Ly-α emission from redshift 3.1 to 4.5, tracing cosmic web filaments on scales of 2.5 − 4 Mpc in the MUSE Extremely Deep Field, which partially overlaps the Hubble Ultra Deep Field. Also, Daddi et al. (2021) report on a plausible observation of cold gas accretion streams towards the massive potential well of the protocluster RO-1001 at z = 2.91 and on the detection of three filaments apparently connecting to the associated 300 kpc-wide Ly-α nebula. In cosmological simulations, the relationship between SMGs and the network of cosmic filaments of dark matter is debated (Miller et al. 2015).

Nowadays, a region of the sky several hundreds of arcseconds in size, corresponding to several cMpc (comoving) at redshift z > 1, can be efficiently mapped at millimetre wavelengths using a modern camera such as NIKA2 installed on the IRAM 30-metre telescope (Adam et al. 2018). In maps covering such a scale, the spatial distribution of SMGs can be compared to the cosmic web in cosmological simulations. We present NIKA2 data that are a serendipitous contribution to this endeavour.

In 2007, as part of a survey designed to search for new circumstellar discs around 50 stars, the field in the direction of the nearby star GJ526 was mapped in the 1.2-mm continuum waveband using the MAMBO camera installed on the IRAM 30-metre telescope and a quasi-alignment of five sources along with this star was found in projection on the sky. At the time, it was not possible to conclude whether these sources were background SMGs or part of a large, clumpy circumstellar disc seen edge-on and centred on this star (Lestrade et al. 2009). In 2017, we again observed the same field with the NIKA2 camera to test whether these sources had remained fixed and were in the background, or whether they had moved at the same rate and in the same direction as the star, and were clumps in a circumstellar disc. Based on the proper motion of the star, this possible displacement was estimated to be as large as 22.6″ in the south-west direction after ten years.

In this paper, we report our finding with NIKA2 that the MAMBO sources did not move over this ten-year period and so they must be SMGs in the distant background. In Sect. 2 we present the NIKA2 observations, and in Sect. 3 we present the NIKA2 maps at 1.15 mm and 2.0 mm, and discuss the overdensity of sources in the field. In Sect. 4 we discuss the nature of the sources characterized by their spectral energy distributions (SEDs) constrained with our NIKA2 photometry and Herschel supplemental data. Finally, in Sect. 5 we discuss the remarkable spatial distribution of these sources in the 1.15-mm NIKA2 map along a filament-like structure and a possible association with the cosmic web of the cosmological simulations.

2. Observations

The observations were conducted with the NIKA2 dual-band millimetre camera installed on the IRAM 30-metre telescope at Pico Veleta (altitude 2850 m) in Spain (Monfardini et al. 2011; Catalano et al. 2014; Calvo et al. 2016; Bourrion et al. 2016; Adam et al. 2018; Perotto et al. 2020). NIKA2 operates with kinetic inductance detectors (KIDs) arranged into three arrays : two with a maximum transmission around 260 GHz (1.15 mm) with 1140 detectors in each polarization, and one at 150 GHz (2.0 mm) with 616 detectors at a single polarization. The transmission bands of the detectors are broad with a ∼50 GHz full width at half maximum (FWHM) transmission at both wavelengths. The combination of its large field of view (6.5 arcmin) and high sensitivity provides the NIKA2 camera with a high mapping speed (111 and 1388 arcmin² mJy⁻² h⁻¹ at 260 GHz and 150 GHz, respectively, at zero atmospheric opacity; Perotto et al. 2020). The angular resolutions (main beam FWHM) are 17.7 and 11.2 arcsec at 150 and 260 GHz, respectively.

The data presented in this study were obtained during the science verification observations of the NIKA2 camera on 24 April 2017, on 14, 17, 18, 19, and 20 February 2018, and on 14 March 2018, in average weather conditions (atmospheric opacity between 0.10 and 0.25 at 2 mm, and between 0.20 and 0.45 at 1 mm) and at elevations between 33° and 67°. The observing strategy consisted of a series of 9′×4.5′ on-the-fly scans with a scanning speed of 40″ s⁻¹ along two orthogonal orientations (−60°, +30°) in the AZ-EL coordinate system of the telescope to minimize residual stripping patterns in the maps. We defined the map centre at the J2000 coordinates of the star GJ526 as in the initial MAMBO map (RA: 13h45m44.52s and Dec: +14 ° 53′20.6″). These observations are split in 167 scans that correspond to a total of 9.3 h of observations on source. For calibration purposes, skydips were conducted to determine the coupling of each KID to the sky emission. An opacity correction could then be computed along the line of sight of each observation scan.

The time-ordered information (TOI) of each detector was processed in order to remove the contribution of atmospheric and electronic correlated noises and to obtain the maps at both wavelengths following the procedure initially described in Adam et al. (2014), and modified in Ponthieu et al. (in prep.). In short, an iterative process between maps and TOI was carried out until the maps did not change anymore (less than 0.1σ) after three iterations. For each sub-scan, the first four harmonics at 2 mm (five at 1 mm), each made of one sine and one cosine, were fit to the TOI with a least-square method and removed. In more detail, there was no masking of the source but rather a downweighting of large signal-to-noise excursions. This downweighting becomes more accurate when the map is close to the final map after three iterations. On the first iteration, the initial map was set to zero. The process includes some inhomogeneous weighing of the data samples and so leads to some large spatial scale filtering. Thus, the procedure is designed for optimized detection of point-sources and not extended sources. Nonetheless, the more aggressive filtering and the avoidance of masks in this procedure induce a filtering effect that has to be determined by simulations of point sources injected into the data. Thus, we established that a photometric correction of 25% that is constant for the range of S/Ns displayed by the detected sources in our study must be applied. It should be noted that this iterative part of the pipeline used to process the data is an addition to the standard pipeline used by Perotto et al. (2020). Finally, we add that the noise of each sample was propagated to the map of each scan with a standard inverse variance weighing. A global rescaling of the noise for each map was done to make the standard deviation of S/Ns in noise map unity (see Sect. 3), in order to measure an accurate S/N for point sources.

Point sources were detected by fitting the map with a Gaussian point spread function, in other words effectively applying a noise-weighted matched-filtering of the map with a fixed width Gaussian. Position accuracy is better than 3 arcsec (for each coordinate) for sources detected above 4σ. The photometric correction of 25% described above was applied to the fit amplitude. The flux density scale was based on specific observations of the primary calibrator Uranus to calibrate all individual detectors of each array and on frequent standard observations of the secondary calibrator MWC349 to check the stability of this calibration over a NIKA2 run of several days. The resulting absolute photometric accuracy is of the order of 15%. Although the weather conditions were not optimal, the expected NIKA2 sensitivity at zero opacity (30 ± 3 and 9 ± 1 mJy$.\sqrt{s}$ at 1.15 mm and 2.0 mm, respectively; Perotto et al. 2020) was recovered once the opacity was taken into account.

3. NIKA2 maps

The NIKA2 maps of the GJ526 field 560″ × 560″ in size are shown in Fig. 1 at 1.15 mm and in Fig. 2 at 2.0 mm. Millijansky point-like sources are apparent at both wavelengths. In order to determine the detection threshold, we analysed the noise distribution in computing S/Ns at all pixels of the jack-knifed map at each wavelength. In such a map of the noise, the sources were eliminated after subtraction between individual maps of consecutive pairs of scans with the two orthogonal scanning orientations and finally the co-addition of all pairs. At 1.15 mm, we found that the S/N distribution is Gaussian, centred almost on zero and has a standard deviation of almost unity, as expected for normally distributed noise (see Fig. 3). However, at 2 mm, although the distribution is also centred close to zero and has a standard deviation close to unity, there is an excess at the positive tail that makes it deviate slightly from Gaussian, despite careful reprocessing (see Fig. 4). This was unexpected since the atmospheric fluctuations are less severe at this wavelength. Close inspection of the data showed that electronic correlated noise in the 2.0 mm sampling modules of the instrument was higher than usual during the observations. With this caveat in mind, source extraction and photometry were carried out with caution.

Fig. 1. GJ526 field, mapped with NIKA2 in the 1.15-mm continuum waveband (colour scale: mJy beam−1, contours: S/N). Source IDs are those in Table 1. The curve traces the filament-like structure identified across the map. S/N contours are −3, +3, +4, +5, +6, .... (the negative contours are dashed white lines and can be seen better when the map is zoomed in the electronic version. The S/N = −4 level was processed in plotting the map but no such negative feature was apparent. For the noise rms per pixel, we refer the reader to Fig. 5. Pixel size is 2 arcsec. J2000 coordinates of the map centre are 13h45m44.52s and +14 ° 53′20.6″. The filled black circle in the bottom right corner shows the HPBW.

Fig. 2. Same as for Fig. 1, but observed at 2.0 mm.

Fig. 3. S/N histogram of the jack-knifed map at 1.15 mm. The mean μ and standard deviation σ, as well as the minimum and maximum of the S/Ns, are indicated in the figure. In red, a Gaussian not fit but set with zero mean, standard deviation unity, and its integral summing up to the total number of S/Ns, is overplotted to show that the theoretical distribution of S/Ns for Gaussain noise satisfactorily matches the histogram of the data in the map.

Fig. 4. Same as for Fig. 3, but observed at 2 mm. The slight excess found at the positive tail is discussed in the text.

To set the source detection threshold, we searched for any possible detections in the jack-knifed maps using the same fitting procedure of a Gaussian point spread function as described in Sect. 2. Clearly, any statistically significant detection in this noise-only map would be spurious. At 1.15 mm, we find that there are nine spurious detections with S/N > +3 and ten with S/N < −3, evenly distributed over the jack-knifed map, and none with |S/N|≥4. Similarly, at 2.0 mm, there are eight spurious detections with S/N > +3 and nine with S/N < −3, and none with |S/N|≥4. We conclude that the threshold S/N = 4 is a reliable criterion for a robust source detection at both wavelengths.

In Table 1 we list the 21 sources with S/N ≥ 3 found in the 1.15-mm NIKA2 map, among which ten are unambiguously detected according to our criterion (S/N ≥ 4) and have flux densities between ∼2 and ∼5 mJy. Statistically, about two additional sources should be real among the remaining eleven sources (3 ≤ S/N < 4) since only nine positive detections are expected to be spurious above +3σ, as explained above. Evidently, these potential sources at 1.15 mm cannot be identified unless additional inputs are included. For instance, ID 8 is a 3.6σ source in the 1.15-mm NIKA2 map but can also be identified in the MAMBO map at slightly above 3σ (see map in Fig. 1 of Lestrade et al. 2009); thus ID 8 is likely a real source at the level of 4.7σ in quadratically combining the two S/Ns. Among the unambiguously detected sources at 1.15 mm, only two are also unambiguously detected at 2 mm (ID 16 and ID 9), and four additional sources have S/Ns of between 3 and 4 at 2.0 mm (see Table 1). This deficit of detections at 2 mm can be explained by the dust emissivity index β of SMGs that is between 1 and 2, or sometimes steeper (see also Sect. 4). Finally, there is a source which is slightly above 4σ at 2.0 mm (north-east of map and see Table 1) but has surprisingly no 1.15 mm counterpart. As already mentioned, electronic correlated noise higher than standard during the observations must explain this deviation from normal statistics in the 2.0 mm map, which we consider not as reliable as the 1.15-mm map.

Finally, we address flux boosting that alters the measured flux densities provided in Table 1. In a flux-limited survey, it is well known that sources with a low S/N are detected at flux densities systematically higher than their true flux densities and this effect is exacerbated when the source population increases in number with decreasing flux (e.g., Coppin et al. 2005). We estimated flux boosting with 1000 simulations of the sky in our NIKA2 maps following the prescriptions by Scott et al. (2008) with their semi-Bayesian approach. We find that it is statistically less than 15% for a S/N above 4, and steadly increases up to 50% for a S/N of 3, similarly at 1 mm and 2 mm. We chose not to correct downwards the quoted flux densities in Table 1.

3.1. Nature of the MAMBO sources

As already mentioned, the field around GJ526 was mapped in 2007 with MAMBO. Five sources were found quasi-aligned and within 100″ from the star GJ526 coordinates, and they were reported in Table 2 of Lestrade et al. (2009). They were hypothesized to be the clumps of a large, broken disc seen almost edge-on surrounding this star. With our new NIKA2 observations, four of these five MAMBO sources were detected anew in the NIKA2 map at 1.15 mm, and are IDs 7, 9, 10, and 16 in Table 1 herein. The fifth MAMBO source − MM134543+145317, a 4.3σ MAMBO detection in Table 2 of Lestrade et al. (2009) − was also identified in the 1.15-mm NIKA2 map but only at the level of 2.8σ and so was not retained in our final NIKA2 source list. Such a difference in detection levels would arise from a S/N underestimate and overestimate of 0.75σ by NIKA2 and MAMBO, respectively.

Astrometric and photometric comparisons between the MAMBO map and the NIKA2 map, which is 40% deeper, are given in Table 2 herein. The coordinate differences between the four MAMBO sources at epoch 2007 (Table 2 of Lestrade et al. 2009) and their NIKA2 counterparts at epoch 2017 in Table 1 herein are only a few arcseconds, and so much less than the displacement of 22.6″ in the south-east direction inferred from the star motion. Rather, these coordinates differences are consistent with telescope pointing errors (rms < 3″, see Table 18 of Perotto et al. 2020) combined with astrometric errors for sources with low S/Ns. Thus, we conclude that the MAMBO sources have not moved in concert with the star and are not clumps of a disc associated with GJ526. They must instead be SMGs in the background, probably unlensed owing to their low flux density level.

The MAMBO flux densities in Table 2 of Lestrade et al. (2009) are comparable with the NIKA2 flux densities at 1.15 mm in Table 1 herein to better than 1.5 times the quadratically combined uncertainties of the two measurements, except for source ID 10 for which the difference is 3.1 times the combined uncertainty. This source would have faded from 6.3 ± 1 mJy at epoch 2007 to 2.97 ± 0.62 mJy ten years later if this difference were real. This is physically unrealistic for an SMG and we favour calibration difficulties with MAMBO and/or NIKA2, inherent to low S/N sources. It should be noted that the differences in flux densities due to differences in reference frequencies of the two cameras, 250 GHz for MAMBO and 260 GHz for NIKA2, are negligible.

3.2. Counterparts in the optical, infrared, and radio domains

To search for counterparts to the NIKA2 sources in the optical, mid-infrared, far-infrared, and radio domains, we inspected images of the Sloan Digital Sky Survey (SDSS), the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the Herschel Observatory (Griffin et al. 2010; Poglitsch et al. 2010) and the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). Owing to the NIKA2 position uncertainties, we found no convincing counterpart in SDSS, WISE, or NVSS. In NVSS, with a 3σ limit of 1.5 mJy at 1.4 GHz, this is not surprising unless it is from an AGN.

More interesting were our findings in the Herschel archive. Seven NIKA2 sources have supplemental photometric data from SPIRE (250 μm, 350 μm and 500 μm) and two sources have PACS data (100 μm and 160 μm). For SPIRE, the flux densities of three sources (IDs 7, 10, and 18) are taken from the Herschel/SPIRE Point Source Catalogue (HSPSC; Schulz et al. 2017)¹ and are reported in Table 3 herein, along with their HSPSC identifiers. We also uploaded the SPIRE observation ID_obs 1342234791 of the Herschel Open-Time Key project DEBRIS (Matthews et al. 2010) to determine the flux densities of three additional sources in Table 3 (IDs 14, 17, and 21). Their coordinate differences between SPIRE and NIKA2 are, at most, 1.5 times the quadratically combined astrometric uncertainties of the two instruments (∼5″). Additionally, we determined SPIRE upper limits for the other NIKA2 sources since they are within the SPIRE maps of this programme and not close to the map border of lesser quality. We measured a representative rms in an area where there is no apparent source in the maps to estimate the SPIRE 3σ upper limits of 34 mJy, 33 mJy, and 31 mJy at 250 μm, 350 μm, and 500 μm, respectively. For PACS, we uploaded the observation ID_obs 1342213079 of the same Open-Time Key project and found that ID 7 is detected and ID 10 is not. We estimated the flux densities of ID 7 using aperture photometry (radius = 30″) and found S_ν = 10 ± 2 mJy at 100 μm and S_ν = 52 ± 3 mJy at 160 μm. The 3σ upper limits of ID 10 from this observation are 6 mJy at 100 μm and 9 mJy at 160 μm. The other NIKA2 sources are not part of the PACS maps (ID 18 is at the border of the maps and the data are unusable).

3.3. Background source number count in the GJ526 field

In Sect. 3.1, we concluded that there are ten sources unambiguously detected in the 1.15-mm NIKA2 map with flux densities larger than 2 mJy. This can be usefully compared to the cumulative source number counts in other fields at a similar wavelength. The COSMOS field and the Lockman Hole North field have been surveyed at 1.2 mm with MAMBO/IRAM 30m (Greve et al. 2004; Bertoldi et al. 2007; Lindner et al. 2011) and at 1.1 mm with AzTEC/ASTE (Scott et al. 2010; Hatsukade et al. 2011). Also, part of the GOODS field has been surveyed at 1.1mm with ALMA (Franco et al. 2018).

The noise rms in the NIKA2 maps is not uniform, and increases radially from the centre, as shown in Fig. 5, because of the scanning strategy used for the observations. However, as it can be seen in this figure, the noise rms increases by less than 10% within a radial distance of 150″ from the map centre and so we first consider this central part of the map as the region of interest with a quasi-uniform noise rms. Practically, the statistical uncertainty of the flux density of a source in this central region is ∼0.50 mJy beam⁻¹ at 1.15 mm. Thus, the detection threshold of 4σ corresponds to a flux density of 2 mJy at this wavelength. Lindner et al. (2011) combined their observations with other single dish studies to report the cumulative source number count N(S_1.2 mm ≥ 2.0 mJy) = 500 ± 200 sources deg⁻². Franco et al. (2018) report the lower formal value $N(S_{1.1\mathrm{mm}}\ge 1.8\mathrm{mJy})={209}_{-119}^{+178}$ sources deg⁻² with ALMA at a similar wavelength. It is known that the number counts derived from interferometric data are lower than the number counts derived from single-dish data because of source blending (Karim et al. 2013; Béthermin et al. 2017). Therefore, consistently with our NIKA2 observations, we adopted the single dish count of Lindner et al. (2011) for our probability calculation. There are six sources (IDs 7, 9, 10, 13, 16, and 18) that are unambiguously detected (S/N > 4) in the central region of the 1.15-mm map with quasi-uniform rms (r < 150″), while 2.7 sources are expected with the density 500 sources deg⁻². The Poisson probability of exactly six sources with S_1.2 mm > 2 mJy being found in this region, given this mean surface density, is as low as 3.7%. It may be argued that the likelihood of finding an excess of sources in the GJ526 field is actually enhanced because this field is one in a collection of 50 other similar-sized fields that were mapped in the initial MAMBO survey (Lestrade et al. 2006, 2009). Although it is relevant to recall this context, a clear-cut conclusion cannot be drawn because the other fields were observed with significantly less integration times than GJ526, and so are less deep. We discuss now the more intriguing spatial distribution of the sources of this overdensity along a filament-like structure.

Fig. 5. Radial variation in the noise rms per pixel in the 1.15-mm map (pixel size = 2″). At 2mm, the variation is the same but reduced by a factor of roughly three. The variation is isotropic to a good approximation over the maps.

3.4. Filament-like structure

Now using the whole map, it is remarkable that seven sources (IDs 3, 6, 9, 10, 16, 18, and 25) are spatially distributed along an arc in projection on the sky across the whole map, as highlighted in Fig. 1. The first six sources are unambiguously detected with flux densities larger than 2 mJy at 1.15 mm, while the last one (ID 25) is a 3.3σ possible detection with a flux density of 2.61 ± 0.79 mJy at this wavelength (Table 1). We can estimate the probability that such a quasi-alignment of sources occurs by chance using small, consecutive sectors all along the arc with the small opening angle 2δ, as sketched in Fig. 6, and in which the averaged source number surface density μ is N(S_1.2 mm ≥ 2 mJy) = 500 sources deg⁻², as referenced in Sect. 3.3. For any triplet of sources along the arc, the probability that the third source does not deviate more than ±δ in degrees from alignment with the first two sources is the Poisson probability μ × πr² × 2δ/360 ° )exp(μ × πr² × 2δ/360° that exactly one source is found within the sector centred on the second source, and having a radius r, which is the distance between the second and third sources. In applying this formula for the first triplet in the NIKA2 map (source IDs 3, 6, and 9), the probability of a chance occurrence of quasi-alignment within ±10° is ∼500 sources deg⁻² × π(70″/3600″)² × 20 ° /360 ° = 3.3% with r₁ = 70″ for the distance between source IDs 6 and 9. Similarly, the probability is 0.6% for the second triplet (IDs 6, 9, and 10) with r₂ = 30″ between source IDs 9 and 10, it is 10.1% for the third triplet (IDs 9, 10, and 16) with r₃ = 130″, it is 4.1% for the fourth triplet (IDs 10, 16, and 18) with r₄ = 80″, and it is 13.0% for the fifth triplet (IDs 16, 18, and 25) with r₅ = 150″. These triplets defined five consecutive sectors (only two are represented in Fig. 6) that are independant, and so the resulting probability is the product of the five probabilities estimated above, that is $\sim \prod _{i=1}^5\mu \pi r_i^{2}2\delta /360°$ (the exponential is approximatively unity in our application), and this combined probability is as low as 1.1 × 10⁻⁷. Therefore, it is very unlikely that this filament-like structure has occurred by chance. This is true whatever is our concern that 50 similar-sized fields have been observed in the initial MAMBO survey. In fact, to account for this, the probability can be boosted 50 times and is 5.5 × 10⁻⁶ (i.e. still very low).

Fig. 6. Estimation of the probability for the chance alignment of sources.

4. Spectral energy distribution

Physical quantities of an SMG can be estimated from a model of its SED fit to photometric data. In an SMG, copious amounts of diffuse small dust grains absorb the intense UV radiation field of massive stars and re-radiate in the far-infrared. Its intrinsic SED peaks at a wavelength corresponding to dust temperature typically between 20 K and 60 K, and has a slope steeper than the Rayleigh-Jeans regime because of the dust emissivity in the (sub-)millimetre domain. For an SMG with bolometric luminosity L_bol at redshift z, its intrinsic SED f_ν is shifted towards longer wavelengths and the detectable flux density at an observed frequency ν is (Blain et al. 2002, Eq. (3)):

$$\begin{array}{c}\hfill S_{\nu }=\frac{1+z}{4\pi D_L^2}{L}_{\mathrm{bol}}\frac{f_{\nu (1+z)}}{\int f_{{\nu }^{\prime }}\mathrm{d}{\nu }^{\prime }},\end{array}$$(1)

where D_L is the luminosity distance to redshift z (Weinberg 1972, Eq. (15.3.24)) and can be computed using the on-line calculator by Wright (2006)² in standard ΛCDM cosmology³. In this equation, the intrinsic SED f_ν is from a modified black-body to account for dust emissivity in the far-infrared and sub-millimetre. In other words, it is a Planck function at dust temperature T_dust multiplied by the attenuation factor 1 − exp(−(ν/ν_c)^β) with the critical frequency ν_c at which the dust opacity reaches unity. In our study, we fixed ν_c to the standard value 3000 GHz (100 μm; Blain et al. 2003).

For our sources, photometry data were limited to the flux densities from NIKA2 and Herschel in Tables 1 and 3, and only to the SPIRE upper limits for some of them. Fitting the SED of each source to such a limited amount of data cannot yield a reliable photometric redshift because of the known degeneracy between redshift and dust temperature when only far-infrared and sub-millimetre data are available (e.g., Sect. 2.10.1 in Blain et al. 2002). Instead, we made the most of the data at our disposal in testing whether or not they are compatible with the hypothesis that all sources are at the same redshift. To this end, we constrained T_dust to the range from 20 K to 60 K (e.g., Chapman et al. 2005; Strandet et al. 2016), and constrained all sources to be at the same median redshift z = 2.5 of the phenomenological model of SMGs by Béthermin et al. (2015) for a survey depth similar to our NIKA2 observations at 1.15 mm. Figure 7 shows that the photometric data are compatible with the hypothesis that all sources are at the same redshift, and thus are in the same cosmic structure. This remains true whether or not flux boosting is applied. Evidently, it is only when their spectroscopic redshifts are measured that this hypothesis could be confirmed. As far as the three parameters varied in the model, we provide their acceptable ranges for each source in Fig. 7; it is from 26 K to 60 K for T_dust (keeping in mind our initial constraint), it is from 2.5 × 10¹² L_⊙ to 23 × 10¹² L_⊙ for L_bol, and it is from 1.0 to 2.0 for β. These luminosities are in the upper range of star-forming ultra-luminous infrared galaxies (ULIRG) in the local universe (10¹² ≤ L_IR/L_⊙ ≤ 10¹³ L_⊙). The brightest source, ID 18, is even formally in the range of hyper-luminous infrared galaxies (HyLIG). In addition, in Fig. 7 we provide the star formation rates (SFRs) between 330 and 2255 M_⊙ yr⁻¹, as derived from the luminosity L_bol using SFR(M_⊙ yr⁻¹) = 1.1 × 10⁻¹⁰ × L_bol(L_⊙) (Hayward et al. 2014).

Fig. 7. SEDs of all NIKA2 sources of interest. Each SED is from a single-component modified black-body adjusted to the Herschel and NIKA2 data by varying the dust temperature Tdust, the bolometric luminosity Lbol, and the dust emissivity index β, and keeping the redshift fixed at z = 2.5 for all sources. Empty markers mark the upper limits. The seven sources in the candidate cosmic filament discussed in this study are labelled in red. ${\chi }_{\nu }^2$ is reduced with the number of degrees of freedom ν (see details in Sect. 4).

5. Discussion

The distinctive quasi-alignment of the seven sources, IDs 3, 6, 9, 10, 16, 18, and 25, across the 1.15-mm NIKA2 map is qualitatively reminiscent of the imprint of a cosmic filament in cosmological hydrodynamical simulations (e.g., Illustris (Springel et al. 2018), SIMBA (Davé et al. 2019), NewHorizon (Dubois et al. 2021), and Uchuu (Ishiyama et al. 2021)). Quantitatively, these simulations provide spatial gauges that characterize the distribution of dark matter haloes that host galaxies. These gauges can be usefully compared to our observations. First, it is noticeable that the candidate filament stretches across the whole NIKA2 maps 560″ in size and so extends at least over ∼4 Mpc (comoving) using the scaling factor of 8.3″ kpc⁻¹, which depends only weakly on redshift for 1 < z < 6². Such an observed extent is larger than the size of a typical cluster of galaxies and is comparable to cosmic filaments in cosmological simulations where they extend from 3 Mpc to 15 Mpc (see e.g., NewHorizon snapshot at z = 2 in Fig. 1 of Dubois et al. 2021). Second, it is noticeable that the five angular separations r_i between our NIKA2 sources given in Sect. 3.4 (completing with r₀ = 20″ between IDs 3 and 8) correspond to linear separations between ∼0.25 Mpc and ∼1.25 Mpc, using the same scaling factor. Such observed separations are also comparable to cosmological simulations where haloes hosting galaxies are typically separated by 0.1–1.5 Mpc (see e.g., the NewHorizon snapshot mentioned above). This match of spatial gauges between our NIKA2 observations and cosmological simulations supports our suggestion that the seven SMGs, which are quasi-aligned, bright (L_FIR > 10¹² L_⊙), and possibly at the same redshift, may trace a cosmic filament.

A more detailed comparison between our observations and cosmological simulations is not an easy task. The first simulation that has been specifically motivated to study SMGs is the Gadget-2 simulation by Davé et al. (2010) based on a scenario of gas-rich satellite infall rather than intergalactic gas infall along cosmic filaments. However, the spatial distribution of SMGs of this simulation is not presented. In a more recent study derived from the SIMBA simulation (Davé et al. 2019), Lovell et al. (2021) also focus on SMGs and in their Fig. 4, they provide a 140 cMpc sized map of the sky projection of all simulated galaxies with SFR > 20 M_⊙ yr⁻¹ and within 0.1 < z < 10. This simulated map is as deep as 0.25 mJy at 850 μm (corresponding to SFR = 20 M_⊙ yr⁻¹) and is convolved with the SCUBA2 beam. Their resulting spatial distribution of SMGs is clearly non-uniform and structured along filaments. Further comparison will require a simulation zoomed to match the smaller scale of our NIKA2 maps (∼5 cMpc) and evolved to match our NIKA2 source redshift, when measured.

6. Conclusions

With the NIKA2 dual-band camera at the 30-metre IRAM telescope, we made two deep maps of a relatively large field of ∼90 arcmin² in the direction of the star GJ526, simultaneously in its 1.15-mm and 2.0-mm continuum wavebands with a 1σ sensitivity of 0.45 mJy beam⁻¹ and 0.15 mJy beam⁻¹, respectively, and with 9.3 h of on-source observations owing to the high mapping speed of the NIKA2/IRAM 30-metre system. In total, ten background millijansky sources are unambiguously detected with S/N > 4 and must be SMGs.

In projection on the sky, the spatial distribution of seven of these SMGs along a filamentary structure that completely crosses both maps is remarkable and has only a low probability of chance occurrence (1.1 × 10⁻⁷). We show that all available photometric data, though limited, are compatible with a model of their SEDs in which they are at the same redshift z = 2.5. Of course, this must be confirmed with spectroscopic redshift measurements. At this stage, we can only speculate that this structure could trace a cosmic filament of dark matter, as is predicted by theory and apparent in cosmological hydrodynamical simulations. In assuming the scaling factor of 8.3 kpc/″ at high redshift with standard cosmology for our maps, the extent of this candidate cosmic filament is at least 4 Mpc and the separations between SMGs are a fraction of Mpc comparable to distances between dark matter haloes in high-resolution cosmological simulations. Again, our interpretation is speculative at this stage and a full assessment requires measurements of the spectroscopic redshifts of the SMGs in this candidate filament. Deep spectroscopic follow-up observations of their molecular (CO) or atomic (CI) lines redshifted into the ALMA or NOEMA bands can be envisioned (Inoue et al. 2016; Hashimoto et al. 2018; Solimano et al. 2021; Bing et al. 2022; Birkin et al. 2021). An attempt to optically identify them would require both deeper optical observations and interferometric millimetre-wavelength continuum observations to improve their positions to sub-arcsecond accuracy Finally, the NIKA2 Cosmology Legacy Survey underway at the IRAM 30-metre telescope (Bing et al. 2022), designed to map fields more than 30 times larger than the GJ526 field, will be of great interest to further probe the spatial distribution of SMGs in the millimetre-wavelength domain and address their relationship with the network of cosmic filaments of dark matter in simulations. In particular, establishing this relationship would support the cold gas accretion mechanism to fuel galaxy growth. This endeavour is also pursued with imaging of diffuse Ly α emission along filamentary structures (e.g., Umehata et al. 2019).

Acknowledgments

This work is based on observations carried out under project number 133-17 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We would like to thank the IRAM staff for their support during the campaigns. The NIKA2 dilution cryostat has been designed and built at the Institut Néel. In particular, we acknowledge the crucial contribution of the Cryogenics Group, and in particular Gregory Garde, Henri Rodenas, Jean Paul Leggeri, Philippe Camus. This work has been partially funded by the Foundation Nanoscience Grenoble and the LabEx FOCUS ANR-11-LABX-0013. This work is supported by the French National Research Agency under the contracts “MKIDS”, “NIKA” and ANR-15-CE31-0017 and in the framework of the “Investissements d’avenir” program (ANR-15-IDEX-02). This work has benefited from the support of the European Research Council Advanced Grant ORISTARS under the European Union’s Seventh Framework Program (Grant Agreement no. 291294). This research has made use of Herschel data; Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

All Tables

All Figures

Fig. 1. GJ526 field, mapped with NIKA2 in the 1.15-mm continuum waveband (colour scale: mJy beam−1, contours: S/N). Source IDs are those in Table 1. The curve traces the filament-like structure identified across the map. S/N contours are −3, +3, +4, +5, +6, .... (the negative contours are dashed white lines and can be seen better when the map is zoomed in the electronic version. The S/N = −4 level was processed in plotting the map but no such negative feature was apparent. For the noise rms per pixel, we refer the reader to Fig. 5. Pixel size is 2 arcsec. J2000 coordinates of the map centre are 13h45m44.52s and +14 ° 53′20.6″. The filled black circle in the bottom right corner shows the HPBW.

In the text

	Fig. 2. Same as for Fig. 1, but observed at 2.0 mm.
In the text

Fig. 3. S/N histogram of the jack-knifed map at 1.15 mm. The mean μ and standard deviation σ, as well as the minimum and maximum of the S/Ns, are indicated in the figure. In red, a Gaussian not fit but set with zero mean, standard deviation unity, and its integral summing up to the total number of S/Ns, is overplotted to show that the theoretical distribution of S/Ns for Gaussain noise satisfactorily matches the histogram of the data in the map.

In the text

	Fig. 4. Same as for Fig. 3, but observed at 2 mm. The slight excess found at the positive tail is discussed in the text.
In the text

	Fig. 5. Radial variation in the noise rms per pixel in the 1.15-mm map (pixel size = 2″). At 2mm, the variation is the same but reduced by a factor of roughly three. The variation is isotropic to a good approximation over the maps.
In the text

	Fig. 6. Estimation of the probability for the chance alignment of sources.
In the text

Fig. 7. SEDs of all NIKA2 sources of interest. Each SED is from a single-component modified black-body adjusted to the Herschel and NIKA2 data by varying the dust temperature Tdust, the bolometric luminosity Lbol, and the dust emissivity index β, and keeping the redshift fixed at z = 2.5 for all sources. Empty markers mark the upper limits. The seven sources in the candidate cosmic filament discussed in this study are labelled in red. ${\chi }_{\nu }^2$ is reduced with the number of degrees of freedom ν (see details in Sect. 4).

In the text