Spatial structure of several diffuse interstellar band carriers
MMNRAS , 1–18 (2017) Preprint 4 July 2018 Compiled using MNRAS L A TEX style file v3.0
Spatial structure of several diffuse interstellar band carriers
Janez Kos, Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW 2006, Australia
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Diffuse interstellar bands (DIBs) hold a lot of information about the state and thestructure of the ISM. Structure can most directly be observed by extensive spectro-scopic surveys, including surveys of stars where DIBs are especially important, asthey are conveniently found in all observed bands. Large surveys lack the quality ofspectra to detect weak DIBs, so many spectra from small regions on the sky have tobe combined before a sufficient signal-to-noise ratio (SNR) is achieved. However, theclumpiness of the DIB clouds is unknown, which poses a problem, as the measuredproperties can end up being averaged over a too large area. We use a technique calledGaussian processes to accurately measure profiles of interstellar absorption lines in 145high SNR and high resolution spectra of hot stars. Together with Bayesian MCMCapproach we also get reliable estimates of the uncertainties. We derive scales at whichcolumn densities of 18 DIBs, CH, CH + , Ca I, and Ca II show some spatial correlation.This correlation scale is associated with the size of the ISM clouds. Scales expressed asthe angle on the sky vary significantly from DIB to DIB between ∼ . ◦ for the DIB at5512 ˚A and 3.5 ◦ for the DIB at 6196 ˚A, suggesting that different DIB carriers have dif-ferent clumpiness but occupy the same general space. Our study includes lines-of-sightall over the northern Milky Way, as well as out of the Galactic plane, covering regionswith different physical conditions. The derived correlation scales therefore represent ageneral image of the Galactic ISM on the scales of ∼ pc to pc. Key words: methods: observational – surveys – stars: early types – ISM: abundances– ISM: clouds – ISM: dust, extinction – ISM: molecules – ISM: structure – Galaxy:structure
Diffuse interstellar bands (DIBs) are yet to be understoodabsorption features observed in VIS and NIR spectra of red-dened stars. There are hundreds of known DIBs all shar-ing some common properties: they are resolved (narrow-est ones have a full width at half maximum (FWHM) ofaround 0.5 ˚A), they all correlate at least vaguely with theamount of interstellar extinction in the line-of-sight (LOS),and almost none have a confirmed carrier or origin. DIBswere discovered in the early 1920s (Heger 1922) and theirinterstellar origin has been confirmed soon afterwards (Mer-rill 1936). Despite being among the first molecules observedin the ISM, most of the carriers are still not identified. Abright exception are a few DIBs in the NIR, for which theC + has been confirmed as a carrier by laboratory stud-ies (Campbell et al. 2015; Walker et al. 2015). None ofthe DIBs addressed in this study have confirmed carriers.Among the proposed carriers are polycyclic aromatic hydro-carbons (L´eger & d’Hendecourt 1985; van der Zwet & Alla-mandola 1985; Salama et al. 1999), fullerenes (Kroto 1988;Campbell et al. 2015), and other hydrocarbons (Motylewski et al. 2000; Maier et al. 2011; Kre(cid:32)lowski et al. 2010) from sim-ple molecules with only 4 atoms to long carbon chains. Seereviews Krelowski (1989), Herbig (1995), and Sarre (2006)for a thorough overview of the field.It has been shown (Kre(cid:32)lowski 2014; Baron et al. 2015)that there must be more than just a few molecules producingDIBs. Most of the strongest DIBs are believed to each belongto a different carrier. Correlation between DIB strengths andcorrelation with the interstellar reddening is most commonlyused to identify DIB families and relations between them(Kre(cid:32)lowski & Walker 1987; Westerlund & Kre(cid:32)lowski 1989;Cami et al. 1997; Kos & Zwitter 2013). Correlation is notthe only relation between DIBs we can observe. Less com-mon are studies of radial velocities and spatial correlations.Most direct approach to this are measurements of three-dimensional structure in large spectroscopic surveys (Koset al. 2014; Zasowski et al. 2015). A huge number of ob-served LOS is the minimal requirement for such approach.Despite many observed LOS the spatial resolution of suchmaps has been poor.Studies of DIBs with a fine ( <
100 pc) spatial resolutionare rare, because in order to observe differential strengths of © a r X i v : . [ a s t r o - ph . GA ] M a r J. Kos
DIBs high quality spectra are needed. Large spectroscopicsurveys of stars provide the needed spatial resolution, butrarely include enough spectra with high SNR. Therefore ded-icated observations must be made like Vos et al. (2011), vanLoon et al. (2009), or van Loon et al. (2013). All three studiesshow a highly structured DIBs-bearing ISM in the Galacticplane Vos et al. (2011) and well away from it (van Loonet al. 2009, 2013). Only a few strongest DIBs were studiedin mentioned papers in small regions of the Galaxy (Up-per Scorpius, ω Cen and the Tarantula nebula). The mostdetailed study of small scale variations is (Cordiner et al.2013). They studied DIBs and some other ISM lines in veryhigh SNR spectra of 4 stars in ρ Oph and found variationsof ∼
3% (expressed the same way as variations between LOSin this paper) in LOS as close as 2.8 arc seconds. It is hardto generalize properties of the ISM observed in these regionsover the whole Galaxy. We aim to change this by observingstars over the whole Galactic plane. Poorly sampled but highquality data can be used to statistically constrain the natureof the ISM structure. Instead of producing localized maps,we study small scale variations of DIB strengths throughpairwise correlations in tightly separated LOS.We can hypothesize three different cases for the struc-ture of the DIB bearing ISM: clouds with different clumpi-ness for different DIB species, separate clouds of differentDIB species, or gradients between different species in a sin-gle cloud. From the pairwise correlation between DIB pro-files in tight pairs of LOS we can distinguish between thethree cases. We can put some constraints on respective carri-ers from observed differences in spatial distribution betweendifferent DIB species.We acquired high SNR and medium-high resolutionechelle spectra of 145 hot stars covering a wavelength rangebetween 3700 and 7300 ˚A. Target selection and reductionof the dataset is described in Section 2. Only strong andprominent DIBs are studied in this work. To describe theirprofiles we fit them with asymmetric Gaussians, empower-ing the Gaussian processes and a Bayesian MCMC fittingscheme to calculate uncertainties and deal with stellar con-tinuum and statistical noise. The method is described inSection 3. The core of this paper is the study of close pairsof LOS, the differences between them and how they revealthe structure of the ISM and physics of the DIBs. This ispresented in Section 4. A discussion of the results and theimplications of this study are presented in Section 5.
Spectra were obtained with the 1.82 m telescope of the Ob-servatory of Padua in Asiago, Italy and an echelle-type spec-trograph. Observations were done in 17 observing sessionsbetween January 16 2011 and October 15 2016. The spec-trograph covers wavelengths between 3700 ˚A to 7300 ˚A in34 echelle orders without any gaps. The resolving power ofthe setup used is between 20,000 and 23,000. The variationin the resolving power is due to seeing constraints; the slitwidth was set to 200 µ m in good seeing and to 250 µ m whenthe seeing was worse. Some of the data was also taken inthe shared time with other projects that required the wider slit. Enough signal was collected for each star to have thesignal-to-noise ratio (SNR) of 300 or more per pixel over thewhole wavelength range where DIBs included in this studyare present. A required exposure time to accumulate the re-quired SNR is between 1 minute for brightest stars and 1.5hours for darkest stars in this study, assuming a typical see-ing of ∼ A selection function for program stars was kept simple: amagnitude limit of V=8.5 (darker stars would require in-feasible exposure times) and a color cut at spectral typeB3 was used. Cooler stars have a spectrum populated bymore stellar absorption lines and are unsuitable for the pre-sented study with the approach described further in the text.The color cut was chosen in a way that enough candidatespassed the cut. We showed that the color cut is adequatein Kos & Zwitter (2013). Further candidates were discardedfor not having reliable reddening measurements in the lit-erature. Around 300 candidates remained and we observed145 of them. The color excess of observed stars varies be-tween just measurable level of few hundreds of magnitudeto around E(B-V)=1.2 magnitude. Stars with higher colorexcess would be too faint at the blue end of the wavelengthrange and would not fit the accessible magnitude range. Wetried to keep the distribution of the color excess of the ob-served stars as uniform as possible in order to sample morediverse regions of the Galaxy.Four of the observed stars are colder than type B3 stars.They are all rapid rotators (G(cid:32)l¸ebocki & Gnaci´nski 2005),so their spectrum is appropriate for our analysis becausethe rotationally broadened lines satisfy the same criteria astemperature broadened lines of hotter stars.Distances of targeted stars are taken from Neckel et al.(1980) (photometric distances) and Lindegren, L. et al.(2016) (trigonometric distances). Photometric distance isavailable for almost every star and the trigonometric dis-tance is available for only 60% of the stars. Figure 2 showsthe comparison between the two distances. For most starsboth distances match within the errorbars. There are how-ever a few stars where two distances are several sigmas apart.In the analysis we will use the weighted averages between thephotometric and trigonometric distances wherever it is pos-sible. The exception are the stars that have a negative meanparallax.Table 1 lists all relevant information about the observedstars.
Spectra were reduced by standard IRAF routines within echelle and noao packages. Images were bias corrected and
MNRAS000
MNRAS000 , 1–18 (2017) patial structure of several DIB carriers Figure 1.
Map of the Galaxy with observed stars. Red symbols show LOS where we detect double Na I and Ca II profiles with twodistinct radial velocity components. Rest of the LOS are marked green. More densely populated OB associations are shown in zoomed-ininsets. Red line is the celestial equator. The background picture is the Planck dust map (Planck Collaboration et al. 2015). corrected for cosmic rays. We remove the cosmic rays fromthe raw images, because most often we only recorded twospectra, so the cosmic rays can not be removed by a median,for example. Spectral traces were then optimally extractedfrom flat field image and were flat field corrected. Scatteredlight has been removed in the process. Wavelength calibra-tion has been done against spectra of a Thorium-Argon arclamp taken with every set of science images. Echelle spectrawere normalized and echelle orders were stitched togetherinto a single spectrum. Two or more exposures were madefor each star, depending on the required SNR and weatherconditions. Single spectra were shifted into the barycentricvelocity frame and combined into one spectrum per star.An error spectrum tracing the uncertainties for eachpixel was created along with the reduced spectrum.
Different DIBs have a diverse range of profiles: a few aresymmetrical single lines with no resolved substructure (Sarreet al. 1995), most are bands with a resolved substructure(Galazutdinov et al. 2002), and sometimes they are blendedindividual DIBs (Jenniskens & Desert 1994). Data used inthis work have a resolving power just above R=20,000, whichis not enough to resolve any substructure. The substructureor blended DIBs, however, can distort the profiles into asym-metric shapes. This is detectable in our spectra. We thereforedescribe the DIB profile with a Gaussian: Φ ( λ ; A , λ , σ ) = A √ πσ exp (cid:18) −( λ − λ ) σ (cid:19) (1) where the possible asymmetry ( as y m ) of the profile is intro-duced by making the width σ a function of the wavelength: σ ( λ ; λ , as y m ) = σ + exp ( as y m ( λ − λ )) (2)Such implementation of the asymmetry is just one of thepossibilities and since there is no theoretical prescription forthe asymmetry, our choice is the one that works best in thefitting scheme. See Kos & Zwitter (2013) for a comparisonof two more asymmetry implementations.All DIBs studied here are optically thin, so their equiv-alent widths (EWs) are linearly proportional to the columndensities of the absorbing material. Column densities, how-ever, cannot be calculated, because the oscillator strengthsof DIBs are unknown. The opposite holds for atomic andmolecular absorption lines studied in this work. They canbe optically thick, so the column density must be used todescribe the amount of the ISM in the LOS. DIB strengthsare therefore expressed as equivalent widths and atomic andmolecular abundances as column densities. The same schemais used for fitting DIBs, molecular, and atomic absorptionprofiles. We will describe the absorption lines with an ab-sorption line shape I ( λ ) : I ( λ ) = exp (− τ ( λ )) , (3)where τ ( λ ) = N π e m e c f lu Φ ( λ ) , (4)where N is the column density of the absorber, and f lu is theoscillator strength of the observed transition. The strength MNRAS , 1–18 (2017)
J. Kos
Name l b
SpectralType
V E ( B − V ) PhotometricDistance n Parallax g Observingsession Peak S/N ( ◦ ) ( ◦ ) (mag) (mag) (kpc) (mas)HD108 117.9268 1.2498 O4 7.40 0.49 a ± a a , v s a ± a ± a ± m a ± a ± n , v , a ± v ± a ± a ± a a ± a , v ± a a ± a ± a ± a ± s ? 1.15 ± m v , l m ± l a ± a , v a , v a , v ± m ± v , s ± s , a ± m ? ? 1 430HD25940 153.6543 -3.0450 B3V 4.00 0.173 s a a , s , v v , a ± a ± n , s n , s ± v , l , a ± a ± a ± n a a s v , a s a ± v , a a ± a ± a , v a a a , v Table 1.
Basic properties of observed stars. Observing sessions: 1 (1. 16. 2011), 2 (9. 15. 2011 – 9. 16. 2011), 3 (2. 3. 2012 – 2. 6. 2012),4 (5. 3. 2012 – 5. 5.2012), 5 (8. 29. 2012 – 9. 5. 2012), 6 (24. 11. 2012 – 26. 11. 2012), 7 (24. 4. 2013 – 25. 4. 2013), 8 (24. 6. 2013 – 26.6. 2013), 9 (21. 9. 2013 – 23. 9. 2013), 10 (18, 10, 2013 – 18. 10. 2013), 11 (12. 1. 2014 – 13. 1. 2014), 12 (19. 3. 2014 – 21. 3. 2014),13 (10. 4. 2014 – 13. 4. 2014), 14 (3. 1. 2015 – 5. 1. 2015), 15 (31. 1. 2015 – 3. 2. 2015), 16 (1. 4. 2015 – 3. 4. 2015), 17 (15. 10. 2016).References: v: Valencic et al. (2004), s: Snow et al. (1977), a: Savage et al. (1985), m: Maeder (1975), l: S(cid:32)lyk et al. (2006), n: Neckel et al.(1980), f: Fitzpatrick & Massa (2007), g: Lindegren, L. et al. (2016). Coordinates, magnitudes and spectral types are taken from Simbad.Uncertainties of photometric distances are not given explicitly for each star but are estimated to be between 16 and 25%(Neckel et al.1980). Typical uncertainty for E ( B − V ) is 0.03 mag. Peak S/N gives the signal-to-noise ratio in the region with highest photon counts,usually around 5800 ˚A MNRAS000
Basic properties of observed stars. Observing sessions: 1 (1. 16. 2011), 2 (9. 15. 2011 – 9. 16. 2011), 3 (2. 3. 2012 – 2. 6. 2012),4 (5. 3. 2012 – 5. 5.2012), 5 (8. 29. 2012 – 9. 5. 2012), 6 (24. 11. 2012 – 26. 11. 2012), 7 (24. 4. 2013 – 25. 4. 2013), 8 (24. 6. 2013 – 26.6. 2013), 9 (21. 9. 2013 – 23. 9. 2013), 10 (18, 10, 2013 – 18. 10. 2013), 11 (12. 1. 2014 – 13. 1. 2014), 12 (19. 3. 2014 – 21. 3. 2014),13 (10. 4. 2014 – 13. 4. 2014), 14 (3. 1. 2015 – 5. 1. 2015), 15 (31. 1. 2015 – 3. 2. 2015), 16 (1. 4. 2015 – 3. 4. 2015), 17 (15. 10. 2016).References: v: Valencic et al. (2004), s: Snow et al. (1977), a: Savage et al. (1985), m: Maeder (1975), l: S(cid:32)lyk et al. (2006), n: Neckel et al.(1980), f: Fitzpatrick & Massa (2007), g: Lindegren, L. et al. (2016). Coordinates, magnitudes and spectral types are taken from Simbad.Uncertainties of photometric distances are not given explicitly for each star but are estimated to be between 16 and 25%(Neckel et al.1980). Typical uncertainty for E ( B − V ) is 0.03 mag. Peak S/N gives the signal-to-noise ratio in the region with highest photon counts,usually around 5800 ˚A MNRAS000 , 1–18 (2017) patial structure of several DIB carriers Name l b
SpectralType
V E ( B − V ) PhotometricDistance n Parallax g Observingsession Peak S/N ( ◦ ) ( ◦ ) (mag) (mag) (kpc) (mas)HD42352 196.1555 -2.5736 B1III 6.93 0.26 a ± v , a ± a , v , l ± m v , s , a , n ± v , s , a , n ± a ± a , v ± a ± a , v ± a , v f ? 0.01 ± a , v a , v a a , v a ± a a , v a a , v ± a ± a ± a ± a ± a , n a , n a , n ± a ± a , v , n n ± a , n ± a ± n ± a , n , v a , n a , n , v ± a , n a , n a , n ± a , n ± a , n a , n , v ± a , n ± s ± a , n ± a , v , n ± v , a ± s ± a , n , v a a , n ± a a ± a , v , n ± a , n , v ± a , n s ± a , n , v , l ± a , s s ± a a , n , v , l n ± s ± a , n , v ± s ± Table 1.
Continued.MNRAS , 1–18 (2017)
J. Kos
Name l b
SpectralType
V E ( B − V ) PhotometricDistance n Parallax g Observingsession Peak S/N ( ◦ ) ( ◦ ) (mag) (mag) (kpc) (mas)HD209975 104.8707 5.3906 O9.5I 5.11 0.35 a , v ± a a , n ± s a , n a , n ± a ± a , s , n v , a , s , n ± s , a a ± a ± a ± a ± a , n ± v , a a , n ± n ± a , n Table 1.
Continued
Photometric distance / kPc P a r a ll a x − / k P c Figure 2.
Comparison of photometric and trigonometric dis-tances for stars with both distances available. Four stars havetrigonometric distances larger than 3 kpc or negative parallaxand are missing from the plot. Their errorbars are still shown ifthey reach into the plotted range. 1:1 relation is shown with adashed line. of a spectral line is measured by the equivalent width (as-suming a normalized spectrum): EW = ∫ ∞ ( − I ( λ )) d λ. (5) Some of the observed spectra show interstellar absorptionlines with two components with different radial velocities. This means that the LOS passes two distinct clouds of theISM. We detected 30 such LOS. It is possible that some ofthe undetected ones also pass more than one distinct cloudof the ISM but remain undetected, if the clouds have theradial velocity with a difference lower than our resolution.In no cases we found more than two components.To detect the multi-component LOS we use Na I andCa II lines (there are two of each in our wavelength range)which are detectable in every spectrum. Atomic lines, unlikeDIBs, are unresolved, so they are as narrow as the resolu-tion allows. It is possible that Na I and Ca II doublets arenaturally broadened by cloud rotation or turbulence, but itis highly unlikely that any of these velocities would reach15 km s − , which is approximately our resolution.We fitted the Na I and Ca II doublets with a profilethat allowed the second component within a plausible veloc-ity range and strength ratio. Na and Ca profiles were fittedseparately, but both lines in each doublet were fitted at thesame time. All LOS with strength ratios larger than 0.1 andline separation larger than 15 km s − were marked possiblemulti-component. If the fitted values for the Ca doublet andthe Na doublet agreed to 10 km s − or better in line sep-aration and 50% or better for the strength ratio, the LOSwas considered multi-component. A mild criterion for thestrength ratio is due to the fact that different clouds canhave different abundances of each element.As it is evident from Kos & Zwitter (2013) or later inthis paper, some DIBs have a very poor correlation betweentheir EWs and strengths of Ca and Na lines. Therefore wecan not use the double profiles we calculate from Na andCa lines as templates for DIBs. Some DIBs like DIB 5780 that are the key to the conclusions of this paper can hardly We name DIBs by their wavelengths rounded to the nearest ˚A.A DIB at 5780 ˚A is therefore called DIB 5780. We use the samename for the molecule that produces the DIB (the carrier), asthe actual molecules are usually not known. ”DIB 5780” mighttherefore refer to the molecule and not to the absorption feature,depending on the context. MNRAS000
Comparison of photometric and trigonometric dis-tances for stars with both distances available. Four stars havetrigonometric distances larger than 3 kpc or negative parallaxand are missing from the plot. Their errorbars are still shown ifthey reach into the plotted range. 1:1 relation is shown with adashed line. of a spectral line is measured by the equivalent width (as-suming a normalized spectrum): EW = ∫ ∞ ( − I ( λ )) d λ. (5) Some of the observed spectra show interstellar absorptionlines with two components with different radial velocities. This means that the LOS passes two distinct clouds of theISM. We detected 30 such LOS. It is possible that some ofthe undetected ones also pass more than one distinct cloudof the ISM but remain undetected, if the clouds have theradial velocity with a difference lower than our resolution.In no cases we found more than two components.To detect the multi-component LOS we use Na I andCa II lines (there are two of each in our wavelength range)which are detectable in every spectrum. Atomic lines, unlikeDIBs, are unresolved, so they are as narrow as the resolu-tion allows. It is possible that Na I and Ca II doublets arenaturally broadened by cloud rotation or turbulence, but itis highly unlikely that any of these velocities would reach15 km s − , which is approximately our resolution.We fitted the Na I and Ca II doublets with a profilethat allowed the second component within a plausible veloc-ity range and strength ratio. Na and Ca profiles were fittedseparately, but both lines in each doublet were fitted at thesame time. All LOS with strength ratios larger than 0.1 andline separation larger than 15 km s − were marked possiblemulti-component. If the fitted values for the Ca doublet andthe Na doublet agreed to 10 km s − or better in line sep-aration and 50% or better for the strength ratio, the LOSwas considered multi-component. A mild criterion for thestrength ratio is due to the fact that different clouds canhave different abundances of each element.As it is evident from Kos & Zwitter (2013) or later inthis paper, some DIBs have a very poor correlation betweentheir EWs and strengths of Ca and Na lines. Therefore wecan not use the double profiles we calculate from Na andCa lines as templates for DIBs. Some DIBs like DIB 5780 that are the key to the conclusions of this paper can hardly We name DIBs by their wavelengths rounded to the nearest ˚A.A DIB at 5780 ˚A is therefore called DIB 5780. We use the samename for the molecule that produces the DIB (the carrier), asthe actual molecules are usually not known. ”DIB 5780” mighttherefore refer to the molecule and not to the absorption feature,depending on the context. MNRAS000 , 1–18 (2017) patial structure of several DIB carriers be fitted with multiple components, because they are muchbroader than the separation between the components. Suchfits would be degenerate and would not provide any moreinformation than a single component fit. We will thereforeuse LOS with multiple components qualitatively only. DIBs are observed in spectra of hot stars with few spec-tral lines of their own. Spectral lines of hot stars are wideand resolved, but still pose a problem when weak DIBs areobserved with amplitudes of an order of 1% below the con-tinuum in regions populated by stellar lines. Using any ofthe libraries of synthetic spectra to remove the stellar com-ponent from our spectra does not solve the problem, as nomodel of stellar atmospheres is able to predict the spectrumwith sub-one-percent accuracy.Accounting for stellar spectral features is trivial in somecases, as they are wider than DIBs and a local continuumaround the DIB can be estimated. This approach requires acertain level of human input and the estimates for the contin-uum are subjective. It can also introduce systematic errorsand complicates the estimation of uncertainties. In a casewhere the stellar spectral features are of the same width asDIBs, such approach can produce ambiguous measurementsof DIB profiles.In addition to stellar spectral features the signal is af-fected by a correlated noise originating from the steps wedo in the data reduction. The correlation varies with wave-length and it is also not smooth, as the final spectrum is pro-duced from combined and overlapping echelle orders. Suchcorrelated noise brings additional complications to calcula-tion of the uncertainties of measured quantities.We solve both of the above problems (spectra pollutedby stellar spectral features and correlated noise) in one stepby using Gaussian processes to fit stellar spectral featuresand correlated noise at the same time as we fit a DIB profileto the spectra.Only a 20 ˚A wide region around each DIB or an ISMabsorption line is used to fit the profile, so any spectral fea-tures with a scale > ˚A are irrelevant. The Gaussian process is a statistical model, where obser-vations occur in a continuous domain. Each data-point isa normal distribution and each collection of data-points isa multivariate normal distribution (Rasmussen & Williams2005). The Gaussian process can be used as a machine learn-ing algorithm to predict values where no data is observedand also predict the uncertainties. We use Gaussian pro-cesses as a regression engine.We want to fit a function f Θ ( x ) to N data-points ( x n , y n ) with uncertainties σ n , where Θ are free parameters. We canwrite the ln-likelihood as: ln p ( y n | x n , σ n , Θ ) = − N (cid:213) n = ( y n − f Θ ( x n )) σ n + C , (6)which can be rewritten as a matrix equation: ln p ( y n | x n , σ n , Θ ) = − r T K − r −
12 ln det ( K ) − N π, (7) where r is the residual vector r = ( y n − f Θ ( x n )) of length N and K is the covariance matrix of size N × N : K = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) σ · · · σ · · · ... ... . . . ... · · · σ N (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (8)The covariance matrix above is diagonal, so taking into theaccount white noise only. If we want to introduce correlatednoise, we simply populate off-diagonal elements with non-zero values. In our case we want to fit the correlated noise,which would be impossible if every element in the covariancematrix was a free parameter. Instead we model the covari-ance matrix with a covariance function: K ij = σ i δ ij + k ( x i , x j ) , (9)where δ ij is the Kronecker delta and k ( x i , x j ) is a covari-ance function. In our case the covariance functions will befunctions of only two parameters.We are left with maximizing the ln-likelihood to get theoptimal parameters Θ of the f Θ ( x ) function. The parametersof the covariance function can also be fitted, if the noisecharacteristics are not known. We use two different covariance functions (or kernels)to describe two components of the correlated noise. Theexponential-squared covariance function models the widestellar absorption lines: k se ( x , x (cid:48) ) = a exp (cid:18) || x − x (cid:48) || l (cid:19) , (10)and a Mat´ern / covariance function models the correlatednoise: k m / ( x , x (cid:48) ) = a (cid:32) + √ || x − x (cid:48) || l (cid:33) exp (cid:32) − √ || x − x (cid:48) || l (cid:33) , (11)where l is a characteristic width of each covariance functionand a is the scaling. Many other covariance functions arecommonly used and some are compared in Rasmussen &Williams (2005).Figure 3 shows exponential-squared and Mat´ern / co-variance functions and a few samples of random functions f ∼ GP( , k ( x , x (cid:48) )) constructed from each covariance func-tion. a equals to 1 in the given examples.Scalings a in front of each covariance function in equa-tions 10 and 11 can be estimated immediately and does nothave to be fitted. a scales as the variance of the noise wewant to describe with the given covariance function. Oneused to model the stellar component is scaled to the vari-ance of the spectrum in the wavelength range that we areinterested in, and the second one, used to model the noise, isscaled to the squared average uncertainty in the wavelengthrange of interest. The final covariance function entering themaximization of the ln-likelihood is a sum of a Mat´ern / and exponential-squared covariance functions.Characteristic widths l remain to be fitted, as they varywith the stellar type and noise characteristics. We fit thefree parameters with a Bayesian MCMC scheme (Foreman-Mackey et al. 2013). MNRAS , 1–18 (2017)
J. Kos || x − x || Exponential-squared Covariance Function l =0 . l =0 . l =1
202 20210 12 14 16 18 20 x || x − x || Matern 3/2 Covariance Function l =0 . l =0 . l =1
202 20210 12 14 16 18 20 x Figure 3.
Top: Exponential-squared covariance functions withthree different values of the characteristic width l are plottedin the left panel. Three random functions drawn from a Gaus-sian process with the exponential-squared covariance function areplotted on the right for a different l in each of the three panels.Bottom: Mat´ern / covariance functions with three different val-ues of the characteristic width l are plotted in the left panel.Three random functions drawn from a Gaussian process with theMat´ern / covariance function are plotted on the right for a dif-ferent l in each of the three panels. Our spectra are populated by many absorption lines andin order to fit the right one, we have to set some initialconditions that will constrain the fit.We fit six parameters of which four describe the DIBand two are used to describe the noise. Initial parametersfor the latter are simply the average characteristic widthsof the two covariance functions. The initial parameters forthe DIB profile are measured from the spectrum by makinga preliminary fit of the DIB profile. The same model as inthe final fit is used, but noise or uneven continuum is disre-garded. The profile is fitted with the Levenberg-Marquardmethod and the fitted parameters are used as initial param-eters for the final Bayesian fit. Initial parameters for eachof 128 walkers are perturbed, so that different walkers startat different initial conditions. Initial parameters for the co-variance function widths, DIB profile amplitude, width, andposition are now defined by a normal distribution aroundthe original value with a standard deviation of 10% and theasymmetry defined by a normal distribution around the orig-inal value with a standard deviation of 0.05. Perturbation ofthe initial conditions is unnecessary (but harmless) for fit-ting prominent DIBs, as posteriors from every walker willconverge toward the same distribution. For very weak DIBs,however, there is a possibility that a carefully selected ran-dom function can fit the DIB better than the given modelor vice versa. In this case we want to be sure that a large
Complete Fit
Wavelength / Å N o r m a li z e d f l u x DIB Models f ∼ GP (1 . , k se ( x, x )) f ∼ GP (1 . , k m / ( x, x )) Figure 4.
DIB 5850 in a spectrum of star HD 18326 with a fittedDIB and noise profile. Top: The complete fit. Black points showthe measured flux and measured uncertainties. In blue are tenrandom samples from the posterior distribution. Bottom: Decom-position of the blue curves from the top panel. In green are ran-dom functions sampled from the Gaussian process (top to bottomare the random functions sampled from a Mat´ern / covariancefunction, squared-exponential covariance function, and from thecombination of both.) In red are ten models for the DIB profile. parameter space is explored by MCMC in order to get thecorrect posterior distribution.128 Walkers are progressed for 250 steps and after thatthe parameters of the best fit are used as the initial parame-ters for the final chain. 1000 steps are made with 128 walkersand last 750 steps from each walker are used to sample theposterior.Whenever an interstellar absorption line that is not aDIB is fitted, the asymmetry parameter is set to 0, is notperturbed and is not fitted.The values of the initial fit are also used to construct aprior. The reason is that the all the initial fits were inspectedby eye and provide additional knowledge that should be usedin the prior. This is done mostly to prevent the wrong ab-sorption line be fitted when there are several prominent ab-sorption features around the DIB of interest. We use priors to prevent the fitted parameters to take val-ues that are not physically possible or outside the expectedranges. Such rigours priors as presented below are usuallynot needed, as strong DIB profiles in regions without any in-terfering stellar lines are well constrained even if we use flatpriors. But as soon a stellar line of a comparable strength isblended with a DIB, priors assure us that the DIB profile isfitted and not the stellar feature. The amount of interferenceis then only reflected in the measured uncertainties.
MNRAS000
MNRAS000 , 1–18 (2017) patial structure of several DIB carriers The characteristics of the stellar spectrum componentis hard to guess, as it changes between different regions ofthe wavelength range. If we are fitting a DIB close to a Hy-drogen line, for example, one very broad feature is expected,so a large value for l se must be used. In a region populatedonly by narrower metallic lines, the l se will be smaller. Thecorrect value for l se is therefore best to be fitted. The up-per limit for l se is not given, while the lower limit is at l se = . . The limit is somewhat arbitrary, because it onlydistinguishes the stellar spectrum noise from the correlatednoise. All the noise with the characteristic length above thisvalue is considered of stellar origin and all the noise withthe characteristic length below this value is considered to berandom correlated noise. If a star happens to have featuresnarrower than that in its spectrum, they will be treatedas correlated random noise. This also limits the correlationlength for the correlated noise to just under 3 pixels, whichmakes sense. The prior for the l m / therefore has an upperlimit of 0.22 ˚A. It also has a lower limit of 0.08 ˚A, which isthe pixel size. Because the DIBs widths fall into the samerange as the l se , we make the prior for l se lower at valuesthat are smaller than the FWHM of the DIB.We first fitted DIBs in most spectra without using anypriors in order to get the average profile shapes. The aver-age values are then used to construct priors for the remain-ing four parameters. Prior for the asymmetry is constructedfrom the average value only, priors for the FWHM and cen-tral wavelength are constructed from the averaged valuesand the initial fit, and prior for the EW is constructed solelyfrom the initial fit.Prior for the asymmetry is simply the average asymme-try measured for each DIB individually. Prior has a Gaussianshape with the FWHM of 0.15 and is limited into the rangebetween -0.75 and 0.75.Prior for the EW is a Gaussian centred on the EW weget from the initial fit. It has a lower limit of EW=0 and anupper limit which is three times the initial EW. Shape of theprior is such that the ln P drops to 0.1 at the upper limit.Because we visually inspected all initial fits, we are confidentthat the EW of three times the fitted value is impossible.Prior for the FWHM is constructed from the averagedvalue for the FWHM of each DIB and the value from theinitial fit.Prior for the central wavelength is a wide peak centredat the average value and a narrower peak centred at thevalue from the initial fit. The central wavelength is also lim-ited to ± There should always exist LOS that are close enough to eachother that the light travels through the same clouds of theISM. If two LOS also pass the same region of the cloud, weexpect the absorption lines produced in the cloud to be thesame in both LOS. Assuming some structure in the ISM wewill only see differences in LOS when the distance betweenthem approaches or exceeds the size of the structure.The two closest LOS we observed are 18 arc seconds l se / Å -1 l n P FWHM =2 . Å FWHM =0 . Å l m / / Å -1 l n P EW / Å -1 l n P EW init =0 . Å EW init =0 . Å Asymmetry -1 l n P FWHM / Å -1 l n P Inconsistent
FWHM init
Consistent
FWHM init
Central wavelength / Å -1 l n P Rough init estimateGood init estimate
Figure 5.
Examples of priors for all six fitted variables. Top-left: Prior for l se is flat at larger values but drops to −∞ at l se = . . At lower values the squared exponential covariancefunction would fall into the domain of the correlated noise. Theprior drops to −∞ gradually, starting at the value of the FWHMof the DIB, preventing the DIB to be fitted by the random func-tion. Top-right: A flat prior is used for l m / between 0.08 ˚A (oursampling and pixel size) and 0.22 (the lower limit for l se ). Middle-left: Prior for the equivalent width is a Gaussian centred on theinitial value. It is scaled as such that the ln P drops to 0.1 atthree times the initial value. We chose the equivalent width to belimited between zero (physical limit) and three times the initialequivalent width (the initial value is precise enough that largervalues should be impossible). Middle-right: Asymmetry prior iscalculated for each DIB from the distribution of asymmetries wemeasured without using this prior. A normal distribution is as-sumed for each DIB but it is limited to − . < asym < . .Bottom-left: Width prior is calculated the same way as the asym-metry prior, except that the initial value is added to the prior.This is not possible for the asymmetry prior, as the asymmetryis less constrained after the initial fit. Two examples are plotted.In black is a prior where calculated and initial widths are simi-lar (FWHM=2.05 and FWHM=2.00) and in gray is a case wherethe initial width falls into the wing of the calculated distribution(FWHM=3.1 and FWHM=2.6). Bottom-right: Prior for the cen-tral wavelength is composed of two components. The wide one iscentred at an average value we calculated without using the prior.The narrow component comes from the initial value. The widthis proportional to the uncertainty we get from the covariance ma-trix after the initial fit. Priors after a well constrained (grey) andnot so well constrained (black) initial fits are shown.MNRAS , 1–18 (2017) J. Kos apart and the next one 3.9 arc minutes. We cover LOS any-where between 3.9 arc minutes and 180 ◦ without any majorgaps. There are 21 pairs of LOSs with the separation lessthan 30 arc minutes. Assuming a distance of 1 kpc and ex-cluding the difference in distance of the two stars, two LOSwith separation of 30 arc minutes will never be more than8.7 pc apart. There are three possible explanations for differences we ob-serve between DIBs at small scales. First case is one wheredifferent DIB carriers are present in the same cloud of theISM, but the distribution within the cloud is more clumpyfor one DIB carrier than the other. Such state is illustratedin Figure 6. In the second case different DIB carriers or dif-ferent groups of DIB carriers exist in two or more clouds,where each carrier or group dominates one cloud. One LOSmight therefore pass one cloud but not the other. Such stateis illustrated in Figure 7. The third option is that there is agradient within the cloud, so the two LOS pass the samecloud but pass the regions with different composition. Asource of the gradient can be a nearby hot star, a shockwave or other external radiation. This case is illustrated inFigure 8. This case can be seen as a special case of the firstmodel, but can be tested independently from the first model,so we treat it separately. It is also possible that the actualstate of the ISM is a combination of the three options. Ourobservations are able to detect any combination of the abovemodels as well.Case 2 can be best confirmed by observing DIBs withdifferent radial velocities. If there exist two separate cloudsof gas there is no reason that they should always have thesame radial velocity and in as many pairs as we observed,we expect some pairs to have different radial velocities evenat small angle-of-separation.Case 3 can be confirmed if we observe several LOS thatsample the same gradient. A smooth transition between theratios of different DIBs confirms a smooth gradient. If theratios do not vary smoothly, the ISM is in the state fromcase 1.
An excellent example of a big difference in DIB strengths atsmall angles-of-separation are pairs HD 25443 and a doublestar HD 25638 + HD 25639. HD 25443 is 0.3 ◦ away fromthe other two. Figures 9 and 10 show the pairs HD25443,HD25638, and HD25638, HD25639, respectively. Each fig-ure shows all relevant DIBs, atomic and molecular lines. Alsoshown in each figure is the map with positions of stars on thesky together with dust map and its polarisation, CO emis-sion map, synchrotron emission and its polarisation (PlanckCollaboration et al. 2015), and near UV diffuse radiation(Murthy 2014). In the panel showing each absorption lineis written the ratio between strengths of both lines in thepair and difference in radial velocity. Figures for other pairswith the angle-of-separation smaller than 1 ◦ are available inAppendix B.In Figure 9 DIB 5780 is equally strong in both spectra. Case 1: Different species have different clumpiness
Figure 6.
Two different species can be in the same cloud ofthe ISM, but one of them is more clumpy than the other. Thedifference in the equivalent width between both LOS is thereforemore probable to be larger for the more clumpy species. DIBs5780 and 5797 are used as an example of two different species ofthe ISM.Case 2: Different species are in separate clouds
Figure 7.
There exist two clouds in the LOS. In each cloud onespecies dominates.Case 3: Different species create a gradient due to an externalinfluence
Figure 8.
A cloud of the ISM in a quiescent state would havesimilar distribution of both species. If an external source createsa gradient in the physical properties within the cloud, a chemicalgradient follows. An external source is not necessarily a star (star3 in figure), it can be a shockwave, cosmic ray radiation, or otherforms of radiation. MNRAS000
A cloud of the ISM in a quiescent state would havesimilar distribution of both species. If an external source createsa gradient in the physical properties within the cloud, a chemicalgradient follows. An external source is not necessarily a star (star3 in figure), it can be a shockwave, cosmic ray radiation, or otherforms of radiation. MNRAS000 , 1–18 (2017) patial structure of several DIB carriers Figure 9.
Maps of the region around the pair, DIBs, atomic, and molecular lines for HD 25443 (red) and HD 25638 (blue). Some detailsabout the two LOS are given in the top-right corner.Note the difference in some DIBs.
On contrary, the DIB 5797 is weaker toward HD 25443. Themost extreme ratio is seen in DIB 5512, which is almostabsent in HD 25443 LOS. DIB 6196 is the only one strongerin HD 25443 LOS. Interesting is also the CH to CH + ratiobeing higher in HD 25638/25639 LOS. Note that the colorexcess of the pair is not the same either. In this case thecolor excess ratio corresponds to the ratio of the DIB 5797,but this is not the rule in other pairs.In Figure 10, where the two LOS are much closer toeach other, there are no significant differences between anyabsorption lines.If we try to test the three models for the structure ofthe ISM clouds based on the DIBs observed in the men-tioned two pairs, we see that the observations support thefirst model and do not exclude the third model, but thesecond model is only plausible if all clouds have the same radial velocity. Figures 6, 7, and 8 have been produced withthe pair HD25443, HD2563 in mind, trying to explain theobserved strengths of DIBs 5780 and 5797. To describe the differences between EWs between pairs ina more general way, without concentrating on selected pairslike in the previous section, we make a diagram where weplot the difference for every possible pair as it depends onthe angle-of-separation.The difference d between measured EWs of DIBs in LOS and is defined as: d = || EW − EW || EW + EW (12) MNRAS , 1–18 (2017) J. Kos
Figure 10.
Maps of the region around the pair, DIBs, atomic, and molecular lines for HD 25638 (red) and HD 25639 (blue). Somedetails about the two LOS are given in the top-right corner.Note that all DIBs are the same in both LOS. or in the case of color excess, the difference is d = || E ( B − V ) − E ( B − V ) || E ( B − V ) + E ( B − V ) . (13)The difference is expected to be close to zero at zeroangle-of-separation. It is not expected to be exactly zero, asthe two stars on the same LOS can be positioned at differentdistances. It is most likely, however, that the same cloud ofthe ISM is in the way, so assuming a very small differenceis safe. At large separations the expected difference is large,but not necessary one, due to the selection function of ourstars. We observe reddened stars, so it is unlikely that anyDIB will have a zero EW. More important, the difference isexpected to be constant regardless the separation angle, aslong as the LOS separation is larger than the angle at whichthe correlation between different LOS exists. The highestpossible difference of one is possible for weaker DIBs where the measured EW is consistent with zero. Thanks to our rig-orous approach of estimating the probability density func-tion for each measured EW, such measurements are com-pletely consistent with measurements of stronger DIBs. Onecan see from Figure 11 and Appendix A that a higher max-imum difference is indeed more probable for weaker DIBs.Maximum difference is therefore one of the parameters todescribe the shape of the pairwise correlation diagram.In addition to the difference at large angles we intro-duce two more parameters by fitting a broken power lawto the 95th percentiles of the pairwise correlation diagram.The two parameters are the position of a knee (the break inthe power law) and the slope at lower angles. In contrast tothe difference at large angles, the latter two parameters havean actual physical meaning. The reason we chose the 95thpercentile to describe the maximum difference at each angle- MNRAS000
Maps of the region around the pair, DIBs, atomic, and molecular lines for HD 25638 (red) and HD 25639 (blue). Somedetails about the two LOS are given in the top-right corner.Note that all DIBs are the same in both LOS. or in the case of color excess, the difference is d = || E ( B − V ) − E ( B − V ) || E ( B − V ) + E ( B − V ) . (13)The difference is expected to be close to zero at zeroangle-of-separation. It is not expected to be exactly zero, asthe two stars on the same LOS can be positioned at differentdistances. It is most likely, however, that the same cloud ofthe ISM is in the way, so assuming a very small differenceis safe. At large separations the expected difference is large,but not necessary one, due to the selection function of ourstars. We observe reddened stars, so it is unlikely that anyDIB will have a zero EW. More important, the difference isexpected to be constant regardless the separation angle, aslong as the LOS separation is larger than the angle at whichthe correlation between different LOS exists. The highestpossible difference of one is possible for weaker DIBs where the measured EW is consistent with zero. Thanks to our rig-orous approach of estimating the probability density func-tion for each measured EW, such measurements are com-pletely consistent with measurements of stronger DIBs. Onecan see from Figure 11 and Appendix A that a higher max-imum difference is indeed more probable for weaker DIBs.Maximum difference is therefore one of the parameters todescribe the shape of the pairwise correlation diagram.In addition to the difference at large angles we intro-duce two more parameters by fitting a broken power lawto the 95th percentiles of the pairwise correlation diagram.The two parameters are the position of a knee (the break inthe power law) and the slope at lower angles. In contrast tothe difference at large angles, the latter two parameters havean actual physical meaning. The reason we chose the 95thpercentile to describe the maximum difference at each angle- MNRAS000 , 1–18 (2017) patial structure of several DIB carriers Figure 11.
Top left: Pairwise correlation between the angle-of-separation and difference in the color excess. To take the uncertaintiesof the color excess into the account, a normal distribution around the mean value (tabulated in Table 1) is taken. 100 samples are takenfrom each distribution which are used to calculate the difference. To plot the distribution, we simply plot the samples. We also scatterthe samples around the exact angle separation for the plot to be easier to read. The 95th percentiles of the samples are calculated withinlog-sized bins (shown in red). A broken power law is fitted and shown in blue. Dashed lines show the uncertainties of the fit. The fittedposition of the knee and the slope are given in the top-left corner. Top right: Same as in previous panel, but for the pairwise correlationbetween the angle-of-separation and difference in the EW of the DIB 5780. The 100 samples are taken from the measured probabilitydensity function of the EW. Bottom left: Same as previous panel but for DIB 5797. Bottom right: Same as previous panel but for DIB5512. of-separation is purely arbitrary. The percentile chosen mustbe less than 100 to avoid any outliers (statistical or system-atic) and choosing any other large percentile gives consistentresults. Figure 11 shows four pairwise correlation diagramsused to illustrate the main differences between species. Therest of the diagrams are collected in Appendix A.The diagrams for three DIBs in Figure 11 show thatthe distribution of DIB 5780 is smoother than that for DIBs5797 and DIB 5512. In the case of DIB 5780 the knee isat the angle-of-separation of more than 3 ◦ . The knee is atobviously lower value for DIB 5797, and for DIB 5512 itis practically undetectable as only three data-points lie atangles-of-separation smaller than the position of the knee.In the latter case more data at low angle separation shouldbe taken to actually determine the position of the knee withany significance. The sole data-point at the lowest angle-of-separation is more or less irrelevant to the fitted brokenpower law, as it has a big uncertainty. The differences atthe level of a few percent carry a huge uncertainty, as themeasured EWs are only precise to a few percent even forthe strongest DIBs. The difference is therefore almost alwaysconsistent with zero. It is also consistent with measurements in Cordiner et al. (2013) made at the same or smaller scaleas ours. Results presented in the previous section are in strong agree-ment with the first model, where different DIBs have differ-ent clumpiness.We were unable to find any tight pairs of LOS whereDIBs would have different radial velocities and the angle-of-separation is smaller than the break in the pairwise corre-lation diagram. We also found no such LOS where differentDIBs would have different radial velocities. With the amountof observed LOS it is safe to conclude that the second model,where different DIBs dominate different clouds, can be re-jected.We also did not find any regions with a smooth gradi-ent in a ratio between any two DIBs. We conclude that thegradients only exist at smaller scales then the scales we canreach with our approach or they are completely dominatedby the first model.
MNRAS , 1–18 (2017) J. Kos
Line Position ofthe knee Line Position ofthe knee( ◦ ) ( ◦ )DIB 5512 . ± . DIB 4780 . ± . Ca I . ± . DIB 5850 . ± . DIB 5705 . ± . DIB 6090 . ± . DIB 6445 . ± . Ca II . ± . DIB 6379 . ± . E(B-V) . ± . CH . ± . DIB 6613 . ± . CH + . ± . DIB 4964 . ± . DIB 6270 . ± . DIB 5797 . ± . DIB 4762 . ± . DIB 6202 . ± . DIB 6660 . ± . DIB 5780 . ± . DIB 5545 . ± . DIB 6196 . ± . DIB 4726 . ± . Table 2.
Positions of the knee are collected for all measured DIBsand ISM species in an ascending order.
We explored three simple but meaningful models for the spa-tial distribution of the DIB species. Only one can be con-firmed with our data, showing that the DIB species havedistinctly different clumpiness. This was showed by plottingpairwise correlations between DIB strengths in single-linedLOS and fitting a two parameter model showing that thecorrelation breaks at different scales for different DIBs (seetable 2 for a condensed list of results). Assuming a typicaldistance to the observed stars to be 1 kpc, the correlationfor DIB 5512 breaks down at 4 pc and for DIB 6196 at 62 pc.Since the actual clouds are closer than the stars, these num-bers are only the upper limits. This range is in agreementwith scales for some atomic absorption lines calculated in asimilar way as ours (e.g. Smoker et al. 2015).We did not explore the second parameter of the fittedmodel (the slope of the correlation at small scales) as thedata is rarely sufficient to constrain it. It can easily be imag-ined that the slope depends on the fractal structure of theISM, a relation between big and small clumps of the DIB-bearing medium. More of the closest pairs should have beenobserved to constrain the slope, but there are not many inthe magnitude and color range of the targeted stars. Accu-racy of the measured DIB strengths is another concern whenmeasuring the slope, because small differences expected attight angles-of-separation are hard to measure. More pre-cise measurements require better resolution, higher SNR andability to remove the stellar absorption lines from the spec-tra. These are beyond the limitations of our small survey.Observing finer structure than we targeted in this workis a sensible goal for large spectroscopic surveys where a hugenumber of close pairs can be observed. This work has evenfurther application to large spectroscopic surveys of starsand the study of the ISM within them. DIBs are commonlyused as tracers of the ISM in large spectroscopic surveys.Unlike other ISM absorption lines in the spectra of stars,DIBs are found in all bands covered by contemporary VISand NIR surveys. Together with the distance to the observedstars they can give a three-dimensional picture of the inter-stellar medium (Kos et al. 2014; Zasowski et al. 2015). DIBsalso serve as tracers of physical conditions that other absorp-tion lines are not sensitive to (Bailey et al. 2016). But oftenthe spectra in a small region have to be combined to achievea SNR required to observe DIBs (e.g. Kos et al. 2014; Za-
Rank P o s i t i o n o f t h e k n ee / ◦ C a I C H + C H C a II E ( B - V ) EWE(B − V) = 0 . Å magEWE(B − V) = 0 . Å magEWE(B − V) = 1 . Å mag Å P o s i t i o n o f t h e k n ee / ◦ Figure 12.
Position of the knee for every observed line. Size of thesymbol is proportional to an average line strength normalized tocolor excess. DIBs are plotted in blue and other interstellar linesin red. Two insets in the top-left show the relations between theposition of the knee and the remaining two parameters describingthe line shape (width and asymmetry). sowski et al. 2015; Puspitarini et al. 2015). We derived thesize of a region in which we can expect the LOS to pene-trate the same cloud of the ISM. It turns out that for someDIBs this size is considerably larger than for others. Thesevalues (properly modified due to the selection function) cantherefore be used in large spectroscopic surveys to amplifythe SNR when studying the ISM and select most suitableDIBs for a particular case. The fine structure can then beexplicitly measured, instead of inferred like it is in this work.A significant diversity of the clumpiness of the DIB-bearing ISM reveals the nature of the DIB carriers. We ex-pect the DIB carriers with a smoother (less clumpy) distri-bution to be ionised. Ionised molecules are more likely to befound in the outer layer of the clouds, and are more likelyto smooth out the finer structure due to higher tempera-ture. At least DIB 5780 has been attributed to a chargedmolecule in independens studies (Milisavljevic et al. 2014).However, we can not divide DIB carriers into ionised andneutral molecules solely from the measured clumpiness, asthe structure depends on molecular mass and chemistry ofthe ISM. DIBs 6202, 5780, and 6196 however stand out fromthe rest (see Figure 12), as their correlation breaks at an-gles larger than 3 ◦ , which is considerably more than for anyother DIB. This suggests that these are absorption lines ofionized molecules. We do not observe any significant corre-lation between the position of the knee and either strength,wavelength, width, or asymmetry of DIBs.Pairwise correlation as presented in this paper is a sim-ple description of the shape of the ISM clouds and dependsonly on a couple parameters related to the physical proper-ties of the actual carriers. It is therefore very convenient fortesting the theoretical models for the DIB carriers, shouldany such models appear in the future. MNRAS000
Position of the knee for every observed line. Size of thesymbol is proportional to an average line strength normalized tocolor excess. DIBs are plotted in blue and other interstellar linesin red. Two insets in the top-left show the relations between theposition of the knee and the remaining two parameters describingthe line shape (width and asymmetry). sowski et al. 2015; Puspitarini et al. 2015). We derived thesize of a region in which we can expect the LOS to pene-trate the same cloud of the ISM. It turns out that for someDIBs this size is considerably larger than for others. Thesevalues (properly modified due to the selection function) cantherefore be used in large spectroscopic surveys to amplifythe SNR when studying the ISM and select most suitableDIBs for a particular case. The fine structure can then beexplicitly measured, instead of inferred like it is in this work.A significant diversity of the clumpiness of the DIB-bearing ISM reveals the nature of the DIB carriers. We ex-pect the DIB carriers with a smoother (less clumpy) distri-bution to be ionised. Ionised molecules are more likely to befound in the outer layer of the clouds, and are more likelyto smooth out the finer structure due to higher tempera-ture. At least DIB 5780 has been attributed to a chargedmolecule in independens studies (Milisavljevic et al. 2014).However, we can not divide DIB carriers into ionised andneutral molecules solely from the measured clumpiness, asthe structure depends on molecular mass and chemistry ofthe ISM. DIBs 6202, 5780, and 6196 however stand out fromthe rest (see Figure 12), as their correlation breaks at an-gles larger than 3 ◦ , which is considerably more than for anyother DIB. This suggests that these are absorption lines ofionized molecules. We do not observe any significant corre-lation between the position of the knee and either strength,wavelength, width, or asymmetry of DIBs.Pairwise correlation as presented in this paper is a sim-ple description of the shape of the ISM clouds and dependsonly on a couple parameters related to the physical proper-ties of the actual carriers. It is therefore very convenient fortesting the theoretical models for the DIB carriers, shouldany such models appear in the future. MNRAS000 , 1–18 (2017) patial structure of several DIB carriers ACKNOWLEDGEMENTS
JK is grateful to the Observatory of Padua TAC for a gener-ous amount of observing time and to Maruˇska ˇZerjal, TomaˇzZwitter, and Gregor Traven for assisting with observations.This work is based on observations collected at Copernicotelescope (Asiago, Italy) of the INAF - Osservatorio Astro-nomico di Padova. JK is supported by a Discovery Projectgrant from the Australian Research Council (DP150104667)awarded to J. Bland-Hawthorn and T. Bedding. JK is grate-ful to them both for helpful conversations during the courseof this research.
REFERENCES
Bailey M., van Loon J. T., Farhang A., Javadi A., KhosroshahiH. G., Sarre P. J., Smith K. T., 2016, A&A, 585, A12Baron D., Poznanski D., Watson D., Yao Y., Cox N. L. J.,Prochaska J. X., 2015, MNRAS, 451, 332Cami J., Sonnentrucker P., Ehrenfreund P., Foing B. H., 1997,A&A, 326, 822Campbell E. K., Holz M., Gerlich D., Maier J. P., 2015, Nature,523, 322Cordiner M. A., Fossey S. J., Smith A. M., Sarre P. J., 2013, TheAstrophysical Journal Letters, 764, L10Fitzpatrick E. L., Massa D., 2007, ApJ, 663, 320Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Galazutdinov G., Moutou C., Musaev F., Kre(cid:32)lowski J., 2002,A&A, 384, 215G(cid:32)l¸ebocki R., Gnaci´nski P., 2005, in Favata F., Hussain G. A. J.,Battrick B., eds, ESA Special Publication Vol. 560, 13th Cam-bridge Workshop on Cool Stars, Stellar Systems and the Sun.p. 571Heger M. L., 1922, Lick Observatory Bulletin, 10, 146Herbig G. H., 1995, ARA&A, 33, 19Jenniskens P., Desert F.-X., 1994, A&AS, 106, 39Kos J., Zwitter T., 2013, ApJ, 774, 72Kos J., et al., 2014, Science, 345, 791Krelowski J., 1989, Astronomische Nachrichten, 310, 255Kre(cid:32)lowski J., 2014, in Cami J., Cox N. L. J., eds, IAU Sym-posium Vol. 297, The Diffuse Interstellar Bands. pp 23–33,doi:10.1017/S1743921313015536Kre(cid:32)lowski J., Walker G. A. H., 1987, ApJ, 312, 860Kre(cid:32)lowski J., Beletsky Y., Galazutdinov G. A., Ko(cid:32)los R.,Gronowski M., LoCurto G., 2010, ApJ, 714, L64Kroto H., 1988, Science, 242, 1139L´eger A., d’Hendecourt L., 1985, A&A, 146, 81Lindegren, L. Lammers, U. Bastian, U. Hern ˜A (cid:44) andez, J. Klioner,S. Hobbs, D. Bombrun, A. Michalik, D. 2016, A&AMaeder A., 1975, A&A, 42, 471Maier J. P., Walker G. A. H., Bohlender D. A., Mazzotti F. J.,Raghunandan R., Fulara J., Garkusha I., Nagy A., 2011, ApJ,726, 41Merrill P. W., 1936, ApJ, 83, 126Milisavljevic D., et al., 2014, ApJ, 782, L5Motylewski T., et al., 2000, ApJ, 531, 312Murthy J., 2014, ApJS, 213, 32Neckel T., Klare G., Sarcander M., 1980, A&AS, 42, 251Planck Collaboration et al., 2015, preprint, ( arXiv:1502.01588 )Puspitarini L., et al., 2015, A&A, 573, A35Rasmussen C. E., Williams C. K. I., 2005, Gaussian Processesfor Machine Learning (Adaptive Computation and MachineLearning). The MIT PressSalama F., Galazutdinov G. A., Kre(cid:32)lowski J., Allamandola L. J.,Musaev F. A., 1999, ApJ, 526, 265 Sarre P. J., 2006, Journal of Molecular Spectroscopy, 238, 1Sarre P. J., Miles J. R., Kerr T. H., Hibbins R. E., Fossey S. J.,Somerville W. B., 1995, MNRAS, 277, L41Savage B. D., Massa D., Meade M., Wesselius P. R., 1985, ApJS,59, 397S(cid:32)lyk K., Galazutdinov G. A., Musaev F. A., Bondar A. V.,Schmidt M. R., Kre(cid:32)lowski J., 2006, A&A, 448, 221Smoker J. V., Keenan F. P., Fox A. J., 2015, A&A, 582, A59Snow Jr. T. P., York D. G., Welty D. E., 1977, AJ, 82, 113Valencic L. A., Clayton G. C., Gordon K. D., 2004, ApJ, 616, 912Vos D. A. I., Cox N. L. J., Kaper L., Spaans M., Ehrenfreund P.,2011, A&A, 533, A129Walker G. A. H., Bohlender D. A., Maier J. P., Campbell E. K.,2015, ApJ, 812, L8Westerlund B. E., Kre(cid:32)lowski J., 1989, A&A, 218, 216Zasowski G., et al., 2015, ApJ, 798, 35van Loon J. T., Smith K. T., McDonald I., Sarre P. J., FosseyS. J., Sharp R. G., 2009, MNRAS, 399, 195van Loon J. T., et al., 2013, A&A, 550, A108van der Zwet G. P., Allamandola L. J., 1985, A&A, 146, 76MNRAS , 1–18 (2017) J. Kos
APPENDIX A: PAIRWISE CORRELATIONDIAGRAMS
Figures in this section show pairwise correlation between theangle-of-separation and difference in the EW of DIBs. Totake the uncertainties of the EW into the account, 100 sam-ples are taken from each measured probability distributionfor the EW of the DIB in question. To plot the distribution,we simply plot the samples. We also scatter the samplesaround the exact angle separation for the plot to be easierto read. The 95th percentiles of the samples are calculatedwithin log-sized bins (shown in red). A broken power law isfitted and plotted in blue (uncertainties are plotted with adashed line). The fitted position of the knee and the slopeare given in the top-left corner of each diagram.
MNRAS000
MNRAS000 , 1–18 (2017) patial structure of several DIB carriers APPENDIX B: DIAGNOSTIC PLOTS FORPAIRS OF LOS
Plots in this appendix show all the analysed ISM lines in allpairs that are 1 ◦ or closer in separation. Each figure showsall relevant DIBs, atomic and molecular lines. Also shownin each figure is the map with positions of stars on the skytogether with dust map and its polarisation, CO emissionmap, synchrotron emission and its polarisation, and nearUV diffuse radiation. In the panel showing each absorptionline is written the ratio between strengths of both lines inthe pair and difference in radial velocity. MNRAS , 1–18 (2017) J. Kos
This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000