Let $\tau_0:=\frac{1}{s}\tau_1$, then $\tau_i=s^{i-1}\tau_1$ and $v_i=D(\tau_{i+1}+\tau_i)=s^i v_0$ $$\int\limits_{0}^{\infty} G_{3D}(2D\tau,r)e^{-\mu_a\tau}d\tau \approx \sum\limits_{i=0}^{k-1}\int\limits_{\tau_i}^{\tau_{i+1}} G_{3D}(2D\tau,r)e^{-\mu_a\tau}d\tau = \sum\limits_{i=0}^{k-1}(\tau_{i+1}-\tau_i) G_{3D}(D(\tau_{i+1}+\tau_i),r)e^{-\mu_a\frac{\tau_{i+1}+\tau_i}{2}} = \frac{s-1}{s+1}\sum\limits_{i=0}^{k-1}(\tau_{i+1}+\tau_i) G_{3D}(D(\tau_{i+1}+\tau_i),r)e^{-\mu_a\frac{\tau_{i+1}+\tau_i}{2}} = \frac{s-1}{s+1}\sum\limits_{i=0}^{k-1}\frac{s^i v_0}{D} G_{3D}(s^i v_0,r)e^{-\mu_a\frac{s^i v_0}{2D}}$$

Given $s=\frac{1+\sqrt{5}}{2}$, $ \frac{s-1}{s+1} \approx 0.236068$ rather than $0.240606$ in the paper. Can someone tell me how to get that magic number? And why $v_0$ is specified in terms of $\mu_a$? Thanks for your kindness.