The goal of this project is to compute the dust echo from a jetted transient source. We have developed a code (in Python 3.7) that calculates the sublimation radius for dust grains of different sizes and the extinction in an iterative way, and then the volumetric emissivity at the observer's wavelength is calculated from the dust temperatures for all grains that have not sublimated. Finally, we account for light-travel delays and obtain the lightcurve of dust echo from a given observer’s viewing angle and observing wavelength.
Below is a guide to calculate the dust echo lightcurve for arbitrary viewing angle (wrt. the jet axis) and observer's wavelength.
The first step is to change the two directories in "fpath.py": 'codedir' should be the directory containing all the codes in the current project, and 'savedir' should be the directory where the user wants to save the results.
The second step is to run "generate_Td_asub.py", which computes the dust temperature T_d(r, t) and the critical grain size below which dust has already sublimated a_sub(r, t), both as 2-dimensional functions of radius r (distance to the source) and retarded time t at that radius (since the arrival of the first photon from the source). This is the most computational intensive part and should take somewhere between 10 seconds up to few minutes, depending on the grid resolution adopted. The results are stored in files 'nH%.1e_Td.txt' % nH0 and 'nH%.1e_asub.txt' % nH0 in the directory 'savedir' specified in "fpath.py". There are a number of adjustable parameters in "generate_Td_asub.py", the main ones being: nH0 (hydrogen number density of the CSM), tmax (the duration of the UV-optical pulse from the source), LUV (the isotropic equivalent luminosity of the source). The current version of the code assumes that the CSM is uniform and that the source lightcurve is a step-function for duration tmax. Within the current code architecture, it takes very little additional effort to drop these two assumptions and add more general functional forms of nH0(r) and LUV(t), because the dust temperature is computed on a 2-dimensional grid of (r, t) anyway.
The third and final step is to run "run_Ldnu_obs.py", which computes the isotropic equivalent specific luminosity L_nu(tobs) and projected position of the flux centroid xcentr(tobs), both as functions of the observer's time tobs. The results are stored in a file 'nH%.1e_lam%.2fum_theobs%.2f_Ldnu_xcentr.txt' % (nH0, lamobs, theobs) in the directory 'savedir' specified in "fpath.py". The adjustable parameters in "run_Ldnu_obs.py" are: nH0 (hydrogen number density in cm^-3), lamobs (observer's wavelength in micro-meter), and theobs (observer's viewing angle in radian). Each of these three parameters could be a list of many values, and the code will compute the results for all combinations of the listed values. One has to make sure for all choices of nH0, the corresponding dust temperature profiles have already been computed by "generate_Td_asub.py" in the earlier step. There is another important hidden adjustable parameter in "Ldnu_obs.py" (which is called by "run_Ldnu_obs.py") and that is the half opening angle of the jet: thej (in radian). The default value is 4 degrees = 0.07 rad. The jet opening angle determines how far away from the jet axis dust grains are heated by the jet emission. The current version of the code only allows for a "top-hat" jet angular structure, meaning that there is no jet emission beyond the angle thej.
The resulting txt files from "run_Ldnu_obs.py" are structured as follows. There are four rows (separated by '\n'). The first row shows the info for the observer's time grid, which is taken to be logarithmic grid between 'tobsmin' and 'tobsmax' with 'Ntobs' number of grid points and the right boundary is not included. For Python users, the time grid is given by
tobs = numpy.logspace(tobsmin, tobsmax, Ntobs)
The second line shows the meaning of the later two rows of numbers: 'Ldnu' is the specific luminosity [in erg/s/Hz] at the wavelength specificed in the filename; 'xcentr' is the projected distance [in pc] between the flux centroid and the center of explosion on the sky. To obtain the projected angualr position of the flux centroid, one needs to divide xcentr by the angular diameter distance to the source. The third row shows the value of 'Ldnu' and the last row shows 'xcentr', at each time tobs. When Ldnu = 0 (dust echo hasn't arrived yet), we have xcentr = numpy.nan (does not exist).
Note that the dust echo from the counter jet (which propagates away from the observer) is not included, and we do not expect the counter-jet echo to contribute significantly to the total flux in the first decade after the explosion. It takes very little additional effort to compute the counter-jet echo, if needed.
The biggest uncertainty is the isotropic equivalent UV luminosity 'LUV' and duration 'tmax' from the source. In the case of short GRBs, dust heating is dominated by the UV-optical photons from the reverse shock, so LUV and tmax will depend on the CSM density at a distance of ~0.1pc from the source, initial jet Lorentz factor Gamma_0, and the isotropic equivalent energy of the jet E_j. For typical parameters (E_j ~ 1e52 erg and Gamma_0 ~ 100), we expect the total energy E ~ LUV * tmax to be in the range (1e48, 3e51) ergs. The default choice is LUV = 3e47 erg/s and tmax = 300 sec (which is for a reasonably high CSM density of ~0.1 cm^-3). Lower CSM density will decrease the total radiation energy of the source.
See arXiv:2108.04243 for a detailed description of the method used in this project.