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?