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
same level dirs
FITS
v42.3

• GALPROP
• (c fortran)

v42.3
nuclei (fits) file (RRsun) nuclei_full (fits)
file (whole galaxy) ?-ray emissivities (fits)
files (brems, IC, p0) ?-ray skymaps (fits) files
(rings brems, IC, p0)
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 -ICFITSIO/includ
  e
• 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 _at_ 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
• includeltiostreamgt
• includeltcstdlibgt
• includeltstringgt
• includeltcctypegt
• includeltfstreamgt
• includeltcmathgt

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 maxA64
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, po-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 p0-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
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 (p0-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 b40 -using integrated H I
  distribution

11
Dark Matter in GALPROP

• DM annihilation products
• ?0?0-gt 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 NetworkCross Sections
Secondary, radioactive 1 Myr K-capture isotopes
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)
sources (SNR, nuclear reactions)
diffusion
convection
diffusive reacceleration
convection
E-loss
radioactive decay
fragmentation

• ?(r,p,t) density per total momentum

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 p0-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 -gt 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
Distributed Stochastic Reacceleration
Simon et al. 1986 Seo Ptuskin 1994
Scattering on magnetic turbulences
Dpp p2Va2/D D vR1/3 - Kolmogorov spectrum
Icr
1/3
?E
strong reacceleration

• Fermi 2-nd order mechanism

Dxx 5.2x1028 (R/3 GV)1/3cm-2 s-1 Va 36 km
s-1 ? R-d, d1.8/2.4 below/above 4 GV
weak reacceleration
E
32
Convection
Galactic wind
Jones 1979
DR0.6

• Escape length

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-d, d2.46/2.16 below/above 20 GV
33
Damping of Interstellar Turbulence
Kolmogorov cascade
Iroshnikov-Kraichnan cascade
nonlinear cascade
Mean free path
W(k)

• Simplified case
• 1-D diffusion
• No energy losses

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 pp-gtp,K-gte (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-gt pbar scaling using Gaisser Schaeffer
  1992 or Simon et al. 1998 results similar
• Total inelastic x-section (p TN83, A Moiseev
  Ormes 1997)
• ppbar annihilation x-section (ppbar)tot
  (pp)tot (LE TN83, HE PDG00 Regge
  parameterization)

36
Elemental Abundances CR vs. Solar System
Volatility
CR abundances ACE
output
O
Si
Na
Fe
S
CNO
Al
CrMn
Cl
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 cant be done without modern nuclear codes.
• Use LANL nuclear database and modern computer
  codes.

40
LiBeB Major Production Channels

• Well defined (65)
• C12, O16 -gtLiBeB
• N14 -gt Be7
• (see Moskalenko Mashnik 28 ICRC, 2003)
• Few measurements
• C13,N -gt LiBeB
• B -gt BeB
• Unknown
• LiBeB,C13,N -gt LiBeB
• Tertiary reactions also important! -35

Propagated Abundance Cross-section
Li6
O
C
16
12
Li
B
N
Be
7
13
A
10
9
15
14
11
41
Effect of Cross Sections Radioactive Secondaries
Different size from different ratios
T1/2?
Zhalo,kpc

• Errors in CR measurements (HE LE)
• Errors in production cross sections
• Errors in the lifetime estimates
• Different origin of elements (Local Bubble ?)

Ek, MeV/nucleon
42
How It Is Really Done Nucleons

• Calculated for pA reactions and scaled for aA
  (Ferrando et al. 1988)
• Calculation of total nuclear cross sections
• Letaw et al. 1983
• Wellisch Axen 1996 (corrected), Zgt5
• Barashenkov Polanski 2001
• Calculation of isotopic production cross sections
• Webber et al. 1993 (non-renormalized,
  renormalized) Egt200 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
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 Digel05
46
Gas Distribution
Molecular hydrogen H2 is traced using J1-0
transition of 12CO, concentrated mostly in the
plane (z70 pc, Rlt10 kpc) Atomic
hydrogen H I (radio 21 cm) has a wider
distribution (z1 kpc, R30 kpc) Ionized
hydrogen H II (visible, UV, X) small proportion,
but exists even in halo (z1 kpc)
Z0,0.1,0.2 kpc
Sun
47
Gas Rings HI (Inner Outer Galaxy)
Seth Digel05
48
Gas Rings HI (Our Neighborhood)
Seth Digel05
49
How It Is Really Done Gammas

• Bremsstrahlung (Koch Motz 1959, SMR2000) many
  different regimes
• LE 0.01 lt Ekin lt 0.07 MeV nonrelativistic
  non-screened brems
• Intermediate 0.07 lt Ekin lt2 MeV
• HE Ekin gt 2 MeV arbitrary screening unshielded
  charge, 1-, 2-electron atoms (form factors,
  Hylleraas, Hertree-Fock wave functions)
• Fano-Sauter limit k-gtEkin
• 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/Z1 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!

_