### External module function mqcalc - used by mcdiff

''going beyond'' dipolar approximation is a desirable feature of an accurate calculation of magnetic neutron scattering intensity and can be performed using mcphase. The internal module so1ion does this in an excellent way for the ground state multiplet of rare earth atoms (see section 8.5) using equations (49) for the scattering operator . A more general formula for this scattering operator operator relates this operator to the Fourier transform of the magnetisation density of the unfilled shell of a specific ion [26]:

 (75)

Note that the magnetisation density consists of two contributions , the spin and orbital contribution. The orbital contribution is not uniquely defined due to a gauge freedom ( , the curl of the magnetisation must give the current density, so any gradient of a potential may be added to without changing the result). However, this does not matter, because adding a gradient of a potential to will just give a contribution to the equation (78), thus the neutron scattering cross section is not sensitive to the chosen gauge, which is an important feature of the theory.

Technically, ''going beyond'' dipolar approximation in the program mcdiff can be done with module functions mq and estates. The output of mq is the Fourier transform of the magnetization density in units of , which is equivalent to two times the scattering operator for a given orientation of the scattering vector. mq is called many times, for every scattering vector. In order to do an efficient calculation the eigenstates should be calculated only once, this is the task of function estates.

The format to be used is:

extern "C" void mqcalc(ComplexVector & mq,double & th,double & ph,double J0,
double & J2, double & J4, double & J6,ComplexMatrix & est)


The meaning of the symbols is as follows:

on input
th     polar angle theta of the scattering vector Q (angle with the axb axis=c axis) in rad
ph     polar angle phi of the scattering vector Q (angle with bx(axb)=a in the
projection into the  bx(axb),b plane = angle with a in the projection into the
J0,J2,J4,J6     form factor functions <jn(Q)>
est             eigenstate matrix (as calculated by estates),
it should also contain population numbers of the states
on output
mq(1..3)        Kartesian components of the scattering operator mq=<M(Q)>/muB=-2<Q>_TH
Note that:
mq(1)=<Mbx(axb)(Q)>/muB
mq(2)=<Mb(Q)>/muB
mq(3)=<Maxb(Q)>/muB
according to Lovesey Neutron Scattering equation 6.87b the
scattering operator is given in  spherical coordinates
Q-1,Q0,Q+1 (introduced as described above on input of th and ph)
these are related to cartesion coordinates by 11.123
thus at Q=0  <M(Q)>=2<S>+<L>


Martin Rotter 2017-01-10