GALPROP tutorial

1. GALPROP tutorial Igor Moskalenko (Stanford U.)

2. Obtaining GALPROP The link: http://www.mpe.mpg.de/~aws/dlp/zxc/kty/v42.3p/ Contains c++, fortran routines & input files (dat-files, compass gas maps, and isrf) Dedicated GALPROP Web-site will go online soon! • Controlled changes in GALPROP: tests + documentation +… • New version(s) + archive versions • Post relevant information: best models, gas maps, ISRF, nuclear cross sections… • Allow for communication with users • Ability to run GALPROP on-line…

3. I/O galdef-file gas: COR, HIR, IAQ-files RFComposite –(fits) file dat-files (xsec, nuc.network) GALDEF GALPROP (c++ & fortran) same level dirs FITS v42.3 v42.3 nuclei –(fits) file (R=Rsun) nuclei_full –(fits) file (whole galaxy) γ-ray emissivities –(fits) files (brems, IC, π0) γ-ray skymaps –(fits) files (rings; brems, IC, π0) FITS Heliospheric modulation: on-the-fly in a plotting routine

4. Example Makefile #CXX = g++-2.95 #FC = g77-2.95 CFITSIO = ${GLAST_EXT}/cfitsio/v2470 CPPFLAGS = -O3 -Wno-deprecated -I${CFITSIO}/include FFLAGS = LIBS = -lm -lg2c FITSLIB = -L$(CFITSIO)/lib -lcfitsio -Wl, rpath,$(CFITSIO)/lib LDFLAGS = $(FITSLIB) $(LIBS) FOBJS := $(patsubst %.f,%.o,$(wildcard *.f)) CCOBJS := $(patsubst %.cc,%.o,$(wildcard *.cc)) galprop: ${FOBJS} ${CCOBJS} $(CXX) *.o -o $@ ${LDFLAGS}

5. Some Editing Tested for g++ v2.9x compiler. New g++ compiler v3.x is more strict –routines require some editing: using namespace std; #include<iostream> #include<cstdlib> #include<string> #include<cctype> #include<fstream> #include<cmath>

6. GALPROP Input: galdef-files GALPROP is parameter-driven (user can specify everything!) Grids • 2D/3D –options; symmetry options (full 3D, 1/8 -quadrants) • Spatial, energy/momentum, latitude & longitude grids • Ranges: energy, R, x, y, z, latitude & longitude • Time steps Propagation parameters • Dxx, VA, VC & injection spectra (p,e) • X-factors (including R-dependence) Sources • Parameterized distributions • Known SNRs • Random SNRs (with given/random spectra), time dependent eq. Other • Source isotopic abundances, secondary particles (pbar , e±, γ, synchro), anisotropic IC, energy losses, nuclear production cross sections…

7. Algorithm primary source functions (p, He, C .... Ni) source abundances, spectra primary propagation -starting from maxA=64 source functions (Be, B...., e+,e-, pbars) using primaries and gas distributions secondary propagation tertiary source functions tertiary propagation (i) CR –fixing propagation -rays (IC, bremsstrahlung, πo-decay) radio: synchrotron (ii) γ-rays

8. GALPROP Output/FITS files Provides literally everything: • All nuclei and particle spectra in every grid point (x,y,R,z,E) -FITS files Separately for π0-decay, IC, bremsstrahlung: • Emissivities in every grid point (x,y,R,z,E,process) • Skymaps with a given resolution (l,b,E,process) • Output of maps separated into HI, H2, and rings to allow fitting X, metallicity gradient etc.

9. Z 1 6 1 5 0 2 4 7 R,x,y 3 -1 Spatial Grids Typical grid steps (can be arbitrary!) Δz = 0.1 kpc, ΔΔz = 0.01 kpc (gas averaging) ΔR = 1 kpc ΔE = x1.2 (log-grid)

10. GALPROP Calculations Constraints • Bin size (x,y,z) depends on the computer speed, RAM; final run can be done on a very fine grid ! • No other constraints ! –any required process/formalism can be implemented Calculations (γ -ray related) • Vectorization options • Now 64 bit to allow unlimited arrays • Heliospheric modulation: routinely force-field, more sophisticated model ? • For a given propagation parameters: propagate p, e, nuclei, secondaries (currently in 2D) • The propagated distributions are stored • With propagated spectra: calculate the emissivities (π0-decay, IC, bremss) in every grid point • Integrate the emissivities over the line of sight: • GALPROP has a full 3D grid, but currently only 2D gas maps (H2, H I, H II) • Using actual annular maps (column density) at the final step • High latitudes above b=40˚ -using integrated H I distribution

11. Dark Matter in GALPROP DM annihilation products: χ0χ0−> p, pbar, e+, e−, γ A set of routines (gen_DM_source.cc) to assign • The DM density profile (NFW, isothermal etc.) • Source functions for p, pbar, e+, e− • Source function for γ’s • A set of user-defined parameters (10 int, 10 double precision) in galdef-file • DM annihilation products particles are propagated in the same model as CR particles. • Calculation of skymaps for DMγ-rays

12. galdef_44_599278 -I

13. galdef_44_599278 -II

14. galdef_44_599278 -III

15. galdef_44_599278 -IV

16. Nuclear Reaction Network+Cross Sections Secondary, radioactive ~1 Myr & K-captureisotopes Co57 Fe55 Mn54 Cr51 V49 Ca41 Ar37 Cl36 β-, n Al26 p,EC,β+ Be7 Be10 Plus some dozens of more complicated reactions. But many cross sections are not well known…

17. Nuclear Reaction Network I II III IV V nuc_package.cc

18. nuc_package.cc: Stable & Long-lived Isotopes

19. nuc_package.cc: Long-Lived Isotopes

20. nuc_package.cc: “Boundary” Nuclei

21. isotope_cs_eval.dat

22. Transport Equations ~90 (no. of CR species) ψ(r,p,t) – density per total momentum sources (SNR, nuclear reactions…) diffusion convection diffusive reacceleration convection E-loss radioactive decay fragmentation

23. Finite Differencing

24. Finite Differencing: Example

25. Tri-Diagonal Matrix

26. Coefficients for the Crank-Nicholson Method

27. Near Future Developments Full 3D Galactic structure: • 3D gas maps (from S.Digel, S.Hunter and/or smbd else) • 3D interstellar radiation & magnetic fields (A.Strong & T.Porter) Cross sections: • Blattnig et al. formalism for π0-production • Diffractive dissociation with scaling violation (T.Kamae –param.) • Isotopic cross sections (with S.Mashnik, LANL; try to motivate BNL, JENDL-Japan, other Nuc. Data Centers) Modeling the local structure: • Local SNRs with known positions and ages • Local Bubble, local clouds –may be done at the final calculation step (grid bin size ??) Energy range: • Extend toward sub-MeV range to compare with INTEGRAL diffuse emission (continuum; 511 keV line) Heliospheric modulation: • Implementing a modern formalism (Potgieter, Zank etc.) Visualization tool (started) using the classes of CERN ROOT package: images, profiles, and spectra from GALPROP to be directly compared with data Improving the GALPROP module structure (for DM studies)

28. More developments • Point sources: develop algorithm(s) for modeling the background and interface to the rest of GLAST software • Instrumental response: how to implement • Diffuse emission analysis has to include point source catalog! • At least, two diffuse models: with/without the “excess” • Develop test case(s) to test the accuracy of the numerical model (simple gas distribution, no energy losses, uniform ISRF etc.) • Complete C++ package: rewrite several fortran routines in C++ • Develop a fitting procedure to make automatic fitting to B/C ratio, CR spectra and abundances • Develop a dedicated Web-site: • Controlled changes in GALPROP: tests +documentation +… • Allow for communication with users • Post relevant information: best models, gas maps, ISRF, nuclear cross sections… • Ability to run GALPROP on-line…

29. Fixing Propagation Parameters: Standard Way • Using secondary/primary nuclei ratio: • Diffusion coefficient and its index • Propagation mode and its parameters (e.g., reacceleration VA, convection Vz) B/C Interstellar Be10/Be9 Ek, MeV/nucleon Radioactive isotopes: Galactic halo size Zh Zh increase Ek, MeV/nucleon

30. Peak in the Secondary/Primary Ratio • Leaky-box model: fitting path-length distribution -> free function • Diffusion models: • Diffusive reacceleration • Convection • Damping of interstellar turbulence • Etc. B/C Ek, MeV/nucleon too sharp max? Accurate measurements in a wide energy range may help to distinguish between the models

31. B Distributed Stochastic Reacceleration Simon et al. 1986 Seo & Ptuskin 1994 Scattering on magnetic turbulences Fermi 2-nd order mechanism Dpp~ p2Va2/D D ~ vR1/3 - Kolmogorov spectrum Icr 1/3 ΔE strong reacceleration Dxx = 5.2x1028 (R/3 GV)1/3cm-2 s-1 Va = 36 km s-1 γ ~ R-δ, δ=1.8/2.4 below/above 4 GV weak reacceleration E

32. Convection Escape length Galactic wind Jones 1979 D~R0.6 Xe v R-0.6 wind or turbulent diffusion resonant diffusion E problem:too broad sec/prim peak Dxx = 2.5x1028 (R/4 GV)0.6cm-2 s-1 dV/dz = 10 km s-1 kpc-1 γ ~ R-δ, δ=2.46/2.16 below/above 20 GV

33. Damping of Interstellar Turbulence Kolmogorov cascade: Simplified case: • 1-D diffusion • No energy losses Iroshnikov-Kraichnan cascade: nonlinear cascade Mean free path W(k) dissipation k 1/1020cm 1/1012cm Ptuskin et al. 2003, 2005

34. Dxx – Diffusion Coefficient ~R0.6 Reacceleration with damping Plain diffusion ~β-3 Diffusive reacceleration

35. How It Is Really Done: Secondary Particles • Positrons/electrons p+p->π,K->e± (MS 1998) • Dermer 1986 method: LE –Stecker Δ-isobar model (isotropic decay), HE –scaling (inv. x-section: Stephens & Badhwar 1981), plus interpolation in between • Pion decay “includes” polarization of muons • Kaon decay – scaling (Stephens & Badhwar 1981) • Antiprotons (M et al. 2002) • pp Inclusive production x-section (Tan & Ng 1983) • pA, AA-> pbar –scaling using Gaisser & Schaeffer 1992 or Simon et al. 1998 –results similar • Total inelastic x-section (p: TN’83, A: Moiseev & Ormes 1997) • p+pbar annihilation x-section = (p+pbar)tot – (p+p)tot (LE: TN’83, HE: PDG’00 –Regge parameterization)

36. Volatility Elemental Abundances: CR vs. Solar System CR abundances: ACE “output” O Si Na Fe S CNO Al Cl CrMn LiBeB F ScTiV Solar system abundances “input”

37. Fitting to Measured CR Abundances (ACE vs HEAO-3) Fitting to measured CR abundances in the wide energy range (~0.1 – 30 GeV) is problematic: May indicate systematic or cross-calibration errors

38. Total Nuclear Cross Sections Ekin, MeV/nucleon Wellisch & Axen 1996

39. Isotopic Production Cross Sections of LiBeB Semi-empirical systematics are not always correct. Results obtained by different groups are often inconsistent and hard to test. Very limited number of nuclear measurements: Evaluating the cross section is very laborious and can’t be done without modern nuclear codes. Use LANL nuclear database and modern computer codes.

40. LiBeB: Major Production Channels Propagated Abundance * Cross-section • Well defined (65%): C12, O16 ->LiBeB N14 -> Be7 (see Moskalenko & Mashnik 28 ICRC, 2003) • Few measurements: C13,N -> LiBeB B -> BeB • Unknown: LiBeB,C13,N -> LiBeB • “Tertiary” reactions also important! -35% Li6 O C 16 12 Li B N Be 13 7 A= 10 9 15 14 11

41. natSi+p26Al ST W 27Al+p26Al W ST Effect of Cross Sections: Radioactive Secondaries Different size from different ratios… • Errors in CR measurements (HE & LE) • Errors in production cross sections • Errors in the lifetime estimates • Different origin of elements (Local Bubble ?) T1/2=? Zhalo,kpc Ek, MeV/nucleon

42. How It Is Really Done: Nucleons • Calculated for p+Areactions and scaled for α+A (Ferrando et al. 1988) • Calculation of total nuclear cross sections • Letaw et al. 1983 • Wellisch & Axen 1996 (corrected), Z>5 • Barashenkov & Polanski 2001 • Calculation of isotopic production cross sections • Webber et al. 1993 (non-renormalized, renormalized); E>200 MeV/nucleon, essentially flat • Silberberg & Tsao 2000 (non-renormalized, renormalized); claim that it works at all energies, but is problematic sometimes • Fits to the available data (LANL, Webber et al., etc.) in the form of a function or a table (see .dat files), but data may not be always available • Use the best of all three, but very time consuming work • Nuclear reaction network • Nuclear Data Tables (includes several decay channels + branching) • Standard β± -decay, emission of p, n • K-capture isotopes can be treated separately

43. 20 0 -20 -40 Galactic Latitude 220 200 180 160 140 120 100 80 60 40 20 0 Galactic Longitude CfA 1.2m Interstellar Gas • Extend CO surveys to high latitudes • newly-found small molecular clouds will otherwise be interpreted as unidentified sources, and clearly limit dark matter studies • C18O observations (optically thin tracer) of special directions (e.g. Galactic center, arm tangents) • assess whether velocity crowding is affecting calculations of molecular column density, and for carefully pinning down the diffuse emission Dame, Hartmann, & Thaddeus (2001) Dame & Thaddeus (2004)

44. Distribution of Interstellar Gas • Neutral interstellar medium • 21-cm H I & 2.6-mm CO (stand in for H2) • Near-far distance ambiguity, rotation curve, optical depth effects, … (25°, 0°) CO 25° Dame et al. (1987) G.C. H I Hartmann & Burton (1997) W. Keel Seth Digel

45. Seth Digel’05

46. Gas Distribution Molecular hydrogenH2is traced using J=1-0 transition of 12CO,concentrated mostly in the plane (z~70 pc, R<10 kpc) Atomic hydrogen H I(radio 21 cm)has a wider distribution (z~1 kpc, R~30 kpc) Ionized hydrogen H II(visible, UV, X)small proportion, but exists even in halo (z~1 kpc) Z=0,0.1,0.2 kpc Sun

47. Gas Rings: HI(Inner & Outer Galaxy) Seth Digel’05

48. Gas Rings: HI (Our Neighborhood) Seth Digel’05

49. How It Is Really Done: Gammas • Bremsstrahlung (Koch & Motz 1959, SMR2000) –many different regimes: • LE: 0.01 < Ekin < 0.07 MeV –nonrelativistic non-screened brems • Intermediate: 0.07 < Ekin <2 MeV • HE: Ekin > 2 MeV arbitrary screening: unshielded charge, 1-, 2-electron atoms (form factors, Hylleraas, Hertree-Fock wave functions) • Fano-Sauter limit: k->Ekin • “Anisotropic” IC (MS2000) • Takes into account the anisotropic angular distribution of background photons • Neutral pion decay (see “secondary positrons/electrons”) • Synchrotron radiation (Ginzburg 1979, Ghisellini et al. 1988) • Averaging over pitch angle • Uses total magnetic field (regular + random) • Emissivities: uses “real” H2, H I gas column densities (“rings”) • Skymap calculations: integration over the line of sight

50. Uncertainties in the Propagation Models • Production cross sections of isotopes and pbars • Typical errors ~50% • Fitting to B/C ratio may introduce errors in Dxx • Propagation models and parameters • Gas distribution in the Galaxy • Ambient spectrum of CR (solar modulation, GeV excess in γ’s !) • Current knowledge of CR diffusion • Heliospheric modulation • Depends on rigidity • Similar for all nuclei A/Z~ +2 • Different for protons A/Z=+1 and pbars A/Z=-1 • Systematic errors of measurements • Difficult to account for Simultaneous measurements of Li, Be, B, C, secondary e+, p in a wide energy range ~100 MeV/nucleon – 100 GeV/nucleon are needed to understand CR propagation and distinguish between the models: looking forward to Pamela launch! _