The SOPHIA Monte Carlo Code
Monte Carlo simulations of photohadronic processes in astrophysics
- Anita Mücke (now: Reimer)
- Ralph Engel
- Jorg P. Rachen
- Raymond J. Protheroe
- Todor Stanev
Computer Physics Communications, 124, 290-314
The cosmic ray spectrum extends to extremely high energies. Giant air showers have been observed with energy exceeding ~1011 GeV. Energy losses due to interactions with ambient photons can become important, even dominant for such energetic nucleons, above the threshold for pion production. Photoproduction of hadrons is expected to cause a distortion of the ultra-high energy cosmic ray (CR) spectrum by interactions of the nucleons with the microwave background (the Greisen-Zatsepin-Kuzmin cutoff), but it may also be relevant to the observed high energy gamma ray emission from jets of Active Galactic Nuclei (AGN) or Gamma-Ray Bursts (GRB). Moreover, it is the major source process for the predicted fluxes of very high energy cosmic neutrinos.
The photohadronic cross section at low interaction energies is dominated by the Delta(1232) resonance. Since the low energy region of the cross section is emphasized in many astrophysical applications, the cross section and decay properties of the prominent Delta-resonance have often been used as an approximation for photopion production, and the subsequent production of gamma rays and neutrinos. As discussed in e.g. Mücke et al. 1999, this approximation is only valid for a restricted number of cases, and does not describe sufficiently well the whole energy range of photohadronic interactions. A more sophisticated photoproduction simulation code is needed to cover the center-of-mass energy range of about s1/2~1 - 103 GeV, which is important in many astrophysical applications.
This was the motivation for developing the Monte-Carlo event generator SOPHIA (Simulations Of Photo Hadronic Interactions in Astrophysics), which we wrote as a tool for solving problems connected to photohadronic processes in astrophysical environments, but can also be used for radiation and background studies at high energy colliders such as LEP2 and HERA, as well as for simulations of photon induced air showers. The philosophy of the development of SOPHIA has been to implement well established phenomenological models, symmetries of hadronic interactions in a way that describes correctly the available exclusive and inclusive photohadronic cross section data obtained at fixed target and collider experiments.
We plan to continue investigating a variety of astrophysics topics related to photomeson production using the SOPHIA code. Past studies by members of the SOPHIA collaboration include:
- hadronic blazar emission models (Synchrotron-Proton blazar model)
- cosmic ray propagation
- ...
How to download the SOPHIA package
The original SOPHIA code (version 2.00) can be downloaded from the CPC program library.
The latest version (2.01) may be derived upon request from A. Reimer anita.reimer@uibk.ac.at.
Please cite the following reference in all publications which make use of the SOPHIA code:
A. Mücke, Ralph Engel, J.P. Rachen, R.J. Protheroe, and Todor Stanev, 2000, Comp. Phys. Commun., 124, 290.
If you wish to get notified in case of any code changes, bug fixes, etc. please send an email to anita.reimer@uibk.ac.at
General Information
SOPHIA is designed to simulate minimum bias photo-production by interactions of unpolarized nucleons with unpolarized photons from the particle production threshold through the resonance region up to s1/2~103 in the multiparticle production region. The implemented interaction processes include: incoherent interaction of photons with the virtual structure of the nucleon at low energies near threshold; resonance excitation/decay (9 resonances with masses 1232-1950 MeV); diffractive scattering (vector meson production) and multipion production on the basis of QCD string fragmentation at high energies s1/2>2GeV. A modified JETSET version 7.4 is used for the latter high energy event generation. SOPHIA simulates individual events for given nucleon and photon energies where photon energies are sampled from various distributions. Corresponding incident ambient photon fields are limited to power law and blackbody spectra in the public code version. A detailed description of all processes and the SOPHIA code can be found in A. Mücke, Ralph Engel, J.P. Rachen, R.J. Protheroe, and Todor Stanev, 2000, Comp. Phys. Commun., 124, 290.
Program structure
Coding:
SOPHIA is coded in FORTRAN 77 (operating systems : UNIX, Linux, Open-VMS) and has been tested thoroughly on DEC-Alpha and Intel-Pentium based workstations. No tests were done for center-of-mass energies s1/2>1000GeV. Memory required to execute is <1 megabyte. 10000 events at a center-of-mass energy of 1.5 GeV require a typical CPU time of about 75 seconds. Other programs used in SOPHIA in a modified form are: Rndm (a processor independent random number generator based on Marsaglia & Zaman 1987, preprint FSU-SCRI-87-70), Jetset 7.4 (a Lund Monte Carlo for jet fragmentation; Sjostrand 1994, Comp.Phys.Commun., 82, 74), DECSIB (Sibyll-routine which decays unstable particles; Fletcher et al., 1994, Phys.Rev.D 50, 5710).
The SOPHIA source code consists of several files which contain a number of routines. SOPHIA20.f main program containing the routines which organize the various tasks to be performed. Furthermore, the input is handled here and some kinematic transformations needed as input to several routines are performed. initial.finitialization routine for parameter settings. sampling.f collection of routines/functions needed for sampling the CMF energy squared s and the photon energy eps in the observer frame. eventgen.f event generator for photomeson production in proton-photon and neutron-photon collisions. This is the heart of SOPHIA.
The output of the final states is organized in the common block S_PLIST /P(2000,5), LLIST(2000), NP, Ideb/the array P(i,j) contains the 4-momenta and rest mass of the final state particle in cartesian coordinates (P(i,1) = P_x, P(i,2) = P_y, P(i,3) = P_z, P(i,4) = energy, P(i,5) = rest mass).LLIST()gives the code numbers of all final state particles and NPis the number of stable final state particles. Note that before the initial call of eventgen the initialization routine has to be called. output.fcontains output routines.
Details on each routine can be found in A. Mücke et al., 2000, Comp. Phys. Commun., 124, 290.
The particle identity is given by LLIST with the following numbering scheme (identical to the DECSIB particle identity list):
| code | particle | mass |
|---|---|---|
| 1 | gam | 0.0000 |
| 2 | e+ | 0.0005 |
| 3 | e- | 0.0005 |
| 4 | mu+ | 0.1057 |
| 5 | mu- | 0.1057 |
| 6 | pi0 | 0.1350 |
| 7 | pi+ | 0.1396 |
| 8 | pi- | 0.1396 |
| 9 | k+ | 0.4937 |
| 10 | k- | 0.4937 |
| 11 | k0l | 0.4977 |
| 12 | k0s | 0.4977 |
| 13 | p | 0.9383 |
| 14 | n | 0.9396 |
| 15 | nue | 0.0000 |
| 16 | nueb | 0.0000 |
| 17 | num | 0.0000 |
| 18 | numb | 0.0000 |
| 21 | k0 | 0.4977 |
| 22 | k0b | 0.4977 |
| 23 | eta | 0.5488 |
| 24 | etap | 0.9576 |
| 25 | rho+ | 0.7714 |
| 26 | rho- | 0.7714 |
| 27 | rho0 | 0.7717 |
| 28 | k*+ | 0.8921 |
| 29 | k*- | 0.8921 |
| 30 | k*0 | 0.8965 |
| 31 | k*0b | 0.8965 |
| 32 | omeg | 0.7826 |
| 33 | phi | 1.1020 |
| 34 | SIG+ | 1.1894 |
| 35 | SIG0 | 1.1925 |
| 36 | SIG- | 1.1973 |
| 37 | XI0 | 1.3149 |
| 38 | XI- | 1.3213 |
| 39 | LAM | 1.1156 |
| 40 | DELT++ | 1.2300 |
| 41 | DELT+ | 1.2310 |
| 42 | DELT0 | 1.2320 |
| 43 | DELT- | 1.2330 |
| 44 | SIG*+ | 1.3828 |
| 45 | SIG*0 | 1.3837 |
| 46 | SIG*- | 1.3872 |
| 47 | XI*0 | 1.5318 |
| 48 | XI*- | 1.5350 |
| 49 | OME*- | 1.6724 |
Antibaryons have negative ID numbers correspondingly. Decayed particles are marked by adding 10000 to their ID code.The decay of a particle with code ID can be turned off viaIDB(ID) = -ABS(IDB(ID))and turned on viaIDB(ID) = ABS(IDB(ID))This has to be done only once at the beginning of event generation or might be changed on an event-by-event basis.
Example Programs
The following are example files of input parameter sets and resulting SOPHIA output files:
- input and output for straight power law target photon field
- input and output for broken power law target photon field
- input and output for black body radiation target photon field
History of Updates
Version 2.01 (April 12, 2007)
Last modification on March 27, 2007:
An error in the routine eventgen has been corrected (pointed out by P. Lipari): The problem comes from inconsistent sign conventions for the cos(theta) specified in the parameters. The definition of the explicit particle momenta is only used in the Lorentz transformation. The relevant quantity affected is the boost parameter along the z-axis.
Since the photon momentum is much smaller than the nucleon momentum in "standard" applications in astrophysics, the bug has no effect on simulation results. For a typical UHECR propagation configuration of E_p = 1011GeV, E_ph = 10-10GeV, the relative error would be of the order of 10-20. This applies to isotropic or anisotropic radiation fields.
The situation is different for configurations in which the photon and nucleon energies are of the same order. In this case the total energy would be still fine but the total z-momentum would be wrong.
If you find any bugs in the SOPHIA code, please inform:
Anita Reimer anita.reimer@uibk.ac.at
or
Ralph Engel ralph.engel@ik.fzk.de
Publications where the SOPHIA package has been used:
A. Mücke, Ralph Engel, J.P. Rachen, R.J. Protheroe, and Todor Stanev, 2000, Comp. Phys. Commun., 124, 290:
Monte Carlo simulations of photohadronic processes in astrophysics
A. Mücke & R.J. Protheroe, 2000, AIP, 515, 149:
Modeling the April 1997 Flare of Mkn 501
T. Stanev, R. Engel, A. Mücke, R.J. Protheroe, J.P. Rachen, 2000, Phys. Rev. D, 62, 093005:
Propagation of ultrahigh energy protons in the nearby universe
R.J. Protheroe & A. Mücke, 2001, AIP, 558, 700:
Application of the Synchrotron Proton Blazar Model to BL Lac Objects
R.J. Protheroe & A. Mücke, 2001, ASPC, 250, 113:
Estimating jet power in proton blazar models
A. Mücke & R.J. Protheroe, 2001, ICRC, 3, 1153:
Neutrino Emission from HBLs and LBLs
A. Mücke & R.J. Protheroe,2001, APh, 15, 121:
A proton synchrotron blazar model for flaring in Markarian 501
A. Mücke, R.J. Protheroe, R. Engel, J. Rachen, T. Stanev, 2003, APh, 18, 593:
BL Lac objects in the synchrotron proton blazar model
...