Survey of near-infrared diffuse interstellar bands in $Y$ and $J$ bands. I. Newly identified bands

We searched for diffuse interstellar bands (DIBs) in the 0.91$<\lambda<$1.33 $\mu$m region by analyzing the near-infrared (NIR) high-resolution ($R=20,000$ and 28,000) spectra of 31 reddened early-type stars ($0.04<E(B-V)<4.58$) and an unreddened reference star. The spectra were collected using the WINERED spectrograph, which was mounted on the 1.3 m Araki telescope at Koyama Astronomical Observatory, Japan, in 2012--2016, and on the 3.58 m New Technology Telescope at La Silla Observatory, Chile, in 2017--2018. We detected 54 DIBs -- 25 of which are newly detected by this study -- eight DIB candidates. Using this updated list, the DIB distributions over a wide wavelength range from optical to NIR are investigated. The FWHM values of the NIR DIBs are found to be narrower than those of the optical DIBs on average, which suggests that the DIBs at longer wavelengths tend to be caused by larger molecules. Assuming that the larger carriers are responsible for the DIBs at longer wavelengths and have the larger oscillator strengths, we found that the total column densities of the DIB carriers tend to decrease with increasing DIB wavelength. The candidate molecules and ions for the NIR DIBs are also discussed.


INTRODUCTION
Corresponding author: Satoshi Hamano satoshi.hamano@nao.ac.jp Diffuse interstellar bands (DIBs) are ubiquitous absorption features that are detected in the visible to nearinfrared (NIR) spectra of reddened stars. Their absorbing matter is considered to be gas-phase carbonaceous molecules such as carbon chain molecules, polycyclic aromatic hydrocarbons (PAHs), and fullerenes.
However, there is no critical evidence to confirm this. The ionized buckminsterfullerene (C + 60 ) was almost confirmed as the carrier of the five NIR DIBs via a comparison between the astronomical spectra and the spectrum of the gas-phase C + 60 obtained from laboratory experiments Walker et al. 2016;Campbell & Maier 2018;Cordiner et al. 2019), although the identification is still under discussion, because the intensity ratios of the two strong bands were quite variable after the correction of the contaminated stellar line and the detection of the weak bands was not robust (Galazutdinov et al. 2017a(Galazutdinov et al. , 2021. The gas-phase absorption spectra of some neutral PAH molecules have also been measured and compared to DIBs (Gredel et al. 2011;Salama et al. 2011). Although this comparison did not result in the detection of corresponding absorption bands in the astronomical spectra, Gredel et al. (2011) and Salama et al. (2011) succeeded in setting stringent upper limits for the abundances of specific neutral PAHs in some translucent clouds. Further laboratory experiments are required to test the so-called PAH-DIB hypothesis (Leger & D'Hendecourt 1985;van der Zwet & Allamandola 1985).
Spectroscopic surveys of DIBs with high signal-tonoise (S/N) ratios over a wide range of wavelengths are essential for the identification of DIB carriers via a comparison of the DIB spectra with experimentally obtained absorption spectra of the candidate molecules. Some spectroscopic surveys in the optical wavelength range have successfully detected hundreds of DIBs within a detection limit of a few mÅ in equivalent widths (EWs) (Jenniskens & Desert 1994b;Weselak et al. 2000;Tuairisg et al. 2000;Galazutdinov et al. 2000;Hobbs et al. 2008Hobbs et al. , 2009Fan et al. 2019). However, DIB observations in the NIR range have been limited, because the performance of NIR high-resolution spectrographs (in terms of spectral resolution and sensitivity) has been lower than that of optical spectrographs. In addition, many telluric absorption bands in the NIR range prevent the detection of weak absorption features. Recent progress in NIR spectroscopy has enabled searches for DIBs in the NIR spectra over the last decade (Geballe et al. 2011;Cox et al. 2014;Hamano et al. 2015;Elyajouri et al. 2017;Galazutdinov et al. 2017b;Lallement et al. 2018).
We have conducted a comprehensive survey of NIR DIBs using the high-sensitivity high-resolution NIR WINERED spectrograph (Ikeda et al. 2016). In our first study (Hamano et al. 2015, hereafter, H15), we observed 25 stars and successfully identified 15 new DIBs in the range of 0.91 < λ < 1.32 µm. Moreover, we found that the EWs of some NIR DIBs are highly correlated with each other, which suggests that their carriers have similar molecular properties. In our subsequent study (Hamano et al. 2016, hereafter, H16), we investigated the environmental dependence of NIR DIB strengths using the high-quality spectra of seven bright stars in the Cygnus OB2 association, which is one of the most massive clusters or OB associations veiled with large interstellar extinction (the extinction of the most reddened member, No. 12, reaches A V = 10.2 mag; Wright et al. 2015;Whittet 2015). Owing to the large extinction and high flux density in the NIR wavelength range, we could detect even the weakest DIBs reported in H15 with high precision. We found that the NIR DIBs are not correlated with the column densities of the line-of-sight C 2 molecules, which suggests that the DIB carriers are not distributed in the molecular clouds traced with C 2 molecules. We also found that the carrier of DIB λ10504 would be destroyed by the strong UV radiation in the Cyg OB2 association, whereas the DIB carriers of other strong DIBs survive in this environment, which suggests differences in carrier properties, such as the ionization potential and dissociation energy.
In this series, we explore the properties of the DIBs in the Y and J bands, using the large quantity of highquality WINERED data. The most important DIBs, the C + 60 bands, which are covered by the spectrograph, will be comprehensively investigated in order to reveal the properties of C + 60 in the interstellar medium. It is also of great interest to explore the relation of other DIBs with C + 60 . We will also study the correlations between NIR DIBs and optical DIBs, which have been extensively investigated. Through this series, we aim to constrain the carriers of DIBs in the Y and J bands.
In this first study of the series, we updated the catalog of DIBs in the 0.91-1.33 µm range, using a larger sample and higher quality of spectra than those of H15. We analyzed the high-S/N spectra of 32 objects, comprising of an unreddened reference star (β Ori) and 31 reddened stars in the range 0.04 < E(B − V ) < 4.58. To find very weak DIBs, we included large-extinction lines of sight of the Cyg OB2 association (2.2 < A V < 10.2 mag; Wright et al. 2015) and Westerlund 1 (Wd 1) cluster (8.5 < A V < 17 mag; Damineli et al. 2016). In comparison with H15, the sensitivity and wavelength stability of WINERED were improved by the upgrade (Ikeda et al. 2022). The accuracy of the removal of telluric absorption lines is also improved (Sameshima et al. 2018). With these improvements, we can detect much weaker DIBs than those in H15. The updated catalog of DIBs in the Y and J bands will be the basis of this series of papers.
The remainder of this paper is organized as follows. In Section 2, we describe our observations and targets. In Section 3, our data reduction procedures are described. In Section 4, we describe the newly detected DIBs, and in Section 5, we discuss the distributions of DIBs from the optical wavelength to the NIR range and the carriers of the newly found DIBs. Finally, we present a summary of the study in Section 6.

Observations
The data were collected using the high-resolution NIR echelle spectrograph, WINERED (Ikeda et al. 2016(Ikeda et al. , 2022, which uses a 1.7 µm cutoff 2048 × 2048 HAWAII-2RG IR array. We used the WIDE mode, with a 100 µm slit, which includes the wavelength range of 0.91-1.35 µm, with a spectral resolving power of R ≡ λ/∆λ = 28,000 or ∆v = 11 km s −1 . WINERED was mounted on the F/10 Nasmyth focus of the 1.3-m Araki telescope at Koyama Astronomical Observatory, Kyoto Sangyo University, Japan (Yoshikawa et al. 2012), from 2012 to 2016 December. From 2017 to 2018, WINERED was mounted on the ESO 3.58 m New Technology Telescope (NTT) at the La Silla Observatory, Chile. The pixel scales of WINERED were 0 .8 pixel −1 and 0 .27 pixel −1 for the Araki telescope and the NTT, respectively. Most of the data were obtained with the telescope dithered by 30 arcsec (Araki telescope) or 10 arcsec (NTT) -the so-called "ABBA" sequence. In the few cases in which the seeing was not good, we alternately obtained the spectra of the object and sky (the "OSO" sequence). The telluric-standard A0V-A2V type stars were observed at airmasses that were similar to those of the targets.
The observations made with the Araki telescope were primarily conducted in 2014 January and 2014 August-October. Only the observation of β Ori was conducted on 2016 February 11. The typical seeing size at the Koyama Astronomical Observatory is approximately 3-5 arcsec. Some of the data have already been published in H15 and H16. Immediately prior to conducting the observations in 2014 August, we installed an H-band blocking filter to decrease the thermal leak in the 1.7-1.8 µm range. The new filter successfully decreased the background noise, and it yielded spectra with high-S/N ratios. However, the filter was bent by the tight mechanical mounting, which caused a slight off-focus of the light on the IR array. Consequently, the spectral resolving power was reduced to R = 20, 000, which is lower than the nominal value, for the observations made during 2014 August-October (since then, the problem has been solved).
The observations using the NTT were conducted in 2017 July (ESO program ID: 099.C-0850(B)). The typ- ical seeing size was approximately 1 arcsec. Compared with the spectra obtained using the Araki telescope, the telluric absorption lines of the NTT spectra were considerably weaker, owing to the lower humidity and higher elevation of the La Silla Observatory. Figure 1 shows the comparison of the spectra of the telluric-standard stars. For the spectra obtained using the Araki telescope, the telluric absorption lines strongly depend on the season in which the observations were made. Very strong absorption lines of atmospheric water vapor can be seen in the Araki spectra obtained in August, which is the season of high humidity and temperatures in Japan. The difference in the intensities of the telluric absorption lines affected on the DIB search. The wavelength ranges in which the water vapor lines are strong, such as 0.93 < λ < 0.95 µm, 1.11 < λ < 1.17 µm and 1.33 < λ < 1.35 µm, were not available for the DIB search in the Araki spectra obtained in August.

Targets
The targets and observation information are shown in Table 1. We observed 31 reddened early-type stars and one 1 reference star without any interstellar extinction (β Ori). Figure 2 shows the distribution of E(B − V ) and the spectral types of our targets. The reddened stars cover a wide range of spectral types, from A3 to O4.5; however, our observations lack of the spectral-type coverage of the reference star: only β Ori (B8Iab). To avoid the misidentification of the stellar absorption lines as the DIBs, we also utilized the model spectrum synthesized using SPTOOL (Y. Takeda, private communication) as a reference for the stellar spectra. Our targets include a wide range of color excesses, 0.04 < E(B − V ) < 4.58, which enable us to find the DIBs that originate from interstellar clouds. We comment on some targets below.
HD183143 is the most extensively observed star for DIB studies, owing to its large flux density and DIB strength, which is a result of its large extinction (E(B -V ) = 1.27). The DIB properties of this star are well known, owing to previous surveys of DIBs in optical (e.g., Herbig 1975;Herbig & Leka 1991;Jenniskens & Desert 1994b;Tuairisg et al. 2000;Galazutdinov et al. 2000;Hobbs et al. 2009;Cox et al. 2014) and NIR ranges (Cox et al. 2014). The distributions and properties of DIBs over a wide wavelength range are discussed in Section 5.1 using the data of HD183143.
The Cyg OB2 association is one of the most massive clusters in our galaxy. We examined the spectra of six stars from Cyg OB2, which were also investigated in H16. The distance to the Cyg OB2 is estimated as 1669 ± 6 pc based on Gaia Data Release 2 data (Orellana et al. 2021), and it has a very large extinction (Wright et al. 2015). In H16, we detected DIBs with very large EWs in the spectra of the Cyg OB2 members. In addition, various gaseous environments can be traced using this target, because the gas clouds in the lines of sight of Cyg OB2 have a complex structure, consisting of a diffuse component and a patchy dense component (Whittet 2015). Therefore, Cyg OB2 is an ideal target for finding very weak DIBs.
Wd 1 is a massive young stellar cluster with a mass of ∼ 10 5 M (Negueruela et al. 2010) located at a distance of d = 4.12 +0.66 −0.33 kpc from the Sun (Beasor et al. 2021). The typical extinction of the cluster members is A V = 11.4 ± 1.2 (Damineli et al. 2016). W7 and W33, which are our targets, are very luminous B supergiants in the cluster (Clark et al. 2005). The color excesses of these two stars were calculated by using the stellar intrinsic colors of Cox (2000). Owing to these color excesses being the highest among our targets, the DIBs are the strongest; thus, these stars of Wd 1 are the best targets for detecting new weak DIBs.  Notea Telluric-standard stars used for the correction of the telluric absorption lines. b The S/N ratio per pixel after the division by a telluric-standard spectrum at the center of the Y band.

Data reduction
The obtained raw data were processed using the pipeline that was developed to reduce the WINERED data.
The pipeline automatically produces the wavelength-calibrated one-dimensional spectra from the raw data. The pipeline conducts the subtraction of sky frames, the subtraction of scattered light, flat-fielding, bad pixel interpolation, the correction of echellogram distortion, spectrum extraction, wavelength calibration, and continuum normalization. A brief description of the pipeline software is provided in Ikeda et al. (2022). A more detailed description will be given by S. Hamano et al., in preparation. Since most of the optical components of the WINERED spectrograph are maintained at room temperature, the wavelength of the spectrum of each frame shifted slightly with the changes in ambient temperature during the observations. The relative wavelength shifts among the multiple frames of a target are corrected in the pipeline by aligning the wavelengths of the input frames with that of the frame with the highest count. However, the absolute wavelength shift from the dispersion solution was left in the pipeline-reduced spectra. In this study, the shift was measured using the cross-correlation between the model spectra of telluric absorption and the pipeline-reduced spectra. Because we obtained high-S/N spectra of early-type stars, for which the stellar lines are fewer than for late-type stars, we could measure the wavelength shifts with a high accuracy. The wavelengths of the spectra were recalibrated with the shift values.
The telluric lines of the target spectrum were removed using the method detailed in Sameshima et al. (2018). In the method, we made synthetic telluric spectra using molecfit Kausch et al. 2015). Using these spectra as a reference, we identified stellar features in the standard spectrum and fitted them with multiple Gaussian curves, which were synthesized to construct a stellar absorption-line spectrum. Because the telluricstandard stars are A0V-A2V dwarf stars, the metal lines of telluric-standard stars are not complex and can be fitted well with multiple Gaussian curves. The observed standard spectra were divided by the synthesized spectrum, to cancel out intrinsic stellar features, and the resulting spectra were used to correct the telluric absorption lines of the targets. A detailed description of this procedure is provided in Sameshima et al. (2018).

DIB measurement
Before the search for new DIBs, we measured the parameters of the DIBs that were found using the origi-nally developed software. The spectrum files, the restframe DIB wavelengths, and the guess velocities of the DIBs were inputted to the software. First, the software cut out the input spectrum around the DIBs in the range of ±500 km s −1 and normalized it by fitting a Legendre function to the surrounding continuum region. Then, the S/N ratio per pixel was calculated from the standard deviation of the continuum region. Using the S/N ratio and the normalized spectrum of the telluricstandard star, the uncertainty of each pixel of the input spectrum was calculated. The absorption regions, for which the depths by the peak are lower than the continuum level (=1) by 2σ, are then automatically searched for in the normalized spectrum. The region in which the center velocity is the closest to the input velocity is picked up as the DIB region. If the absorption region cannot be found within a ±50 km s −1 range from the input velocity, the DIB is judged as a nondetection.
For the detected DIBs, the EWs and their uncertainties, the central wavelengths and corresponding heliocentric velocities, the FWHM in a wavelength scale, the wavelengths of the absorption peak, and the depth were measured. For the undetected DIBs, the 3σ upper limits of the EWs are calculated. The DIB EWs (W ) were measured from the normalized spectra as follows: where x i is the i th pixel, N is the total number of summed pixels, I n (x i ) is the normalized flux at x i , and ∆λ is the wavelength width per pixel. The EW uncertainty was then calculated using the following equation: where σ stat and σ cont are the statistical uncertainty and the systematic uncertainty from continuum fitting, respectively. σ stat was calculated as follows: where δI n (x i ) is the uncertainty of the normalized flux at x i . σ cont was estimated with the rms shift method, using the statistical uncertainties (Sembach & Savage 1992). In the case of a nondetection, the upper limit of the EW is calculated using 3σ W with Eq. (2), by setting the integration range from the typical width of the DIBs and the input velocity.
The central wavelength of each band, λ c , was measured as the weighted average for the overall band as follows: where τ (x i ) is the optical depth at x i . A formulation was adopted from a preceding survey of DIBs in the optical wavelength range in Fan et al. (2019). For a DIB with an asymmetric profile, the central wavelength measured using Eq. 4 differs from the wavelengths at the absorption peak. The FWHMs of each DIB were determined by calculating the difference between the wavelengths where the depths are half of the peak depth.
All the results of the continuum normalization, the search for DIBs, and the DIB measurements were checked by eye, and if some of the procedures had failed, owing to the systematic noise by telluric absorption lines and the blending of the other features, the fitting parameters of the continuum normalization and/or the integration range were adjusted appropriately.
The velocities of the line-of-sight interstellar clouds are necessary for the search for the DIBs in the spectrum. We do not have reliable interstellar features in the 0.91 < λ < 1.33 µm range to determine the lineof-sight velocities of the clouds, such as the K I line at λ = 7698.9645Å, which has frequently been used to investigate the velocity profiles of the interstellar absorption. In this study, we use the DIB λ10780, which is relatively narrow and strong in this wavelength range; therefore, it is likely to be best suited for the velocity measurement. First, the rest-frame wavelength of DIB λ10780 is determined by comparing the velocities calculated by DIB λ10780 with those measured using the K I line. Among our 31 targets, the K I line profiles of 10 targets have been obtained in previous studies. We used the average of the velocities of the components, weighted by the column densities for the targets toward which multiple-velocity components are detected using a K I line. Table 2 shows the velocities and references of the K I velocities. We determined the restframe wavelength of DIB λ10780 as λ rest = 10780.6Å by minimizing the difference between the K I velocities and the λ10780 velocities. In the calculation, we used only six targets, for which the FWHMs of DIB λ10780 are lower than 1.5Å, in order to reduce the systematic uncertainty in the rest-frame wavelength, owing to the blending of multiple-velocity components. The FWHMs of DIB λ10780 are also listed in Table 2. Figure 3 shows the comparison between the spectra of DIB λ10780 and K I 7699 for the six targets.
Our value, 10780.6Å, is slightly larger than the previous values of 10780.3 ± 0.2Å (Cox et al. 2014) and 10780.46 ± 0.10Å (H15), which were determined using a Gaussian fitting. The difference in the rest-frame wavelength is caused by the difference in the method of the DIB wavelength measurements because the λ10780 profile is asymmetric. We compared the wavelengths mea-sured using our method (Eq. 4) with a Gaussian fitting. It was found that the central wavelengths of DIB λ10780 measured using a Gaussian fitting are smaller than those calculated with Eq. (4) by ∼ 0.2Å, which is comparable to the difference between the rest-frame wavelength of λ10780 in this study and those in the previous studies.
We further checked the robustness of the rest-frame wavelength of DIB λ10780 by comparing the velocities measured using a C + 60 band, DIB λ9577. Although DIB λ9577 is not suited for a velocity measurement, owing to the contamination of the strong telluric absorption lines and the intrinsically broad profile, the rest-frame wavelength is measured in the laboratory experiments as 9577.0 ± 0.2Å (Campbell & Maier 2018). Galazutdinov et al. (2017a) measured the rest-frame wavelengths of C + 60 bands λλ9577 and 9632 for 19 reddened stars, and showed that the mean wavelengths, 9577.0 and 9632.2 A, matched with the wavelengths measured in the laboratory experiments but with a large scatter. The right panel of Figure 4 shows the comparison between the velocities. It was found that the velocities measured using λ10780 and λ9577 are consistent, which suggests that the rest-frame wavelength of λ10780, 10780.6Å, is robust for determining the line-of-sight velocities of the interstellar clouds of our targets.

DIB search
We used the reduced spectra of all the targets to search for new DIBs in the range of 0.91-1.33 µm. Our strategy for distinguishing DIBs from stellar absorption lines was to compare the heavily reddened targets to reference stars without any extinction. Initial candidate DIBs were acquired by comparing the spectra of the most heavily reddened stars -HD183143, HD168625, HD168607, and the stars of Cyg OB2 and Wd 1-to those of the reference star, β Ori. Since the spectral types of some highly reddened stars are earlier than that of β Ori, the strengths of several lines change from object to object. The clearest difference between late-B and early-B/O stars is the line-strength ratio of He II to He I. Thus, we compared the DIB candidates detected in the Cyg OB2 members, but not in HD183143 against the spectra of HD 210191 (B2III), to confirm that they were not stellar absorption lines. Although there are intervening clouds toward HD 210191 (E(B -V ) = 0.04), this star can be used as a semi-standard star, owing to the weakness of the strongest DIBs in the range (e.g., λ10780 < 5 mÅ). Note that the spectral resolution difference between the reference stars (R = 28, 000 or ∆v = 11 km s −1 ) and a subset of the reddened stars (R ∼ 20, 000 or ∆v ∼ 15 km s −1 ) did not affect the  . Comparison between the velocities determined by the DIB λ10780 (λrest = 10780.6Å) and those with a K I line (left) and the C + 60 DIB λ9577 (right). The stars with a FWHM10780 < 1.5Å are shown with red diamonds, and those with a FWHM10780 < 1.5Å are shown with black circles. DIB search, because most of the DIBs and stellar lines are broader than 30 km s −1 , which is broader than the velocity widths of the instrumental profiles . To find new DIB candidates, we first shifted the spectra of the heavily reddened stars into alignment using the DIB λ10780 velocities (Table 2), and we then extracted candidates by comparing the spectra of the reddened stars with those of the reference stars. The candidate search was conducted by eye, while the spectra of telluric-standard stars were simultaneously checked, to avoid misidentification, owing to residual features from the strong telluric absorption lines. The rest-frame wavelengths of the DIB candidates were measured using the shifted spectra of the reddened stars.
After this initial search, we searched for absorption in the wavelengths of the DIB candidates toward all the targets within a range of −30 < v 10780 < 30 km s −1 using the software described in §3.2. If any absorption bands with similar wavelengths to those of the DIB candidates were detected with EWs at the 5σ level, we measured the parameters of the DIBs as in §3.2; in the case where no absorption bands were detected, we calculated the upper limits. In cases where stellar absorption lines, residuals of telluric absorption lines, or spurious features significantly affected the DIB candidates, we evaluated neither the EW nor the upper limit. After checking the measurements of the DIB candidates for all targets, the parameters of each detected DIB candidate were checked. Compared with the velocities of DIB λ10780, the DIB candidate absorptions for which the velocities deviated from v 10780 at more than a 3σ level were rejected. Additionally, an FWHM histogram was produced for each DIB candidate, and the DIB candidate absorptions that were peculiarly broad or narrow were also rejected. The DIB candidates that significantly blended with other features were also removed. In the case of a slight blending with other features, the integration range was arbitrarily adjusted. In that case, therefore, the systematic uncertainties caused by the blending should be included in the parameters of the DIB candidates.
The candidates that satisfy the following requirements are regarded as DIBs: (1) the candidates are detected toward more than 10 reddened stars; and (2) the correlation coefficients between the candidates' EWs and E(B − V ) are larger than zero, with p < 0.05. The multipeaked DIBs were treated as a single feature, because we could not determine whether these multiple peaks originated from the intrinsic profile of an identical DIB carrier or the blended profile of multiple DIBs.   Notea The references that confirmed the DIBs for the first time: (1) Joblin et al. (1990); (2)  As a result of the search for new DIBs in §3.3, we detected 34 new DIBs, nine of which were independently reported by Ebenbichler et al. (2022). In addition, we also detected eight candidates. Figure 5 shows the spectra of 54 detected DIBs, 25 of which were newly detected in this study. The spectra of two stars with no or low reddening, β Ori and HD 210191, and the stars with the largest reddening in the Cyg OB2 association (No. 12) and Wd 1 (W 33) are plotted. In Figure 5, the DIBs detected toward W 33 have broader profiles compared to those of Cyg OB2 No. 12, on average, probably because of the Doppler broadening. Tables 3 and 4 summarize the results of the measurements for the detected DIBs and candidate DIBs, respectively. Tables 3 and 4 show the rest-frame wavelengths (column 1), wavenumbers (column 2), numbers of detections (column 3), EWs per a unit of E(B − V ) (column 4), the correlation coefficients with E(B − V ) (column 5), FWHMs in a wavelength scale (column 6), and comments on the blending features (column 7). The DIB rest-frame wavelengths were determined by minimizing the squared difference between the velocities of DIB λ10780 (Table 2) and the DIB velocities, which depend on the rest-frame wavelengths of the DIBs. The EWs per a unit of Tables 3 and 4 are the average of the FWHMs of the detected DIBs with a sigma clipping. To avoid the effect of Doppler broadening in the FWHM calculation, as far as possible, the FWHMs larger than the 1σ level were clipped, whereas the FWHMs smaller than the 2σ level are clipped. In the following subsections, we comment on some specific DIBs and candidates.

C + 60 DIBs
All five absorption bands identified as the electronic band of C + 60 are located in the wavelength range of our observations . Two main bands, λ9577 and λ9632, which were contaminated with the strong water vapor lines, could be detected toward many targets. The latter band is known to overlap with a pair of stellar Mg II lines at 9631.9 and 9632.4Å (Jenniskens et al. 1997;Galazutdinov et al. 2017a;Walker et al. 2017). We did not try to remove the effect of the overlapping Mg II lines because that is outside the scope of this study. Therefore, we note that the measured parameters of λ9632 in Table 3 are affected by Mg II lines. Three weak subbands, λ9348, λ9365, and λ9428, are contaminated with stronger water vapor lines. Therefore, these bands could not be detected at all in the spectra obtained using the Araki telescope, owing to the strong water vapor lines (Figure 1). We tried to analyze the three subbands in the NTT spectra; however, we could not detect them clearly, owing to the intrinsic weakness of the bands, the systematic uncertainties caused by the removal of the telluric lines, and the blending stellar lines. Therefore, we did not include these three bands in Table 3 in this study. The properties of C + 60 in the interstellar clouds and their relations with the other interstellar features, including the newly found NIR DIBs, are of great interest. An analysis dedicated to these C + 60 bands will be conducted in the subsequent study.

λ11691, λ11695, and λ11698
Three DIBs, λ11691, λ11695, and λ11698, were detected in the wavelength range where the telluric lines are strong. Because the profiles of these DIBs are clearly resolved in the spectra of most of the targets, we have treated these features as independent DIBs in this study. However, the EWs of these DIBs are well correlated with each other (r ∼ 0.9). Considering that these DIBs are affected by strong telluric lines, which increase the scatter in the correlation, it is possible that these DIBs originate from the same carrier, and the separations and the relative intensities of the three peaks may be important for constraining the molecular properties of the carrier molecule.

Note-
The electronic bands of C 2 and CN are located in the range of 0.91-1.33 µm (Hamano et al. 2019). In particular, the C 2 bands have many rotational lines in a wide wavelength range, because C 2 can be rotationally excited to higher levels. The rotational lines of the J = 0 − 20 of the C 2 (1,0) and (0,0) bands range between 10130-10300Å and 12070-12332Å, respectively. It was difficult to find weak DIBs in these wavelength ranges of the C 2 bands, because the C 2 bands were strong in the spectra of most of the stars in the Cyg OB2 association and the Wd 1 cluster, which were used for the initial search for new DIBs. Four DIBs-λ12194, λ12200, λ12294, and λ12313-and a DIB candidate λ12110, were detected in the region of the C 2 (0,0) band, while we could not find any new DIBs in the region of the C 2 (1,0) band. In particular, λ12110 is blended with a Q(8) line, which is relatively strong among the C 2 lines. We did not deblend the Q(8) line from λ12110 by simulating the C 2 line profile. Therefore, λ12110 was measured only in the spectra of the stars toward which the C 2 band was not detected or very weak.

DIBs detected in the wing of Pa β
DIBs λ12799, λ12838, and λ12862 were detected on the wing of the broad Paβ absorption line. We removed this broad wing by fitting a polynomial function for measuring the EWs of these DIBs. 4.2. NIR DIBs in the 0.91-1.33 µm range Figure 6 shows the EWs of the NIR DIBs in the range of 0.91-1.33 µm for Cyg OB2 No. 12, toward which all 54 DIBs were detected, except for DIB λ12110, which was blended with the C 2 band. The EWs of 20 of the DIBs in this object that were detected in previous studies (H15, Cox et al. 2014;Joblin et al. 1990;Foing & Ehrenfreund 1994;Groh et al. 2007) are larger than 30 mÅ, whereas the EWs of the weakest DIBs that are newly detected in this study are as small as 10 mÅ. The lack of DIBs in λ < 9500Å, 11100 < λ < 11500Å, and λ > 13300Å in Fig. 6 is a result of the strong telluric absorption lines. Accordingly, there are fewer DIBs in and near these wavelength ranges in Fig. 7, which shows the distribution of the DIB numbers and EWs for Cyg OB2 No. 12. The peaks of the number density distributions in the Y and J bands are at approximately 10500 and 12500Å, respectively, where the telluric absorption lines are very weak. The weakest DIBs in the Y and J bands are also detected in these bins. These results suggest that the detection sensitivity of the NIR DIBs is strongly influenced by the accuracy of the telluric absorption. Figure 8 shows the overall distribution of the central wavelengths of the DIBs, from the optical wavelength range to the K band, where the longest-wavelength DIBs have been detected. This figure primarily uses the DIB EWs of HD183143 obtained by Fan et al. (2019, for 4000-9000Å), those obtained in this study (for 9100-13300Å) and those from Cox et al. (2014, for 15000-18000Å). For the DIBs in the K band and some weak DIBs in the H band, we adopt the DIB EWs for Cyg OB2 No. 5 from Galazutdinov et al. (2017b), because, to date, these DIBs have not been detected toward HD183143. Although the number of DIBs detected in the NIR region has gradually increased over the course of the past decade, previous surveys could not detect DIBs in the range of W < 50 mÅ, in which the majority of the optical DIBs are distributed. In this study, we were able to detect a number of weak NIR DIBs, which has enabled us to discuss the DIB distribution over a wide wavelength range (see Section 5). It will be possible to investigate the profiles of these DIBs, as well as their correlations with gaseous parameters (e.g., H I column density), and to use them for comparison with the laboratory spectra of candidate molecules for new NIR DIBs. These new DIBs will provide new insights into the    Table 3. The spectra are shifted by the velocities that were measured using the DIB λ10780, to align the DIBs. The regions that are contaminated with the strong telluric absorption lines are interpolated from the surrounding regions, which the strong telluric absorption lines do not contaminate. The interpolated and original spectra are shown using the thick and thin lines, respectively. The spectra are shown with gray lines in cases where the DIBs could not be evaluated, owing to the blending of the other features. The detected DIBs are indicated by the yellow areas. The marks with the black thin lines that are sometimes plotted above the spectra show DIBs and candidates that appear in the panels of the other DIBs.    (Joblin et al. 1990;Foing & Ehrenfreund 1994;Groh et al. 2007). The green diamonds show the DIBs detected in Ebenbichler et al. (2022). The light blue triangles show the candidate DIBs in this study. Lower panel: transmittance spectrum of the atmosphere. The spectrum is synthesized using ATRAN. The locations of the Paschen H I lines and C2 Phillips bands are also shown.
mysterious carriers of DIBs. In particular, in our forthcoming study, we will investigate the DIB-E(B −V ) and DIB-DIB correlations.

Distributions of DIB numbers and strength
In this subsection, we discuss the distributions of the DIB numbers and strength, using the DIBs of HD183143, which is most frequently used for the DIB surveys. The DIBs toward HD183143 have been searched for from the optical wavelength range to the NIR range, at high sensitivities. In Section 5.1.1, we discuss the DIB distributions over a wide wavelength range, from the optical wavelength range to the H band. In Section 5.1.2, we attempt to deduce the intrinsic distribution of the DIBs toward HD183143 by considering the biases among the observations in different wavelength ranges, such as the differences in S/N, the spectral resolutions, and telluric absorption lines. . Number density (top) and EW density (bottom) distributions toward Cyg OB2 No. 12. The bin width is 500 A in all panels. The black, gray, and white bars show DIBs with EW > 100 mÅ, 10 < EW < 100 mÅ, and EW < 10 mÅ, respectively.
length range, we have adopted the most extensive catalog of DIBs, which was reported by Fan et al. (2019), who obtained a high-sensitivity (S/N ∼ 1300 at the median) and high-resolution (R = 38, 000) spectrum of HD183143 and detected 472 bands between 4000 and 9000Å. The EWs of the DIBs of HD183143 in the H band (15,000 < λ < 18,000Å) are adopted from Cox et al. (2014). Figure 9 shows the distributions of number density (top), total EW (middle), and the summation of the product of the oscillator strengths and column densities of the DIB carriers (bottom). The summation is calcu-  (Joblin et al. 1990;Foing & Ehrenfreund 1994;Groh et al. 2007). The green diamonds show the DIBs detected in Ebenbichler et al. (2022). The light blue triangles show the candidate DIBs in this study. Gray and black circles show the DIBs in the range of 4000-9000Å (Fan et al. 2019) and the H-band DIBs (Cox et al. 2014), respectively. For the DIBs in the K band and some weak DIBs in the H band, we adopt the DIB EWs for Cyg OB2 No. 5 from Galazutdinov et al. (2017b), because, to date, these DIBs have not been detected toward HD183143. Note that the points of No. 5 are only plotted to show the DIB distribution at the longest wavelength, and cannot be directly compared with the other data from HD183143, owing to the difference in the DIB EWs between HD183143 and No. 5. The gray shaded areas show the wavelength ranges in which we could not search for DIBs, owing to strong telluric absorption. Lower panel: synthetic spectrum of telluric absorption created by ATRAN (Lord 1992). lated as follows: where i is the DIB index within each wavelength bin, f i and N i are the oscillator strength and column density of the ith DIB carrier, respectively, and W i and λ i are the EW and central wavelength of the ith DIB, respectively. In this calculation, an optically thin condition is assumed.
The peak of the number density distribution (Fig. 9, top panel) was around the bins from 6000 to 7000Å. According to Fan et al. (2019), this peak does not reflect the intrinsic DIB distribution, because a significant fraction of the weak DIBs at λ > 7000Å were not detectable, owing to contamination by strong telluric absorption. The number density of the NIR DIBs at 9000 < λ < 13500Å is considerably lower than that of the optical DIBs, which may be a result of both the intrinsic DIB distributions and the differing spectral qualities in terms of S/N, spectral resolution, and telluric lines. In particular, very few DIBs with W < 10 mÅ are detected in the NIR wavelength range.
The difference in the detection limits among the observations has a lower impact on the EW distribution than on the number density distribution because most of the EW contribution comes from strong DIBs (W > 100 mÅ). For example, the 4000-4500Å bin has a very high value, owing to λ4430, which is the strongest DIB ever found. Although this stochasticity primarily affects the distribution, a trend of decreased EWs at longer wavelengths is seen, which may be attributable to an intrinsic distribution of DIB EWs rather than to observational bias. This trend is more clearly seen in the distribution of i f i N i (the bottom panel of Fig. 9). Figure 10 shows the DIB distributions on a wavenumber scale. The wavelength coverage of our observations corresponds to 7400 < ν < 11100 cm −1 . The wavenumber distribution is probably more appropriate than the wavelength distribution for investigating the intrinsic distributions of the DIB numbers and strengths than the wavelength distribution because the transition energies are simply proportional to the wavenumbers. The trends seen in Fig. 9 can also be seen in these distributions, but the differences between the optical and NIR ranges are diluted, because a bin in the wavenumber units covers a larger wavelength range at a longer wavelength. In particular, the EWs are nearly constant across the wavenumber plot, except for a few bins that contain the strongest DIBs. To understand the intrinsic distribution of these DIB properties, it is necessary to consider the effects on the distributions that arise from differences in observational conditions.

Observational biases
Next, we deduced the intrinsic distribution of the DIBs from the observed distribution by minimizing the effects of observational bias. The observed distributions were affected by observational biases on the strength and density of the telluric absorption and the detection limits.
For the detection limit, we used only those DIBs that were deeper than a fixed peak depth (d p ) threshold for all wavelength ranges, to minimize the effects of different detection limits. Since the S/N ratio and spectral resolution of the optical spectrum of HD183143 obtained by Fan et al. (2019) are both higher than those of our NIR spectrum of HD183143, the threshold from the NIR spectrum was set to be d p = 7 × 10 −3 , which is the 5σ detection limit for our NIR spectrum of HD183143. Figure 11 shows the resultant wavenumber density of the numbers, EWs, and i f i N i of the DIBs in the sampled wavelength ranges. The number density exhibits a peak in the 6600 < λ < 6850Å bin, and a declining trend is seen toward both the shorter-and longer-wavelength ranges. This is consistent with the results shown in Fig. 10, and the peak wavelength corresponds to the 1.8 eV transition energy. In contrast, the EW density is approximately constant over the entire range, with the exception of the 6200 < λ < 6500 A bin, in which the second-strongest DIB, λ6284.3, is located. The almost constant wavenumber density of the DIB EWs results in a rapid decrease of i f i N i toward the longer-wavelength limit, as seen in the bottom panel of Figure 11. Therefore, the carrier abundance of the DIBs tends to be reduced at longer wavelengths, suggesting that the NIR DIBs trace molecular species that are less abundant in interstellar clouds and/or have much smaller oscillator strengths than those in the optical region.

FWHM distribution
In this subsection, we discuss the distributions of the FWHMs of the optical and NIR DIBs for HD147889, using the data of Fan et al. (2019) and this study. HD147889, which is considered to be a single-cloud sightline (Siebenmorgen et al. 2020), is a good target for investigating the DIB profiles. Here, we compare the FWHMs of the DIBs in the optical range and those in the NIR range in order to examine the relation between the DIB width and the wavelength. The intrinsic profiles of the DIBs are important for constraining the DIB carriers (e.g., Huang & Oka 2015). If an intrinsic DIB profile reflects the rotational contour of a molecular electronic transition, the FWHM on a wavenumber scale is a function of the rotational temperature (T rot ) and the rotational constant (B ) (Cami et al. 2004). A difference in the FWHM distribution between the optical and NIR wavelength ranges would imply that the carriers producing the DIBs in each wavelength range have different molecular properties (e.g., electric dipole moment and molecular size), on average. Figure 12 shows the histograms of the DIB FWHMs measured in the wavelengths, normalized by wavelength, and measured in wavenumbers for optical DIBs (5000 < λ < 6500Å; Fan et al. 2019) and NIR DIBs (9100 < λ < 13300Å; this study). The bin size for the NIR DIBs is twice as large as that for the optical DIBs, owing to the difference in their numbers. The spectral resolutions in Fan et al. (2019) and in this study were 8 and 11 km s −1 , respectively. The thresholds corresponding to the spectral resolutions were also plotted with the histogram of the normalized DIB FWHMs. In the FWHM distribution on a wavenumber scale, the DIB width is narrower in the NIR range, on average. Because most of the DIBs have FWHMs broader than the spectral resolutions (see the center panels in Figure 12), the difference of the FWHM distribution on a wavenumber scale is probably not caused by the difference of the spectral resolutions.
The detection limits of the observations also affect the FWHM distributions. Since the broader DIBs have shallower depths, it is relatively difficult to detect broader DIBs at the same EW. Therefore, the detected weak DIBs tend to have narrower widths, and the detection limit can change the FWHM distribution. The FWHM distributions of the DIBs with a central depth larger than 0.009 are also shown with the black bars in Fig.  12. The threshold is sufficiently higher than the detection limits of both observations, without a significant decrease in the DIB numbers. The FWHM distributions of the DIB samples that are limited by depth would reflect the intrinsic DIB properties better than that of a full sample would. Other factors of the observations, such as the contamination of the telluric absorption lines and the wavelength coverages of each echelle order, can bias the FWHM distributions; however, it is difficult to evaluate their effects.
In these depth-limited samples, the difference between the wavenumber-scale FWHM distributions is clear: the medians are 2.1 and 1.5 cm −1 for the optical and NIR wavelength ranges, respectively. To investigate the statistical significance of the difference in the FWHM distributions, we conducted a Mann-Whitney U test, which is a nonparametric statistical test with the null hypothesis that the medians of two samples are the same, under the assumption that the shapes of the distributions are identical. The sample sizes of the optical and NIR DIBs are 29 and 27, respectively, which are sufficient for the statistical test. We rejected the null hypothesis, with a p-value of 0.038, which suggests that the difference in the FWHM distributions between the optical and NIR DIBs is statistically significant. The difference in the FWHM distributions between the optical and NIR wavelength ranges may imply that the DIB width on a wavenumber scale becomes intrinsically narrower at longer wavelengths.

General implications
Herein, we comment on the implications for the carriers of the NIR DIBs, based on our results. In Section 5.2, we suggested that the intrinsic FWHMs of the NIR DIBs in units of wavenumber are, on an average, narrower than those of the optical DIBs. If we assume that the intrinsic widths of the DIBs are determined by the rotational constants of the DIB carriers, the NIR DIB carriers should have smaller moments of inertia (i.e., they should be larger in size) than the optical DIB carriers. This is consistent with the relationship between the transition wavelengths and the molecule sizes of conjugated molecules, such as carbon chains, PAHs, and fullerenes. This proportional relation between the transition wavelength and the number of carbon atoms in molecules has been experimentally confirmed for the π − π electronic transitions of polyacetylene chains (Maier 1998). The relation is also seen for the PAH molecules and ions in the calculation of the wavelengths of the PAH bands that was conducted in Ruiterkamp et al. (2005).
In Section 5.1, we suggested that the summation of the DIB column density ( i f i N i ) decreases with wavelength. We cannot determine whether oscillator strength or column density is the primary contributor to this decrease. If we assume that conjugated molecules are a potential carrier, then their oscillator strengths will be in proportion to the total electron number, i.e., the number of carbon atoms (e.g., Halasinski et al. 2003). If the carriers of the DIBs at longer wavelengths are larger in size, as suggested by their FWHM distribution, then the DIBs at longer wavelengths can be considered to have higher oscillator strengths. In this case, column densities associated with longer-wavelength DIBs decrease more rapidly than is observed for i f i N i with increasing wavelength. It follows that larger DIB carriers have lower column densities in interstellar clouds, because large molecules can have a variety of structural patterns and more atoms are necessary to form bigger molecules. In summary, we suggest that DIBs at longer wavelengths tend to be caused by larger molecules. Ruiterkamp et al. (2005) simulated DIB spectra toward a line-of-sight cloud of HD147889, the physical parameters of which have been well constrained, under the assumption that the DIB absorption in this region is caused by interstellar PAHs. Their results suggest that larger PAH ions tend to have transitions at longer wavelengths, and that the ionization fraction of PAHs strongly depends on their size and structure. Metallicity is another key parameter that determines the DIB distribution (Cox & Spaans 2006). Further high-sensitivity observations in both the optical and NIR wavelength ranges toward various lines of sight would put constraints on the DIB carrier properties.

Possible carriers
DIBs are considered to originate from interstellar carbonaceous molecules (Sarre 2008). Specific candidates include carbon chain molecules, PAHs, and fullerenes. The only positively identified carrier to date is ionized buckminsterfullerene (C + 60 ), which has been identified as the carrier of five DIBs at approximately 0.95 µm Walker et al. 2015;. Other fullerenes, such as C + 70 and C 2+ 70 , were also tested by the same authors , but their transitions were not detected in the astronomical spectra. Omont (2016) reviewed the properties of other fullerenes and fullerene derivatives to explore their ability to produce interstellar DIB features. Fullerene compounds may be candidates for DIB carriers. Tomita et al. (2005) obtained NIR spectra of C − 60 in its gas phase at room temperature (300K) and identified three absorption bands at 9382, 9145, and 10460 cm −1 , which corresponded to air wavelengths of 10655, 10931, and 9557Å, respectively. Near the strongest of the three bands (10655Å), we detected a weak DIB candidate at 10,650Å (Table 4). The EW and FWHM of the λ10650 for Cyg OB2 No. 12 are 27±2 mÅ and 2.5Å, respectively. We could not detect DIBs close to the other two weaker bands at 10931 and 9557Å. The band at 10931Å overlaps the H I line, and the band at 9557Å is heavily contaminated by telluric lines; therefore, we could not even set upper limits. Despite the wavelength difference between the main band at 10655Å and λ10650, a small shift such as this may be induced because the absorption band was obtained at room temperature (Tomita et al. 2005), which can broaden and shift bands. Based on the oscillator strength of f = 0.022 (Strelnikov et al. 2015), the EW toward Cyg OB2 No. 12 (27 mÅ) corresponds to N (C − 60 ) = 1.2 × 10 12 cm −2 . Considering the amount of C + 60 in the line of sight of Cyg OB2 No. 12 (N (C + 60 ) = 2.5 × 10 13 cm −2 estimated from the λ9577 band), C − 60 may be present, and it would favor denser clouds with reduced UV intensity. Although it remains a matter of speculation, DIB λ10650 is a potential candidate for the absorption band of C − 60 . Firm identification requires (1) obtaining the gas-phase spectrum of C − 60 at a low temperature, and (2) detecting the other two subbands. Since the detection of the band at 10931Å is nearly impossible, owing to its overlap with the stellar H I line, detecting the band at 9557Å would be more significant.
Ionized PAHs are another class of DIB carrier candidates in the NIR range. Mattioda et al. (2005) obtained the NIR (0.7-2.5 µm) spectra of cations and anions of 27 PAHs (the largest of which was C 50 H 22 ) using matrix isolation spectroscopy. The Ar matrix that they used can broaden and shift the obtained absorption band. Therefore, it was difficult for them to identify DIBs as the absorption bands of ionized PAHs. They demonstrated that strong and narrow bands exist in the NIR absorption spectra of PAH cations and anions. In H15 and in this study, we have detected DIBs whose wavelengths are close to the absorption bands of PAH cations and anions. If the gas-phase spectra of such ionized PAHs can be obtained, the DIBs detected here can potentially be confirmed.

SUMMARY
We explored weak NIR DIBs the in 0.91-1.33 µm range using the NIR high-resolution (R = 20, 000 and 28,000) spectra of 31 reddened stars that were collected using the WINERED spectrograph. Our findings are summarized as follows: 1. The large DIB EWs toward the heavily reddened lines of sight enabled us to detect 54 DIBs, 25 of which were newly discovered by our observations. We also independently detected nine of 12 DIBs newly detected in Ebenbichler et al. (2022). We succeeded in detecting DIBs as weak as 10 mÅ in the NIR range. We found that, as in the optical range, many weak DIBs populate in the NIR range of 0.91-1.33 µm. The wavelength range of 0.91-1.33 µm that has been explored in this study is of great importance for the study of interstellar molecules, because it contains the absorption bands of both small and large carbon molecules, including C + 60 , C 2 , and CN, and many anonymous DIBs.
2. The FWHMs of the NIR DIBs were found to be narrower than those of the optical DIBs, on average. Assuming that the DIB width is determined by the rotational constant, this difference suggests that the DIBs at longer wavelengths tend to be caused by larger molecules, which is consistent with the properties of conjugated molecules.
3. We investigated the distributions of the DIB numbers, EWs, and column densities, according to wavelength (wavenumber), from the optical to the NIR range. The number density of the DIBs peaks at λ ∼ 6600Å and declines toward both the shorter and longer wavelengths. Additionally, the sum of the DIB column densities decreases with increasing wavelength. Assuming that the DIBs at longer wavelengths tend to originate from larger molecules, as suggested by their FWHM distributions, the oscillator strength can be assumed to be larger for DIBs at longer wavelengths. Thus, we suggest that DIBs at longer wavelengths trace less abundant molecules.
4. The comparison of the DIB catalog that has been compiled in this study with the gas-phase spectra of candidate molecules, such as fullerenes, PAHs, and carbon chains, will contribute to the further identification of DIB carriers. As a trial, we checked the NIR absorption spectra of the gasphase C − 60 that was obtained at room temperature, and we found that the DIB candidate λ10650 was close to the main absorption band of C − 60 at 10655 A. The detection of weaker subbands in astronomical spectra and in the laboratory spectra of gasphase C − 60 at lower temperatures will be necessary to confirm the identification of this ion.
It will be of great interest to investigate the correlations among DIBs, particularly the correlations of all DIBs with the C + 60 DIBs at λ = 9577 and 9632Å, to assess the contribution of fullerenes and their associated compounds to DIBs. Moreover, it will also be important to constrain the sizes and structures of the DIB carriers by fitting the molecular rotational contour to the observed DIB profiles (e.g., Huang & Oka 2015). The new list of DIBs in the NIR range produced by this study represent a valuable input to for further investigations into DIB carriers.
We are grateful to the staffs of Koyama Astronomical Observatory and La Silla Observatory for their support during our observations. This study is based on observations collected at the European Southern Observatory under ESO program 099.C-0850(B). We thank Dr. Mitsunori Araki for the useful comments and suggestions.