\begin{equation}  \label{eqn:system}  \left(\mathbf{B} + \mathbf{D}\right) \mathbf{m} = \mathbf6661616\;,\end{equation}where $\mathbf{m}$ is the vector of extrapolated wavefield, and$\mathbf1111161$ is an appropriate righthand side. The helix transformprovides us with the operator $\mathbf{B}^{-1}$, which we can use toprecondition system (\ref{eqn:system}). Introducing the change ofvariables\begin{equation}  \label{eqn:prec}  \mathbf{m} = \mathbf{B}^{-1} \mathbf{x}\;,\end{equation}we can transform the original system (\ref{eqn:system}) to the form\begin{equation}  \label{eqn:system2}  \mathbf6166111 = \left(\mathbf{B} + \mathbf{D}\right) \mathbf{B}^{-1} \mathbf{x} =   \mathbf{x} + \mathbf{D}\,\mathbf{B}^{-1} \mathbf{x}\;.\end{equation}When the velocity perturbation is small\footnote{In the  linear-algebraic sense, this means that the spectral radius of  operator $\mathbf{D}\,\mathbf{B}^{-1}$ is strictly less than one.}, eventhe simple iteration\begin{eqnarray}  \label{ralg}  \mathbf{x}_0 & = & \mathbf6611661\;; \\  \mathbf{x}_k & = & \mathbf6611161 - \mathbf{D}\,\mathbf{B}^{-1} \mathbf{x}_{k+1}\;\end{eqnarray}will converge rapidly to the desired solution. This interestingpossibility needs thorough testing.\parThe third untested possibility (Papanicolaou, personal communication)is to implement a clever patching in the velocity domain, applying aconstant-velocity filter locally inside each patch. Recently developedfast wavelet transform techniques \cite[]{wavelet}, in particular the\emph{local cosine transform}, provide a formal framework for thatapproach.\section{Conclusions}The feasibility of multidimensional deconvolution, proven by thehelix transform, allows us to revisit the problem of implicitwavefield extrapolation in three dimensions. The attraction ofimplicit finite-difference methods lies in their unconditionalstability, a property invaluable for practical applications.\parWe have shown that at least in the constant coefficient case (that is,laterally invariant velocity), it is possible to implement anextremely efficient implicit extrapolation by a recursive inversefiltering in the helix-transformed computational model. Unfortunately,the case of lateral velocity variations still presents a difficultproblem that may not have an exact solution. We are currentlyexploring different roads to that goal.\section{Acknowledgments}We thank Biondo Biondi for useful discussions and for sharing histhree-dimensional expertise.\bibliographystyle{seg}\bibliography{SEP2,SEG,paper}\append{The 1/6-th trick}\inputdir{Math}Given the filter $D_2 (k)$, defined in formula (\ref{eqn:d2k}), we canconstruct an accurate approximation to the second derivative operator$-k^2$ by considering a filter ratio (another Pad\'{e}-typeapproximation) of the form\begin{equation}  \label{eqn:ratio}  -k^2 \approx \frac{D_2(k)}{1 + \beta D_2 (k)}\;,\end{equation}where $\beta$ is an adjustable constant \cite[]{Claerbout.blackwell.85}.The actual Pad\'{e} coefficient is $\beta=1/12$.  As pointed out byFrancis Muir, the value of $\beta = 1/4 - 1/\pi^2 \approx 1/6.726$gives an exact fit at the Nyquist frequency $k = \pi$. Fitting thederivative operator in the $L_1$ norm yields the value of $\beta\approx 1/8.13$. All these approximations are shown in Figure\ref{fig:sixth}.\sideplot{sixth}{width=2.5in}{The second-derivative operator in the  wavenumber domain and its approximations.}\append{Constructing an ``isotropic'' Laplacian operator}\inputdir{Math}The problem of approximating the Laplacian operator in two dimensionsnot only inherits the inaccuracies of the one-dimensionalfinite-difference approximations, but also raises the issue ofazimuthal asymmetry. For example, the usual five-point filter\begin{equation}\label{eqn:usual}F_5 =\begin{array}{|r|r|r|}\hline0 & 1 & 0 \\ \hline1 & -4 & 1 \\ \hline0 & 1 & 0  \\ \hline \end{array}\end{equation}exhibits a clear difference between the grid directions and thedirections at a 45-degree angle to the grid. To overcome thisunpleasant anisotropy, we can consider a slightly larger filter of theform\begin{equation}\label{eqn:new}F_9 = \begin{array}{|r|r|r|}\hline\alpha & \gamma & \alpha \\ \hline\gamma & -4 (\alpha+\gamma) & \gamma \\ \hline\alpha & \gamma & \alpha  \\ \hline \end{array}\end{equation}where the constants $\alpha$ and $\gamma$ are to be defined. TheFourier-domain representation of filter (\ref{eqn:new}) is\begin{equation}  \label{eqn:l9k}  F_9 (k_x,k_y) = 4\,\alpha\,[\cos{k_x}\,\cos{k_y} - 1] +  2\,\gamma\,[\cos{k_x}+\cos{k_y}-2]\;,\end{equation}and the isotropic filter that we can try to approximate is definedanalogously to its one-dimensional equivalent, as follows:\begin{equation}  \label{eqn:lk}  F (k_x,k_y) = 2 (\cos{k} -1) = 2 (\cos{\sqrt{k_x^2+k_y^2}} - 1)\;.\end{equation}Comparing equations (\ref{eqn:l9k}) and (\ref{eqn:lk}), we notice thatthey match exactly, when either of the wavenumbers $k_x$ or $k_y$ isequal to zero, provided that \begin{equation}  \label{eqn:alpha}  \alpha = \frac{1-\gamma}{2}\;.\end{equation}Therefore, we can reduce the problem to estimating a singlecoefficient $\gamma$. Another way of expressing this conclusion is torepresent filter $F_9$ in equation (\ref{eqn:l9k}) as a linearcombination of filter $F_5$ from equation (\ref{eqn:l9k}) and itsrotated version \cite[]{cole}, as follows:\begin{equation}\label{eqn:new2}F_9 = \gamma\, \begin{array}{|r|r|r|}\hline0 & 1 & 0 \\ \hline1 & -4  & 1 \\ \hline0 & 1 & 0  \\ \hline \end{array}+ (1-\gamma)\,\begin{array}{|r|r|r|}\hline1/2 & 0 & 1/2 \\ \hline0 & -2  & 0 \\ \hline1/2 & 0 & 1/2  \\ \hline \end{array}\end{equation}With the value of $\gamma=0.5$, filter $F_9$ takes the value\begin{equation}\label{eqn:mc}F_9 = \begin{array}{|r|r|r|}\hline1/4 & 1/2 & 1/4 \\ \hline1/2 & -3 & 1/2 \\ \hline1/4 & 1/2 & 1/4  \\ \hline \end{array}\end{equation}and corresponds precisely to the nine-point McClellan filter\cite[]{mc,GEO56-11-17781785}. On the other hand, the value of$\gamma=2/3$ gives the least error in the vicinity of the zerowavenumber $k$. In this case, the filter is\begin{equation}\label{eqn:my}F_9 = \begin{array}{|r|r|r|}\hline1/6 & 2/3 & 1/6 \\ \hline2/3 & -10/3 & 2/3 \\ \hline1/6 & 2/3 & 1/6  \\ \hline \end{array}\end{equation}Errors of different approximations are plotted in Figure\ref{fig:laplace}\footnote{Another way of constructing  circular-symmetric filters is suggested by the rotated McClellan  transform \cite[]{Biondi.sep.77.27}.}\plot{laplace}{width=6in}{The numerical anisotropy error of different  Laplacian approximations. Both the five-point Laplacian (plot a) and  its rotated version (plot b) are accurate along the axes, but  exhibit significant anisotropy in between at large wavenumbers. The  nine-point McClellan filter (plot c) has a reduced error, while the  filter with $\gamma=2/3$ (plot d) has the flattest error around the  origin.}\par\inputdir{laplace}Under the helix transform, a filter of the general form(\ref{eqn:new}) becomes equivalent to a one-dimensional filter withthe $Z$ transform\begin{eqnarray}\label{eqn:f9z}F_9 (Z) & = & \alpha\,Z^{-N_x-1} + \gamma\,Z^{-N_x} + \alpha\,Z^{-N_x+1} + \gamma\,Z^{-1} -4\,(\alpha + \gamma) \nonumber \\ & & + \gamma\,Z + \alpha\,Z^{N_x-1} + \gamma\,Z^{N_x} + \alpha\,Z^{N_x+1}\;,\end{eqnarray}where $N_x$ is the helix period (the number of grid points in the $x$dimension). To find the inverse of a convolution with filter(\ref{eqn:f9z}), we factorize the filter into the causal minimum-phasecomponent and its adjoint:\begin{equation}\label{eqn:factor}F_9 (Z) = P (Z) P (1/Z)\;.\end{equation}To find the coefficients of the filter $P$, any one-dimensionalspectral factorization method can be applied. It is important to pointout that the result of factorization (neglecting the numerical errors)does not depend on $N_x$. Another approach is to define a residualerror vector for the coefficients of Z in equation (\ref{eqn:factor})and minimize it for some particular norm. For example, minimizing the$\mathcal{L}_1$ norm when $F_9$ is the McClellan filter(\ref{eqn:mc}), we discover that the filter $P$, after transformingback to two dimensions, takes the form\begin{equation}\label{eqn:l1long}\begin{array}{|r|r|r|r|r|r|r|r|}\hline     &         &           &        &     -1.6094 & 0.4293 & 0.05157 & 0.017406 \\ \hline0.01428 & 0.033513 & 0.0808 & 0.2543 & 0.3521  & 0.1553 &   &          \\ \hline\end{array}\end{equation}The results of applying a recursive deconvolution with filter(\ref{eqn:l1long}) are shown in Figure \ref{fig:inv-laplace}.  Anessentially similar procedure, only with a different set of filters,works for implicit wavefield extrapolation.  \plot{inv-laplace}{width=4in,height=4in}{Inverting the Laplacian  operator by a helix deconvolution.  The top left plot shows the  input, which contains a single spike and the causal minimum-phase  filter $P$. The top right plot is the result of inverse filtering.  As expected, the filter is deconvolved into a spike, and the spike  turns into a smooth one-sided impulse. After the second run, in the  backward (adjoint) direction, we obtain a numerical solution of  Laplace's equation! In the two bottom plots, the solution is shown  with grayscale and contours.}%\activeplot{name}{width=6in,height=}{}{caption}%\activesideplot{name}{height=1.5in,width=}{}{caption}%%\begin{equation}%\label{eqn:}%\end{equation}%%\ref{fig:}%(\ref{eqn:})