Investigation of cosmic ray penetration with wavelet crosscorrelation analysis
MaxPlanckInstitut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany
email: ryang@mpihd.mpg.de
Received: 15 March 2016
Accepted: 20 May 2016
Aims. We use FermiLATand Planck data to calculate the cross correlation between γray signal and gas distribution in different scales in giant molecular clouds (GMCs). Then we investigate the cosmic rays (CRs) penetration in GMCs with these informations.
Methods. We use the wavelet technique to decompose both the γray and dust opacity maps in different scales, then we calculate the wavelet cross correlation functions in these scales. We also define wavelet response as an analog to the impulsive response in Fourier transform and calculate that in different scales down to FermiLATangular resolution.
Results. The γray maps above 2 GeV show strong correlation with the dust opacity maps, the correlation coefficient is larger than 0.9 above a scale of 0.4 degree.The derived wavelet response is uniform in different scales.
Conclusions. We argue that the CR above 10 GeV can penetrate the GMC freely and the CRs distributions in the same energy range are uniform down to parsec scale.
Key words: gamma rays: ISM / cosmic rays / ISM: clouds
© ESO, 2016
1. Introduction
Cosmic rays (CRs) play an important role in the determination of ionization level in molecular clouds (MC), and thus have an impact on the dynamical process and star formation (for a review, see Dalgarno 2006). The penetration of CRs inside the MCs is a key issue to evaluate the CR density therein. On the other hand, MCs are also regarded as the CR barometers (Aharonian 1991, 2001; Casanova et al. 2010) and thus are important targets in the γray astronomy. The exclusion of CRs from MCs may reduce the γray brightness and alter the spectral shape. The issue of the penetration or exclusion of CRs from MCs has been investigated by several groups (Skilling & Strong 1976; Dogiel & Sharov 1990; Gabici et al. 2007; Morlino & Gabici 2015). However, in these papers quite different results have been derived, from almost freepenetration to exclusion of CRs up to 100 GeV. Generally, the propagation of CRs depends strongly on the turbulence spectrum inside MCs. The high gas density in MCs cause high damping rate, which reveal that the propagation of CRs inside MCs are dominated by convection rather than diffusion. In such a scenario the CRs can freely penetrate into the dense core of the MCs except the low energy particles, which may suffer from the ionization losses (Morlino & Gabici 2015). On the other hand Istomin & Kiselev (2013) argue that the magnetic turbulence can also be generated in the weak ionized gas by the kinematic turbulence in the neutral gas. Thus the CRs exclusion can occur due to the slower diffusion induced by the higher turbulence level (see e.g. Gabici et al. 2007).
Fig. 1 Upper panel: γray counts map above 2 GeV after point source subtraction in the region Orion A (left) and Taurus (right). Lower panel: τ_{353} maps for Orion A (left) and Taurus (right). 

Open with DEXTER 
Various observational efforts have also been made on the interaction of CRs and MCs, both on γrays (see e.g. Ackermann et al. 2012a,b; Yang et al. 2014, 2015; Planck Collaboration Int. XXVIII 2015) and ionizations (Padovani et al. 2009; Nath et al. 2012). However, the γray analysis treated MCs as a whole rather than investigated the CRs penetration in different scales, while the ionization study can only account for low energy part of CRs distributions and the contamination from electrons is hard to exclude. Direct observation evidence is still lacking in determination of the different scenarios mentioned above. In this paper we analyze the FermiLATand Planck dust opacity data in Taurus and Orion A region and use a wavelet crosscorrelation technique (see e.g. Frick et al. 2001; Arshakian & Ossenkopf 2016) to study correlation between γray maps and dust opacity maps in different scales down to the instrument angular resolution. The correlation factor and wavelet response derived can give us unambiguous information on the CRs density in different scales. We also compare the derived cross correlation and wavelet response with the simulation results to test the technique.
Fig. 2 Wavelet decomposition of τ_{353} map of Taurus with the scale of 0.2, 0.5, 1.0, 1.5 degrees (from top left to bottom right). 

Open with DEXTER 
This paper is organized as follows. In Sect. 2 we describe the data reduction and maps used in this study, in Sect. 3 we briefly describe the wavelet decomposition technique and calculate the wavelet cross correlation and wavelet impulsive response at different scales, in Sect. 4 we perform simulations with different scenarios and compare them with the results from the data, and in Sect. 5 we conclude and discuss the implication of our results.
2. γray and dust column maps
We chose Orion A and Taurus region as our targets. Both of them are famous nearby giant molecular clouds (GMCs) and significantly detected in γrays (Yang et al. 2014). Orion A is located 500 pc away from the Sun and has a mass of about 1.2 × 10^{5}M_{⊙}, while Taurus is 140 pc away and has a mass of about 0.3 × 10^{4}M_{⊙}. The diffuse emission induced mainly by neutral pion decay in the interaction of CRs and ambient gas dominates the contributions from the point sources in both regions, which makes it a ideal place to investigate the interaction of CRs and ambient gas.
For FermiLATdata we have selected observations with over the period which covers seven years exposure time (MET 239 557 417 – MET 451 533 077). For the data reduction, we used the standard LAT analysis software package v9r33p0^{1}. We want to use the data with high angular resolution to investigate the very inner core of MCs. The FermiLATangular resolution improves with the rising energy. The 68% containment angles are about 5.0°,0.8°,0.5°, and 0.2° at 0.1, 1, 2 and 10 GeV, respectively^{2}. Thus the higher angular resolution can be only achieved at higher energy bands, where the statistics are much worse. In wavelet study the shooting noise may dominate the signal in small scales with lower significance. On the other hand, the front converted photon events have a better angular resolution, the 68% containment angle at 2 GeV is about 0.35°, rather than 0.5° degree for all events. Thus we selected front converted events with energies exceeding 2000 MeV to satisfy the requirement on both statistics and angular resolution.
The regionofinterest (ROI) was defined to be a 10° × 10° square centered on the position of Taurus and Orion A, respectively. In order to reduce the effect of the background related to the Earth’s albedo, we excluded from the analysis the time intervals when the Earth was in the fieldofview (more specifically, when the centre of the fieldofview was 52° above the zenith), as well as the time intervals when parts of the ROI were observed at zenith angles >100°. Although point sources only contribute small part of the total emission in this regions, to minimize these contaminations we subtract the best fit point sources’ flux in both regions. The P8R2_v6 version of the postlaunch instrument response functions (IRFs) was adopted in the likelihood fitting. We also use the 3FGL catalog (Acero et al. 2015) and standard diffuse emission templates^{3}. The spectral parameter of the point sources with ROI and the diffuse components are left free in the fitting. Then we subtract the best fitting point sources model maps generated by gtmodel tool from the raw counts map. The point sources subtracted counts maps can be found in Fig. 1.
The Planck dust opacity maps are used to model the gas distribution. The detailed description of the dust opacity maps can be found in Planck Collaboration XIX (2011). Although Planck Collaboration XI (2014) argue that the radiance may be a more proper tracer of gas in the diffuse interstellar medium, we still chose τ_{353} as the tracer since the column in Taurus and Orion A is relatively high (~10^{22} cm^{2}). The τ_{353} map can also be found in Fig. 1. And we use the relation(1)where cm^{2} (Planck Collaboration XIX 2011) to convert dust opacity into gas column.
3. Wavelet transform and wavelet cross correlation
Wavelet transform of function f(x) is defined as(2)Here x is the physical coordinates, ψ(x) is the analysing wavelet (real or complex, * indicates the complex conjugation), a is the scale parameter, and κ is a normalization parameter.
Wavelet transform can decompose images into different scales (Frick et al. 2001). Here we use the Mexican hat wavelet transform to decompose the images, which can be expressed as(3)As an example we show the decomposition of τ_{353} map of Taurus in different scales in Fig. 2. It is clear that on high resolution maps only the dense cores emerge and in low resolution ones there are only large scale structures.
Fig. 3 Wavelet crosscorrelation coefficient of Orion A (left panel) and Taurus (right panel), respectively. The gray shaded areas show the simulation results assuming a uniform CR density at all scales, while the red areas show the results when assuming a smaller diffusion coefficient inside MC. χ is defined in Eq. (7). 

Open with DEXTER 
The wavelet cross correlation coefficient of two maps can be defined as(4)where(5)is the wavelet spectrum.
To calculate the cross correlation coefficients between the γray maps and dust opacity maps we first smooth both map into the same angular resolution. This is done by following the relation that , where σ_{f} and σ_{o} are the final and origin map resolution and σ_{s} is the width of the smoothing kernel. The angular resolution for FermiLATand beam width for Planck maps can be found in the calibration database (CALDB) files in Fermi Science tools and Planck Collaboration VII (2016), respectively. If the point spread function has a Gaussian form, the 68% containment angle is identical to σ of the Gaussian function by definition. Although the PSF of FermiLATis not a perfectly Gaussian^{4}, the tail is not relevant for the analysis described here. Thus we assume for FermiLATmaps σ_{o} = 0.35°. To make full use of the FermiLATangular resolution we chose σ_{f} to be 0.4° in this work.
The results of the crosscorrelation coefficients are shown in Fig. 3. The results show a strong correlation above 0.4 degrees (correlation coefficient larger than 0.9). The error bars are derived using jackknife method. In jackknife we divided the ROI into 100 groups, 1° × 1° each, then we masked one of the subgroups each time and calculated the corresponding wavelet cross correlation coefficients in the masked map. The final errors can be estimated as , where σ_{cc} is the standard deviation of the wavelet crosscorrelation coefficients in all the jackknife samples.
Fig. 4 The wavelet response of Orion A (left panel) and Taurus (right panel), respectively. Gray shaded areas show the simulation results assuming a uniform CR density at all scales, while the red areas show the results when assuming a smaller diffusion coefficient inside MC. 

Open with DEXTER 
Fig. 5 Wavelet crosscorrelation coefficient (left panel) and wavelet response (right panel) of 100 realization of null test simulations. 

Open with DEXTER 
To further investigate the relations between γray and dust opacity maps we also define a wavelet response as (6)which is an analog of the impulsive response in Fourier transform. If we assume that subscript 1 represents the dust opacity maps in Eq. (3), the wavelet response should be proportional to the γray emissivity per Hatom, that is, the CR density in different scales. The derived wavelet response are shown in Fig .4. The flatness shows a uniform CR density in all scales.
To test the method we perform a null test on the random simulated maps. At each time step we sample two random maps with poisson distribution and with the same average value as the γray map and dust column map, respectively. Then we calculate the wavelet cross correlation coefficient and wavelet response as described above. We perform 100 realizations of the simulations and calculate the mean and standard deviation of the wavelet cross correlation coefficient and wavelet response. The results are shown in Fig. 5, which are perfectly consistent with the zero results.
4. CR penetrations and exclusions
To test the power of our method we also calculate the wavelet cross correlation coefficient and wavelet response by using the simulated γray maps, for both the CR freepenetration case and CR exclusion case.
At first we simulate the FermiLATγray maps by assuming a uniform CR density, which is identical to the local interstellar spectrum (LIS) of CRs (Casandjian 2015), inside MC. To do this we first calculate the γray emissivities per Hydrogen atom by using the LIS spectrum and the neutral pion decay cross section parametrized by Kafexhiu et al. (2014). We then multiply this value with the gas column to derive the γray emissivities at each pixel of the τ_{353} map. We then multiply the exposure map to get the model map. Finally we draw samples at each pixel of the “model map” from a Poisson distribution. We then calculate the wavelet cross correlation coefficient and wavelet response with the simulated γray map and dust opacity maps. The results of 100 realizations of such simulations are shown in gray shaded area in Figs. 3 and 4, which are in very good agreement with the data. We neglect the shooting noise in the Planck opacity maps in our simulations because of the relatively higher sensitivity of Planck.
We also simulate the results for CR exclusion cases. The estimation of CR exclusions from MC is done in the same framework as described in Gabici et al. (2007), where a parameter χ was adopted to describe the ratio between the diffusion coefficient inside MC and in the Galaxy. Thus the diffusion coefficient inside MC can be parametrized as(7)The CR inside MCs is characterized by two time scales, say, diffusion time scale τ_{diff} and cooling time scale τ_{cool}, where (8)and(9)By equating τ_{diff} and τ_{cool}, and taking into account the correlation between magnetic field strength and gas density(10)we get(11)where σ is the gas column. Below E^{∗} CR cannot penetrate into the core of MC. It is clear for χ = 1 case E^{∗} is so small that all the CRs can penetrate into the cloud freely. In this case the CR density inside MC should be uniform. On the other hand, for χ = 0.01, E^{∗} ~ 10 GeV, below this value the CRs cannot penetrate the cloud with column σ ~ 10^{22} cm^{2}. We chose χ = 0.01 as our fiducial model for CR exclusion case. Remarkably, our rule of thumb estimation came to a similar result with the numerical solution of the transport equations in Gabici et al. (2007) (see Fig. 1 and text therein). Thus to calculate the γray emissivity per H atom in this case we simply assume a sharp cutoff below E^{∗} in the CR spectrum. Then we calculate E^{∗} for each pixel in the τ_{353} map and then calculate the corresponding γray emissivities. The following procedure is the same as the simulations for the uniform CR case. We note that if E^{∗} is 10 GeV the derived emissivities should be similar to the uniform CR case since in our analysis we only consider γray maps above 2 GeV. But E^{∗} scales as σ^{2.5} and in the ROI of Taurus and Orion A dense regions with σ> 10^{22} cm^{2} exit. Thus the γray emissivities in the χ = 0.01 case should be significantly different from that in the CR uniform case.
These results are shown in the red area in Figs. 3 and 4. Due to the CR exclusion the correlation between the γray map and dust opacity map are no longer perfect, especially at smaller scales, which correspond to dense cores. In this region the CRs cannot penetrate effectively. The absolute value of wavelet response is also significantly smaller, which means a lower CR density inside MC. The rising of wavelet response with the scales in the red shaded area in Fig. 4 also illustrate the severe CRs exclusion in dense cores.
The simulation results show a very good agreement between the data and simulations of the CR freepenetration case (uniform CR density), both in shape and absolute values. we note that the agreement in the absolute value of the the wavelet response reveal the same CR density in the MCs as the LIS. This is consistent with the analysis for the GMCs in the Gould belt in Yang et al. (2014). Thus the wavelet crosscorrelation method can be used to investigate not only the CR penetration problem but also the CR densities in different structures. On the other hand, the simulation results for CR exclusion cases deviate significantly with the data and the simulations for CR freepenetration case, which show that the wavelet method can clearly distinguish these two different cases, and our result is robust.
5. Conclusion and discussion
The results show that the CR can penetrate the MC cores freely and the CRs density is uniform down to the scale of the FermiLATangular resolution, which corresponds to 0.7 pc for Taurus and 2.6 pc for Orion A. The current analysis was done for the γray with energy above 2 GeV, which relate to CRs above 10 GeV. The results reveal that there is no suppression of the diffusion coefficient in GMCs, or propagation inside the GMCs is dominated by advection. The results imply that the assumption of MC as CR barometers is correct. Furthermore, the comparison of the simulation results and data also reveal that the CR density inside the GMCs are the same as the LIS, which is consistent with the results in Yang et al. (2014).
The results show that wavelet cross correlation is a unique and powerful tool to investigate CR propagation. In the current study, this method was applied to nearby passive clouds, with γray energy above 2 GeV. It would be also interesting to extend the analysis to lower energy to investigate ionization loss of CRs in the MCs (Morlino & Gabici 2015). However the poor PSF of FermiLATat lower energies makes the extension difficult. It would also be interesting to investigate the energydependent wavelet cross correlation to check weather there are energydependent effects of CR propagations. This would be extremely important near the young CR accelerators, for example, supernova remnants. However the Poisson noise would also limit the application of the wavelet method. In the current study for Taurus and Orion A the Poisson noise dominates the smallscale structures for counts map above 5 GeV. For the young accelerators which are more distant than the local GMCs the statistic may be even worse. But the much larger collected areas for air Cherenkov telescope arrays (ACTA) such as H.E.S.S., MAGIC, VERITAS and the forthcoming CTA make it possible to perform this study by using the VHE image near young accelerators.
gll_iem_v06.fit and iso_P8R2_SOURCE_V6_v06.txt, available at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html
References
 Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2012a, ApJ, 756, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2012b, ApJ, 755, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F. A. 1991, Ap&SS, 180, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonian, F. A. 2001, Space Sci. Rev., 99, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Arshakian, T. G., & Ossenkopf, V. 2016, A&A, 585, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casandjian, J.M. 2015, ApJ, 806, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Casanova, S., Aharonian, F. A., Fukui, Y., et al. 2010, PASJ, 62, 769 [NASA ADS] [Google Scholar]
 Dalgarno, A. 2006, Proc. Nat. Acad. Sci., 103, 12269 [NASA ADS] [CrossRef] [Google Scholar]
 Dogiel, A. V., & Sharov, S. G. 1990, International Cosmic Ray Conf., 4, 109 [NASA ADS] [Google Scholar]
 Frick, P., Beck, R., Berkhuijsen, E. M., & Patrickeyev, I. 2001, MNRAS, 327, 1145 [NASA ADS] [CrossRef] [Google Scholar]
 Gabici, S., Aharonian, F. A., & Blasi, P. 2007, Ap&SS, 309, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Istomin, Y. N., & Kiselev, A. 2013, MNRAS, 436, 2774 [NASA ADS] [CrossRef] [Google Scholar]
 Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014 [NASA ADS] [CrossRef] [Google Scholar]
 Morlino, G., & Gabici, S. 2015, MNRAS, 451, L100 [NASA ADS] [CrossRef] [Google Scholar]
 Nath, B. B., Gupta, N., & Biermann, P. L. 2012, MNRAS, 425, L86 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2011, A&A, 536, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXVIII. 2015, A&A, 582, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2016, A&A, in press, DOI: 10.1051/00046361/201525844 [Google Scholar]
 Skilling, J., & Strong, A. W. 1976, A&A, 53, 253 [NASA ADS] [Google Scholar]
 Yang, R.z., de Oña Wilhelmi, E., & Aharonian, F. 2014, A&A, 566, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yang, R.z., Jones, D. I., & Aharonian, F. 2015, A&A, 580, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Upper panel: γray counts map above 2 GeV after point source subtraction in the region Orion A (left) and Taurus (right). Lower panel: τ_{353} maps for Orion A (left) and Taurus (right). 

Open with DEXTER  
In the text 
Fig. 2 Wavelet decomposition of τ_{353} map of Taurus with the scale of 0.2, 0.5, 1.0, 1.5 degrees (from top left to bottom right). 

Open with DEXTER  
In the text 
Fig. 3 Wavelet crosscorrelation coefficient of Orion A (left panel) and Taurus (right panel), respectively. The gray shaded areas show the simulation results assuming a uniform CR density at all scales, while the red areas show the results when assuming a smaller diffusion coefficient inside MC. χ is defined in Eq. (7). 

Open with DEXTER  
In the text 
Fig. 4 The wavelet response of Orion A (left panel) and Taurus (right panel), respectively. Gray shaded areas show the simulation results assuming a uniform CR density at all scales, while the red areas show the results when assuming a smaller diffusion coefficient inside MC. 

Open with DEXTER  
In the text 
Fig. 5 Wavelet crosscorrelation coefficient (left panel) and wavelet response (right panel) of 100 realization of null test simulations. 

Open with DEXTER  
In the text 