I realise I’ve been neglecting the blog a bit the last few months, with only a few posts that link to other things, but that’s because I’ve been solving a nasty problem for my research and I think it’s now time to update what’s going on with a ‘real’ post.
For the past year or so, we’ve been working on modelling a galaxy using a set of models developed by my supervisor, Larry Widrow, and his collaborator, John Dubinski, who is at Toronto. The models assume a standard exponential disc for spiral galaxies, but the curious thing is that, because discs appear to contain a luminous bulge at the centre (which is generally considered to be a dynamically distinct entity) there is degeneracy in how surface brightness profiles can be broken down: typically, we break them down into a part due to the bulge (which is often assumed to be a de Vaucouleurs-type r^1/4 profile) and a part due to the disc, which is assumed to be exponential. In the centre, the bulge and the disc contribute to the total light seen in a galaxy. There is some theory behind assuming that discs form exponential density distribution, but in principle other functional forms are allowed for disks. One of these is called a Kormendy disc (named after John Kormendy who proposed this in 1977).
A Kormendy disc is basically an exponential disc with a hole in the centre. So the centre of the disc emits no light, because there are no stars in the centre. This sounds weird, but there is nothing preventing us from attributing all the central light to the bulge, and letting the disc take over as the bulge falls off – we can find fits using the Kormendy functional form as well as the standard exponential. In fact, a siginificant minority of galaxies appear to be better fit by a bulge + Kormendy disc than by a bulge + exponential disc.
The presence of a hole is highly counterintuitive, and you may wonder if ‘seeing’ a hole is merely a measurement artifact or a real measurement due to some other effect. Neither of these appears to be the case; measurement artifacts have been ruled out, and effects such as central dust (which could obscure disc light) are unlikely to be the cause. But, there are good reasons to think it reflects a real physical effect; many of the galaxies that display holes are known to be barred (i.e., they have a central bar of stars out of which typically emanate spiral arms), and bars could be responsible for the effect, but even some that don’t have bars are well fit by holes.
So what has this got to do with my research? Well, the galaxy we’re looking at might in fact be fit better by a central hole than by the exponential disc in our models. Originally we accounted for this by pretending that the real mass distribution was exponential while the light distribution had a Kormendy hole. This is reasonable, because the galaxy appears to have some dust at the centre, although, as I mentioned above, holes are not likely due to dust. So we decided to fundamentally change the structure of the models to generate a Kormendy disc instead, and that’s what I’ve been doing for the past two months… and it’s a pain, because the code is written in Fortran (the F-word again) 77, and Fortran 77 is just plain annoying.
Fortunately, the code is structured in such a way that the radial profile can be changed separately without altering anything else, so it was in principle just a matter of changing the density wherever it was called in the code. Unfortunately, the code itself has been changed several times over the past 15 years and the syntax and variable names used across different subroutines are inconsistent. Moreover, just changing the density profile at first did not work, partly because there is a mathematical (though not a physical) singularity at the centre; to rectify this, I just redefined the central density (and the first two derivatives) to be zero. Debugging was annoying because outputting important quantities in subroutines ran into fortran’s absurd I/O recursion error, which really seems to be an annoying bug that has just never been fixed (as evidenced by the fact that it would pop up sometimes and not at other times when I had made no changes to the code). Finally, the code has a bothersome habit of giving different results for the same input when compiled on different machines… I’ve quietly decided to ignore that since it seems to work on my department machine (but on my home desktop or my laptop). Regardless, after two months of fighting through oceans of NANs and INFs, outputting dozens of different quantities at all points of the program, and banging my head against the wall (not literally, but it wasn’t going to stay not literal for long), I tried an arbitrary but well-behaved radial profile and found that it worked perfectly, which indicated that the main problem was likely a coding bug or something having to do with the centre of the Kormendy disc. I finally traced the problem to the calculation of certain spline quantities at the centre… where the density was zero. Instead of fighting my way through tedious, possibly buggy numerical subroutines, I changed the density (and the first two derivatives) to a very small nonzero quantity… and presto! We have self-consistent Kormendy discs.
So now we have real Kormendy discs, and we can model galaxies with either a Kormendy profile or an exponential profile, and the sun is shining, and life is grand (well, sort of, but that’s another story :P).