live_or_die doesn't agree with extract for spherical and supernova models

[jhmatthews]

I've been working with @smangham a little trying to get the reverberation mapping stuff to work, and have come across an issue in which live or die and extract modes produce different results.

To reproduce, run ~/Dropbox/Python/bug_pf_files/sphere_out.pf in extract and live or die.

the results are as follows, for one observer - top is wavelength and bottom is velocity. vertical lines mark line centre in the top panel, and the min and max velocities in the wind on the bottom panel.

This doesn't necessarily mean there is something majorly wrong with extract mode or anything but we need to understand the differences as Sam was planning on using live or die for the moment.

I've run this as far back as Python 77 with identical results to the dev version, and reproduced problems in supernovae and spherical wind models.

[kslong]

It does mean that something has changed, and I regard this as quite serious. I ran various tests like this previously. Does it fail in this dramatic fashion with a SV wind?

[lazygun37]

A couple of thoughts:

-- I think the velocity plot is "wrong" in that the x-axis is inverted -- blue shifts should correspond to negative velocities, i.e. P Cyg profiles should lok the same (absorption on the left) regardless of whether you plot in wavelength or velocity space.

-- is this a pure spherical outflow? no rotation? and is it illuminated by a disk or by the central star? does it contain a disk? is thermal wind emission included? If its a purely spherical, radial outflow without rotation, without wind emission and without any disk in the system, then both modes are giving the wrong answer. In that situation, all sightlines are equivalent, so energy/photon number conservation demands that the integrated flux in absorption and emission must be the same, since the input spectrum contained no line. To my eye, it doesn't look like the run contained rotation (which would tend to smooth the profile). But it may be illuminated by a disk and/or include wind emission. THe simple model described above is a useful test case.

[kslong]

I have run several checks at this point:

First, with the existing dev branch version of python and a standard SV model, there are significant differences in the live or die and normal extraction of the spectrum.

Second, with py68e, the oldest version I could easily compile, the spectra for the two types of extractions are very similar.

My tentative conclusion is that we have broken live or die somewhere along the way.

It will be interesting to see how things fall out with the first version in Git.

[jhmatthews]

Christian-

• The model should be spherical outflow, no emission, no rotation, no disk etc, and so I agree that the lack of symmetry in the P-Cyg is worrying
• yes I guess I should have plotted velocity the other way around technically.

The SV wind model produces differences in the CIV profile, for example (this is a viewing angle of 45 degrees):

I'll upload a better plot tomorrow when I'm in the office, but basically I think you're just seeing subtle differences in the absorption profile, and it strikes me as rather similar to issue #39 - something I've tried to solve several times but with no answer as yet.

This will be my priority - it's essential for CV models and Sam's work. I'm worried that breaking live or die may mean we have broken radiative transfer / scattering in general...

[jhmatthews]

I've updated the above plot to a better one. Here's the full SV spectrum in both modes:

The differences are pretty subtle, but they are there, and it really does remind me of #39 - I will test with some older models. We need to sort this pretty sharpish.

[jhmatthews]

I've also run the spherical model with two different angles- 45 and 80 degrees. The results for both modes are not good:

It is my understanding that all four of the lines plotted here should look the same and be classic P-Cygs, absorption and emission summing to zero with no angle dependence. Clearly the bottom live or die plot has something badly wrong as well.

This seems to be a much more exaggerated problem than in the SV model, so there could be a problem in general with our spherical model.

[jhmatthews]

Ok - I've traced most of this problem. Between Python 70f and 70g this line was deleted from line ~1387 of resonate.c:

doppler (&pold, p, v, *nres);   /* Get the final frequency of the photon */

I think this was done accidentally because it was embedded in some logging statements.

What should happen, I believe is the following (correct me if I'm wrong):

• 1. if a resonance line scatter occurs, it does so at a point such that nu = ( 1 - v/c) nu_line, where v here is the projected velocity along the photon direction.
• 2. after the scatter, we choose a new random direction, depending on the scattering mode
• 3. we should then adjust the frequency of the photon to nu_out = nu_line / ( 1 - v/c), where v is now the projected velocity along the new direction

Bottom line: I think what is happening is that the last step was accidentally deleted in resonate, meaning that for the last 2+ years we have not been shifting the frequencies of photons in ionization cycles or live or die mode.

The reason we haven't noticed is because extract has this final step, and probably we never thought to check live or die mode. This bug implies we should, and I also suggest that duplicating code here is a bad idea (i.e. there should be a governing routine for frequency shifts or something).

Adding the line containing the doppler function in the appropriate place in resonate.c produces the following results:

There is still an angle dependence here, which must be wrong. With that addition the live or die and extract modes agree.

[kslong]

Good detective work, James. By “there is still an angle dependence”, do you mean that for a spherical model we are still seeing differences in the output spectrum. If so, this was also something I thought I had checked in the distant past. If not, could you clarify what you meant by this comment, and with what .pf file.

[jhmatthews]

pf file for 2 angles, spherical model: ~/Dropbox/Python/bug_pf_files/sphere_out_2angles.pf

What I mean is that for both extract and live or die the spectra are different at 45 and 80 degrees. I believe they should look identical for a spherical outflow from a star.

I realised this pf file is in spherical coordinates, so I will check in cylindrical too when I get a minute.

[jhmatthews]

The above commit re-implements the correct line of doppler shifting in resonate.c. This issue should remain open until we understand why line profiles change with angle in the spherical model, and also until the tests on the fiducial and Proga models have been done.

[Higginbottom]

I have rerun the baseline Proga model used to produce the almost featureless spectra in figure 6. The new run is essentially identical (I won't mention the significant irritation of rewriting my plotting routines to take account of a new field "Created" in the spectral files - Oh wait - I didn't...)
I shall rerun the other runs, and just remake all the figures, and I think if they look the same by eye, then there is no rewrite needed!!!

[jhmatthews]

Thanks Nick. I'm closing this issue then. The results of mine and Nick's tests can be seen at:

https://sites.google.com/site/agnwinds/current-work-in-progress/fiducial-qso-model-correct-doppler-shift-and-scattering

https://sites.google.com/site/agnwinds/current-work-in-progress/higginbottom-proga-runs-with-d

and for SN models:

https://sites.google.com/site/agnwinds/current-work-in-progress/tardis-v-python/spectral-comparisons