1 / 50

GALPROP tutorial

GALPROP tutorial. Igor Moskalenko (Stanford U.). 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!

roz
Download Presentation

GALPROP tutorial

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  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! _

More Related