A Principal Component Analysis of the Diffuse Interstellar Bands

Our goal in this paper is to determine the parameters that drive the variations in the DIB spectrum. Thus, we need to select a sample of sightlines whose physical conditions are reasonably well determined, and that represents the overall observed variations of the DIB spectrum. The first requirement implies that we should restrict ourselves as much as possible to single-cloud lines of sight to avoid having to deal with ill-defined averages throughout multiple intervening clouds. The second requirement stresses the importance of including lines of sight that are as observationally different as possible. Finally, in order to be able to relate any changes to known observables, we need to pick lines of sight for which auxiliary data (e.g., hydrogen column densities, $E(B-V)$, extinction properties, etc) are known and available. In this pilot study, we chose to restrict ourselves to include only a limited number of DIBs, and to lines of sight for which we can find high-resolution spectra that allow us to exclude possible blends with stellar lines.

We started our selection of targets from the thorough and detailed study of elemental depletion in the lines of sight toward 243 stars published by Jenkins (2009). This study critically reviews available literature data for $E(B-V)$, N(H i), N(H₂), N(H), and introduces a depletion strength factor, F_⋆, describing the collective level of elemental depletion in a line of sight. Starting from this sample thus ensures a consistent treatment of the required auxiliary line-of-sight data and allows us to ensure a wide coverage of environmental conditions (to the extent that they can be traced by any of these parameters).

We searched the VLT/UVES and ELODIE archives to find good-quality, high-resolution spectra of these targets. UVES (the Ultraviolet and Visual Echelle Spectrograph) is a high-resolution instrument on the VLT covering wavelength ranges from 3000–4000 Å and 4200–11,000 Å, with maximum resolutions of 80,000 and 110,000, respectively (Dekker et al. 2000). ELODIE is an echelle spectrograph on the 1.93 m telescope at the Observatoire de Haute-Provence in France. ELODIE has a resolution of 42,000 and covers the wavelength range from 3906 to 6801 Å (Moultaka et al. 2004). We found appropriate data for 91 of the Jenkins targets: 43 targets in the UVES database; the remaining 48 from the ELODIE database. In a few rare cases, parts of the spectrum would be of too low quality or simply missing from the data, and in those cases we supplemented our data with archival spectra from the ESPaDOns (Echelle spectropolarimetric device for the observation of stars) instrument on the Canada–France–Hawaii Telescope (CFHT)—a bimodal instrument with maximum resolutions of 68,000 and 81,000 for its spectropolarimetric and non-polarimetric modes, respectively, covering a wavelength range of 3700–10,000 Å (Donati 2003).

To further select only single-cloud lines of sight, we examined the interstellar Na i D lines at 5890 and 5895 Å and the K i lines at 7665 and 7699 Å (see, e.g., Bhatt & Cami 2015, illustrated in Figure 1). For our current study, we consider a sightline to be a single cloud if these interstellar lines show only one dominant component at the spectral resolution of UVES or ELODIE. Thus, a target was still considered to be a single cloud if there were multiple radial velocity components, but one component had significantly stronger features than the others. For example, HD 149757 is known to have two strong radial velocity components at approximately −27 and −15 km s⁻¹, but the one at −15 km s⁻¹ has much larger column densities (Herbig 1968, see also Figure 1). For the targets with only ELODIE spectra, the K i lines are outside the available wavelength range, and thus we could only use the Na i D lines to inspect the number of radial velocity components. An obvious and inherent drawback of using these lines is that they are easily saturated. We searched for Ca i and CH lines too, but in most cases, these were too weak to be seen.

Zoom In Zoom Out Reset image size

With this exercise, we established that 33 of our targets can be identified as single-cloud lines of sight. However, for three of them, the UVES archival spectra have a gap in the wavelength coverage from approximately 5760 Å to 5830 Å and thus two important DIBs, λ5780 and λ5797, are missing. We therefore excluded these targets from our data set. After these considerations, we were left with a sample of 30 single-cloud lines of sight. These targets are listed in Table 1 along with their line-of-sight parameters we will use in this paper (see below).

Table 1.  Basic Target Data

Target	Alt.	R.A.	Decl.	V	$E(B-V)$	N(H i)	N(H2)	f(H2)	F⋆	$\tfrac{{\rm{W}}(\lambda 5797)}{{\rm{W}}(\lambda 5780)}$	vISMa	Referencesa	Data
	Name	(J2000)	(J2000)			(1021 cm−2)	(1020 cm−2)				(km s−1)		Source
HD 15137	⋯	02 27 59.81	+52 32 57.6	7.86	0.24	${1.29}_{-0.40}^{+0.57}$	${1.86}_{-0.12}^{+0.26}$	${0.22}_{-0.06}^{+0.09}$	0.37 ± 0.09	0.30 ± 0.02	−9.58	1	ELODIE
HD 22951	40 Per	03 42 22.65	+33 57 54.1	4.98	0.19	${1.10}_{-0.32}^{+0.35}$	${2.88}_{-0.98}^{+1.48}$	${0.35}_{-0.18}^{+0.27}$	0.73 ± 0.05	0.35 ± 0.02	12.47	1	ELODIE
HD 23180	o Per	03 44 19.13	+32 17 17.7	3.86	0.22	${0.76}_{-0.23}^{+0.26}$	${3.98}_{-1.16}^{+1.64}$	${0.51}_{-0.24}^{+0.33}$	0.84 ± 0.06	0.65 ± 0.04	13.45	2	ELODIE
HD 23630	η Tau	03 47 29.08	+24 06 18.5	2.87	0.05	${0.22}_{-0.07}^{+0.10}$	${0.35}_{-0.12}^{+0.18}$	${0.28}_{-0.15}^{+0.23}$	0.89 ± 0.10	0.16 ± 0.04	16.76	2	ELODIE
HD 24398	ζ Per	03 54 07.92	+31 53 01.1	2.88	0.27	${0.63}_{-0.07}^{+0.06}$	${4.68}_{-1.59}^{+2.40}$	${0.59}_{-0.31}^{+0.46}$	0.88 ± 0.05	0.55 ± 0.02	14.54	2	ELODIE
HD 24534	X Per	03 55 23.08	+31 02 45.0	6.10	0.31	${0.54}_{-0.07}^{+0.08}$	${8.32}_{-0.73}^{+0.80}$	${0.76}_{-0.11}^{+0.13}$	0.90 ± 0.06	0.62 ± 0.04	14.5	5	ELODIE
HD 24760	Per	03 57 51.23	+40 00 36.8	2.90	0.07	${0.25}_{-0.05}^{+0.05}$	${0.33}_{-0.15}^{+0.27}$	${0.21}_{-0.14}^{+0.25}$	0.68 ± 0.04	0.18 ± 0.02	7.06	2	ELODIE
HD 24912	ξ Per	03 58 57.90	+35 47 27.7	4.04	0.26	${1.29}_{-0.24}^{+0.26}$	${3.39}_{-0.99}^{+1.40}$	${0.35}_{-0.15}^{+0.21}$	0.83 ± 0.02	0.26 ± 0.01	11.2	1	ELODIE
HD 27778	62 Tau	04 23 59.76	+24 18 03.6	6.33	0.34	${0.22}_{-0.22}^{+0.55}$	${5.25}_{-0.88}^{+1.06}$	${0.82}_{-0.27}^{+0.45}$	1.19 ± 0.07	0.43 ± 0.03	15.22	2	ELODIE
HD 35149	23 Ori	05 22 50.00	+03 32 40.0	5.00	0.08	${0.43}_{-0.13}^{+0.12}$	${0.03}_{-0.03}^{+0.00}$	${0.02}_{-0.02}^{+0.00}$	0.54 ± 0.11	0.20 ± 0.04	24.09	2	UVES
HD 35715	Ψ Ori	05 26 50.23	+03 05 44.4	4.60	0.03	${0.31}_{-0.13}^{+0.13}$	6 ± 2 × 10−6	4 ± 2 × 10−6	0.66 ± 0.11	0.10 ± 0.04	25.2	1	ELODIE
HD 36822	ϕ1 Ori	05 34 49.24	+09 29 22.5	4.40	0.07	${0.65}_{-0.12}^{+0.13}$	${0.21}_{-0.06}^{+0.09}$	${0.06}_{-0.03}^{+0.04}$	0.74 ± 0.08	0.19 ± 0.04	25.53	1	ELODIE
HD 36861	λ Ori A	05 35 08.28	+09 56 03.0	3.30	0.10	${0.60}_{-0.16}^{+0.16}$	${0.13}_{-0.05}^{+0.08}$	${0.04}_{-0.02}^{+0.04}$	0.57 ± 0.04	0.48 ± 0.04	25.2	3	ELODIE
HD 40111	139 Tau	05 57 59.66	+25 57 14.1	4.82	0.10	${0.79}_{-0.15}^{+0.16}$	${0.54}_{-0.20}^{+0.31}$	${0.12}_{-0.07}^{+0.10}$	0.49 ± 0.04	0.20 ± 0.04	15.29	2	ELODIE
HD 110432	BZ Cru	12 42 50.27	−63 03 31.0	5.32	0.39	${0.71}_{-0.21}^{+0.29}$	${4.37}_{-0.38}^{+0.42}$	${0.55}_{-0.11}^{+0.13}$	1.17 ± 0.11	0.25 ± 0.01	6.8	3	UVES
HD 143275	δ Sco	16 00 20.01	−22 37 18.1	2.29	0.00	${1.41}_{-0.29}^{+0.29}$	${0.26}_{-0.09}^{+0.15}$	${0.03}_{-0.02}^{+0.03}$	0.90 ± 0.03	0.19 ± 0.02	−10.90	2	UVES
HD 144217	β1 Sco	16 05 26.23	−19 48 19.6	2.62	0.18	${1.23}_{-0.11}^{+0.12}$	${0.68}_{-0.09}^{+0.10}$	${0.10}_{-0.02}^{+0.02}$	0.81 ± 0.02	0.11 ± 0.01	−8.95	2	UVES
HD 145502	ν Sco	16 11 59.74	−19 27 38.5	4.13	0.20	${1.17}_{-0.59}^{+0.56}$	${0.78}_{-0.23}^{+0.32}$	${0.12}_{-0.07}^{+0.08}$	0.80 ± 0.11	0.18 ± 0.01	−8.49	2	ELODIE
HD 147165	σ Sco	16 21 11.32	−25 35 34.0	2.91	0.31	${2.19}_{-0.87}^{+0.90}$	${0.62}_{-0.18}^{+0.25}$	${0.05}_{-0.03}^{+0.04}$	0.76 ± 0.06	0.13 ± 0.01	−6.26	2	UVES
HD 147933	ρ Oph A	16 25 35.10	−23 26 48.7	5.02	0.37	${4.27}_{-0.80}^{+0.98}$	${3.72}_{-1.09}^{+1.53}$	${0.15}_{-0.07}^{+0.09}$	1.09 ± 0.08	0.27 ± 0.03	−8.02	2	UVES
HD 149757	ζ Oph	16 37 09.54	−10 34 01.5	2.58	0.29	${0.52}_{-0.04}^{+0.02}$	${4.47}_{-0.75}^{+0.90}$	${0.63}_{-0.17}^{+0.20}$	1.05 ± 0.02	0.50 ± 0.04	−14.98	2	UVES
HD 164284	66 Oph	18 00 15.80	+04 22 07.0	4.78	0.11	${0.42}_{-0.39}^{+0.23}$	${0.71}_{-0.21}^{+0.29}$	${0.25}_{-0.20}^{+0.18}$	0.89 ± 0.18	0.15 ± 0.02	−15.32	1	ELODIE
HD 170740	⋯	18 31 25.69	−10 47 45.0	5.76	0.38	${1.07}_{-0.47}^{+0.59}$	${7.24}_{-1.22}^{+1.47}$	${0.58}_{-0.18}^{+0.22}$	1.02 ± 0.11	0.26 ± 0.01	−12.9	6	UVES
HD 198478	55 Cyg	20 48 56.29	+46 06 50.9	4.86	0.43	${2.04}_{-0.63}^{+0.84}$	${7.41}_{-2.17}^{+3.06}$	${0.42}_{-0.20}^{+0.27}$	0.81 ± 0.05	0.24 ± 0.01	−10.04	2	ELODIE
HD 202904	υ Cyg	21 17 55.08	+34 53 48.8	4.43	0.09	${0.23}_{-0.23}^{+0.21}$	${0.14}_{-0.05}^{+0.07}$	${0.11}_{-0.10}^{+0.12}$	0.39 ± 0.11	0.13 ± 0.05	−12.90	4	ELODIE
HD 207198	⋯	21 44 53.28	+62 27 38.0	5.96	0.47	${3.39}_{-0.50}^{+0.59}$	${6.76}_{-0.59}^{+0.65}$	${0.28}_{-0.05}^{+0.05}$	0.90 ± 0.03	0.53 ± 0.01	−15.28	2	ELODIE
HD 209975	19 Cep	22 05 08.79	+62 16 47.3	5.11	0.27	${1.29}_{-0.38}^{+0.41}$	${1.20}_{-0.41}^{+0.62}$	${0.16}_{-0.09}^{+0.12}$	0.57 ± 0.26	0.31 ± 0.01	−11.39	2	ELODIE
HD 214680	10 Lac	22 39 15.68	+39 03 01.0	4.88	0.08	${0.50}_{-0.15}^{+0.14}$	${0.17}_{-0.04}^{+0.05}$	${0.06}_{-0.03}^{+0.03}$	0.50 ± 0.06	0.34 ± 0.02	−9.2	1	ELODIE
HD 214993	12 Lac	22 41 28.65	+40 13 31.6	5.23	0.06	${0.58}_{-0.18}^{+0.20}$	${0.43}_{-0.14}^{+0.22}$	${0.13}_{-0.07}^{+0.10}$	0.68 ± 0.10	0.17 ± 0.02	−9.44	1	ELODIE
HD 218376	1 Cas	23 06 36.82	+59 25 11.1	4.84	0.16	${0.89}_{-0.26}^{+0.28}$	${1.41}_{-0.48}^{+0.73}$	${0.24}_{-0.13}^{+0.19}$	0.60 ± 0.06	0.28 ± 0.01	−12.65	1	ELODIE

Note. R.A. and decl. are taken from SIMBAD. V, $E(B-V)$, N(H i), N(H₂), and F_⋆ values are from Jenkins (2009), where we assume an uncertainty of ±0.02 mag on the $E(B-V)$ values. $f({{\rm{H}}}_{2})=2N({{\rm{H}}}_{2})/[N({\rm{H}})+2N({{\rm{H}}}_{2})]$ is the fraction of molecular hydrogen.

^aValue and reference refer to the velocity of the dominant interstellar component.

References. (1) This paper, (2) Welty & Hobbs (2001), (3) Bhatt & Cami (2015), (4) Welty et al. (1994), (5) Sonnentrucker et al. (2007), (6) Adams (1949).

Download table as:  ASCIITypeset image

It should be noted that six of the targets in our sample—HD 23630, HD 24534, HD 110432, HD 149757, HD 164284, and HD 202904—are Be stars. Such objects have an intrinsic $E(B-V)$ value relative to non-Be stars of the same spectral type (Schild 1978; Sigut & Patel 2013); the hot, circumstellar gas produces excess emission in the V filter and therefore, their $E(B-V)$ values are overestimated. Furthermore, any dust existing in the circumstellar shell (CS) can contribute to $E(B-V)$, while there is no evidence to suggest that DIBs exist in CS environments (Snow & Wallerstein 1972; Snow 1973; Krełowski & Sneden 1995). Hence, we expect weaker-than-normal DIB strengths relative to $E(B-V)$ for these six targets.

We applied a heliocentric correction to all targets, and then shifted all spectra to their interstellar rest frames. Interstellar velocity components were obtained from the literature, or measured from a known interstellar feature, if not available. These velocities are listed in column 14 of Table 1.

For our purposes, we only wanted to include those DIBs for which EWs can be confidently measured, i.e., with small relative errors. This limits our selection to fairly strong and often narrow DIBs that are as much as possible free of contamination from stellar lines. The only exceptions we made were to include four of the C₂ DIBs—λλ4964, 5513, 5546, and 5769—despite the latter three being quite weak. We thus include a sample of 16 DIBs in our analysis: λλ4428, 4964, 5494, 5513, 5545, 5546, 5769, 5780, 5797, 5850, 6196, 6270, 6284, 6376, 6379, and 6614.

We measured the EWs for these 16 DIBs in all lines of sight. Most of these were straightforward to measure, as they are not heavily contaminated by stellar and/or telluric features. The largest uncertainty regarding these measurements stems from establishing the continuum. To obtain a good estimate of these uncertainties, we used the following Monte-Carlo approach to vary the continuum level and perform direct integration of the spectra (similar to that in Bhatt & Cami (2015); see Figure 2 for illustration). First, we measured the standard deviation of the flux values over a specified, featureless range of data in the vicinity of the feature to estimate the uncertainty on the flux values—i.e., we measured the signal-to-noise ratio (S/N). Next, we selected a point on each side of the feature and defined a continuum baseline by adopting a linear continuum between those two points. We also selected points as our integration limits; note that for consistency, we used the same integration limits for the same feature in all lines of sight. To simulate the process of determining the continuum, we varied the two selected continuum points 1000 times by adding a random number to the flux values selected from a normal distribution with a mean of zero and standard deviation equal to one-third of the measured standard deviation in the featureless continuum; we believe that this represents well how accurately one can position the continuum (which corresponds to determining the mean flux in the adjacent continuum), and we found that this produces reasonable continuum estimates (see Figure 2). Using one full standard deviation produces many continuum points, which are clearly too high or too low and thus result in unrealistic continuum levels. In essence, we thus simulated the entire process of determining a continuum line 1000 times. For each continuum, we then measured the EW. The EW we use in this paper is then the mean of these 1000 measurements, and the standard deviation of these measurements provides the uncertainty.

Zoom In Zoom Out Reset image size

We kept a few precautions in mind when using this method. For instance, the sharp and narrow λ5797 is known to be blended with the broader and shallower λ5795. To avoid measuring a contribution from λ5795, we measured the λ5797 while treating the λ5795 as continuum (see Figure 2). If bad pixels were found within the integration range, that data point was replaced with the average of the neighboring points. For HD 23180, the λ5494 feature was strongly contaminated by a stellar line. We could not resolve the two features to obtain a proper measurement, so instead we adopted the value obtained from higher resolution observations by Bondar (2012).

The two broad DIBs in our sample—λ4428 and λ6284—cannot be measured using the methods described above, because they are heavily contaminated by stellar and telluric features, respectively. For the λ6284 DIB, we first applied a telluric correction using molecfit (version 1.1.0; Kausch et al. 2015; Smette et al. 2015) and then proceeded as for the other DIBs.

The λ4428 DIB is extremely broad, and the region spanned by this DIB is plagued by stellar features. Snow et al. (2002) showed that the intrinsic profile of the band is Lorentzian, and thus, rather than numerically integrating the DIB profile, we preferred to fit a Lorentzian profile to the observations and determine the EWs from the fitted parameters. In addition to our best fit, we also determined Lorentzians that represent the upper and lower envelope of the observed profiles; these then have a different full width at half maximum (FWHM) and central depth (CD) value (while we kept the same continuum). We found the difference between the best-fit values and upper or lower envelope values and determined the uncertainties on EW through error propagation.

We compared our measurements to values found in the literature whenever possible and found a generally good agreement. For instance, 23 out of our 30 lines of sight were also studied by Friedman et al. (2011) and we found a very good correlation between their reported EW values and our own (see, e.g., Figure 3).

Zoom In Zoom Out Reset image size

All of our EW measurements are shown in Table 5 in the Appendix of this paper.

To have a better chance of understanding what drives the variations in the DIB strengths, we need to include other parameters that offer some description of the lines of sight. The Jenkins (2009) study provides a critically reviewed set of measurements that describe the lines of sight we study here and we reproduce some of their values in Table 1.

There are two quantities that are related to the dust in the line of sight. The amount of dust is traditionally characterized by using the color excess, $E(B-V)$, and is often used as a normalization factor for DIB strengths. We adopt an uncertainty of ±0.02 mag for our $E(B-V)$ measurements. Jenkins (2009) furthermore showed that the depletion of different elements in a line of sight can be described by a single parameter, F_⋆, that is a measure for the "total depletion"; individual elemental depletion factors scale with F_⋆. Since depletion may play a role in the formation and/or destruction of the DIB carriers, we include it here in our analysis.

Information on the gas in the line of sight stems from measurements of hydrogen, in particular, the column densities of neutral and molecular hydrogen, N(H i) and N(H₂), respectively. Note that total hydrogen column densities N(H) are not listed, but can easily be calculated as N(H) = N(H i) + 2N(H₂). For HD 23630, a reliable N(H i) value could not be obtained by Jenkins; consequently, he uses synthetically derived N(H) and F_⋆ values, which we adopt for this paper. We use this synthetic N(H) value and the measured N(H₂) value to derive a synthetic N(H i) from N(H) − 2N(H₂). For two other targets, HD 27778 and HD 202904, Jenkins provides only upper limits and best values for N(H i). We take the lower limit to be zero in both cases. For N(H), he lists only upper and lower limits. We adopted N(H)_best = N(H i)_best + 2 N(H₂)_best for these targets, as well as synthetic F_⋆ values. For HD 35149, only an upper limit for N(H₂) is provided. In this case, N(H₂) is quite small, so Jenkins takes N(H) = N(H i) and calculates F_⋆ normally. Because N(H₂) is so small in this case, we take Jenkins' value to be the best value, set the lower limit equal to zero, and set the upper limit equal to the best value. It has been suggested that f(H₂) can be used as an indicator for the amount of interstellar UV radiation that can penetrate since H₂ is dissociated when not shielded (see, e.g., Cami et al. 1997; Sonnentrucker et al. 1997). This could be an important parameter in our study, and we therefore calculated the fraction of molecular hydrogen, f(H₂), for each line of sight from N(H₂) and N(H).

Similarly, it has long been known that W(λ5797)/W(λ5780) is somehow related to the physical conditions in a region of space (Cami et al. 1997; Krelowski et al. 1997); more specifically, it is also suggested to trace UV exposure, where larger values of the W(λ5797)/W(λ5780) ratio correspond to more sheltered (ζ-type) environments. We thus also include the W(λ5797)/W(λ5780) ratio as a line-of-sight parameter in our study.

Our sample is quite diverse in terms of f(H₂): reported f(H₂) values in diffuse clouds run from 0 to ∼0.8 (Snow & McCall 2006); our sample covers a similar range from 0 to 0.824. Our sample also covers a broad range of physical conditions because it includes both archetypal σ and ζ sightlines. For comparison, our W(λ5797)/W(λ5780) values range from 0.10 to 0.65 whereas those calculated from EWs in a large sample of 133 stars by Friedman et al. (2011) range from 0.15 to 0.85. Furthermore, the boundary between σ and ζ clouds is usually placed around W(λ5797)/W(λ5780) = 0.3 and therefore, we include a large sample of both cloud types. For F_⋆, we similarly cover much of the possible range. Jenkins defines F_⋆ such that a line of sight with minimal depletions has a value of 0, and the −15 km s⁻¹ component of HD 149757—the archetype for large depletions—has a value of 1.0. Our values for F_⋆ are found between 0.37 and 1.19. For $E(B-V)$, the situation is slightly more complicated. Because we restrict ourselves to single clouds, we are biased toward objects of lower reddening. Indeed, our $E(B-V)$ values range from 0.00 to 0.47, while it is not uncommon for multiple-cloud sightlines to exceed 1. In particular, $E(B-V)$ values below 0.08 can be problematic; below this threshold, H₂ does not exist in appreciable amounts (Savage et al. 1977) and many DIBs similarly fall below the limit of detection. Within our sample there are six of these low-$E(B-V)$ targets: HD 23630, HD 24760, HD 35715, HD 36822, HD 143275, and HD 214680. Studies by Kumar et al. (1982), Federman et al. (1984), and Kumar (1986) confirm that strong DIBs such as λ5780, λ5797, λ6284, and λ6614 can be detected toward such objects. Moreover, Galazutdinov et al. (1998) confirm that W(λ5797)/W(λ5780) is still sufficient to differentiate between σ and ζ environments, despite the weak reddening. Some weaker DIBs may not appear in these poorly reddened objects; in such cases, sightlines having larger $E(B-V)$ values will still be observable and will dominate the overall trends.

The starting point of our analysis is a set of n variables for which we have m measurements x_i; in our case, we use n = 23 variables (representing DIB EWs and line-of-sight parameters) that are measured for m = 30 lines of sight. Our measurements thus span an n-dimensional parameter space—each dimension corresponding to a different variable—and a full set of measurements for a single line of sight can be represented by an n × 1 vector in this parameter space. Starting from this set of n possibly correlated (and thus not necessarily orthogonal) variables, a PCA finds a coordinate transformation that casts the variables in terms of a new, orthogonal, n-dimensional reference frame:

$$\begin{eqnarray}&&{y}_{i}\,{\hat{{\boldsymbol{n}}}}_{i}^{\prime }={a}_{i,1}{x}_{1}{\hat{{\boldsymbol{n}}}}_{1}+{a}_{i,2}{x}_{2}{\hat{{\boldsymbol{n}}}}_{2}\,+\,...\,+\,{a}_{i,n}{x}_{n}{\hat{{\boldsymbol{n}}}}_{n}\end{eqnarray} \tag{ 1 }$$

where the ${\hat{{\boldsymbol{n}}}}_{{\boldsymbol{i}}}$ represent unit vectors in the original parameter space and ${\hat{{\boldsymbol{n}}}}_{i}^{\prime }$ unit vectors in the new reference frame—these are the so-called principal components (PCs). The coefficients, ${a}_{i,j}$, are constrained such that their squared sums equal one:

$$\begin{eqnarray}&&{a}_{i,1}^{2}+{a}_{i,2}^{2}\,+\,...\,+\,{a}_{i,n}^{2}=1.\end{eqnarray} \tag{ 2 }$$

The coordinate transformation is chosen such that when representing the measurements in the new reference frame, the largest amount of variance in the data occurs along the axis defined by the first PC, and the second PC accounts for as much of the remaining variation in the data as possible, with the constraint that it is orthogonal to PC₁. Each successive PC has the same properties: they are linear combinations of the original variables, and each accounts for as much of the remaining variation as possible, while being orthogonal to (and thus, uncorrelated with) all previous PCs.

The entire coordinate transformation can be written in matrix notation:

$$\begin{eqnarray}&&Y={AX}\end{eqnarray} \tag{ 3 }$$

where X is an n × m matrix containing the original set of data; A is the n × n transformation matrix; and Y is an n × m matrix containing the transformed data points in the new reference frame defined by the set of PCs.³

The PCA obtains the transformation matrix A from the eigenvalues (λ₁, λ₂,..., λ_n) of the covariance matrix of the original variables. The rows of this matrix A are the corresponding eigenvectors; each eigenvector describing the projections of the original variables onto the new coordinate system defined by the PCs. The eigenvalues are furthermore ordered such that λ₁ is the largest and λ_n is the smallest with ${\sum }_{i}^{n}{\lambda }_{i}=n$. Each eigenvalue indicates how many original variables a single PC can account for, i.e., an eigenvalue of two indicates that the corresponding PC carries the same weight as two of the original variables.

One of the key features of a PCA is the ability to reduce the number of dimensions in a data set. Indeed, it is often the case that we can use the results of a PCA to accurately describe the variation in a data set with a reduced number of variables. As such, PCs that account for very little variation are often ignored, and only the leading p components are kept. Selecting a value for p is somewhat arbitrary, but the decision is usually based on one of three criteria: (1) once a certain amount of variation is accounted for (e.g., 90%), all further PCs are ignored; (2) only PCs with eigenvalues greater than one are kept (representing variables that can replace more than one of the original variables); or, (3) if an "elbow" (a sharp, sudden drop-off) is noticed in a so-called "screeplot" (i.e., a plot of eigenvalue versus component number, see, e.g., Figure 6), all PCs beyond the elbow are ignored (Jolliffe 2014). If p components are kept, then the last n − p rows of A are eliminated, and Y becomes a p × m matrix. At that point, we effectively express the original n variables with only p new variables, while keeping as much of the variance in the data as desired.

The variables x_i that we want to include in our analysis are listed in column 1 of Table 2. Since PCA is concerned with linear combinations of variables, it is redundant to include N(H) in addition to N(H i) and N(H₂); however, whereas N(H i) and N(H₂) measure the column densities of individual species, N(H) approximates the total amount of gas in the line of sight. Recognizing this distinction, we choose to keep N(H) in addition to N(H i) and N(H₂).

Table 2.  Input Variables Used in the PCA, and Their Mean Values and Standard Deviations

Variable	Mean	Standard
Name		Deviation
(xi)	(${\bar{x}}_{i}$)	(${s}_{{x}_{i}}$)
x1 = W(λ4428)	645.9	350.9
x2 = W(λ4964)	6.8	5.7
x3 = W(λ5494)	5.6	4.1
x4 = W(λ5513)	3.9	3.9
x5 = W(λ5545)	5.9	4.3
x6 = W(λ5546)	2.8	2.2
x7 = W(λ5769)	2.6	2.2
x8 = W(λ5780)	131.5	77.0
x9 = W(λ5797)	37.7	27.8
x10 = W(λ5850)	15.6	13.6
x11 = W(λ6196)	13.5	8.3
x12 = W(λ6270)	20.7	13.9
x13 = W(λ6284)	149.4	84.1
x14 = W(λ6376)	9.3	7.4
x15 = W(λ6379)	24.7	18.0
x16 = W(λ6614)	52.5	35.2
x17 = $E(B-V)$	0.20	0.13
x18 = N(H i)	1.0 × 1021	9.2 × 1020
x19 = N(H2)	2.4 × 1021	2.6 × 1020
x20 = N(H)	1.5 × 1021	1.2 × 1021
x21 = f(H2)	0.27	0.23
x22 = F⋆	0.75	0.21
x23 = $\tfrac{{\rm{W}}(\lambda 5797)}{{\rm{W}}(\lambda 5780)}$	0.2907	0.1567

Note. Entries for i = 1, 2,...16 refer to equivalent widths of named features in mÅ. Units for the remaining entries are those specified in Table 1. Standardized variables z_i correspond to the x_i, e.g., z₈ is the standardized equivalent width of the λ5780 DIB.

Download table as:  ASCIITypeset image

There are a few underlying assumptions associated with PCA, which we will discuss. The first is that raw data is comparable in units and magnitude; otherwise, the results of PCA will be more strongly influenced by variables that have larger variances (i.e., variations among N(H) values that are on the order of 10²¹ will completely dominate those of $E(B-V)$, which are on the order of 10⁻¹). To address this problem, we standardize each variable prior to performing PCA:

$$\begin{eqnarray}&&{z}_{i,j}=\displaystyle \frac{{x}_{i,j}-{\bar{x}}_{i}}{{s}_{{x}_{i}}}\end{eqnarray} \tag{ 4 }$$

where i refers to the ith variable and j refers to the jth observation. Essentially, we subtract the mean and divide the residuals by the standard deviation, and use the resulting variables ${z}_{i,j}$ as input to the actual PCA. Columns 2 and 3 of Table 2 show the mean and standard deviation of each of our variables, respectively. After this process, the standardized variables (z₁, z₂,... z₂₃) each have a mean of zero and a standard deviation of one.

A second assumption of PCA is that all values are accurately known. This assumption is challenged by the fact that our measurements have uncertainties. Since PCA does not consider uncertainties, we need a way to quantify how reliable our results are. Since the PCs represent vectors in a multi-dimensional parameter space, this is not straightforward.

We therefore used a Monte-Carlo (MC) simulation to estimate the uncertainties on our results. We first performed PCA using our measured values. These results provide the PCs that we discuss in the sections that follow. Next, we generated 1000 perturbed data sets, by adding random noise to each observation, selected from a normal distribution with a mean of zero and a standard deviation corresponding to the error bar on each measurement. By definition, none of our measured quantities can be negative, so if a perturbed measurement was negative after the addition of random noise, the value was replaced with zero. For each of these new data sets, we standardized the variables⁴ and performed a separate PCA. The result is 1000 unique transformation matrices. For each entry in the transformation matrix, we found the mean (${\bar{a}}_{i,j}$) and standard deviation (${s}_{{a}_{i,j}}$) among the 1000 unique matrices. Then we computed the upper and lower limit on each value in the matrix A according to ${\bar{a}}_{i,j}\pm {s}_{{a}_{i,j}}$. We compared these upper and lower limits to the original, unperturbed results and ultimately, were left with positive and negative errors on each entry in the transformation matrix (i.e., each component of each eigenvector). We applied the same methods to obtain uncertainty measurements on each entry in Y, the matrix of transformed data points.

A third assumption associated with PCA is that each variable follows a normal distribution. We tested our variables for normality and found that, in some cases, the distributions were not normal. The assumption of normality, however, is not critical and the results of PCA should still be valid if this condition is not met, so we decided to continue with the analysis.

Finally, the results of PCA can be sensitive to outliers. Despite including some unusual targets (e.g., HD 147933 is found within a dark cloud, and HD 36822 and HD 36861 are found near H ii regions), we do not remove any targets from our sample. In order to probe a wide range of interstellar conditions, it is important to include these abnormal lines of sight: any conclusions made about the DIBs through this analysis should extend to these sightlines as well. The only potential complication is for the six Be stars discussed in Section 2.1. The atypical $E(B-V)$ values arise from circumstellar material and are therefore not a part of the interstellar DIB environment. We acknowledge that they could influence our results and therefore we differentiate these points in the plots and analyses that follow.

It is insightful to illustrate our methodology by performing a PCA on a simple, two-dimensional data set; we therefore performed a PCA analysis using only the measurements for the variables $E(B-V)$ and N(H). Using the notation provided in Table 2, we denote these variables x₁₇ and x₂₀ and their corresponding standardized forms as z₁₇ and z₂₀. The 2 × 30 matrix of (standardized) measurements X from Equation (3) is then

$$\begin{eqnarray}&&X=\left[\begin{array}{c}{z}_{17}\\ {z}_{20}\end{array}\right]=\left[\begin{array}{cccc}0.273 & -\,0.106 & ... & -\,0.333\\ 0.121 & 0.121 & ... & -\,0.285\end{array}\right].\end{eqnarray} \tag{ 5 }$$

From this, the PCA results in the following 2 × 2 transformation matrix A,

$$\begin{eqnarray}A=\left[\begin{array}{cc}{a}_{\mathrm{1,1}} & {a}_{\mathrm{1,2}}\\ {a}_{\mathrm{2,1}} & {a}_{\mathrm{2,2}}\end{array}\right]=\left[\begin{array}{cc}0.707 & 0.707\\ 0.707 & -\,0.707\end{array}\right],\end{eqnarray} \tag{ 6 }$$

as well as the 2 × 30 matrix, Y, containing the transformed data points:

$$\begin{eqnarray}&&Y=\left[\begin{array}{c}{\mathrm{PC}}_{1}\\ {\mathrm{PC}}_{2}\end{array}\right]=\left[\begin{array}{cccc}0.279 & 0.011 & ... & -\,0.437\\ 0.107 & -\,0.161 & ... & -\,0.034\end{array}\right].\end{eqnarray} \tag{ 7 }$$

Note that we use the notation PC_i to describe the measurement values in the new coordinate frame described by the PCs. That is, PC_i = y_i in Equation (1). The results pertaining to the PCs for this example are summarized in Table 3. Column 1 lists the PC number, while Column 2 shows the eigenvalue associated with that PC. Since there are two original variables, the sum of the eigenvalues is two. The fraction of the total variation in the data for which each PC is responsible is shown in Column 3, followed by the cumulative fraction in Column 4. The (unit length) eigenvectors are shown in Column 5.

Figure 4 illustrates these results in two ways. In the top panel, we show the original data set along with the PC vectors (i.e., the eigenvectors corresponding to PC₁ and PC₂). From this figure, it is clear that the vector PC₁ lies in the direction of maximum variance. From Table 3, we see that in fact almost all of the variation in the data set (>90%) is accounted for by PC₁. This indicates that a single variable can adequately describe variations among $E(B-V)$ and N(H)—an interstellar finding that is well known (e.g., Bohlin et al. 1978, see also the discussion below).

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 4, we illustrate the results in the form of a so-called "biplot." Here, we show the transformed data points (i.e., those of matrix Y) and the projections of the original variables onto the PC₁–PC₂ plane. This format is much more convenient for interpreting the PCs. For instance, we can examine the projections of $E(B-V)$ and N(H) onto PC₁ (i.e., only the horizontal components of the vectors) to see how PC₁ influences each original variable. Both projections onto PC₁ are positive; thus, an increase of PC₁ implies an increase of both $E(B-V)$ and N(H). An obvious interpretation for PC₁ is therefore the amount of material in the line of sight: more material will result in larger $E(B-V)$ and N(H) values at the same time.

Doing the same for PC₂, we see that the vertical component of the $E(B-V)$ vector is positive, while that of N(H) is negative. PC₂ therefore represents some difference between $E(B-V)$ and N(H). A possible interpretation is therefore the dust-to-gas ratio: an increase in the dust-to-gas ratio can result in an increase in the $E(B-V)$ or a decrease in the N(H).

Using the eigenvectors, we can construct equations for PC₁ and PC₂. Since the PCs were constructed from the standardized variables, the resulting equations will be in terms of z₁₇ and z₂₀; using Equation (4) and the mean and standard deviations in Table 2, however, we can express them in terms of the original variables, $E(B-V)$ and N(H).

$$\begin{eqnarray}{\mathrm{PC}}_{1} & = & 0.707({z}_{17})+0.707({z}_{20})\\ & = & 5.35(E(B-V))+(5.93\times {10}^{-22})N({\rm{H}})-1.99,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}{\mathrm{PC}}_{2} & = & 0.707({z}_{17})-0.707({z}_{20})\\ & = & 5.35(E(B-V))-(5.93\times {10}^{-22})N({\rm{H}})-0.19.\end{eqnarray} \tag{ 9 }$$

We can also use our example to demonstrate how to reconstruct the observed data with reduced dimensionality. In our 2D example, PC₁ accounts for most of the variation in the data, and one could argue that the variations in PC₂ are mostly insignificant. In other words, we can attribute variations in PC₂ to noise in the data, and thus we can describe our data by just using one variable: PC₁. Mathematically, we thus eliminate the final column of matrix A, forming a 1 × 30 matrix, Y; in doing so, we leave Equation (8) as is and set Equation (9) equal to zero. Solving for N(H) then reveals

$$\begin{eqnarray}&&N({\rm{H}})=9.0\times {10}^{21}E(B-V)-0.3\times {10}^{21}\end{eqnarray} \tag{ 10 }$$

while a least-squares fit to the original data set yields

$$\begin{eqnarray}&&N({\rm{H}})=10.0\times {10}^{21}E(B-V)-0.50\times {10}^{21}.\end{eqnarray} \tag{ 11 }$$

Both results are comparable (see the blue and green lines in Figure 5), illustrating that a PCA can be used to reliably reconstruct relationships between variables that are present in the original data set. Note, however, that this is somewhat different from the well-known finding by Bohlin et al. (1978):

$$\begin{eqnarray}&&N({\rm{H}})=5.8\times {10}^{21}E(B-V).\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

This relation is also shown in Figure 5 (red line). Two important factors contribute to the difference between Bohlin et al. (1978) and our results. First, our sample of sightlines is biased to contain only single cloud lines of sight, with diverse physical conditions whereas Bohlin et al. (1978) used a large sample of much more reddened lines. Second, the Bohlin et al. (1978) relation was forced to go through the origin; the PCA and least-squares fit were not. However, a discussion of such differences is not relevant for our work.

This example shows the steps involved in carrying out a PCA and using its results. For what follows, we will be using multi-dimensional data; the only complication is then that, for the discussion and interpretation, we need to work with projections onto a 2D space for proper visualization.

Biplots for PC₁ and PC₂ are shown in Figure 7; the top panel shows the full PCA results while the other panels show the results of the additional PCA's using only line-of-sight parameters and using only DIBs. As was the case for the 2D example (see Section 3.3), the biplot shows the measurements for all lines of sight transformed into the PC reference frame, and only shows the values for the first two PCs—i.e., the two parameters that together carry 80% of the variance in the data set. The biplot also shows the projection of the original variables onto the plane defined by these two PCs. Such projections are useful to find trends and correlations: the larger the projection of an original variable is onto an axis, the better the original variable correlates with the PC. Furthermore, original variables whose projections onto the PC plane are similar, must necessarily represent tight correlations as well. Note as well the rectangular areas at the vector heads (best seen in the close-up in the top panel to the right) that represent the uncertainties on the projections as derived from an MC simulation.

Since PC₁ and PC₂ capture most of the variation in the data, the PC1–PC2 biplot is particularly insightful for identifying trends and correlations. For example, W(λ6614) and W(λ6196) align almost perfectly with one another within their respective uncertainty ranges in PC₁. This implies that the two DIBs must correlate very tightly—a well-known fact for these DIBs (see, e.g., McCall et al. 2010). Similarly, W(λ5780) is known to correlate well with W(λ6284), W(λ6270), and W(λ4428) (Krelowski & Walker 1987; Herbig 1993) and these appear to be grouped together on the PC₁–PC₂ plane. Many more known relations between DIB strengths (see, e.g., Lan et al. 2015, for an overview) can be recovered from these plots.

To infer clues about what PC₁ physically represents, we can look at the projections of the original variables onto the PC₁ axis. These projections can be easily determined from Equations (24) and (25). For example, the vector that corresponds to the projection of the λ5780 DIB (i.e., z₈, see Table 2) onto the PC1–PC2 plane is given by the coefficients in front of z₈ in Equations (24) and (25), i.e., the corresponding vector is (0.207, 0.306). This vector is shown (albeit scaled up for visualization) in Figure 7. As PC₁ increases, so too does each of the original variables. Therefore, we interpret PC₁ as measuring the amount of material in the line of sight; as the amount of material increases, so does the strength of each DIB along with variables such as color excess and the various column densities. Intuitively, this makes sense: The factor that most strongly contributes to the column density of any species in a given line of sight is the amount of material.

We can elaborate further by examining the other two PCA cases. The bottom right panel of Figure 7 shows the PCA results when only DIBs are considered. Again, each variable increases in the same direction, so PC₁ must once again be some proxy for the amount of material in the line of sight. However, since we only consider DIBs in this case, PC₁ is specifically tracing the amount of material that contributes to DIBs. If we compare these results to those in the top panel, we find the same ordering of DIBs—i.e., W(λ5797) still has the strongest PC₁ projection and W(λ5769) still has the weakest. This confirms that PC₁ traces the amount of DIB-producing material in the top panel as well. This conclusion is further established by the bottom left panel, illustrating the PCA results when only the line-of-sight parameters are considered. Here, PC₁ appears to be quite different: both N(H₂) and $E(B-V)$ have stronger PC₁ projections than N(H), whereas in the top panel, the N(H) and $E(B-V)$ projections on PC₁ are comparable, while that of N(H₂) is much shorter. This means that PC₁ is tracing something slightly different in each case. In the former, it seems to trace something more closely related to the amount of dust (given the large $E(B-V)$ projection), whereas in the latter, the amount of gas plays a stronger role. Since the DIB carriers are widely accepted as begin gaseous species, it seems reasonable that N(H) would more strongly influence the strength of DIBs than $E(B-V)$ as illustrated in our full PCA results, so this too is consistent with our interpretation that PC₁ traces the amount of DIB-producing material. The fact that these cases are so different while the just-DIB results are nearly identical to those obtained using the full data set implies that our results are driven by the DIBs, rather than the line-of-sight parameters.

Interestingly, the largest projection onto PC₁ comes from the λ5797 DIB, meaning that PC₁ has the strongest correlation with W(λ5797) (r = 0.957—see Figure 8). This result also carries over to the just-DIB case (r = 0.963—not shown). Hence, the best measure for the amount of DIB-producing material in the line of sight is W(λ5797) and not N(H) or $E(B-V)$. Therefore, W(λ5797) is conceivably a more appropriate normalization factor for DIB strengths than $E(B-V)$. This is discussed further in Section 5.

Zoom In Zoom Out Reset image size

After establishing that PC₁ traces the amount of DIB-producing material along the line of sight, we can define a new parameter, N_DIB, describing this value. The transformed data points derived from PCA are, by definition, centered at zero and therefore approximately half of our sightlines have negative PC₁ values. For N_DIB to be physically meaningful, we need to have only positive values. Since our PCs were calculated from standardized variables, we have a way of adjusting our values. Using Equation (4), we subtracted the mean divided by the standard deviation from each variable prior to performing PCA. To correct for this, we must add these values to each term in our PC₁ equation (see Equation (24)). Thus, we calculate a constant, C₁, from the weighted fractions of $\bar{{x}_{i}}/{{s}_{x}}_{i}$ for each variable, x_i. Explicitly,

$$\begin{eqnarray}&&{C}_{1}=\displaystyle \sum _{i}\,\displaystyle \frac{{a}_{i,1}(\bar{{x}_{i}})}{{s}_{{x}_{i}}}\end{eqnarray} \tag{ 13 }$$

and thus, PC₁ can be expressed in its uncentered form:

$$\begin{eqnarray}{\mathrm{PC}}_{1,\mathrm{uncentered}} & = & {\mathrm{PC}}_{1}+{C}_{1}\\ & = & {\mathrm{PC}}_{1}+6.616.\end{eqnarray} \tag{ 14 }$$

It is also worth noting that, since the PCs define directions in space, the scaling is arbitrary. We can use the tight correlation with W(λ5797) to rescale PC₁ such that as the strength of λ5797 doubles, so too does N_DIB. The result is Equation (15).

$$\begin{eqnarray}{N}_{\mathrm{DIB}} & = & 0.136({\mathrm{PC}}_{1,\mathrm{uncentered}})\\ & = & 0.136({\mathrm{PC}}_{1})+0.900.\end{eqnarray} \tag{ 15 }$$

Since PC₁ is not directly measurable, we provide an alternate expression in terms of W(λ5797).

$$\begin{eqnarray}&&{N}_{\mathrm{DIB}}\approx 0.0185W(\lambda 5797)+0.19.\end{eqnarray} \tag{ 16 }$$

The constant in Equation (16) indicates that some DIB carriers are present before the λ5797 carriers start to form. This result is further discussed in Section 5.

PC₂ represents the second largest source of variations in the DIB strengths and line-of-sight parameters in our sample, and must be a parameter that is uncorrelated with PC₁. The largest contribution to PC₂ comes from W(λ5797)/W(λ5780) (see Figure 7 and Equation (25)), which correlates strongly (but negatively) with PC₂. The second largest projection onto PC₂ is W(λ6284), in the positive direction. Note that some of our original parameters are hardly influenced by PC₂, most notably W(λ5797), W(λ6379), N(H), and $E(B-V)$.

These correlations (and lack thereof) are valuable clues to infer the physical quantity that is represented by PC₂. We start by noting that W(λ5797)/W(λ5780) is a variable that is clearly related to physical conditions, and has, in particular, been linked to the amount of exposure to the ambient UV radiation field (see, e.g., Herbig 1995; Krełowski & Sneden 1995; Cami et al. 1997; Sonnentrucker et al. 1997; Krełowski et al. 1999; Vuong & Foing 2000; Cox et al. 2006; Friedman et al. 2011; Vos et al. 2011). More specifically, these studies have shown that the λ5780 DIB is weak in shielded environments, but becomes significantly stronger with increasing UV exposure in diffuse clouds; the λ5797 DIB, on the other hand, is already strong in shielded environments, and is largely indifferent to UV exposure in most diffuse clouds. In the harshest environments (e.g., Orion), both DIBs are very weak or absent. This different dependence on the radiation field then results in the more sheltered ζ-type clouds having a high W(λ5797)/W(λ5780) ratio and the more exposed σ-type cloud environments having a low ratio. This dependence on the radiation field is also responsible for the so-called "skin effect"—the fact that DIB carriers seem to be concentrated near the surfaces of clouds (see, e.g., Adamson et al. 1991; Herbig 1995; Cami et al. 1997; Vos et al. 2011).

One way such behavior can be understood is by considering that the DIB carriers are molecules in a specific ionization state, and the changing DIB strengths reflect the molecules undergoing ionization with the carrier of the λ5797 having a lower ionization potential than the λ5780 DIB carrier (Cami et al. 1997; Sonnentrucker et al. 1997). Note that the charge balance is really determined by the interplay between ionization and recombination, and thus not only depends on the radiation field, but also on the density. Similarly, the DIB behavior can be explained by the DIB carriers being a specific hydrogenation or protonation state (Vuong & Foing 2000; Le Page et al. 2001, 2003); also in this case, the density plays a role.

Thus, the W(λ5797)/W(λ5780) ratio certainly traces cloud depth, which encompasses UV exposure, density, and possibly also temperature. Given the strong correlation with PC₂, PC₂ could thus trace directly the strength of the ambient radiation field G₀,⁵  or alternatively it could be a measure for G₀/n_e and reflect the balance between ionization and recombination (or similarly between hydrogenation and de-hydrogenation).

We dismiss this second interpretation for two reasons. First, Gnaciński et al. (2007) showed that several of the DIBs included in our sample are independent of n_e, contrary to the large variations we see among their PC₂ values. Furthermore, Kos & Zwitter (2013) found no differences in n_e values between σ and ζ sightlines, despite W(λ5797)/W(λ5780) having the largest projection on PC₂. Thus, we conclude that PC₂ must be tracing changes in G₀ only, though of course this in turn can be determined by the cloud depth.

If we dichotomize sightlines into σ and ζ clouds, this becomes more clear. Referring again to Figure 7, the negative projections onto PC₂ are all indicators of a sheltered ζ environment: we see a strong W(λ5797)/W(λ5780), hydrogen in molecular form (i.e., large f(H₂)), and the existence of C₂ (demonstrated by the four C₂-DIBs—λ5769, λ5546, λ5513, and λ4964; Thorburn et al. 2003). In contrast, positive PC₂-projections are consistent with more exposed conditions: we observe hydrogen in neutral form, as well as a strong projection from λλ5780, 4428, and 6284—DIBs known to correlate well with H i (Herbig 1993). Finally, both N(H) and $E(B-V)$ contribute almost nothing to PC₂. Again, this is consistent with our interpretation since both variables contribute to the amount of material but are generally unaffected by changes in G₀. Note that λ5797 and λ6379 exhibit the same behavior, as was found previously (Cami et al. 1997; Kos & Zwitter 2013; Lan et al. 2015). The case of λ5797 is especially interesting. Although we use W(λ5797)/W(λ5780) to differentiate between regions of high and low UV exposure, these changes almost exclusively reflect the λ5780 carrier. Thus, while the carrier of λ5780 is sensitive to changes in G₀, λ5797 survives at a steady strength per unit reddening over a wide range of conditions. This indifference to physical conditions is essentially why this DIB is a good tracer for the amount of DIB carrier material.

Comparing the three separate panels of Figure 7, we once again see that the PCA results from the line-of-sight parameters (bottom left) are quite different. N(H), in particular, has a strong, positive projection on PC₂, making it clear that PC₂ traces something other than G₀ in this case. This is not surprising since PC₂ necessarily depends on PC₁, and PC₁ was different between the full PCA and the PCA for the line-of-sight parameters only. The PCA results using only the DIBs, however, are nearly identical to those obtained with all 23 variables. This suggests that changes in PC₂—that is to say, changes in G₀—are directly observable through the relative strengths of DIBs.

For the sake of clarity, we further formalize our association of PC₂ with the radiation field by defining a new parameter G_DIB. Ideally, we would like to scale this parameter to represent actual G₀ values; however, without knowledge of specific G₀ values for our targets, we cannot do this. Therefore, we leave PC₂ as it is, and simply assign the name G_DIB:

$$\begin{eqnarray}&&{G}_{\mathrm{DIB}}={\mathrm{PC}}_{2}.\end{eqnarray} \tag{ 17 }$$

Since W(λ5797)/W(λ5780) correlates most strongly with PC₂, we will approximate G_DIB as a function of W(λ5797)/W(λ5780). Figure 9 shows G_DIB as a function of the W(λ5797) to W(λ5780) strength ratio. Fitting a straight line to the data gives the following relation with r = −0.716 (shown in green in Figure 9):

$$\begin{eqnarray}&&{G}_{\mathrm{DIB}}\approx -9.07[{\rm{W}}(\lambda 5797)/{\rm{W}}(\lambda 5780)]+2.55.\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

Arguably, a linear fit is not appropriate since there are three sightlines that strongly deviate from this relation—HD 15137, HD 198478, and HD 209975. However, the interstellar features for these sightlines are quite broad, suggesting that there may be more than one cloud in the line of sight that we are unable to resolve. If this is the case, then we cannot expect them to follow the same trend as true single clouds. Ignoring these three sightlines in our linear fit, we obtain (with r = −0.832; shown in blue in Figure 9):

$$\begin{eqnarray}&&{G}_{\mathrm{DIB}}=-9.98[{\rm{W}}(\lambda 5797)/{\rm{W}}(\lambda 5780)]+2.67.\end{eqnarray} \tag{ 19 }$$

Alternatively, these sightlines are in fact single clouds and represent a true trend in the data. With this interpretation, the linear dependence breaks down for positive G_DIB values (i.e., σ-type clouds) and G_DIB increases independently of W(λ5797)/W(λ5780). As previously discussed, W(λ5797) is independent of PC₂, while W(λ5780) shows a strong dependence; hence, this trend is most easily explained as an ionization effect of the λ5780 carrier. We propose that carrier of λ5780 is a cation or a dehydrogenated species (see, e.g., Cami et al. 1997; Sonnentrucker et al. 1997). As G₀ increases, more of these molecules become ionized/dehydrogenated so W(λ5780) becomes stronger. Finally, when a certain G₀ is reached, we enter the σ-cloud regime, where the cation/dehydrogenated state dominates. Hence, W(λ5780) is no longer dependent on the environment, and λ5797 and λ5780 exist in approximately equal proportions. We can model this using a piece-wise function. For σ clouds, W(λ5797)/W(λ5780) has a mean value of 0.198 ± 0.064, with a maximum W(λ5797)/W(λ5780) value around 0.3. For ζ clouds, the linear relation holds. We fit a line to the negative G_DIB values, forcing the intercept through the σ cloud mean using the MPFITEXY⁶ routine (Williams et al. 2010) based on the MPFIT package (Markwardt 2009). We obtain a correlation coefficient of −0.75 for ζ clouds. The resulting equation is

$$\begin{eqnarray}&&{G}_{\mathrm{DIB}}\\ \quad \approx \,\left\{\begin{array}{ll}\mathrm{undefined}, & \mathrm{if}\ \tfrac{{\rm{W}}(\lambda 5797)}{{\rm{W}}(\lambda 5780)}\lesssim 0.3\\ -8.44\left(\tfrac{{\rm{W}}(\lambda 5797)}{{\rm{W}}(\lambda 5780)}\right)+1.94, & \mathrm{if}\ \tfrac{{\rm{W}}(\lambda 5797)}{{\rm{W}}(\lambda 5780)}\gtrsim 0.3.\end{array}\right.\end{eqnarray} \tag{ 20 }$$

From this, it is clear that W(λ5797)/W(λ5780) is a suitable measure for G₀ within ζ clouds, whereas for σ clouds, it may or may not be appropriate.

Figure 10 shows the PC₁–PC₃ and PC₂–PC₃ biplots; note that we no longer consider the PCA results from just DIBs because the uncertainties are too large to properly interpret the results. Instead, we focus only on the full PCA results.

Zoom In Zoom Out Reset image size

Because PC₃ is necessarily uncorrelated with PC₁ and PC₂, it must trace something that is fairly independent of both the amount of material and the ambient UV field. Some clues about PC₃ are revealed in Figure 10. For negative PC₃ values, we see strong projections from f(H₂) and N(H₂). For positive values, we see strong projections from N(H i) and the C₂-DIBs. We may therefore interpret PC₃ as tracing some difference between the amount of H₂ and the amount of both C₂ and H i. van Dishoeck & Black (1982) showed that C₂ tends to form in cooler, denser regions than H₂ so one could interpret PC₃ to trace the temperature and density. Note that in that interpretation, colder and denser environments correspond to higher values for PC₃. We test this hypothesis by considering the projections on PC₃. Although the vector for F_⋆ has large uncertainties, the projection is clearly negative. This would then imply that more depletion occurs in warmer, more diffuse regions. This result seems counterintuitive, and it is therefore unlikely that PC₃ traces temperature and density.

A careful analysis of the biplots reveals more clues about PC₃. There is one clear outlier, HD 147933, which has a PC₃ value of 4.42. PC₃ therefore traces some characteristic that sets HD 147933 apart from the other lines of sight included in this study. This target, found within the ρ Ophiuchi dark cloud, has been well-studied. A notable characteristic is its unusually large dust grains (Carrasco et al. 1973), resulting in a ratio of hydrogen to color excess 2.7 times greater than the mean interstellar value (Bohlin et al. 1978). Two plausible explanations are therefore grain size or the gas-to-dust ratio.

We first consider the interpretation that PC₃ traces grain size, with positive PC₃ projections corresponding to larger grains and negative projections corresponding to smaller grains. Since small dust grains dominate the total dust grain surface area, and the only effective mechanism for producing interstellar H₂ is on the surface of dust grains (Gould & Salpeter 1963), it makes sense that the strongest negative PC₃ projections come from f(H₂) and N(H₂). Moreover, $E(B-V)$ has a negative projection since small grains dominate reddening as well. Finally, F_⋆ has a negative PC₃ projection since a larger surface area of dust grains implies that there is more surface onto which elements may deplete. In terms of positive projections, we see a large projection from N(H i). This too is consistent since, with larger dust grains and a reduced surface area, H i cannot effectively be converted into H₂.

To further test this idea, we plot the total-to-selective extinction, R_v, as a function of PC₃ in Figure 11. R_v tends to increase with grain size, so if our interpretation is correct, we would expect to see R_v increasing with PC₃. Although a weak trend is apparent, HD 110432 (shown as a red triangle) is hard to reconcile with the rest of the data. Part of this discrepancy may be due to the fact that HD 110432 is a Be star. Rachford et al. (2001) corrected for the CS contribution and found that an R_v value of 3.3 is more appropriate, though this still implies an above-average grain size.

Zoom In Zoom Out Reset image size

Alternatively, PC₃ may trace the gas-to-dust ratio. With this interpretation, positive projections indicate larger gas-to-dust ratios while negative projections indicate larger dust-to-gas ratios. Both f(H₂) and N(H₂) have large, negative projections, which again is consistent with a larger dust grain surface area and therefore a higher dust-to-gas ratio. $E(B-V)$ also increases in this direction, with similar reasoning. Finally, we see a large projection from F_⋆. Noting that the large grains toward HD 147933 come from grain coagulation rather than elemental depletions (Jura 1980), this too is consistent: if there is more dust, then there is more material onto which gas may deplete. Both N(H) and N(H i) have positive projections, consistent with higher ratios of gas. N(H i), in particular, has a large projection since a large quantity of dust is required to transform H i into H₂.

Figure 12 shows the ratio of hydrogen to color excess, N(H)/$E(B-V)$, plotted against PC₃. The general trend is apparent, and can be represented by (using IDL's FITEXY routine):

$$\begin{eqnarray}&&\displaystyle \frac{N({\rm{H}})\times {10}^{-21}}{E(B-V)}=(7.29\pm 0.23)+{\mathrm{PC}}_{3}\times (1.63\pm 0.26).\end{eqnarray} \tag{ 21 }$$

However, the trendline appears to be influenced by the two outliers—HD 147933 and HD 207198—and we therfore performed another fit excluding these two data points; this yields similar results:

$$\begin{eqnarray}&&\displaystyle \frac{N({\rm{H}})\times {10}^{-21}}{E(B-V)}=(7.31\pm 0.27)+{\mathrm{PC}}_{3}\times (1.89\pm 0.26).\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

There is quite a bit of scatter around this trend though. This is somewhat expected since N(H)/$E(B-V)$ does not measure the gas-to-dust ratio directly, as $E(B-V)$ measures the degree of reddening but not the dust mass. Moreover, our values for N(H) do not include any hydrogen that may exist in ionized form, while both HD 36822 and HD 36861 are known to lie near H ii regions. Finally, $E(B-V)$ values tend to be overestimated for Be stars (see Section 2.1), resulting in smaller-than-expected N(H)/$E(B-V)$ values. Be stars are differentiated in Figure 12 by square data points.

Although we cannot conclusively state what precisely PC₃ is tracing, it seems clear that it is tracing a property of the dust and its relation with the gas. We thus favor the gas-to-dust ratio interpretation here.

Once again, we define a new parameter, GTD, describing the observed changes along PC₃. Bohlin et al. (1978) showed that the mean interstellar value of N(H)/$E(B-V)$ is 5.8 × 10²¹ atoms cm⁻² mag⁻¹. We rescale PC₃ such that GTD gives the ratio of gas to dust with respect to the mean interstellar value: A GTD value of 1 corresponds to the mean interstellar value of N(H)/$E(B-V)$, and HD 147933 has a GTD value of 2.7. GTD is therefore expressed as

$$\begin{eqnarray}&&\mathrm{GTD}=0.318({\mathrm{PC}}_{3})+1.29.\end{eqnarray} \tag{ 23 }$$

Figure 13 shows the PC₄ biplots with PC₁ (top), PC₂ (middle), and PC₃ (bottom). It is difficult to physically interpret PC₄ due to the large uncertainties on the projections; however, it is clear that the strongest projections come from F_⋆ and W(λ5797)/W(λ5780), and thus, PC₄ traces some difference between these two variables. In the previous PCs, F_⋆, and W(λ5797)/W(λ5780) increased in the same direction, suggesting that depletions generally increase with cloud depth. In terms of PC₄, however, these two variables increase in opposite directions. Consequently, a possible interpretation is that PC₄ traces the deviations from the mean relation between these two variables, though the physical cause for these deviations is unknown. Figure 14 shows the ratio of these variables plotted as a function of PC₄.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size