Transcript GALPROP tutorial - Stanford University

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!
– 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…
Igor V. Moskalenko 2
November 17, 2005
GALPROP manual/Rome2, Italy
I/O
same
level
dirs
GALDEF
FITS
v42.3
galdef-file
gas: COR, HIR, IAQ-files
RFComposite –(fits) file
dat-files (xsec, nuc.network)
GALPROP
(c++ & fortran)
v42.3
FITS
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)
Heliospheric modulation: on-the-fly in a plotting routine
Igor V. Moskalenko 3
November 17, 2005
GALPROP manual/Rome2, Italy
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}
Igor V. Moskalenko 4
November 17, 2005
GALPROP manual/Rome2, Italy
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>
Igor V. Moskalenko 5
November 17, 2005
GALPROP manual/Rome2, Italy
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…
Igor V. Moskalenko 6
November 17, 2005
GALPROP manual/Rome2, Italy
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
(i) CR –fixing propagation
(ii) γ-rays
Igor V. Moskalenko 7
tertiary source functions
tertiary propagation
-rays (IC, bremsstrahlung, πo-decay)
radio: synchrotron
November 17, 2005
GALPROP manual/Rome2, Italy
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.
Igor V. Moskalenko 8
November 17, 2005
GALPROP manual/Rome2, Italy
Spatial Grids
Z
1
0
1
2
3
4
5
6
7
R,x,y
-1
Typical grid steps (can be arbitrary!)
Δz = 0.1 kpc, ΔΔz = 0.01 kpc (gas averaging)
ΔR = 1 kpc
ΔE = x1.2 (log-grid)
Igor V. Moskalenko 9
November 17, 2005
GALPROP manual/Rome2, Italy
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 ?
1. For a given propagation parameters: propagate p, e, nuclei, secondaries
(currently in 2D)
2. The propagated distributions are stored
3. With propagated spectra: calculate the emissivities (π0-decay, IC, bremss)
in every grid point
4. 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
Igor V. Moskalenko 10
November 17, 2005
GALPROP manual/Rome2, Italy
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
Igor V. Moskalenko 11
November 17, 2005
GALPROP manual/Rome2, Italy
galdef_44_599278 -I
Igor V. Moskalenko 12
November 17, 2005
GALPROP manual/Rome2, Italy
galdef_44_599278 -II
Igor V. Moskalenko 13
November 17, 2005
GALPROP manual/Rome2, Italy
galdef_44_599278 -III
Igor V. Moskalenko 14
November 17, 2005
GALPROP manual/Rome2, Italy
galdef_44_599278 -IV
Igor V. Moskalenko 15
November 17, 2005
GALPROP manual/Rome2, Italy
Nuclear Reaction Network+Cross Sections
Secondary,
radioactive ~1 Myr
& K-capture isotopes
55
V
Fe
Mn54
Cr51
49
Co57
Ca41
Ar37
36
Cl
26
Al
β-, n
p,EC,β+
Be7 Be10
Plus some dozens of more complicated reactions.
But many cross sections are not well known…
Igor V. Moskalenko 16
November 17, 2005
GALPROP manual/Rome2, Italy
Nuclear Reaction Network
I
II
III
IV
V
nuc_package.cc
Igor V. Moskalenko 17
November 17, 2005
GALPROP manual/Rome2, Italy
nuc_package.cc: Stable & Long-lived Isotopes
Igor V. Moskalenko 18
November 17, 2005
GALPROP manual/Rome2, Italy
nuc_package.cc: Long-Lived Isotopes
Igor V. Moskalenko 19
November 17, 2005
GALPROP manual/Rome2, Italy
nuc_package.cc: “Boundary” Nuclei
Igor V. Moskalenko 20
November 17, 2005
GALPROP manual/Rome2, Italy
isotope_cs_eval.dat
Igor V. Moskalenko 21
November 17, 2005
GALPROP manual/Rome2, Italy
Transport Equations ~90 (no. of CR species)


 ( r , p, t )
 q( r , p ) sources (SNR, nuclear reactions…)
t
diffusion

   [D
xx


  V  ]
convection
  2
  
p D

diffusive reacceleration 
pp p p 2 
p 


  dp
1   
  p  V  convection
E-loss 

p  dt
3

fragmentation 





f
radioactive decay
d
ψ(r,p,t) – density per total momentum
Igor V. Moskalenko 22
November 17, 2005
GALPROP manual/Rome2, Italy
Finite Differencing
Igor V. Moskalenko 23
November 17, 2005
GALPROP manual/Rome2, Italy
Finite Differencing: Example
Igor V. Moskalenko 24
November 17, 2005
GALPROP manual/Rome2, Italy
Tri-Diagonal Matrix
Igor V. Moskalenko 25
November 17, 2005
GALPROP manual/Rome2, Italy
Coefficients for the Crank-Nicholson Method
Igor V. Moskalenko 26
November 17, 2005
GALPROP manual/Rome2, Italy
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)
Igor V. Moskalenko 27
November 17, 2005
GALPROP manual/Rome2, Italy
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…
Igor V. Moskalenko 28
November 17, 2005
GALPROP manual/Rome2, Italy
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
Be10/Be9
Ek, MeV/nucleon
Radioactive isotopes:
Galactic halo size Zh
Zh increase
Ek, MeV/nucleon
Igor V. Moskalenko 29
November 17, 2005
GALPROP manual/Rome2, Italy
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
Igor V. Moskalenko 30
November 17, 2005
GALPROP manual/Rome2, Italy
Distributed Stochastic Reacceleration
Simon et al. 1986
Seo & Ptuskin 1994
Scattering on
magnetic turbulences
Dpp~ p2Va2/D
D ~ vR1/3 - Kolmogorov spectrum
Icr
B
1/3
ΔE
Fermi 2-nd order mechanism
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
Igor V. Moskalenko 31
November 17, 2005
strong
reacceleration
weak
reacceleration
E
GALPROP manual/Rome2, Italy
Convection
Galactic wind
Jones 1979
Escape length
Xe
v
wind or
turbulent
diffusion
D~R0.6
R-0.6
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
Igor V. Moskalenko 32
November 17, 2005
GALPROP manual/Rome2, Italy
Damping of Interstellar Turbulence
Kolmogorov cascade:
Iroshnikov-Kraichnan cascade:
nonlinear
cascade
W(k)
1/1012cm
•
•
l(p)
Simplified case:
dissipation
1/1020cm
Mean free path
1-D diffusion
No energy losses
k
p
Ptuskin et al. 2003, 2005
Igor V. Moskalenko 33
November 17, 2005
GALPROP manual/Rome2, Italy
Dxx – Diffusion Coefficient
~β-3
Plain
diffusion
Reacceleration
with damping
~R0.6
Diffusive
reacceleration
Igor V. Moskalenko 34
November 17, 2005
GALPROP manual/Rome2, Italy
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)
Igor V. Moskalenko 35
November 17, 2005
GALPROP manual/Rome2, Italy
Elemental Abundances: CR vs. Solar System
CR abundances: ACE
“output”
O
Si
Na
Fe
S
CNO
Al
Cl
LiBeB
CrMn
F
ScTiV
Solar system abundances
“input”
Igor V. Moskalenko 36
November 17, 2005
GALPROP manual/Rome2, Italy
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
Igor V. Moskalenko 37
November 17, 2005
GALPROP manual/Rome2, Italy
Total Nuclear Cross Sections
Ekin, MeV/nucleon
Wellisch & Axen 1996
Igor V. Moskalenko 38
November 17, 2005
GALPROP manual/Rome2, Italy
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.
Igor V. Moskalenko 39
November 17, 2005
GALPROP manual/Rome2, Italy
LiBeB: Major Production Channels
Propagated Abundance * Cross-section
6
Li
O
C
Li
Be
B
N
A=
•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%
Igor V. Moskalenko 40
November 17, 2005
GALPROP manual/Rome2, Italy
Effect of Cross Sections: Radioactive Secondaries
Different size from different ratios…
27Al+p26Al
T1/2=?
W
ST
ST
natSi+p26Al
W
Zhalo,kp
c
Ek, MeV/nucleon
Igor V. Moskalenko 41
• Errors in CR measurements (HE & LE)
• Errors in production cross sections
• Errors in the lifetime estimates
• Different origin of elements (Local Bubble ?)
November 17, 2005
GALPROP manual/Rome2, Italy
How It Is Really Done: Nucleons
•
•
Calculated for p+A reactions 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
Igor V. Moskalenko 42
November 17, 2005
GALPROP manual/Rome2, Italy
Interstellar Gas
►
Extend CO surveys to high latitudes
– newly-found small molecular clouds will otherwise be interpreted as
20
0
-20
-40
220
200
180
160
140
120
100
80
Galactic Longitude
►
60
40
20
0
CfA 1.2m
Dame, Hartmann, & Thaddeus (2001)
Dame & Thaddeus (2004)
Galactic Latitude
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
Igor V. Moskalenko 43
November 17, 2005
GALPROP manual/Rome2, Italy
Distribution of Interstellar Gas
• Neutral interstellar medium
– 21-cm H I & 2.6-mm CO (stand in for H2)
(25°, 0°)
CO
Dame et al.
(1987)
G.C.
HI
Hartmann &
Burton (1997)
W. Keel
• Near-far distance ambiguity, rotation curve, optical depth
effects, …
Seth Digel
Igor V. Moskalenko 44
November 17, 2005
GALPROP manual/Rome2, Italy
Seth Digel’05
Igor V. Moskalenko 45
November 17, 2005
GALPROP manual/Rome2, Italy
Gas Distribution
Z=0,0.1,0.2 kpc
Molecular hydrogen H2
is 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)
Sun
Igor V. Moskalenko 46
November 17, 2005
Ionized hydrogen H II
(visible, UV, X) small
proportion, but exists
even in halo (z~1 kpc)
GALPROP manual/Rome2, Italy
Gas Rings: HI (Inner & Outer Galaxy)
Seth Digel’05
Igor V. Moskalenko 47
November 17, 2005
GALPROP manual/Rome2, Italy
Gas Rings: HI (Our Neighborhood)
Seth Digel’05
Igor V. Moskalenko 48
November 17, 2005
GALPROP manual/Rome2, Italy
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-, 2electron 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
Igor V. Moskalenko 49
November 17, 2005
GALPROP manual/Rome2, Italy
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
e +,
_
Simultaneous measurements of Li, Be, B, C, secondary
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!
Igor V. Moskalenko 50
November 17, 2005
GALPROP manual/Rome2, Italy