README File for the Code Galaxy =========================================== First version: 12 Aug 1998 (CL) Claus Leitherer, Carmelle Robert, Daniel Schaerer, Jeff Goldader, Rosa Gonzalez-Delgado, & Duilia de Mello Last update: 21 July 2002 (CL) Claus Leitherer & Richard Norris Disclaimer: This code is distributed freely to the community. The user accepts sole responsibility for the results produced by the code. Although every effort has been made to identify and eliminate errors, we accept no responsibility for erroneous model predictions. The code and the models should be referenced as Leitherer et al. (1999, ApJS, 123, 3) 1. A bit of history As so often, this code started out as a small project in 1987. Back then, Claus Leitherer wanted to have a little computer program to calculate ionizing fluxes for individual O stars. The program evolved from handling a single star to a whole population. At the same time, there was the desire to assign positions in the HRD for comparison with observations. Results from this rather simple code were reported in ApJS, 73, 1 (1990). By 1990 the program had grown to compute the output of the mechanical luminosity of winds and of supernovae. A paper on that is in ApJ, 401, 498 (1992). The "real" spectral synthesis capability was implemented in 1991, when we added Kurucz and Schmutz model atmospheres from the far-UV to the near-IR. This allowed us to calculate spectra and colors for young populations. Carmelle Robert included an IUE spectral library in the code in order to compute synthetic lines at 0.75 A resolution (ApJ, 418, 749 [1993], ApJS, 99, 173 [1995]). By 1994, we had computed so many models, synthesizing all kinds of galaxy parameters that the suggestion was made to combine everything, plot it in a homogeneous way, and publish it (ApJS, 96, 9 [1995]). Who would have predicted that quite a few astronomers found this useful? Later we began to add a few more specialized routines, like the computation of O VI line profiles by Rosa Gonzalez-Delgado (ApJ, 489, 601 [1997]). The most recent subroutine was written by Jeff Goldader. This part calculates the strengths of selected stellar features in the near-IR. After having served us so well over 10 years, it become clear that the code needed a major upgrade. Daniel Schaerer replaced the ancient Maeder (1990) tracks by the most recent (1992-94) model series of Geneva. He also implemented the capability to perform isochrone synthesis in addition to classical evolutionary synthesis. One of the subroutines used for this was written by Georges Meynet. We also decided to take advantage of the homogeneous atmosphere grid compiled by Lejeune et al. (1997) and replaced the original Kurucz models by the new set. The updated code was made available to the community via a spiffy website called Starburst99 (Leitherer et al. 1999, ApJ, 123, 3). Since then, quite a few important improvements were made. First, Duilia de Mello added a high-resolution spectral library of B stars in fall 1999. A description is in ApJ, 530, 251 (2000). This library replaced the low-resolution B-star library we used before. (O stars have always been from a high-dispersion library.) Alessandra Aloisi added supernova yields to the code. Before that, only yields from stellar winds were taken into account, and the supernova yields of those elements which undergo no nuclear processing during a supernova explosion. Now the models are up-to-date with respect to the nucleosynthesis in type II supernovae (and only those -- we don't do type I's yet). This change was made in August 2000. Claus Leitherer and summer student Joao Leao Souza added a spectral library of LMC and SMC stars in August 2000 (Leitherer et al. 2001, ApJ, 550, 724). This allows calculations of models for UV spectra with sub-solar metallicities. Only O-stars are used for SMC/LMC metallicity, everything else is solar. This is not perfect -- but it is a starting point. Miguel Cervino, in collaboration with Daniel Schaerer greatly improved several numerical aspects of the code. In particular he wrote new code to handle the calculation of the supernova rate in isochrone synthesis. Previously, all SN related parameters showed rather ugly discontinuities due to numerical instabilities (see Leitherer et al. 1999). This was fixed in October 2000. Anne Pellerin and Claus Leitherer added a high-resolution spectral library of hot stars observed with FUSE. The library stars are in the Galaxy, the LMC, and the SMC. The spectral range is from 1000 to 1183 A, and the structure corresponds to that of the IUE/HST stars used at longer wavelengths. The library is discussed in Robert et al. (2002). A major shortcoming of the original Starburst99 models was removed by Richard Norris, Linda Smith, and Paul Crowther. They replaced the unblanketed WR atmospheres of Schmutz et al (1992) by a fully blanketed set calculated with John Hillier's code. The updated Starburst99 code, released as v4.0 in July 2002, gives more realistic EUV fluxes during WR dominated phaes. More to come in the future --- stay tuned! 2. File structure The package comes with quite a few files. Some are essential whereas others are not, and one can do without them. 2.1 The README file -- you should read this first. 2.2 galaxy.f -- this is the code. It is standard Fortran 77. Nothing fancy, and some software wizards may consider the coding pretty pedestrian. Keep in mind the code was written for ourselves and not as shareware for the community. The main goal was to make it easily understandable and to make the structure be driven by astrophysical requirements. Elegance and speed were not our prime concern. 2.3 Makefile -- this file is used to compile the code. You may want to tailor this to your need. 2.4 go_galaxy -- this is a script to run the code. Its most important purpose is to assign the directories where the auxiliary files reside. You may also want to modify this depending on your directory structure. 2.5 save_output -- this is a script to save the output files and give them reasonable names. (The input file with the model parameters is also included here.) To get started, we suggest to keep the names as they are in this file. 2.6 mod****.dat -- 10 files containing the Geneva evolutionary models for standard ("c") and enhanced ("e") mass-loss rates. The chemical composition is 2x solar (040), solar (020), 0.4x solar (008), 0.2x solar (004), and 0.05x solar (001). 2.7 lcb97_***.flu -- 5 files containing the model atmospheres of Lejeune et al. They match the metallicities of the evolutionary tracks. Note that some models were interpolated from the original Lejeune set since the required metallicities were not available. 2.8 wr_beta*.fluxes -- 2 files with the WR model atmospheres of Schmutz et al. (1992). 2.9 CMFGEN*.dat -- 10 files containing the WR models computed by the UCL group (5 metallicities, 2 WR types). 2.10 WMbasic*.dat -- 5 files with the Pauldrach O-star atmospheres for 5 metallicities. 2.11 sp.dat -- IUE spectral library of O, B, and WR stars used in the subroutine linesyn. This library is also available (in a more readable form) via the CD accompanying the article in PASP, 108, 996 (1996), or from our web page. 2.12 sp_low.dat -- FOS and STIS spectral library of LMC/SMC O stars used in the subroutine linesyn. The structure is identical to that of sp.dat. All stars other than O stars are the same as in sp.dat. (Leitherer et al. 2001; ApJ, 550, 724.) 2.13 fuse_high.dat -- FUSE spectral library of Galactic O and early B stars used in the subroutine fusesyn. The structure is identical to that of sp.dat. All stars other than O and early B have their continua set to 1. 2.14 fuse_low.dat -- same as fuse_high.dat but for LMC and SMC stars. Molecular hydrogen lines were removed. The reference for both FUSE libraries is Robert et al. (2002). 2.15 schkal.dat -- spectral-type calibration of Schmidt-Kaler (1982). The table simply contains a list of log L (in solar luminosities) and log Teff (in K). Each line corresponds to a certain spectral type. 2.16 irfeatures.dat -- contains data for the near-IR CO features at 1.62 and 2.29 microns, and the silicon feature at 1.59 microns from Origlia, Moorwood, and Oliva (1993) A&A, 280, 536. 2.17 standard.input1 -- this is the input file specifying the model parameters. It is explained in more detail below. 2.18 standard.* -- results of a standard model run. These files are obtained by using the parameters in the galaxy.input file. They can be used for test runs when implementing the code. 3. How to run the code Get organized first. Create a directory for the code and the unix scripts (files 2.1 -- 2.5). Create separate directories for the evolutionary tracks (2.6), for the model atmospheres (2.7 -- 2.10), and for the auxiliary files (2.11 -- 2.16). A directory for the input (2.17) and for the output (2.18) should also be created. The directory structure should correspond to what you have declared in 2.4 and 2.5. The output directory will be populated by files produced during a successful run. The source code should be pretty much machine independent, except for the names of the auxiliary files which are called in the code. Search for the string "*** THE FILE NAME IS LOCATION DEPENDENT ***" and modify accordingly. If you extract the code and all the files from the gzipped tar file, you will obtain the same filenames and folder structure as at STScI, and little renaming should be required. 3.1 The input Once all the files are in place and the declarations are complete, the input parameters need to be specified. If this is your first attempt, we suggest to leave the parameters as they are. They produce reasonable results and should give you a first impression of what is in store. Once you have gained more experience and have become more adventurous, you can modify the parameters to suit your needs. These are the parameters to play with: MODEL DESIGNATION: [NAME] standard --- any identifier you want to assign to the model. You will find it in the header of each output file. CONTINUOUS STAR FORMATION (>0) OR FIXED MASS (<=0): [ISF] -1 -- if this is a negative integer, star formation is instantaneous, otherwise it is continuous. TOTAL STELLAR MASS [10^6 SOLAR MASSES] IF 'FIXED MASS' IS CHOSEN: [TOMA] 1. -- this is the total stellar mass (spread out between the upper and lower cut-off masses). It is only used if an instantaneous burst is specified. SFR [SOLAR MASSES PER YEAR] IF 'CONT. SF' IS CHOSEN: [SFR] 1. -- the star formation rate (only used for a continuous rate). The total accumulated mass is spread out between the upper and lower cut-off masses. IMF EXPONENT (2.35 = SALPETER): [ALPHA] 2.35 -- IMF exponent. A power law is assumed. UPPER MASS LIMIT FOR IMF [SOLAR MASSES]: [UPMA] 100. -- upper mass limit for the IMF. LOWER MASS LIMIT FOR IMF [SOLAR MASSES]: [DOMA] 1. -- lower mass limit for the IMF. SUPERNOVA CUT-OFF MASS [SOLAR MASSES]: [SNCUT] 8.0 -- stars with ZAMS masses of 8 M and higher form supernovae. This is the suggested standard value but can be modified if desired. BLACK HOLE CUT-OFF MASS [SOLAR MASSES]: [BHCUT] 120. -- stars with ZAMS masses of 120 M and lower form supernovae. An alternative scenario would be to let stars above a certain threshold form a black hole. For instance, BHCUT=40. results in SNe only from the mass range 40 to 8 M. METALLICITY + TRACKS: [IZ] 92-94 STD. MASS-LOSS: 11=0.001; 12=0.004; 13=0.008; 14=0.020; 15=0.040 92-94 HIGH MASS-LOSS: 21=0.001; 22=0.004; 23=0.008; 24=0.020; 25=0.040 24 -- this integer indicates the evolutionary tracks to be used. The choices are 11-15 or 21-25, where the former select one of five metallicities with the standard tracks, and the latter do the same for the high mass-loss tracks. Example: "23" selects 40% solar metallicity with the high mass loss tracks. WIND MODEL (0=MAEDER; 1=EMP.; 2=THEOR.; 3=ELSON; 4=UCL): [IWIND] 2 -- this selects the wind model to be used for the calculation of the wind power. The four models are discussed in ApJ, 401, 498 (1992) and in Smith et al. (2002). "4" is the suggested default parameter. INITIAL TIME [1.E6 YEARS]: [TIME1] 0.01 -- the epoch of the onset of the star formation. In almost all cases you want this to be close to 0. It should not be exactly 0 for numerical reasons. 0.01 (i.e. 10e4 yr) is a good number. TIME STEP [1.e6 YEARS]: [STEP] 0.1 -- this is the timestep used for the calculations. It is a very important parameter. On the one hand, the computing time scales with STEP, so you want to avoid too high resolution, but on the other, short evolutionary phases can be missed. 0.1 (i.e. 10e5 yr) is a good value if you use full isochrone synthesis. If full isochrone synthesis is not used, 0.1 or larger is suggested only for tests --- be aware that WR or RSG numbers are no longer properly calculated for a STEP of 0.1 unless full isochrone synthesis is selected! LAST GRID POINT [1.e6 YEARS]: [TMAX] 20. -- the oldest age of the model. SMALL (=0) OR LARGE (=1) MASS GRID; ISOCHRONE ON LARGE GRID (=2) OR FULL ISOCHRONE (=3): [JMG] 3 -- these are four options for the interpolation in mass. They are explained in the code. Shortly: 0 -- evolutionary synthesis with a mass resolution of 5 M (only recommended for tests); 1 -- same as 0, but with a resolution of 1 M. This method was used in Leitherer & Heckman (1995); 2 -- isochrone synthesis with a fixed mass resolution of 1 M; 3 -- isochrone synthesis with a variable mass grid. This is the fanciest method and is the recommended mode. LMIN, LMAX (ALL=0): [LMIN,LMAX] 0 -- LMIN and LMAX are the indices of the evolutionary tracks, sorted by mass. Normally you do not want to mess with the variable and leave it at 0. However, if you want to track down some peculiarity of the output, you may want to compute the parameters for only one track. For instance, specifying 21,21 indicates that only a 100 M star should by used, and everything else is suppressed. The cross-ID's between index and mass are at the bottom of the input file. The example here refers to JMG=1 or 2. For JMG=0, you would have chosen 5,5. This does not apply to JMG=3 since the mass grid is variable and LMIN, LMAX are not used. TIME STEP TO PRINT THE CONTINUOUS AND LINE SPECTRA [1.e6YR]: [TDEL] 1.0 -- the file containing the output spectrum can be pretty big. This parameter controls the time step to print out the spectrum. This is independent of the time resolution -- only the print out is affected! 1 Myr is usually a good value but if you compute the starburst up to 100 Myr, you may prefer TDEL=5 Myr unless you have many Gb of disk space. ATMOSPHERE FOR SYNTHETIC SPECTRUM: 1=PLA,2=LEJ,3=LEJ+SCH,4=LEJ+HIL,5=PAU+HIL [IATMOS] 5 -- this is the choice of the model atmosphere. 1 is a bare-bone version with black bodies, good only for tests. 2 uses the Kurucz models as compiled by Lejeune for all stars. 3 uses Lejeune for stars with plane-parallel atmospheres and Schmutz for stars with strong winds. 4 uses Lejeune, but replaces the Schmutz by the Hillier atmospheres. 5 is like 4, except for the O atmospheres, for which we use the Pauldrach models. 5 is the recommended value. METALLICITY OF THE UV LINE SPECTRUM: (1=SOLAR, 2=LMC/SMC) [ILINE] 1 -- a switch for the choice of the UV spectral library. This switch applies to both the FUSE and the HST/IUE libraries. It is independent of the metallicity of the tracks/atmospheres. Normally one would use ILINE=1 with IZ=24 and ILINE=2 with IZ=22. RSG FEATURE: MICROTURB. VEL (1-6), SOL/NON-SOL ABUND (0,1) [IVT,IRSG] 3,0 -- atmospheric parameters used for the spectral features in the near-IR. Detailed explanations are in the sp-feature subroutine. Defaults are 3,0, i.e. microturbulent velocities of 3 km/sec and solar abundance ratios for alpha-element/Fe OUTPUT FILES (NO<0, YES>=0) [IO1,...] +1,+1,+1,+1,+1,+1,+1,+1,+1,+1,+1,+1,+1,+1 These are options to generate various outputs. We recommend to set all flags to "yes", at least until you become more familiar with the code. Some of the subroutines are interrelated. If you choose such a subroutine but not the other, required one, a warning will be issued. The 14 output flags are explained in the next section. 3.2. The output (1) Computation of the number of ionizing photons. (7) must be set to "yes" since the spectrum below 912 A is needed. Output is the number of ionizing photons in the HI, HeI, and HeII continuum, their fractions relative to the total luminosity, and the total luminosity. (2) Calculation of the supernova rate and the mechanical luminosities (winds and supernovae). It requires (4) to obtain the stellar wind luminosities. Otherwise it is independent of other subroutines. (3) HRD with a few evolutionary tracks. This is mostly useful for test purposes. This part is independent of all other subroutines and can be turned on/off without doing any harm. The HRD cannot be produced in Isochrone Synthesis mode. (4) Mechanical luminosity and related quantities due to stellar winds. (No supernovae!) This does not depend on any other subroutine since no information on the energy distribution is needed. (5) Two output files containing the stellar spectral types during each time step and the relative numbers of WR stars. The spectral types follow the scheme by Schmidt-Kaler, oversampled by a factor of 2. For instance, there are 18 entries for spectral type B. They are the number of stars for types B0, B0.5, B1,...B9.5 (total of 18). Schmidt-Kaler's table has B0, B1,....B9 (total of 9). The spectral types are printed out only every TDEL. Otherwise it is too bulky. (6) The mass in individual elements released via stellar winds and supernovae. No other subroutines are needed. The supernova yields for type II supernovae are taken into account. (7) The spectrum of the stellar population for each time step. The columns are time, wavelength, stellar+nebular, stellar only, and nebular only fluxes. (1) is needed in order to calculate the nebular continuum. (8) The ultraviolet line spectrum at 0.75 A resolution from 1200 to 1600 A (LMC/SMC library) or to 1800 A (Milky Way library). The subroutine needs (7) to compute the stellar continuum and (1) for the nebular continuum. If (1) is turned off, the nebular contribution can not be added (it is often small, though). The columns have time, wavelength, absolute luminosity, and rectified (continuum=1) luminosity. (9) Calculation of colors and magnitudes. The subroutine needs (7) to compute the stellar continuum and (1) for the nebular continuum. If (1) is turned off, the nebular contribution can not be added and the computed colors are for stars only (this may sometimes be desirable). The filter system is defined in the code. (10) Calculation of the strengths of H_alpha, H_beta, Pa_beta, and Br_gamma. For each line we give the continuum luminosity, the line luminosity, and the equivalent width (everything logarithmic). The subroutine needs (7) to compute the stellar continuum and (1) for the number of ionizing photons. (11) Calculation of the strengths of various IR spectral features. First is the CO index as computed by Doyon et al. (1994, ApJ, 421, 101). (Please note that this calculation has no metallicity dependence. A later version of this routine will compute the CO index using the model atmospheres themselves and give metallicity-dependent results.) Next are two computations of the CaII IR triplet using the relations of Diaz et al. (1989, MNRAS, 239, 325). The relations from Diaz et al. have no temperature dependence; the first calculation has the feature present in stars of all temperatures; the second has the index set to zero strength for stars with T>7200K (spectral type A or earlier). Next come the 1.62 and 2.29 micron CO features, and the 1.59 micron Si feature, which were modeled for individual stars by Origlia et al. (1993, A&A, 280, 536.) The indices can be computed for solar [Si/Fe] and [C/Fe], or a model with enhanced [Si/Fe] and depleted [C/Fe] (as for young systems enriched primarily by Type II SNe), and for stellar atmospheric microturbulent velocities (MTVs) of 1-6 km/s. (Note that the changes to the abundance ratios and MTVs are self-contained in this routine and have no effect upon the other outputs, e.g., colors, of the code.) (12) This subroutine is equivalent to (8), but it computes the spectral region between 1000 and 1180 A. (13) This subroutine is still under construction. Eventually it will compute a high-resolution line spectrum in the optical. Stay tuned! (14) Calculation of the most important WR emission lines. These are only those lines originating in WR winds --- not the nebular lines in the HII region. Quantities given are the line fluxes and the equivalent widths. The subroutine needs (7) to compute the stellar continuum and (1) for the nebular continuum. If (1) is turned off, the nebular contribution cannot be added. 4. Need Help? Feel free to ask us if you run into trouble or have questions. But please understand that we have no official help desk, and we try to help you in our spare time. There is no guarantee that we can trouble-shoot your problems.