Details on Detailed Balancing

[henrikjacobsenfys]

General

Detailed Balancing, briefly mention in a comment in #2 is going to cause some headaches, and needs to be handled and documented carefully.

The core is simple: it is easier to create excitations than to absorb them at low temperature, and therefore, the scattering intensity is vanishing on the neutron energy gain/sample energy loss (E<0) side of the spectrum.
We define:

• $\omega$ is the energy transfer,
• $\beta = 1 / (k_B T)$ is the inverse temperature,
• $n(\omega) = \frac{1}{e^{\beta \omega}-1}$ is the Bose occupation factor.
• $n(\omega)+1 = \frac{1}{1 - e^{-\beta \omega}}$ is a useful relation to know

I have set $\hbar=1$ for simplicity.

Mathematically, detailed balancing is captured by the following relation of the scattering function (also known as the dynamical structure factor)

$$S(-Q,\omega)=\exp(-\beta \omega) S (Q,\omega).$$

It is related to the generalized susceptibiity $\chi''(\omega)$ by

$$S(\omega) = \frac{1}{\pi} \cdot (n(\omega)+1) \cdot \chi''(\omega).$$

Another way to write this important relation is (notice the difference between $\chi'$ and $\chi''$ )

$$S(\omega) = \omega \cdot (n(\omega)+1) \cdot \chi'(Q,\omega=0)\cdot F(\omega),$$

where $F(\omega)$ is the spectral weight function, which has area 1. More on that in a moment. Importantly, $\chi'(Q,\omega=0)$ does not depend on $\omega$.

$\chi'$ and $\chi''$ is typically what is calculated theoretically, while $S$ is typically what is measured. When fitting, we usually allow fitting of the area of the peak in the case DBF=1 (usually in arbitrary units), and is proportional to $\chi'(Q,\omega=0)$.

In other words, the fit function, if the data is $S(Q,\omega)$, is

$$I=DBF\cdot A \cdot F(\omega),$$

where
$DBF=\omega*(n(\omega)+1)$ is called the Detailed Balance Factor
$A$ is the "area" of the peak (or would be, if DBF=1)
$F$ is the shape of the peak, usually a Lorentzian, Gaussian, pairs of Lorentzians or Gaussians, or Damped Harmonic Oscillator

The main points are the following:

1. Many people are unaware that they are doing this, and ignore/misuse the DBF, i.e. they fit simply $I=A\cdot F(\omega)$ (some people fit amplitude instead of area as well). This may give wrong results. For example, in the figure below, the T=3 Lorentzian with DBF looks pretty much like a Lorentzian on a sloping background
2. Many people are unaware that the DBF is temperature dependent
3. Because of1) and 2), $A$ is very different from what one would naively expect the "amplitude" of the peak to be. This is illustrated in the image below, which shows 4 Lorentzians with amplitude 1, 3 of which have been multiplied by the DBF for various temperatures. This gives two problems
   • Finding the right start guess is difficult for users, as the "amplitude" doesn't match what they see. This is exacerbated by the resolution function smearing out peaks, lowering the amplitude further
   • Interpreting the fitted amplitude parameter is error prone - an apparent change in amplitude with temperature may simply be the DBF

Proposal/Things to discuss

• I will implement the DBF as described here, as it's the correct thing to do. Care must be taken at $T=0$, and at low $T$, but as shown in the figure, it's already solved.
• In addition to this text, when it's time to make tutorials, I will make a detailed (haha) tutorial on fitting using detailed balancing, with references, to inform users.
• It makes sense to work with $A$ and $F$ as defined here, which means $A$ is closer to an area than amplitude. Since area is conserved under convolutions, this also makes a lot of sense. However,
  • Area and width are strongly correlated when fitting. We should therefore convert $A$ to amplitude and fit that, then convert back. EDIT: as shown below, amplitude and width are also correlated due to the convolution. We therefore save nothing by using amplitude.
  • It's more natural for users to work with apparent amplitude, i.e. the height of the peak as they see it in their data. It may therefore be a good idea to allow users to input an apparent amplitude parameter, to make it easier to figure out start guesses.
  • This means we have "$A$", "amplitude", and "apparent amplitude". It's unclear to me how many of these 3 amplitude/area parameters we need to work with - at least 2, but perhaps all 3?

References: Principles of Neutron Scattering from Condensed Matter by Andrew Boothroyd

[henrikjacobsenfys]

To illustrate the issue we discussed today, here's what a (simple) data set might look like:

Users will want to fit this by giving reasonable start guesses. They'll already have the resolution fixed at this point.

To be clear, with this data, any method will give a good result and the fit will converge to the true values. I'm using this to illustrate the difficulties.

If we use Amplitude

By eye, the data has amplitude 1.0, and a width (full width at half max) of maybe 0.2 meV. There's also a background of 0.05.
Here's how that looks:

This is definitely good enough to get a good fit, but the height of the guess doesn't match the input because of the DBF. If we want to get closer, let's increase the peak height by a factor 2.7:

That's very nice, but it's clearly too wide. Let's adjust:

This actually surprised me - the amplitude changed! This is because of the convolution with the resolution, which conserves area, not amplitude.
This means that amplitude and width are correlated, and there is therefore no benefit to fitting amplitude and width rather than area and width - I was wrong about this.

If we use Area

It's very hard to guess what the area should be. I'll do a reasonable, but not correct, guess of 0.5 meV to illustrate:

This undershoots the peak height by a factor 2.5 - let's adjust:

This is very nice! But a bit too broad, let's make it narrower:

The width is better, but now the peak height is wrong, even though we just got it right. This is very frustrating, especially when data is difficult to fit.
It also clearly illustrates why area and width are correlated. (I wrote the Area section first, then the Amplitude section).

What I believe users want:

They definitely want to get the Area out at the end - it's the physically meaningful parameter. But they will also want to adjust height and width of their start guess to match their data, without worrying about correlated parameters - the multiple steps of "now the height is right, but width is wrong - wait, why is the height not right anymore?!" is very frustrating. We should therefore allow the amplitude after convolution as an input (and make it very clear that it's not a physically meaningful parameter), and calculate the normalisation numerically, simply by dividing by the max value of the function. For now, I don't want to fit the numerically normalised function, as that is very error prone.

As to how to fit it, people absolutely want to be able to let Area/amplitude and width vary simultaneously. Besides, it's clear that this is required, since, as we've illustrated, the height of the peak depends on the width parameter.