The second term in Eq. (4) is Display Formula
$\u222bSn\u22121\u222bSn\u22121\phi *(s^)P(s^,s^\u2032)\phi (s^\u2032)ds^\u2032\u2009ds^,$(21)
which contains the phase function given in Eq. (22) and can be expanded using a Fourier series in powers of $g$^{37}Display Formula$P(s^\xb7s^\u2032;g)=12\pi +1\pi \u2211l=1\u221egl\u2009cos(l\Delta \theta ),$(22)
where $\Delta \theta =arccos(s^\xb7s^\u2032)$. Thus we can write Display Formula$\u222b2\pi \u222b2\pi \phi (s^\u2032)P\theta (s^,s^\u2032)\phi *(s^)ds^\u2009ds^\u2032=\u222b2\pi \u222b2\pi [12\pi a0+1\pi \u2211n=1\u221ean\u2009cos(n\theta \u2032)+1\pi \u2211n=1\u221ebn\u2009sin(n\theta \u2032)]\u2062{12\pi +1\pi \u2211l=0\u221egl\u2009cos[l(\theta \u2212\theta \u2032)]}\u2062[12\pi a0*+1\pi \u2211m=1\u221eam*\u2009cos(m\theta )+1\pi \u2211m=1\u221ebm*\u2009sin(m\theta )]d\theta \u2009d\theta \u2032,$(23)
where $\theta $ and $\theta \u2032$ are the angles between the $z$-axis and $s^$ and $s^\u2032$, respectively. As such, the scattering angle between the previous direction $s^\u2032$ into the new direction $s^$ is given by $(\theta \u2212\theta \u2032)$. It is possible to expand $cos[l(\theta \u2212\theta \u2032)]$ as $cos(l\theta )cos(l\theta \u2032)+sin(l\theta )sin(l\theta \u2032)$ which in turn allows us to employ orthogonality relationships to simplify the above integrals and write Display Formula$\u222bS1\u222bS1\phi (s^\u2032)P\theta (s^,s^\u2032)\phi *(s^)ds^\u2009ds^\u2032=12\pi a0a0*+1\pi \u2211n=1\u221eanan*gn+1\pi \u2211n=1\u221ebnbn*gn.$(24)
Substituting this expression into Eq. (4), we can write the full expression for the functional gradient with respect to the scattering coefficient Display Formula$\u2202\epsilon \u2202\mu s=12\pi a0a0*+1\pi \u2211n=1\u221eanan*+1\pi \u2211n=1\u221ebnbn*\u221212\pi a0a0*+1\pi \u2211n=1\u221eanan*gn+1\pi \u2211n=1\u221ebnbn*gn$(25)
Display Formula$=1\pi \u2211n=1\u221e[anan*+bnbn*](1\u2212gn).$(26)
The ability to calculate these gradients is the first step to finding a computationally efficient way to solve the full QPAT inversion using an MC model of light transport.