Not only can the causality and viscosity features of time domain methods be incorporated in the frequency domain, but the square roots and complex exponentials of the phase shift method can be replaced by complex multiplies and divides. First note that the square root expansion need not begin from a starting guess or from a lower order iteration. It could begin from the square root previously obtained from the preceding or k. (I noticed this when an early version of the program had a bug that made all my 15 calculations look like 90 calculations)! Also, is necessarily small in the phase shift method to accommodate imaging at every time point. So the finite differencing form (60) is probably as good as the complex exponential. Why don't you try it?