Exceptional points of radial potentials\\at positive energies
A numerical investigation is conducted for the exceptional points of the positive energy problem $(-\Delta+q_0)u = Eu$, where $E>0$ and $q_0$ is a real-valued, radially symmetric and compactly supported potential. In this case the scattering transform for positive energy is shown to be real-valued and radially symmetric. A new computational method for the Faddeev Green's function for positive energy is introduced. The numerical results show exceptional points for some potentials and supports the only known theoretical result in our case: the absence of exceptional points for small enough potentials and for large enough spectral parameters.
PE_exceptional_points_v9.tex
\documentclass{amsart}
\usepackage{amsmath}
  \usepackage{paralist}
  \usepackage{graphics} %% add this and next lines if pictures should be in esp format
  \usepackage{epsfig} %For pictures: screened artwork should be set up with an 85 or 100 line screen
 \usepackage[colorlinks=true]{hyperref}
\hypersetup{urlcolor=blue, citecolor=red}

  \textheight=8.2 true in
   \textwidth=5.0 true in
    \topmargin 30pt
     \setcounter{page}{1}

% The next 5 line will be entered by an editorial staff.
\def\currentvolume{X}
 \def\currentissue{X}
  \def\currentyear{200X}
   \def\currentmonth{XX}
    \def\ppages{X--XX}


\newtheorem{theorem}{Theorem}[section]
%\theoremstyle{plain}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
%\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
%\numberwithin{equation}{section}


\newcommand{\sgn}{\mbox{sgn}}
\newcommand{\vct}[1]{{\mathbf #1}}       % vector symbol
\newcommand{\mtx}[1]{{\mathrm #1}}       % matrix symbol
\newcommand{\bmat}[2]{ \begin{array}{#1} #2 \end{array}}
\newcommand{\SSS}{{\mathcal{S}}}
\newcommand{\DtoN}{\Lambda_{\gamma}}
\newcommand{\DtoNs}{\Lambda_{\sigma}}
\newcommand{\bound}{\partial\Omega}
\newcommand{\iprod}[1]{\left\langle {#1} \right\rangle}
\newcommand{\abs}[1]{\left\vert {#1} \right\vert}
\newcommand{\norm}[1]{\left\Vert {#1} \right\Vert}
\newcommand{\diff}{\Omega_2\setminus\overline{\Omega_1}}
\newcommand{\pd}[2]{\frac{\partial {#1}}{\partial {#2}}}
\newcommand{\pdpd}[2]{\frac{\partial^2 {#1}}{\partial {#2}^2}}
\newcommand{\Texp}{\mathbf{t}^{\mbox{{\tiny \rm exp}}}}
\newcommand{\Tapp}{\mathbf{t}^{\mbox{{\tiny \rm app}}}}
\newcommand{\gammaapp}{\gamma^{\mbox{\tiny \rm app}}}
\newcommand{\TexpR}{\mathbf{t}^{\mbox{{\tiny \rm exp}}}_{\mbox{{\tiny R}}}}
\newcommand{\TappR}{\mathbf{t}^{\mbox{{\tiny \rm app}}}_{\mbox{{\tiny R}}}}
\newcommand{\TexpB}{\tilde{\mathbf{t}}^{\mbox{{\tiny \rm exp}}}}
\newcommand{\TB}{\mathbf{t}^{\mbox{{\tiny \rm B}}}}
\newcommand{\TBR}{\mathbf{t}^{\mbox{{\tiny \rm B}}}_{\mbox{{\tiny R}}}}
\newcommand{\TLS}{\mathbf{t}^{\mbox{{\tiny LS}}}}
\newcommand{\muLS}{\mu^{\mbox{{\tiny LS}}}}
\newcommand{\muexpR}{\mu^{\mbox{{\tiny \rm exp}}}_{\mbox{{\tiny R}}}}
\newcommand{\muappR}{\mu^{\mbox{{\tiny \rm app}}}_{\mbox{{\tiny R}}}}
\newcommand{\muBR}{\mu^{\mbox{{\tiny \rm B}}}_{\mbox{{\tiny R}}}}
\newcommand{\muexpRL}{\mu^{\mbox{{\tiny \rm exp}}}_{\mbox{{\tiny R}},{\tiny \lambda}}}
\newcommand{\muexp}{\mu^{\mbox{{\tiny \rm exp}}}}
\newcommand{\gammaexpR}{\gamma^{\mbox{{\tiny \rm exp}}}_{\mbox{{\tiny R}}}}
\newcommand{\gammaBR}{\gamma^{\mbox{{\tiny \rm B}}}_{\mbox{{\tiny R}}}}
\newcommand{\gammBR}{\gamma^{\mbox{{\tiny \rm B}}}_{\mbox{{\tiny R}}}}
\newcommand{\gammaexp}{\gamma^{\mbox{{\tiny \rm exp}}}}
\newcommand{\gammaexpRL}{\gamma^{\mbox{{\tiny \rm exp}}}_{\mbox{{\tiny R}},{\tiny \lambda}}}
\newcommand{\TLSR}{\mathbf{t}^{\mbox{{\tiny LS}}}_{\mbox{{\tiny R}}}}
\newcommand{\psiexp}{\psi^{\mbox{{\tiny \rm exp}}}}
\newcommand{\psiB}{\psi^{\mbox{{\tiny \rm B}}}}
\newcommand{\psiLS}{\psi^{\mbox{{\tiny LS}}}}
\newcommand{\Tpsi}{\mathbf{t}^\psi}
\newcommand{\Tx}{{\mathbf{t^{\infty}_x}}}
\newcommand{\TxR}{{\mathbf{t^R_x}}}
\newcommand{\im}{\mbox{Im\,}}
\newcommand{\re}{\mbox{Re\,}}
\newcommand{\dbar}{\bar{\partial}}
\newcommand{\dkb}{\overline{\partial}_k}
\newcommand{\ra}{\rightarrow}
\newcommand{\dbinv}{\overline{\partial}_k^{-1}}
\newcommand{\CC}{{\mathcal{C}}}
\newcommand{\ABgamma}{A_\gamma^{\mbox{{\tiny \rm B}}}}
\newcommand{\Agamma}{A_\gamma}
\newcommand{\SI}{{\mathcal{S}}}
\newcommand{\F}{{\mathcal{F}}}
\newcommand{\G}{{\mathcal{G}}}
\newcommand{\Q}{{\mathcal{Q}}}
\newcommand{\OO}{{\mathcal{O}}}
\newcommand{\B}{{\mathcal{B}}}
\newcommand{\A}{{\mathcal{A}}}
\newcommand{\D}{{\mathcal{D}}}
\newcommand{\T}{{\mathbf{t}}}
\newcommand{\Om}{\Omega}
\newcommand{\DOm}{\partial\Omega}
\newcommand{\fii}{\varphi}
\newcommand{\eps}{\varepsilon}
\newcommand{\R}{{\mathbb R}}
\newcommand{\RR}{{\mathcal R}}
\newcommand{\FF}{{\mathcal F}}
\newcommand{\Z}{{\mathbb Z}}
\newcommand{\C}{{\mathbb C}}
\renewcommand{\L}{{\mathcal L}}
\newcommand{\ol}{\overline}
\newcommand{\bea}[1]{ \begin{eqnarray*} #1 \end{eqnarray*} }
\newcommand{\matr}[1]{\left[\!\begin{array}{l} #1 \end{array}\!\right]}
\newcommand{\odo}[2]{{  \frac{\partial #1}{\partial #2}  }}
\newcommand{\odt}[2]{{  \partial #1/\partial #2  }}
\newcommand{\kk}[1]{{\overline{#1}}}
\newcommand{\texp}{{ \mathbf{t}^{\mbox{exp}} }}
\newcommand{\texpR}{{ \mathbf{t}_{R_T}^{\mbox{exp}} }}
\newcommand{\Ntrig}{{N_{\mbox{{\tiny \sc trig}}}}}
\newcommand{\Nrefine}{{N_{\mbox{{\tiny \sc refine}}}}}
\newcommand{\stulo}[1]{{\left\langle #1 \right\rangle}}       % Sisätulo
\newcommand{\an}{\abs{n}}
\newcommand{\QQ}{{\mathcal{Q}}}
\newcommand{\TT}{{\mathcal{T}}}
\newcommand{\mupm}{\tilde{\mu}_{\pm}}
%\DeclareMathOperator{\round}{round}                        % rounding
%\DeclareMathOperator{\mfrac}{frac}                          % fraction
\def\res{\mathop{Res}\limits}

\newcommand{\dz}{\partial_z}
\newcommand{\dzb}{\overline{\partial}_{z}}

\newcommand{\QQp}{{\mathcal{Q}}}

\newcommand{\qis}{q}
\newcommand{\qnv}{q^{\mbox{\tiny \rm NV}}}

% set unitlength to mm in picture environment
\setlength{\unitlength}{1mm}

%% Place the running title of the paper with 40 letters or less in []
 %% and the full title of the paper in { }.
\title{Exceptional points of radial potentials\\at positive energies}

% Place all authors' names in [ ] shown as running head;
% No more than 40 letters. Leave { } empty
% Please use `and' to connect the last two names if applicable
\author[Janne P.~Tamminen and Samuli Siltanen]{}

% It is required to enter MSC and Keywords.
\subjclass{}
% Please provide minimum  5 keywords. 
 \keywords{}

% Email address of each of all authors is required.
% You may list email addresses of all other authors, separately.
 \email{janne.tamminen@helsinki.fi}
 \email{samuli.siltanen@helsinki.fi}



\begin{document}

\begin{abstract}
A numerical investigation is conducted for the exceptional points of the positive energy problem $(-\Delta+q_0)u = Eu$, where $E>0$ and $q_0$ is a real-valued, radially symmetric and compactly supported potential. In this case the scattering transform for positive energy is shown to be real-valued and radially symmetric. A new computational method for the Faddeev Green's function for positive energy is introduced. The numerical results show exceptional points for some potentials and  supports the only known theoretical result in our case: the absence of exceptional points for small enough potentials and for large enough spectral parameters.
\end{abstract}

\maketitle

\centerline{\scshape Samuli Siltanen and  Janne P. Tamminen}
\medskip


{\footnotesize
 % please put the address of the second  and third author
 \centerline{Department of Mathematics and Statistics}
   \centerline{00014 University of Helsinki, Finland.}
}

\bigskip



\section{Introduction}

\noindent
Application-specific nonlinear Fourier transforms have been used in the analysis of nonlinear evolution equations since the 1970's \cite{Ablowitz1974,Beals1981} and in inverse boundary-value problems since the 1980's \cite{Sylvester1987,Nachman1988}. The definition of such transforms in dimensions two and higher is based on the use of exponentially behaving {\em complex geometric optics (CGO) solutions} first introduced by Faddeev \cite{Faddeev1966} and later rediscovered in the context of inverse problems in \cite{Sylvester1987}. One of the main technical difficulties in the use of CGO solutions is the possibility of {\em exceptional points}, or parameter values for which the solutions are not uniquely defined. 

We present below the first computational examples of exceptional points at positive energies in  dimension two, complementing the earlier work \cite{Music2013} focusing on the zero-energy case. The central novelty is the first algorithm for the evaluation of Faddeev Green's function for nonzero energies, which is a significant extension of the zero-energy case introduced in \cite{Siltanen1999}. Our numerical examples are restricted to the same radial potentials as in \cite{Music2013}, but our numerical methods do apply to general potentials as well.

We discuss solutions $\psi(z,\zeta)$ of the Schr\"odinger equation
\begin{equation}\label{schrode1}
(-\Delta+q_0)\psi(\,\cdot\,,\zeta) = E\psi(\,\cdot\,,\zeta) \quad \mbox{in }\R^2,
\end{equation}
where $q_0:\R^2\ra\R$ is a real-valued compactly supported potential, $\zeta=[\zeta_1\ \zeta_2]^T\in\C^2$ is a vector parameter satisfying $\zeta\cdot\zeta=E$ and the number $E\in \C$ is called {\em energy}. When $|z|\ra\infty$, The solution behaves asymptotically as $\psi(z,\zeta) \sim e^{i\zeta\cdot z}$ meaning that for a suitable Sobolev space $X$ we have
$$
e^{-i\zeta\cdot z}\psi(z,\zeta)-1\in X.
$$
Here $z=[x_1\ x_2]^T\in\R^2$ and
$$
  \zeta\cdot\zeta=(\zeta_1)^2+(\zeta_2)^2,\qquad \zeta\cdot z = x_1\zeta_1+x_2\zeta_2.
$$
Later we will also identify the plane $\R^2$ as the complex plane by writing $z = x_1+ix_2\in\C$. Also, we denote $q=q_0-E$.

Depending on $q_0$, for some $\zeta\in\C^2$ there exists a unique CGO solution $\psi$ as above. Those $\zeta$ for which this is {\em not} the case are called {\em exceptional points} for $q_0$. This happens when the corresponding integral operator is not invertible as presented in section \ref{CGO}. We discuss exceptional points in two cases: zero-energy ($E=0$) in section \ref{sec:ZEB} and positive energy ($E>0$) in section \ref{sec:PEB} and for the remainder of the paper.

In section \ref{SCAT} we show that for a real-valued and radially symmetric potential the scattering transform of positive energy is also real-valued and radially symmetric. The $\dbar$ -equation for positive energy is presented and used in the numerical part of the paper to check the validity of the computed CGO solutions.

The majority of this paper considers the numerical computation of the Faddeev Green's function for positive energy in section \ref{GZETA}. The resulting MATLAB -code is then used in section \ref{NUMERICS} to compute CGO solutions for symmetric, compactly supported potentials of various size. With this numerical investigation we are able to find exceptional points for some potentials, on which there is no theoretical knowledge.

\subsection{The zero-energy case} \label{sec:ZEB}
In this section $E=0$ implying $q=q_0$. As explained in \cite{Nachman1996}, the vectors $\zeta\in\C^2$ satisfying $\zeta\cdot\zeta=0$ are conveniently parameterized as $\zeta=[k\,\,\ \pm ik]^T$, where $k\in\C$ is called the \emph{spectral parameter}. Furthermore, since we restrict ourselves to real-valued potentials, it suffices to use the vectors $\zeta=[k\ ik]^T$. Then $\exp(i\zeta\cdot z)=\exp(ikz)$ with $kz=(k_1+ik_2)(x_1+ix_2)$. We denote the CGO solutions by $\psi(z,k)$ and define the {\em scattering transform} $\T:\C\ra\C$ by 
\begin{equation}\label{zeroenergytk}
\T(k) := \int_{\R^2} \exp(i\overline{k}\overline{z})\,q_0(z)\,\psi(z,k)\,dz,
\end{equation}
where $dz$ denotes the Lebesgue measure on $\R^2$. It is easy to see that replacing $\psi(z,k)$ by $e^{ikz}$ in the right hand side of (\ref{zeroenergytk}) gives the (dilated) Fourier transform of $q_0$. However, since $\psi$ depends on $q_0$ via equation (\ref{schrode1}), we conclude that $\T(k)$ can be viewed as a nonlinear Fourier transform of $q_0$.

We present two applications of the nonlinear Fourier transform (\ref{zeroenergytk}) at zero energy: electrical impedance tomography (EIT)  and Novikov-Veselov  equation.

The inverse conductivity problem of Calder\'on \cite{Calder'on1980} is the mathematical model of EIT. Let $\Om \subset\R^2$ be the unit disc and let $\sigma: \Om\to (0,\infty)$ be a bounded  function satisfying $\sigma(x)\geq c>0$. Let $u\in H^{1}(\Om)$ be the unique solution to
$$
  \nabla\cdot\sigma\nabla u = 0 \mbox{ in }\Om, \qquad  u\big|_{\partial\Om} = \phi.
$$
Calder\'on's problem is to recover the electric conductivity distribution $\sigma$  from the Dirichlet-to-Neumann (DN), or voltage-to-current, map 
\begin{equation}\label{DN-def}
  \Lambda_{\sigma} : \phi \mapsto \sigma\frac{\partial u}{\partial \nu}\Big|_{\partial \Om}.
\end{equation}
Here $\nu$ is the unit outer normal to the boundary. 

In EIT one attaches electrodes to the skin of a patient, feeds electric currents into the body and measures the resulting voltages at the electrodes. Repeating the measurement with several current patterns yields a current-to-voltage data matrix that can be used to compute an approximation to $\Lambda_{\sigma}$. Since different organs and tissues have different conductivities, recovering $\sigma$  amounts to creating an image of the inner structure of the patient. This technique has  commercial realizations, including breast cancer detection \cite{Assenheimer2001} and monitoring lung function\footnote{http://campaigns.draeger.com/pulmovista500/en/}. See \cite{Mueller2012,Cheney1999} for more information on EIT and its applications.

It was shown in \cite{Nachman1996} that $\sigma$ can be reconstructed in two steps. First, use $\Lambda_\sigma$ to determine $\T(k)$ corresponding to the Schr\"odinger potential $q_0=\sigma^{-1/2}\Delta \sigma^{1/2}$ of \eqref{schrode1}. Second, calculate $\sigma$ as the solution of a so-called D-bar equation. See Figure \ref{fig:EITscheme}. This method was implemented numerically and tested with {\em in vivo} human data in \cite{Siltanen2000,Mueller2003,Isaacson2004,Isaacson2006,Knudsen2009,Mueller2012}. Furthermore, there are two similar EIT methods based on nonlinear frequency-domain techniques: \cite{Brown1997}, implemented  in \cite{Knudsen2003,Knudsen2004a,Hamilton2012}, and \cite{Astala2006a}, implemented  in \cite{Astala2006,Astala2010,Astala2011,Mueller2012}. 

\begin{figure}[t]
\unitlength=.3mm
\begin{picture}(380,155)
\put(17,15){\Large$\sigma(z)$}
\put(228,15){\Large$\sigma(z)$}
\put(110,15){\Large$q_0(z)$}
\put(110,125){\Large $\T(k)$}
\put(20,77){\Large$\Lambda_\sigma$}
\thicklines
\put(119,33){\vector(0,1){84}}
\put(122,63){{\footnotesize nonlinear}}
\put(122,53){{\footnotesize Fourier}}
\put(122,43){{\footnotesize transform}}
\put(210,53){{\footnotesize $\overline{\partial}$ equation}}
\put(35,85){\vector(2,1){70}}
\put(140,120){\vector(1,-1){92}}
\put(27,33){\vector(0,1){40}}
\put(6,53){\footnotesize EIT}
\put(5,43){\footnotesize data}
\thinlines
  \multiput(0,107)(4,0){97}{\line(1,0){2}}
  \multiput(0,39)(4,0){97}{\line(1,0){2}}
\put(220,112){\footnotesize Nonlinear frequency domain ($k$-plane)}
\put(278,27){\footnotesize Spatial domain ($z$-plane)}
\end{picture}
\caption{\label{fig:EITscheme}Using nonlinear Fourier transform in the solution of Calder\'on's inverse conductivity problem as in \cite{Nachman1996}. The potential $q_0=\sigma^{-1/2}\Delta \sigma^{1/2}$ depends on the conductivity $\sigma(z)$. The {\em scattering transform} $\T(k)$ of $q_0$, defined in the frequency domain, can be determined from EIT measurements. Furthermore, $\sigma(z)$ can be recovered from $\T(k)$ using a so-called D-bar equation. In practice, noise-robustness is provided by nonlinear low-pass filtering \cite{Knudsen2009}.}
\end{figure}

A crucial point of \cite{Nachman1996} is to prove that there are no exceptional points for \emph{potentials of conductivity type}, defined as follows:
\begin{definition}\label{def:condtype}
  A compactly supported real-valued potential $q_0\in C_0^\infty(\R^2)$ is of \textsl{conductivity type} if $q_0=\sigma^{-1}\Delta\sigma$ for some real-valued $\sigma\in C^\infty(\R^2)$ satisfying $\sigma(z)\geq c>0$ for all $z$ in a bounded set $\Omega\subset\R^2$ and $\sigma(z)=1$ for all $z\in\R^2\setminus\Omega$.
\end{definition}

Let the potential $q(z,\tau)$ depend on $z\in\R^2$ and time $\tau\geq 0$. The  Novikov-Veselov (NV) equation is the following nonlinear evolution equation:
\begin{equation}  \label{evolution}
  \frac{\partial q}{\partial \tau} =  -\dz^3q -\dzb^3q + \frac{3}{4}\,\dz(q\,v) + \frac{3}{4}\,\dzb(q\,\overline{v}),
\end{equation}
where $v(z,\tau)  =  \dzb^{-1}\dz \,q(z,\tau)$ and $\dzb=\frac 12(\frac{\partial}{\partial x_1} + i\frac{\partial}{\partial x_2})$. It was first introduced in a periodic setting  in \cite{Novikov1986,Veselov1984} as a generalization of the (1+1)-dimensional Korteweg-de Vries (KdV) equation \cite{Korteweg1895} to dimension (2+1).

The study of  (\ref{evolution}) in the non-periodic setting via the inverse scattering method was initiated in \cite{Boiti1986,Boiti1987} and continued in \cite{Tsai1989,Tsai1993,Tsai1994}. They discuss the inverse scattering scheme shown in Figure \ref{fig:NVscheme}, analogously to \cite{Gardner1967} in the case of KdV. However, exceptional points are a central difficulty in these papers.

\begin{figure}[t]
\unitlength=.3mm
\begin{picture}(200,100)
\thinlines
 \put(38,75){\vector(1,0){122}}
 \put(-3,72){$\T(k,0)$}
 \put(165,72){$\T(k,\tau)$}
 \put(55,80){$\exp(i \tau(k^3+\overline{k}^3))\cdot$}
 \put(10,15){\vector(0,1){47}}
 \put(-42,37){{\footnotesize Transform}}
 \put(-3,0){$q_0(z,0)$}
 \put(175,62){\vector(0,-1){47}}
 \put(186,37){{\footnotesize Inverse transform}}
 \put(38,3){\vector(1,0){122}}
 \put(46,7){{\footnotesize Nonlinear evolution (\ref{evolution})}}
 \put(165,0){$q_0(z,\tau)$}
\end{picture}
\caption{\label{fig:NVscheme}Direct and inverse nonlinear Fourier transform in the solution of the Novikov-Veselov equation. The temporal evolution of the scattering data is given simply by point-wise multiplication: $\T(k,\tau)=e^{i \tau(k^3+\overline{k}^3)}\T(k,0)$.}
\end{figure}

Building on the absence of exceptional points for conductivity-type potentials in the sense of Definition \ref{def:condtype}, the following series of papers established the well-posedness of the scheme in Figure \ref{fig:NVscheme}: \cite{Lassas2007,Lassas2012a,Lassas2012b,Perry2012}. The natural further question is: do exceptional points for the initial data $q(z,0)$ correspond to soliton solutions of equation (\ref{evolution})? The study of this question needs information on the presence of exceptional points. There are presently two studies of them: the explicit calculations in \cite{Grinevich2012} and the analysis and numerical computations in \cite{Music2013} based on the concept of {\em criticality} as discussed by Murata in \cite{Murata1986}.


\subsection{The positive energy case}\label{sec:PEB}

In the case $\zeta\cdot\zeta=E>0$ we follow \cite{Novikov1992,Grinevich2000} and parameterize the CGO solutions $\psi(z,\lambda)$ by a complex number $\lambda\in\C$, whose relation to the complex vector $\zeta\in\C^2$ is
\begin{equation}\label{lambda-definition}
\lambda = \frac{\zeta_1+i\zeta_2}{\sqrt{E}},\quad \zeta = \matr{(\lambda+\frac{1}{\lambda})\frac{\sqrt{E}}{2} \\ (\frac{1}{\lambda}-\lambda)\frac{i\sqrt{E}}{2}}. 
\end{equation}
The scattering transform for positive energy is defined by 
\begin{equation}\label{positiveenergytk}
  \vct{t}(\lambda) = \int_{\R^2}\exp\Big(\frac{i\sqrt{E}}{2}(\kk{\lambda}z+\kk{z}/\kk{\lambda})\Big)\,q_0(z)\,\psi(z,\lambda)dz.
\end{equation}
Compare (\ref{positiveenergytk}) with (\ref{zeroenergytk}). In the inverse scattering theory this nonlinear Fourier transform has been used to solve the inverse problem of finding the potential $q_0$ from its scattering data.

An analog of EIT in the positive energy case is acoustic imaging (AT). Let $\Om \subset\R^2$ be again the unit disc and let $\rho: \Om\to (0,\infty)$ be a bounded  function satisfying $\rho(x)\geq c>0$. Consider the acoustic wave equation
\begin{equation}\label{wave_eq} 
  -\kappa(z)\frac{\partial^2}{\partial t^2}f(z,t) = \nabla\cdot(\frac{1}{\rho(z)}\nabla f(z,t)),
\end{equation}
where $f(z,t)$ is acoustic pressure, $\kappa(z)$ is compressibility and $\rho(z)$ is the density. The speed of sound is given by $c(z)=(\kappa(z)\rho(z))^{-1/2}$. Substituting the time-harmonic solution $f(z,t)=p(z)e^{i\omega t}$ into (\ref{wave_eq}) leads to 
$$
 \kappa(z)\omega^2p(z) = \nabla\cdot(\frac{1}{\rho(z)}\nabla p(z))\quad\mbox{ in }\Omega.
$$
Furthermore, the change of variables $p = \rho^{1/2}w$ gives 
\begin{equation}\label{schrode-AT}
  (-\Delta + q)w = 0\quad\mbox{ in }\Omega,
\end{equation}
where
\begin{equation}\label{acoustic_q}
 q = \rho^{1/2}\Delta \rho^{-1/2} - \omega^2\kappa\rho.
\end{equation}
Assuming $\rho\vert_{\DOm}=1$ and $q\vert_{\DOm} = -\omega^2\kappa\rho\vert_{\DOm} = -E$ the related CGO solutions satisfy \eqref{schrode1} and we have the positive energy case. For the values of $E$ for which there exists a unique solution of \eqref{schrode-AT} we can define the Dirichlet-to-Neumann map $\Lambda_q$ similarily to \eqref{DN-def}. The inverse acoustic problem is then to reconstruct the density $\rho$ from the knowledge of $\Lambda_q$.

According to \cite[Section 7.1]{Grinevich2000}, the Novikov-Veselov equation (\ref{evolution}) can be considered for initial data $q(z,0)$ that tends to a negative constant as $|z|\ra\infty$, suggesting a positive energy case like above.

We are now in the position to pose intriguing open question of wide interest: 
\begin{enumerate}
\item Under which assumptions on $\omega$, $\rho$ and $\kappa$ in (\ref{acoustic_q}) is it possible to use a scheme such as in Figure \ref{fig:EITscheme} to solve the acoustic inverse problem? Are there exceptional points and what are the other details of this scheme?
\item What kind of positive-energy initial data $q(z,0)$ allows the use of the diagram in Figure \ref{fig:NVscheme} for the solution of the Novikov-Veselov equation? Do positive-energy exceptional points correspond to soliton solutions?
\end{enumerate}
This paper is the first computational insight into positive-energy exceptional points, paving the way to answering the above questions.

\section{Complex Geometrics Optics solutions and exceptional points}\label{CGO}

Consider special solutions of the Schr\"{o}dinger equation \eqref{schrode1} denoted as 
\begin{equation*}%\label{CGO-definition}
\psi(z,\zeta)=e^{i\zeta\cdot z}\mu(z,\zeta),
\end{equation*}
where $\zeta\in\C^2$ and $\Im(\zeta)\neq\vct{0}$. In the zero-energy case the exponential behaviour required is then 
\begin{equation*}%\label{exp-behaviour}
\mu(z,\zeta)-1\in W^{1,p}(\R^2),\quad p>2.
\end{equation*}
We denote this by $\mu(z,\zeta)\sim 1$ or $\psi(z,\zeta) \sim e^{i\zeta\cdot z}$. In the positive energy case it is not yet known what would be the analog, suitable Sobolev space. As seen in section \ref{absence}, in the theory of inverse scattering we require
\begin{equation*}%\label{exp-behaviour}
\mu(z,\zeta)-1\in C(\R^2).
\end{equation*}
If $\zeta\in\R^2$, we have a setting of a physical scattering of a particle with momentum $\zeta$. Only if $\Im(\zeta)\neq\vct{0}$  we call the solutions $\psi(z,\zeta)$ and $\mu(z,\zeta)$ CGO solutions. They do not represent an actual physical scattering, but their asymptotics can be used as additional scattering data. For large $\abs{z}$ we can write $\exp(i\zeta\cdot z)$ in place of $\psi(z,\zeta)$, and the equation \eqref{schrode1} still has to hold giving us
\begin{eqnarray}\label{E-requirement}
  (-\Delta-0) e^{i\zeta\cdot z} = E e^{i\zeta\cdot z} \quad\Rightarrow\quad\zeta\cdot\zeta = E.
\end{eqnarray}
Given the energy $E>0$ we define the spectral parameter by \eqref{lambda-definition}. The unit circle $\abs{\lambda}=1$ divides the complex plane, since $\zeta\rightarrow\Re(\zeta)$ as $\abs{\lambda}\rightarrow 1$. The CGO solutions will have a jump. In place of $\psi(z,\zeta)$ and $\mu(z,\zeta)$ we write $\psi(z,\lambda)$ and $\mu(z,\lambda)$ respectively, even if the energy is then omitted; in the numerical part we only investigate the case $E=1$.

%By transforming to variable $z' = z/\sqrt{E}$ the equation \eqref{schrode1} transforms to
%\begin{eqnarray*}
%(-\Delta_{z'} +q_0(z')\psi(z',\zeta) &=& E\psi(z',\zeta) \\
%(-E\Delta_{z} +q_0(z/\sqrt{E}))\psi(z/\sqrt{E},\zeta) &=& E\psi(z/\sqrt{E},\zeta) \\
%(-\Delta_{z} +\frac{q_0(z/\sqrt{E})}{E})\psi(z/\sqrt{E},\zeta) &=& \psi(z/\sqrt{E},\zeta) \\
%\end{eqnarray*}
%This means that the CGO solution for $q_0$ and $E$ is the same as for $q_0/E$ and $E=1$. Later in the numerical part we assume $E=1$.  given the assumption $E=1$.

The CGO solution $\mu(z,\zeta)$ satisfies another differential equation; starting from \eqref{schrode1} we get
\begin{eqnarray*}
\begin{array}{rcl}
(-\Delta+q_0-E)e^{i\zeta\cdot z}\mu(z,\zeta) &=& 0 \nonumber\\
e^{i\zeta\cdot z}(-\Delta-2i\zeta\cdot\nabla+\zeta\cdot\zeta+q_0-E)\mu(z,\zeta) &=& 0 \nonumber\\
             (L_\zeta+q_0)\mu(z,\zeta) &=& 0, %\label{Lzmu-equation}
\end{array}
\end{eqnarray*}
where
\begin{equation*}%\label{Lz-definition}
L_\zeta := -\Delta-2i\zeta\cdot\nabla.
\end{equation*}
The corresponding Faddeev Green's function $g_\zeta(z)$ satisfies
\begin{equation}\label{gz-requirement}
L_\zeta g_\zeta(z) = \delta_0.
\end{equation}
Turning to $\lambda$-notation, we have
\begin{equation}\label{Llambda-eq}
  (L_\lambda+q_0)\mu(z,\lambda) = 0, 
\end{equation}
where 
\begin{eqnarray}
  \label{Llambda-def}
L_\lambda &:=& -4\partial_z\dbar_z-2i\sqrt{E}(\lambda\partial_z+\frac{1}{\lambda}\dbar_z)\label{Llambda-discrete}.
%&=& -2(\partial_z(\dbar_z+i\lambda\sqrt{E})+\dbar_z(\partial_z+\frac{1}{\lambda}i\sqrt{E})))\nonumber.
\end{eqnarray}
The Green's function of $L_\lambda$ is denoted as $g_\lambda(z)$. Define the Lippman-Schwinger -type equation
\begin{equation}\label{LS-lambda}
  \mu(z,\lambda)-1 = -g_\lambda(z)\ast(q_0(z)\mu(z,\lambda)).
\end{equation}
If $\mu$ solves \eqref{LS-lambda}, then by applying the operator $L_\lambda$ we see it also solves \eqref{Llambda-eq}, so it is the CGO solution of our interest. Define the operators between suitable Banach-spaces $X$ and $Y$, 
\begin{equation*}
\Q:X\to Y,\quad \G_\lambda:Y\to X,
\end{equation*}
by the equations
\begin{eqnarray}
(\Q f)(z) &=& q_0(z)f(z),\label{Q-def}\\
(\G_{\lambda}f)(z) &=& g_\lambda(z)\ast f(z) = \int_{\C}g_\lambda(z-z')f(z')dz'.\label{G-def}
\end{eqnarray}

\begin{definition}
Let $\lambda\in\C\setminus\{0\}$ and $\abs{\lambda}\neq 1$. The spectral parameter $\lambda$ is an \textit{exceptional point}, if $I+\G_\lambda\Q$ is not invertible.
\end{definition}
If $\norm{\G_\lambda\Q}_{\mathcal{L}(X)}<1$, the operator $I+\G_\lambda\Q$ is invertible using the Neumann argument, and from \eqref{LS-lambda} we can solve $\mu$.

\subsection{Absence of exceptional points for small potentials}\label{absence}
The following originates from the work of Grinevich \cite{Grinevich2000} and Novikov \cite{Novikov1992}. Let $E\in\C$, $\lambda\in\C\setminus\{0\}$, then proposition 3.1 in \cite{Novikov1992} states
  \begin{equation}\label{glambda-smallness}
    \abs{g_\lambda(z)} \leq \frac{\hat{c}}{\sqrt{\abs{z}}\abs{E}^{1/4}\sqrt{\abs{\lambda}+1/\abs{\lambda}}},
  \end{equation}
where
\begin{eqnarray*}
    \hat{c} &=& \frac{1}{4}\max_{0<r<\infty}\abs{\sqrt{r}Y_0(r)}+\frac{1}{4\pi}\max_{{0<r<\infty,0\leq\phi_0\leq\pi}}\abs{\sqrt{r}\int_{\phi_0}^{\phi_0+\pi}e^{ir\sin\phi}d\phi}.
\end{eqnarray*}
We assume $q_0(z)$ is smooth, real valued, and satisfies 
\begin{equation}\label{q0-smallness}
\abs{\partial_z^m\dbar_z^nq_0(z)}<\frac{C_{mn}}{(1+\abs{z})^{2+\eps}},\quad m,n\geq0,
\end{equation}
where $C_{mn}$ and $\eps$ are positive constants. Related to this, we write 
\begin{equation*}
    D_\epsilon = \int_{\R^2}\frac{1}{(1+\abs{x})^{2+\epsilon}}dx_1dx_2 = 2\pi\frac{1}{1+\epsilon}\frac{1}{\epsilon}.
\end{equation*}
We consider the operators $\Q$ and $\G_\lambda$ as
\begin{eqnarray*}
\begin{array}{rrcl}
  \Q:& C(\R^2) &\rightarrow& C(\R^2)\cap L^1(\R^2)\\
  \G_\lambda:& C(\R^2)\cap L^1(\R^2) &\rightarrow& C(\R^2).
\end{array}
\end{eqnarray*}
We have
\begin{equation*}
  \norm{\Q(f)}_{L^\infty(\R^2)} = \norm{q_0(z)f(z)}_{L^\infty(\R^2)}\leq C_{00}\norm{f}_{L^\infty(\R^2)},
\end{equation*}
and using \eqref{glambda-smallness}
\begin{equation}\label{Glambda-limit}
  \norm{\G_\lambda(f)}_{L^\infty(\R^2)}\leq\frac{\hat{c}}{\sqrt{\abs{\lambda}+1/\abs{\lambda}}}\left[\norm{f}_{L^1}+\frac{4\pi}{3}\norm{f}_{L^\infty(\R^2)}\right].
\end{equation}
We have 
\begin{eqnarray*}
  \norm{q_0f}_{L^1(\R^2)} &=& \int_{\C}\abs{q_0(z)f(z)}dz \leq \int_{\C}\frac{C_{00}\abs{f(z)}}{(1+\abs{z})^{2+\epsilon}}dz\\
&\leq& \int_{\C}\frac{C_{00}\norm{f}_{L^\infty(\R^2)}}{(1+\abs{z})^{2+\epsilon}}dz = D_\epsilon C_{00}\norm{f}_{L^\infty(\R^2)}.
\end{eqnarray*}
Thus
\begin{eqnarray*}
  \norm{\G_{\lambda}(\Q f)}_{\L^\infty(\R^2)}&\leq& \frac{\hat{c}}{\sqrt{\abs{\lambda}+1/\abs{\lambda}}}\left[\norm{Qf}_{L^1(\R^2)}+\frac{4\pi}{3}\norm{Qf}_{L^\infty(\R^2)}\right]\\
&\leq& \hat{c}\left[D_\epsilon C_{00}\norm{f}_{L^\infty(\R^2)}+\frac{4\pi}{3}C_{00}\norm{f}_{L^\infty(\R^2)}\right]\\
&=& \hat{c}C_{00}(D_\epsilon+\frac{4\pi}{3})\norm{f}_{L^\infty(\R^2)}.
\end{eqnarray*}
By the Neumann argument this leads to the \emph{norm smallness condition}
\begin{equation}\label{norm-smallness}
C_{00} <\frac{1}{\hat{c}(4\pi/3+D_\epsilon)},
\end{equation}
giving us $\norm{\G_\lambda\Q}_{\mathcal{L}(X)}<1$ and unique solvability for $\mu$. In other words, if the potential is small enough according to \eqref{norm-smallness}, there are no exceptional points. This was first proved in Proposition 4.1 in \cite{Novikov1992}. Also if the term $(\abs{\lambda}+1/\abs{\lambda})^{-1/2}$ is included, we have no exceptional points if $\abs{\lambda}$ is large enough.

\subsection{Non-physical scattering transform}\label{SCAT}
Write
\begin{eqnarray}
e_{\lambda}(z)  &:=& e^{\frac{i\sqrt{E}}{2}(\lambda\kk{z}+\kk{\lambda}z+z/\lambda+\kk{z}/\kk{\lambda})}\label{exp-lambda}\\
e_{-\lambda}(z)  &:=& e^{\frac{-i\sqrt{E}}{2}(\lambda\kk{z}+\kk{\lambda}z+z/\lambda+\kk{z}/\kk{\lambda})}\nonumber.
\end{eqnarray}
The non-physical scattering transform is defined by 
\begin{equation}\label{scattering-transform}
\vct{t}(\lambda) = \int_{\C}e_{\lambda}(z)q_0(z)\mu(z,\lambda)dz.
\end{equation}
As first discovered in \cite{Grinevich1986} the CGO solutions $\mu(z,\lambda)$ satisfy the Dbar-equation
\begin{equation}\label{dbar}
  \dbar_\lambda \mu(z,\lambda) = \sgn(\abs{\lambda}-1)\frac{\vct{t}(\lambda)}{4\pi\kk{\lambda}}e_{-\lambda}(z)\overline{\mu(z,\lambda)}.
\end{equation}
The following result is not found in literature. We repeatedly use the uniqueness of the CGO solution in theorem \ref{scat-symmetry} in the following way; let $\psi_1\sim \exp(i\zeta_1\cdot z)$ and $\psi_2\sim \exp(i\zeta_2\cdot z)$ solve 
$$
(-\Delta+q_0)\psi_j = E\psi_j,\quad j=1,2,
$$
and let $\zeta_1\cdot z = \zeta_2\cdot z$. Then we must have $\psi_1 = \psi_2$.

\begin{theorem}\label{scat-symmetry}
  Assume $q_0$ is radially symmetric, real-valued and satisfies the norm-smallness condition \eqref{norm-smallness}. Then $\vct{t}(\lambda)$ is radially symmetric and real-valued.
\end{theorem}
\begin{proof}
Let $\psi(z,\lambda)$ be the unique solution of $(-\Delta+q_0)u = Eu$ with $\psi(z,\lambda)\sim e^{i\zeta\cdot z}$ for the parameter $\lambda$ and the corresponding $\zeta$ given by \eqref{lambda-definition}. We have 
\begin{equation*}
\zeta\cdot z = \frac{\sqrt{E}}{2}[\lambda \kk{z} + \frac{1}{\lambda}z].
\end{equation*}
Consider $\lambda' = e^{i\theta}\lambda$, $\theta\in\R$, and the corresponding $\zeta'$ defined by \eqref{lambda-definition}. We have
\begin{eqnarray*}
  \zeta'\cdot z &=& \frac{\sqrt{E}}{2}[(\lambda'+\frac{1}{\lambda'})x_1+(\frac{1}{\lambda'}-\lambda')ix_2]\\
&=& \frac{\sqrt{E}}{2}[(e^{i\theta}\lambda+e^{-i\theta}\frac{1}{\lambda})x_1+(e^{-i\theta}\frac{1}{\lambda}-e^{i\theta}\lambda)ix_2]\\
&=& \frac{\sqrt{E}}{2}[\lambda(\kk{e^{-i\theta}z}) + \frac{1}{\lambda}(e^{-i\theta}z)]\\
&=& \zeta\cdot (R^{-1}z),
\end{eqnarray*}
where $R^{-1}z$ is the vector $z$ rotated by $-\theta$. The function $\psi(e^{-i\theta}z,\lambda)$ satisfies $\psi(e^{-i\theta}z,\lambda)\sim e^{i\zeta\cdot (R^{-1}z)} = e^{i\zeta'\cdot z}$. Also, because $\Delta$ and $q_0$ are radially symmetric, it satisfies $(-\Delta+q_0)u = Eu$, so in fact $\psi(e^{-i\theta}z,\lambda) = \psi(z,\lambda') =  \psi(z,e^{i\theta}\lambda)$. Thus $\mu(z,e^{i\theta}\lambda) = \mu(e^{-i\theta}z,\lambda)$, and
\begin{eqnarray*}
  \vct{t}(e^{i\theta}\lambda) &=& \int_{\C}e_{e^{i\theta}\lambda}(z)q_0(z)\mu(z,e^{i\theta}\lambda)dz,\\
&=& \int_{\C}e_{\lambda}(e^{-i\theta}z)q_0(e^{-i\theta}z)\mu(e^{-i\theta}z,\lambda)dz\\
&=& \int_{\C}e_{\lambda}(z')q_0(z')\mu(z',\lambda)dz'\\
&=& \vct{t}(\lambda).
\end{eqnarray*}
This proves that $\vct{t}(\lambda)$ is radially symmetric. Let us write $\lambda = r\in\R$, $r\neq0,1$. Then the exponential function \eqref{exp-lambda} satisfies
\begin{eqnarray}
\kk{e_{r}(z)} &=& e_{r}(-z)\label{e-lambda-conj-sym} \\
e_{1/r}(z)      &=& e_{r}(z)\label{e-lambda-sym}\\
e_{r}(\kk{z})      &=& e_{r}(z)\label{e-lambda-sym2}.
\end{eqnarray}
Let $\zeta = [\zeta_1 \ \zeta_2]^{T}$ be the complex vector corresponding to $r$ through \eqref{lambda-definition}, then for $1/r$ we have $\tilde{\zeta} = [\zeta_1,-\zeta_2]^{T}$ as seen from \eqref{lambda-definition}. The unique solution satisfies $\psi(z,1/r)\sim e^{i\tilde{\zeta}\cdot z} = e^{i\zeta\cdot [x_1 \ -x_2]^{T}}$, so $\psi(z,1/r) = \psi(\kk{z},r)$ and $\mu(z,1/r) = \mu(\kk{z},r)$. Along with \eqref{e-lambda-sym},\eqref{e-lambda-sym2} and the symmetry of $q_0$ this leads to
\begin{equation}\label{scat-inverse-symmetry}
\vct{t}(1/r) = \vct{t}(r).
\end{equation}
Given $\zeta$ and $\tilde{\zeta}$ as above, we have $\kk{i\zeta\cdot z} = i\tilde{\zeta}\cdot(-z)$, giving us $\kk{\psi(z,r)} = \psi(-z,1/r)$ and $\kk{\mu(z,r)} = \mu(-z,1/r)$. With \eqref{e-lambda-conj-sym},\eqref{e-lambda-sym},\eqref{scat-inverse-symmetry} and the symmetry of $q_0$ this leads to 
\begin{equation*}
\kk{\vct{t}(r)} = \vct{t}(1/r) = \vct{t}(r).
\end{equation*}
\end{proof}
%\end{theorem}


\section{Computation of the Faddeev Green's function for positive energy}\label{GZETA}

We aim to compute the CGO solution $\mu$ inside the unit circle from the Lippman-Schwinger equation \eqref{LS-lambda}. For that we need a numerical method for $g_\lambda(z)$ for any $\lambda$, $\abs{\lambda}\neq 1$ and any point $z\in D(0,1)$. We look at the Faddeev Green's function with regards to the parameter $\zeta$ to see useful relations. In the case $\zeta\cdot\zeta = E=0$ and $\zeta^I\neq\vct{0}$ the numerical computation of $g_\zeta(z)$ is presented in \cite{Siltanen1999}. We extend the method for $\zeta\cdot\zeta=E>0$. 

Using the Fourier transform, starting from \eqref{gz-requirement}, we have
\begin{eqnarray*}
\begin{array}{rcl}
(-\triangle-2i\zeta\cdot\nabla) g_\zeta &=& \delta_0 \phantom{\Big(}\\
\FF((-\pdpd{}{x_1}-\pdpd{}{x_2}-2i\zeta_1\pd{}{x_1}-2i\zeta_2\pd{}{x_2})g_\zeta) &=& \FF(\delta)\phantom{\Big(}\\
(-(iy_1)^2-(iy_2)^2-2i\zeta_1(iy_1)-2i\zeta_2(iy_2))\hat{g}_\zeta(y) &=& 1/(2\pi)\phantom{\Big(}\\
(y\cdot y+2\zeta\cdot y)\hat{g}_\zeta(y) &=& 1/(2\pi)
\end{array}
\end{eqnarray*}
Using the inverse Fourier transform we get
\begin{equation}\label{green_zeta}
g_\zeta(z) = \frac{1}{(2\pi)^2}\int_{\R^2}\frac{e^{iy\cdot z}}{y\cdot y+2\zeta\cdot y}dy.
\end{equation}
% Using the $\lambda$-notation the same is
% \begin{equation*}
%   \label{glambda}
%   g_\lambda(z) = \frac{1}{\pi^2}\int_{\C}\frac{e^{i(k\kk{z}+\kk{k}z)}}{4k\kk{k}+2\lambda \kk{k}+2k/\lambda}dk.
% \end{equation*}
The following relations can be seen from \eqref{green_zeta}.
\begin{lemma}
Let $\alpha\in \R\setminus\{0\}$, $R$ a rotational matrix with $\det(R) = 1$ and $R\zeta = R\zeta^R+iR\zeta^I$. Then the Faddeev Green's function $g_\zeta(z)$ satisfies
\begin{eqnarray}
g_\zeta(\alpha z) &=& g_{\alpha\zeta}(z)\label{alpha_relation} \\
g_\zeta(Rz) &=& g_{R^{-1}\zeta}(z)\label{rotation_relation} \\
g_\zeta(\matr{-x_1\\x_2}) &=& g_{\matr{-\zeta_1\\ \zeta_2}}(z)\label{x1_relation} \\
g_\zeta(\matr{x_1\\-x_2}) &=& g_{\matr{\zeta_1\\ -\zeta_2}}(z)\label{x2_relation} \\
g_{\matr{\zeta_1\\ \zeta_2}}(\matr{x_1\\x_2}) &=& g_{\matr{\zeta_2\\ \zeta_1}}(\matr{x_2\\x_1})\label{switch_relation} \\
g_{\kk{\zeta}}(z) &=& \kk{g_{\zeta}(-z)}.\label{conjugate_relation}
\end{eqnarray}
\end{lemma}
Using relations \eqref{conjugate_relation} and \eqref{x2_relation} we have
\begin{equation*}
\kk{g_\zeta(-z)} = g_{\kk{\zeta}}(z) = g_\zeta(\matr{x_1\\-x_2}),
\end{equation*}
which results to a switching relation
\begin{equation}\label{x1-symmetry}
g_\zeta(\matr{-x_1\\x_2}) = \kk{g_\zeta(\matr{x_1\\x_2})}
\end{equation}
We can use the rotation relation \eqref{rotation_relation} to reduce $\zeta$ to the form
\begin{equation}\label{reduced-zeta}
\zeta = \matr{k_1\\0}+i\matr{0\\k_2},\quad \abs{k_1}>k_2>0.
\end{equation}
The strategy for computing $g_\lambda(z)$ is now the following:
\begin{enumerate}
\item Use \eqref{lambda-definition} to compute $\zeta$ from $\lambda$.
\item Find the rotational matrix $R$ that satisfies $R\zeta_I = [0,k_2]^{T}$ for some $k_2>0$; then write
  \begin{equation*}
    g_{\zeta}(z) = g_{R^{-1}R\zeta}(z) = g_{R\zeta}(Rz),
  \end{equation*}
  where $R\zeta$ is in the reduced form \eqref{reduced-zeta}.
\item For very small $\abs{z}$ use a method of single layer potential described in section \ref{singlelayer}.
\item Use relation \eqref{alpha_relation} to scale points outwards (since small $\abs{z}$ are a problem as will be seen later) and relation \eqref{x1-symmetry} to switch from $x_1<0$ to $x_1\geq0$.
\item Use computational domains to compute $g_{\zeta}(z)$ for reduced $\zeta$ and $z$ with $x_1\geq 0$.
\end{enumerate}
It takes some analysis to find suitable computational domains for the last step.

\subsection{Analytical integration, part one}

Assume we have the reduced $\zeta$ of \eqref{reduced-zeta} and the switched $z=x_1+ix_2$ with $x_1\geq0$. Write $t = y_1+k_1$, $a = (y_2+k_2i)^2-E$ and subsequently the denominator of the integrand in \eqref{green_zeta} as
\begin{eqnarray*}
y\cdot y + 2y\cdot\zeta &=& y_1^2 + y_2^2 + 2y_1k_1+2y_2k_2i \\
&=& (y_1+k_1)^2+(y_2+k_2i)^2-E \\
&=& t^2+a \\
&=& (t+i\sqrt{a})(t-i\sqrt{a})
\end{eqnarray*}
We define the square root in the same way MATLAB calculates it by default, that is for a complex number $z=r\exp(i\theta),0\leq\theta<2\pi,r\geq0$ the square root is
\begin{equation*}
\sqrt{z} = \left\{
\begin{array}{cl}
\sqrt{r}e^{i\theta/2}&,\quad 0\leq\theta\leq\pi,\\
\sqrt{r}e^{-i(2\pi-\theta)/2}&,\quad \pi<\theta<2\pi.
\end{array}
\right.
\end{equation*}
This way the square root has the following properties: for any $z\in\C$ we have
\begin{eqnarray}
  \Re(\sqrt{z}) &\geq& 0 \label{sqrt-Re}\\
\sqrt{\kk{z}} &=& \kk{\sqrt{z}} \label{sqrt-kk}.
\end{eqnarray}
The numerator of the integrand in \eqref{green_zeta} numerator becomes
\begin{eqnarray*}
e^{iz\cdot y} &=& e^{i(x_1(t-k_1)+x_2y_2)} \\
&=& e^{i(x_2y_2-x_1k_1)}e^{ix_1t}.
\end{eqnarray*}
The integral in \eqref{green_zeta} is thus transformed into
\begin{equation*}
\int_{\R^2}\frac{e^{iy\cdot z}}{y\cdot y+2\zeta\cdot y}d\vct{y} = \int_{-\infty}^{\infty}e^{i(x_2y_2-x_1k_1)}\left(\int_{-\infty}^{\infty}\frac{e^{ix_1t}}{(t+i\sqrt{a})(t-i\sqrt{a})}dt\right)dy_2.
\end{equation*}
The integral over the real parameter $t$ is complexified with $w=w_R+iw_I\in\C$ and we write
\begin{eqnarray}\label{fw}
f(w) = \frac{e^{ix_1w}}{(w+i\sqrt{a})(w-i\sqrt{a})} = \frac{e^{ix_1w_R}e^{-x_1w_I}}{(w+i\sqrt{a})(w-i\sqrt{a})},
\end{eqnarray}
The poles of the function $f(w)$ are $\pm i\sqrt{a}$, where 
\begin{equation*}
a = (y_2+k_2i)^2-E = y_2^2-E-k_2^2 + 2y_2k_2i. 
\end{equation*}
It follows from our definition of $\sqrt{\cdot}$ that
\begin{itemize}
\item when $y_2>0$, $a$ is in the upper half plane, so $i\sqrt{a}$ is in quadrant 2 and $-i\sqrt{a}$ in quadrant 4 (note, that $k_2>0$), and
\item  when $y_2<0$, $a$ is in the lower half plane, so $i\sqrt{a}$ is in quadrant 1 and $-i\sqrt{a}$ in quadrant 3.
\end{itemize}
When $w_I,x_1\geq0$ we have
\begin{equation*}
\abs{f(w)}\rightarrow 0,\quad \mbox{ as }\abs{w}\rightarrow\infty.
\end{equation*}
We choose the integration path 
\begin{equation*}
\Gamma = [-R,R]\cup\{w | w = R\exp(i\theta),0\leq\theta\leq\pi\},
\end{equation*}
so when  $R\rightarrow\infty$ the pole $w=i\sqrt{a}$ is inside the path. Using Residue calculus we get
\begin{eqnarray*}
\int_{\R\subset\C}f(w)dw &=& \int_{\Gamma}f(w)dw = 2\pi i\res_{w=i\sqrt{a}}f(w) \\
&=& 2\pi i\lim_{w\rightarrow i\sqrt{a}}(w-i\sqrt{a})f(w) = 2\pi i\frac{e^{ix_1(i\sqrt{a})}}{i\sqrt{a}+i\sqrt{a}}\\
&=& \pi\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}
\end{eqnarray*}
Thus
\begin{eqnarray*}
g_\zeta(z) &=& \frac{1}{(2\pi)^2}(\int_{-\infty}^{0}e^{i(x_2y_2-x_1k_1)}\pi \frac{e^{-x_1\sqrt{a}}}{\sqrt{a}} dy_2 + \int_{0}^{\infty}e^{i(x_2y_2-x_1k_1)}\pi \frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}dy_2) \\
&=& \frac{1}{4\pi}e^{-ix_1k_1}(\int_{0}^{\infty}e^{-ix_2y_2}\frac{e^{-x_1\sqrt{\kk{a}}}}{\sqrt{\kk{a}}} dy_2 + \int_{0}^{\infty}e^{ix_2y_2}\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}dy_2).
\end{eqnarray*}
Because of \eqref{sqrt-Re} we have
\begin{equation*}
e^{-ix_2y_2}\frac{e^{-x_1\sqrt{\kk{a}}}}{\sqrt{\kk{a}}} = \kk{e^{ix_2y_2}\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}},
\end{equation*}
and so
\begin{equation}\label{gz1}
g_\zeta(z) = \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_{0}^{\infty}e^{ix_2t}\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}dt),\quad x_1\geq0,
\end{equation}
where $a = (t+k_2i)^2-E$ and we have changed the symbol $y_2$ to $t$.

\subsubsection{Observations}\label{observations}

The starting point in our numerical experiments is the parameter $\lambda$ with $\abs{\lambda}>1$. Using \eqref{lambda-definition} we can see, that as $\abs{\lambda}\rightarrow 1$ we have $k_2\rightarrow0$. When $\abs{\lambda}\rightarrow\infty$, we have $\abs{k_1},k_2\rightarrow\infty$.

We want to numerically compute the integral of \eqref{gz1}. We can only compute up to a finite limit, say from $0$ to $N$. Two problems might occur, the integrand either converges slowly, meaning we have to take $N$ very large, or the integrand might oscillate fast, meaning we have to take a great number of integration points in $[0,N]$.

Taking these considerations into account, we see from \eqref{gz1} that problematic situations arise when $x_2$ is large (oscillation), $x_1$ is small (convergence) and when $k_2$ is large (oscillation). When $x_1$ is large we have oscillation, but also better convergence. Also, using numerical evidence by plotting the integral with several values of $k_2$, it seems some large difficulties arise when $k_2$ is close to zero; this means problems when $\lambda$ is close to one.

This leads to the following enhancements. The different computational domains will use different versions of \eqref{gz1}.

\subsection{Analytical integration, part two}
Write
\begin{equation*}
  g(w) = e^{ix_2 w}\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}},\quad a = (w+k_2i)^2-E,\quad w = w_1+iw_2\in\C,
\end{equation*}
and consider the complexified integral $\int_{\R^+}g(w)dw$ of \eqref{gz1}. We have
\begin{equation*}
  e^{ix_2 w} = e^{ix_2w_1}e^{-x_2w_2},
\end{equation*}
so when $x_2,w_2\geq0$ or $x_2,w_2<0$ we have 
\begin{equation}\label{gz-convergence}
\abs{g(w)}\rightarrow0\mbox{ as }\abs{w}\rightarrow\infty. 
\end{equation}
Note that the numerator $\exp(-x_1\sqrt{a})$ converges, since $x_1\geq0$ and \eqref{sqrt-Re} holds. The branch points of $g(w)$ are $\pm\sqrt{E}-k_2i$. 

\subsubsection{The case $x_2\geq 0$}
We choose the integration path
\begin{equation*}
\Gamma_1 = [0,R]\cup\{w | w = R\exp(i\theta),0\leq\theta\leq\pi/2\}\cup\{iR(1-t),\quad 0\leq t\leq 1\}.
\end{equation*}
We have $\int_{\Gamma_1}g(w)dw = 0$, since $g(w)$ is analytic in the first quadrant. Because of \eqref{gz-convergence} we then have
\begin{equation*}
  \int_{\R^+}g(w)dw = -\int_{\infty}^0g(it)idt = \int_0^\infty e^{-x_2 t}\frac{e^{-x_1 i\sqrt{(t+k_2)^2+E}}}{\sqrt{(t+k_2)^2+E}}dt,
\end{equation*}
and so 
\begin{eqnarray}
%g_\zeta(z) &=& \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_0^\infty e^{-x_2 t}\frac{e^{-x_1 i\sqrt{(t+k_2)^2+E}}}{\sqrt{(t+k_2)^2+E}}dt)\label{pre-gz2}.
&=& \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_0^\infty e^{-x_2 t}\frac{e^{-x_1 i\sqrt{t^2+k_1^2+2tk_2}}}{\sqrt{t^2+k_1^2+2tk_2}}dt)\label{gz2}.
\end{eqnarray}
The integrand in \eqref{gz2} converges quickly for large $x_2$ and has high oscillation for large $x_1$.

\subsubsection{The case $x_2<0$}
The branch point $\sqrt{E}-k_2i$ is avoided by integrating along 
\begin{eqnarray*}
\Gamma_2 &=& [0,\sqrt{E}+1]\cup(\{-iRt:0\leq t\leq 1\}\cup\{w = R\exp(i\theta):\frac{3}{2}\pi\leq\theta\leq 2\pi\}\\
& &\cup\{(\sqrt{E}+1)t+(1-t)R: 0\leq t\leq 1\})\\
&:=& P_1\cup (L_1\cup L_2\cup L_3),
\end{eqnarray*}
where $g(w)$ is analytical inside the loop $L_1\cup L_2\cup L_3$, and $\abs{g(w)}\rightarrow 0$ on the circle $L_2$ as $R\rightarrow\infty$. Thus
\begin{equation*}
  \int_{\R^+}g(w)dw = I_1 + \int_0^\infty g(\sqrt{E}+1-it)(-i)dt,
\end{equation*}
where
\begin{eqnarray*}
\begin{array}{rcl}
  I_1 &=& \int_0^{\sqrt{E}+1}g(t)dt \\
 \int_0^\infty g(\sqrt{E}+1-it)(-i)dt &=& \int_0^\infty e^{i(\sqrt{E}+1)x_2}e^{x_2 t}\frac{e^{-x_1\sqrt{b}}}{\sqrt{b}}(-i)dt,
\end{array}
\end{eqnarray*}
and $b= (\sqrt{E}+1+(k_2-t)i)^2-E$. In total 
\begin{eqnarray}
g_\zeta(z) &=& \frac{1}{2\pi}e^{-ix_1k_1}\Re(I_1-ie^{i(\sqrt{E}+1)x_2}\int_0^\infty e^{x_2 t}\frac{e^{-x_1\sqrt{b}}}{\sqrt{b}}dt)\label{gz3}
\end{eqnarray}
The integrand in \eqref{gz3} converges quickly for large $\abs{x_2}$ and has high oscillation for large $x_1$.

\subsection{Choosing upper limits for the integrals}
Write $g_\zeta^{T_1}$,$g_\zeta^{T_2}$ and $g_\zeta^{T_3}$ for the finite integrals for \eqref{gz1},\eqref{gz2} and \eqref{gz3} respectively. We need to choose the upper limits $T_i$, $i=1,2,3$. There will be numerical error caused by the neglected part of the integral and the numerical integration method used.  It is decided to simply require
\begin{equation}\label{error-requirement}
  \abs{g_\zeta-g_\zeta^{T_i}} < 1E-8.
\end{equation}
The error caused by the numerical integration method is assumed to be not dependant on $\lambda$ or $z$. For $g_\zeta^{T_i}$, the integration range $[0,T_i]$ is divided into $M_i$ points (with $g_\zeta^{T_3}$ there is also the additional integral $I_1$) and the Gaussian quadrature is used. The integers $M_i$ are chosen large enough so that for any integer $M>M_i$ the first 8 digits are not changing in the numerical value of $g_\zeta^{T_i}(z)$. In this test the choice of $z$ has only a minor effect, it is done by choosing the ``worst'' possible point for any given computational domain; for example, for $g_\zeta^{T_1}(z)$ using $z = [1,1]$ the integrand has more oscillation than with the point $z=[1,0]$, and thus needs a larger parameter $M_1$.

Finding $T_i$ is a bit cumbersome and is explained in the following. In \eqref{gz1} we have the term $a$ and
\begin{eqnarray}
  \abs{\sqrt{a}} &=& \abs{\sqrt{(t^2-k_2^2-E)+2k_2ti}} = ((t^2-k_1^2)^2+4k_2^2t^2)^{1/4}\nonumber\\
&=& (t^4-2t^2k_1^2+k_1^4+4k_2^2t^2)^{1/4}\nonumber\\
&\geq& (t^4-2k_1^2t^2)^{1/4} \geq (t^4-1/2t^4)^{1/4}\nonumber\\
&=& \frac{t}{2^{1/4}} \label{gz1-ul},
\end{eqnarray}
when $t\geq 2\abs{k_1}$. Then, writing $\theta$ for the angle $\sqrt{a} = r\exp(i\theta)$,
\begin{eqnarray*}
  \abs{\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}} &\leq& \frac{e^{-x_1\Re(\sqrt{a})}}{t/2^{1/4}} = \frac{e^{-x_1\cos(\theta)\abs{\sqrt{a}}}}{t/2^{1/4}} \\
&\leq& \frac{e^{-x_1\cos(\theta)t/2^{1/4}}}{t/2^{1/4}}.
\end{eqnarray*}
The angle goes to zero as $t\rightarrow\infty$, so $\cos(\theta)\rightarrow 1^-$. Since $t\geq 2\abs{k_1}$ we write $c_1=\cos(\theta_1)$, where the angle of $\sqrt{a}|_{t=2k_1}$ is $\theta_1$, and so we have for the integral
\begin{eqnarray}
 \abs{\int_{T}^\infty e^{ix_2t}\frac{e^{-x_1\sqrt{a}}}{\sqrt{a}}dt} &\leq& \int_{T}^\infty \frac{e^{-x_1c_1t/2^{1/4}}}{t/2^{1/4}}dt \nonumber\\
&=& 2^{1/4}\int_{x_1c_1 T/2^{1/4}}^\infty\frac{e^{-s}}{s}ds = 2^{1/4}E_i(x_1c_1 T/2^{1/4}) \label{remainder}.
\end{eqnarray}
The exponential integral function $E_i$ can be computed in MATLAB with ''expint''. Because of \eqref{error-requirement} we require that the remainder \eqref{remainder} is of the order $2\pi/2^{1/4}\cdot1E-8 \approx 7.47E-8$. We can test with MATLAB that $Ei(14)<6E-8$, so we get
\begin{eqnarray*}
\begin{array}{rcl}
  x_1c_1T/2^{1/4} &=& 14\\
T &=& \frac{14\cdot 2^{1/4}}{x_1c_1}.
\end{array}
\end{eqnarray*}
Finally,
\begin{equation}\label{T1-def}
T_1 = \max\{\frac{14\cdot 2^{1/4}}{x_1c_1},2k_1\}.
\end{equation}

From \eqref{gz2} we easily get
\begin{eqnarray*}
  \abs{\int_{T_2}^\infty e^{-x_2t}\frac{e^{-x_1i\sqrt{t^2+k_1^2+2tk_2}}}{\sqrt{t^2+k_1^2+2tk_2}}dt} &\leq& \int_{T_2}^\infty\frac{e^{-x_2 t}}{t}dt \\
&=& Ei(x_2T_2).
\end{eqnarray*}
Thus the upper limit can be computed, as before, from
\begin{eqnarray}
  x_2T_2 &=& 14\nonumber\\
T_2 &=& \frac{14}{x_2}\label{T2-def}.
\end{eqnarray}

Starting from \eqref{gz3} we have
\begin{eqnarray*}
  \abs{\sqrt{b}} &=& ((k_2-t)^4+2(2\sqrt{E}+1)(k_2-t)^2+(2\sqrt{E}+1)^2\\
& &+4(\sqrt{E}+1)^2(k_2-t)^2)^{1/4}\\
&\geq& ((k_2-t)^4)^{1/4} = \abs{t-k_2}.
\end{eqnarray*}
Using the same argument as preceding \eqref{remainder}, we write $c_2 = \cos(\theta_2)$, where the angle of $\sqrt{b}|_{t=k_2}$ is $\theta_2$. Then for $T_2> k_2$ we have
\begin{eqnarray*}
 \abs{\int_{T_3}^\infty e^{x_2t}\frac{e^{-x_1\sqrt{b}}}{\sqrt{b}}dt} &\leq& \int_{T_3}^\infty \frac{e^{x_2t-x_1c_2(t-k_2)}}{t-k_2}dt \\
&=& e^{x_2k_2}\int_{(c_2x_1-x_2)(T_3-k_2)}^\infty\frac{e^{-s}}{s}ds \\
&=& e^{x_2k_2}Ei((c_2x_1-x_2)(T_3-k_2))\\
&\leq& Ei((c_2x_1-x_2)(T_3-k_2)).
\end{eqnarray*}
As before,
\begin{eqnarray}
  (c_2x_1-x_2)(T_3-k_2) &=& 14\nonumber \\
T_3 &=& \frac{14}{c_2x_1-x_2}+k_2\label{T3-def}.
\end{eqnarray}
Using these choices of upper limits guarantees us
\begin{equation*}
  \abs{g_\zeta(z)-g_\zeta^{T_i}(z)} <1E-8,\quad i=1,2,3.
\end{equation*}

\subsection{Use of single-layer potential for small $z$}\label{singlelayer}
For small values of $z$ there is a problem of slow convergence. We will evade this problem by the use of the single-layer potential for a function that satisfies the Helmholtz equation.

Write 
$$
E = k^2,\quad G_\zeta(z) = \exp(i\zeta\cdot z)g_\zeta(z),\quad G(z) = iH_0^{1}(k\abs{z})/4,
$$
where $H_0^1$ is Hankel's function of the first type. We have
\begin{eqnarray*}
(-\triangle-k^2)G_\zeta(z) &=& \delta_z \\
(-\triangle-k^2)G(z) &=& \delta_z ,
\end{eqnarray*}
so
\begin{equation*}
(-\triangle-k^2)(G_\zeta - G) = 0.
\end{equation*}
Write $H_\zeta:=G_\zeta-G$. For any radius $R$ there exists a single-layer potential $p(z)$, which gives the value of $H_\zeta$ by the integral
\begin{equation}\label{layer-potential}
H_\zeta(z) = \int_{\partial D(0,R)}\frac{i}{4}H_0^{1}(k\abs{z-y})p(y)d\mu(y) := S(p(\cdot))(z).
\end{equation}
Assume we know $H_\zeta(z)$ on the circle $\partial D(0,R)$, where $R$ is large enough so that we don't have the problems of slow convergence. The potential can be recovered by the inverse of the integral operator, $p = S^{-1}(H_\zeta(z))$. Then $H_\zeta(z)$ can be calculated using \eqref{layer-potential} for any $\abs{x}<R$. Finally we have
\begin{equation}\label{gzeta-layerpotential}
g_\zeta(z) = e^{-i\zeta\cdot z}(H_\zeta(z) + G(z)).
\end{equation}

% \subsubsection{Observations}
% The more the integrand of \eqref{gz1} oscillates the larger $x_2$ is, making a simple ''brute force'' approach ineffective. However, when $x_2$ is small and $x_1$ is large, the nominator $\exp(-x_1\sqrt{w})$ converges quickly, since $\Re(\sqrt{w})$ is the dominating term as $t\rightarrow\infty$. In an opposite fashion equations \eqref{pre-x2p} and \eqref{gzeta-x2n} converge quickly, when $\abs{x_2}$ is large, but the more they oscillate the larger $x_1$ is. In all of the equations \eqref{gz1},\eqref{pre-x2p} and \eqref{gzeta-x2n} the convergence is slow if $\abs{x}$ is small.

\subsection{Computational domains and the computation of $g_\zeta(z)$}
The plane $\R^2$ is divided into computational domains as shown in figure \ref{areas}. The disks shown are $D_1 = D(0,0.2)$, $D_2 = D(0,0.5)$, $D_3 = D(0,1)$, and not in scale is pictured $D_4 = D(0,2.5)$. Domain 1a is the disk $D_1$, domain 1b is the annulus $D_2\setminus D_1$, domain 1c is the annulus $D_3\setminus D_2$. Domains 2,3,4,5,6 and 7 form the annulus $D_4\setminus D_3$.

Domains 2 and 3 are divided by the line $x_2 = 0.5x_1$, domains 3 and 4 by $x_1=0$, domains 4 and 5 by $x_2=-0.5x_1$, domains 5 and 6 by $x_2=0.5x_1$, domains 6 and 7 by $x_1=0$ and finally domains 7 and 2 by $x_2=-0.5x_1$.

\begin{figure}[h]
\unitlength=1mm
\begin{picture}(120,140)
\epsfxsize=12cm
\put(5,5){\epsffile{01-areas-v2.eps}}
\put(60,109){$x_2$}
\put(105,65){$x_1$}
\put(65,64){1a}
\put(65,55){1b}
\put(65,41){1c}
\put(112,67){2}
\put(16,67){5}
\put(89,107){3}
\put(89,27){7}
\put(42,107){4}
\put(42,27){6}
\end{picture}
\caption{\label{areas}Computational domains}
\end{figure}
For the reduced $\zeta$ we now have the equations
\begin{equation}
g_\zeta^{T_1}(z) = \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_{0}^{T_1}e^{ix_2t}\frac{e^{-x_1\sqrt{t^2+2tk_2i-k_1^2}}}{\sqrt{t^2+2tk_2i-k_1^2}}dt),\label{gz-T1}
\end{equation}
\begin{equation}
g_\zeta^{T_2}(z) = \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_0^{T_2} e^{-x_2 t}\frac{e^{-x_1 i\sqrt{t^2+2tk_2+k_1^2}}}{\sqrt{t^2+2tk_2+k_1^2}}dt)\label{gz-T2},\quad x_2\geq0,
\end{equation}
\begin{eqnarray}
g_\zeta^{T_3}(z) &=& \frac{1}{2\pi}e^{-ix_1k_1}\Re(\int_0^{\sqrt{E}+1}e^{ix_2t}\frac{e^{-x_1\sqrt{t^2+2tk_2i-k_1^2}}}{\sqrt{t^2+2tk_2i-k_1^2}}dt\label{gz-T3}\\
& &-ie^{i(\sqrt{E}+1)x_2}\int_0^{T_3}e^{x_2 t}\frac{e^{-x_1\sqrt{b}}}{\sqrt{b}}dt),\quad x_2<0,\nonumber
\end{eqnarray}
where $b=(\sqrt{E}+1+(k_2-t)i)^2-E$. In general the point $z$ lies in one of the computational domains;
\begin{itemize}
\item In domain 1a, we use the single-layer potential and the equation \eqref{gzeta-layerpotential}.
\item In domain 1b, we scale the point $z$ to the annulus $D_4\setminus D_3$ using \eqref{alpha_relation}, 
  \begin{equation*}
    g_\zeta(z) = g_\zeta(1/5\cdot5x) = g_{\zeta/5}(5x)
  \end{equation*}
(note that the energy $E$ changes via this scaling transformation).
\item In domain 1c we do the same as above with the scaling factor 2.
\item In domain 2 we use \eqref{gz-T1}, since $\abs{x_2}$ is small and $x_1$ is large. The upper limit $T_1$ is computed from \eqref{T1-def}.
\item In domain 3 we use \eqref{gz-T2}, since $x_1$ is small and $x_2>0$ is large. The upper limit $T_2$ is computed from \eqref{T2-def}.
\item In domain 7 we use \eqref{gz-T3}, since $x_1$ is small, $x_2<0$ and $\abs{x_2}$ is large. The upper limit $T_3$ is computed from \eqref{T3-def}.
\item In domains 4,5,6 we use \eqref{x1-symmetry} to switch them to domains 3,2,7 respectively.
\end{itemize}
A sample of the function $g_\lambda(z)$ is pictured in \ref{fig:gzeta}, in $400\times400$ -grid of points $z$, $\lambda=1+i$, $E=1$.

\begin{figure}
\unitlength=1mm
\begin{picture}(120,60)
\epsfxsize=6cm
\put(0,0){\epsffile{02-gzeta-real.eps}}
\epsfxsize=6cm
\put(62,0){\epsffile{02-gzeta-imag.eps}}
\put(90,49){\small $\Im(g_\lambda(z))$}
\put(28,49){\small $\Re(g_\lambda(z))$}
\end{picture}
\caption{\label{fig:gzeta}The real and imaginary parts of $g_{\lambda}(z)$ in $400\times400$ grid of points $z$, $\lambda=1+i$, $E=1$.}
\end{figure}


\subsubsection{Domain 1a and the single-layer method explained in more detail}
We take two radii $R=1.6$ and $\epsilon=0.1$, and divide $[0,2\pi)$ into 256 equally distanced points forming the vector $\vct{\theta}\in\R^{256}$. We write
\begin{eqnarray*}
  \vct{y_1} &=& Re^{i\vct{\theta}}\in\C^{256}\\
  \vct{y_2} &=& (R+\epsilon)e^{i\vct{\theta}}\in\C^{256},
\end{eqnarray*}
where any complex element $y_i$ refers to the point $[\Re(y_i) \ \Im(y_i)]^{T}$. Using the computational method for domains 2,3,4,5,6,7  we calculate $g_\zeta(\vct{y_1})$, and subsequently $H_\zeta(\vct{y_1}) = G_\zeta(\vct{y_1})-G(\vct{y_1})$. In order to avoid singularities in the integral operator \eqref{layer-potential} we aim to calculate the single-layer potential $p(z)$ on the circle $\partial D(0,R+\epsilon)$. The integral operator can be written as a matrix,
\begin{equation*}
  S_p = d\theta (R+\epsilon)\frac{i}{4}\matr{H_0^1(\sqrt{E}\abs{y_1\vct{1}^T-\vct{y_2}^{T}})\\H_0^1(\sqrt{E}\abs{y_2\vct{1}^T-\vct{y_2}^{T}})\\\vdots\\H_0^1(\sqrt{E}\abs{y_{256}\vct{1}^T-\vct{y_2}^{T}})},
\end{equation*}
where $d\theta(R+\epsilon)$ comes from the numerical integration on the circle $\partial D(0,R+\epsilon)$, $d\theta=2\pi/256$, $y_i$ is the i:th component of $\vct{y_1}$ and $\vct{1}\in\R^{256}$ is a vector of ones. Now, the equation \eqref{layer-potential} in vector-form is
\begin{equation*}
H_\zeta(\vct{y_1}) = S_pp(\vct{y_2}).  
\end{equation*}
We use GMRES to compute
\begin{equation*}
  p(\vct{y_2}) = S_p^{-1}H_\zeta(\vct{y_1}).
\end{equation*}
Let us be given points $z_j$, $1\leq j\leq n$, in domain 1a. Using the complex form $z_j\in\C$ the matrix form of the integral operator of \eqref{layer-potential} is
\begin{equation*}
    S = d\theta (R+\epsilon)\frac{i}{4}\matr{H_0^1(\sqrt{E}\abs{z_1\vct{1}^T-\vct{y_2}^{T}})\\H_0^1(\sqrt{E}\abs{z_2\vct{1}^T-\vct{y_2}^{T}})\\\vdots\\H_0^1(\sqrt{E}\abs{z_{n}\vct{1}^T-\vct{y_2}^{T}})},
\end{equation*}
and so for the complex vector $z = [z_1 \ z_2 \ \cdots \ z_n]^{T}$ we have
\begin{equation*}
H_\zeta(z) = Sp(\vct{y_2})
\end{equation*}
and finally
\begin{equation}\label{gzeta-D1}
g_\zeta(z_j) = e^{-i\zeta\cdot z_j}(H_\zeta(z_j) + G(z_j)),\quad 1\leq j\leq n.
\end{equation}

\section{Numerical results}\label{NUMERICS}

For any $\lambda$, the corresponding CGO solution $\mu(z,\lambda)$ can be computed from the Lippman-Schwinger equation \eqref{LS-lambda} using periodic Sobolev spaces wrapped in a torus \cite{Vainikko2000}. The numerical solver uses $2^M\times 2^M $ grid points of $z$ in which the CGO solution is presented. 

We now fix $E=1$. We have two numerical tests for the new numerical method for $g_\lambda(z)$. First, test the computed CGO solutions by applying a discrete version of the $\dbar$ -equation \eqref{dbar}. Second, compute the scattering transform for various real-valued, radially symmetric potentials in search for exceptional points. As was seen in section \ref{SCAT}, for these type of potentials the scattering transform is also radially symmetric and real-valued. Also we have $\vct{t}(1/\lambda) = \vct{t}(\lambda)$, meaning we only need to take values $\lambda>1$, $\lambda\in \R$.

\subsection{Definition of potentials}
We use exactly the same potentials as in the numerical part of \cite{Perry2012}. Take radii $0<R_1<R_2<1$ and a polynomial $\tilde{p}(t) = 1-10t^3+15t^4-6t^5$. Set for $R_1\leq t \leq R_2$
\begin{equation*}
  p(t) = \tilde{p}(\frac{t-R_1}{R_2-R_1}).
\end{equation*}
Then, the approximate test function 
\begin{equation}
  \varphi(\abs{z}) = 
\left\{
\begin{array}{ccl}
1          & \mbox{ for } &0\leq \abs{z}\leq R_1 \\
p(\abs{z}) & \mbox{ for } &R_1< \abs{z}< R_2 \\
0          & \mbox{ for } &R_2\leq \abs{z}\leq 1,
\end{array}
\right.
\end{equation}
is in $C^2$. The values $R_1=0.8$ and $R_2=0.9$ were used.

We consider the potentials
\begin{eqnarray}
\begin{array}{rcll}
  q_{\alpha}^1 &=& \alpha \varphi,\quad &\alpha\in\R,\label{qalpha1-def}\\
  q_{\alpha}^2 &=& \frac{\Delta\sqrt{\sigma}}{\sqrt{\sigma}}+\alpha\varphi,\quad &\alpha\in\R,\label{qalpha2-def},
\end{array}
\end{eqnarray}
where $\sigma\in C^2(\Omega), \sigma\geq c>0$. The conductivity $\sigma$ and the conductivity-type potential $q_0^2$ are pictured in figure \ref{fig:sigmapoten}.

\begin{figure}
\unitlength=1mm
\begin{picture}(120,150)
 \put(0,75){\includegraphics[width=6cm]{sigma2D.eps}}
\put(0,140){\small Conductivity $\sigma(z)$}
 \put(0,0){\includegraphics[width=6cm]{poten2D.eps}}
\put(0,65){\small Conductivity-type potential $q_0^2(z)$ }
 \put(65,75){\includegraphics[width=5.5cm]{sigma1D.eps}}
\put(65,140){\small Profile of conductivity }
\put(62,132){\small $\sigma(|z|)$}
 \put(65,0){\includegraphics[width=5.5cm]{poten1D.eps}}
\put(65,65){\small Profile of potential }
\put(62,57){\small $q_0^2(|z|)$}
\put(118,-3){\small $|z|$}
\put(118,72){\small $|z|$}
\end{picture}
\caption{\label{fig:sigmapoten}Top row: mesh plot and profile plot of the rotationally symmetric conductivity $\sigma(z)=\sigma(|z|)$. Bottom row: mesh plot and profile plot of the resulting conductivity-type potential $q_0^2(z)=q_0^2(|z|)$.}
\end{figure}

Since $\norm{\varphi}_{L^\infty(\R^2)}=1$, using the norm-smallness condition \eqref{norm-smallness} to $q_\alpha^1$ from \eqref{qalpha1-def}, we have no exceptional points when
\begin{equation*}
\abs{\alpha} <\frac{1}{\hat{c}(4\pi/3+D_\epsilon)}.
\end{equation*}

\subsection{Verification of the computed CGO solutions}

To verify that the CGO solutions are correct, we test the $\dbar$ -equation \eqref{dbar} using the five-point stencil method with the finite difference $d\lambda=0.0001$. Take parameters $\lambda$ from 1.01 to 4.5. Choose the potentials $q_0^2$ and $q_{35}^2$. For each $\lambda=\lambda_1+\lambda_2 i$, compute

\begin{enumerate}
\item The CGO solution $\mu_0$ in the z-grid.
\item The CGO solutions $\mu_1,\mu_2,\mu_3,\mu_4,\mu_5,\mu_6,\mu_7$ and $\mu_8$ using $\lambda+d\lambda$, $\lambda+2d\lambda$, $\lambda-d\lambda$, $\lambda-2d\lambda$, $\lambda+d\lambda i$, $\lambda+2d\lambda i$, $\lambda-d\lambda i$ and $\lambda-2d\lambda i$ respectively.
\item The functions $e_\lambda(z)$ and $e_{-\lambda}(z)$.
\item The potential $q_0(z)$ and the scattering transform $\vct{t}(\lambda)$ with $\mu_0$.
\item The derivatives and the $\dbar$ -operation by
\begin{eqnarray*}
\partial_{\lambda_1}\mu &=& \frac{-\mu_2+8\mu_1-8\mu_3+\mu_4}{12d\lambda}\\
\partial_{\lambda_2}\mu &=& \frac{-\mu_6+8\mu_5-8\mu_7+\mu_8}{12d\lambda}\\
\dbar \mu &=& \frac{1}{2}(\partial_{\lambda_1}+i\partial_{\lambda_2})\mu.
\end{eqnarray*}
\item The error
  \begin{equation}\label{dbar-error}
    \norm{\dbar \mu-\frac{1}{4\pi\kk{\lambda}}\vct{t}(\lambda)e_{-\lambda}(z)\kk{\mu_0}}_{L^2(D(0,1))}.
  \end{equation}
\end{enumerate}
In figure \ref{fig:dbar-error} we see the error \eqref{dbar-error} as a function of $\lambda$ using $q_0^2$ on the left, $q_{35}^2$ on the right. The parameter $M$ is increased from 7 to 9. As expected, the error decreases as $M$ increases as it increases the accuracy of the LS -solver. The smallest values of $\lambda$ were omitted in the pictures, for $\lambda=1.01$ the magnitude of the error was between 3 and 13, for the second smallest $\lambda$ it was between 0.003 and 0.02. For values of $\lambda$ near $\abs{\lambda}=1$ the numerical method of $g_\lambda(z)$ has great error due to very small value of $k_2$, see section \ref{observations}.

\begin{figure}
\begin{picture}(120,60)
\epsfxsize=6cm
\put(1,0){\epsffile{03-Dbar_test_p1.eps}}
\epsfxsize=6cm
\put(63,0){\epsffile{03-Dbar_test_p2.eps}}
\put(123,-3){\small $\lambda$}
\put(0,54){\small $\norm{\cdot}_{L^2}$}
\end{picture}
\caption{\label{fig:dbar-error}Errors in the $\dbar$ -equation for two different potentials and different accuracies of the LS -solver; horizontal axis is the spectral parameter $\lambda$, vertical axis is the norm $\norm{\dbar\mu-\frac{1}{4\pi\kk{\lambda}}\vct{t}(\lambda)e_{-\lambda}(z)\kk{\mu_0}}_{L^2}$. On the left we used $q_0^2$ and on the right $q_{35}^2$. In red using circles is $M=7$, in black using crosses $M=8$ and in blue using squares is $M=9$. Two smallest values for $\lambda$ were omitted, for $\lambda=1.01$ the magnitude of the error was between 3 and 13, for the second smallest it was between 0.003 and 0.02.}
\end{figure}

\subsection{Numerical investigation of exceptional points}
We use 250 discrete points of $\lambda$ and 701 discrete points of $\alpha$,
$$
\lambda = 1.01,..,4.5,\quad \alpha = -35,...,35.
$$
We use $M=8$ for the LS -solver. In figure \ref{q1-plane} we plot the radially symmetric and real-valued scattering transform $\vct{t}(\lambda)$ for the potential $q_\alpha^1 = \alpha \varphi$; the x-axis is the parameter $\alpha$ and the y-axis is the spectral parameter $\lambda$. In figure \ref{q2-plane} we have the same for  $q_\alpha^2 = q_\gamma+\alpha \varphi$. Black color represents very small negative values, and white very large positive values of $\vct{t}(\lambda)$. The lines where it abruptly changes between these colors are exceptional circles that move as the parameter $\alpha$ changes.

In figure \ref{scat_samples} we plot the profile of the scattering transform $\vct{t}(\lambda)$ as a function of $\lambda$, using the potential $q_\alpha^1 = \alpha \psi$, with the values $\alpha=-5,-15,-30$. The exceptional points can be seen as singularities in the profiles.

\begin{figure}
\begin{picture}(120,100)
\epsfxsize=12cm
\put(0,0){\epsffile{04-v3-q1_plane.eps}}
\end{picture}
\caption{Scattering transform for the potential $q_\alpha^1 = \alpha\varphi$; x-axis is $\alpha = -35..35$, y-axis is $\lambda = 1.01..4.5$\label{q1-plane}}
\end{figure}

\begin{figure}
\begin{picture}(120,100)
\epsfxsize=12cm
\put(0,0){\epsffile{04-v3-q2_plane.eps}}
\end{picture}
\caption{Scattering transform for the potential $q_\alpha^2 = q_\gamma+\alpha\varphi$; x-axis is $\alpha = -35..35$, y-axis is $\lambda = 1.01..4.5$\label{q2-plane}}
\end{figure}

\begin{figure}
\unitlength=1mm
\begin{picture}(120,130)
 \put(10,88){\includegraphics[width=6cm]{t_plane_p1.eps}}
 \put(73,90){\includegraphics[width=5cm]{t_vs_lambda_p1.eps}}
\put(0,110){$\alpha=-5$}
 \put(10,43){\includegraphics[width=6cm]{t_plane_p2.eps}}
 \put(73,45){\includegraphics[width=5cm]{t_vs_lambda_p2.eps}}
\put(0,65){$\alpha=-15$}
 \put(10,-2){\includegraphics[width=6cm]{t_plane_p3.eps}}
 \put(73,0){\includegraphics[width=5cm]{t_vs_lambda_p3.eps}}
\put(0,20){$\alpha=-30$}
\put(117,-2){$\lambda$}
\end{picture}
\caption{\label{scat_samples}On the right: the profile of $\vct{t}(\lambda)$ as a function of $\lambda$ using the potential $q_\alpha^1 = \alpha\varphi$ with three different values of $\alpha$. On the left: the plane $\vct{t}(\lambda)$ for all parameters $\alpha$ with an indication of the location of the profile on the right. }
\end{figure}

\section{Conclusions}
The numerical evidence show no exceptional points for small $\alpha$ nor for large $\lambda$, which was expected from the result of section \ref{absence}. Also according to our tests there are no exceptional points for positive $\alpha$. For negative $\alpha$, there are either one or two exceptional circles in the range of parameters investigated. The two types of potentials $q_\alpha^1$ and $q_\alpha^2$ have little difference in their exceptional points, mainly in the second exceptional circle forming at $\alpha=-20$ as $\alpha$ is decreased from zero. The numerical method for $g_\lambda$ is not accurate near $\abs{\lambda}=1$ which resulted in high errors in the test using the $\dbar$ -equation.

\section*{Acknowledgements}
JPT was supported in part by Suomen Kulttuurirahasto and European Research Council (ERC) Advanced Grant.


\bibliographystyle{plain}
\bibliography{Inverse_problems_references}


\end{document}



%
%\begin{thebibliography}{99}
%
%
%\bibitem{Faddeev}
%\newblock L.~D.~Faddeev,
%\newblock \emph{Increasing solutions of the Schrödinger equation},
%\newblock Sov. Phys. Dok \textbf{163} (1966), 1033--1035.
%
%\bibitem{Grinevich2000}
%\newblock P.~G.~Grinevich,
%\newblock \emph{Scattering transformation at fixed non-zero energy for the two-dimensional Schrödinger operator with potential decaying at infinity},
%\newblock Russian Math. Surveys, \textbf{55} (2000), 1015--1083.
%
%\bibitem{GM1986}
%\newblock P.~G.~Grinevich and S.~V.~Manakov,
%\newblock \emph{The inverse problem for the two dimensional Schrödinger operator: The $\dbar$-method and non-linear equations},
%\newblock Funktsional. Anal. i Prilozhen, \textbf{20} (1986), 14--24.
%
%\bibitem{Murata1986}
%\newblock M.~Murata,
%\newblock \emph{Structure of positive solutions to $(-\Delta+V)\xi = 0$},
%\newblock Rn . Duke Math. J., \textbf{53} (1986), 869--943.
%
%\bibitem{MPS2012}
%\newblock M.~Music, P.~A.~Perry and S.~Siltanen,
%\newblock \emph{Exceptional circles of radial potentials},
%\newblock manuscript.
%
%\bibitem{Nachman1996}
%\newblock A.~I.~Nachman,
%\newblock \emph{Global uniqueness for a two-dimensional inverse boundary value problem},
%\newblock Annals of Mathematics \textbf{143} (1996), 71--96.
%
%\bibitem{Novikov1992}
%\newblock R.~G.~Novikov,
%\newblock \emph{The Inverse Scatering Problem on a Fixed Energy Level for the Two-Dimensional Schrödinger Operator},
%\newblock Journal of functional analysis \textbf{103}, No. 2 (1992), 409--463.
%
%\bibitem{samu_vaikkari}
%\newblock S.~Siltanen,
%\newblock \emph{Electrical impedance tomography and Faddeev Green's functions},
%\newblock Ph.D thesis, Helsinki University of Technology (1999).
%
%\bibitem{LS_solver}
%\newblock G.~Vainikko,
%\newblock \emph{Fast Solvers of the Lippmann-Schwinger Equation},
%\newblock Research raport, Helsinki University of Technology (1997).
%
%
%\end{thebibliography}

PE_exceptional_points_v9.bbl
\begin{thebibliography}{10}

\bibitem{Ablowitz1974}
Mark~J. Ablowitz, David~J. Kaup, Alan~C. Newell, and Harvey Segur.
\newblock The inverse scattering transform-{F}ourier analysis for nonlinear
  problems.
\newblock {\em Studies in Appl. Math.}, 53(4):249--315, 1974.

\bibitem{Assenheimer2001}
Michel Assenheimer, Orah Laver-Moskovitz, Dov Malonek, David Manor, Udi
  Nahaliel, Ron Nitzan, and Abraham Saad.
\newblock The {T-SCAN} technology: electrical impedance as a diagnostic tool
  for breast cancer detection.
\newblock {\em Physiological measurement}, 22:1--8, 2001.

\bibitem{Astala2011}
K.~Astala, J.L. Mueller, L.~P{\"a}iv{\"a}rinta, A.~Per{\"a}m{\"a}ki, and
  S.~Siltanen.
\newblock Direct electrical impedance tomography for nonsmooth conductivities.
\newblock {\em Inverse Problems and Imaging}, 5(3):531--549, 2011.

\bibitem{Astala2010}
K.~Astala, J.L. Mueller, L.~P{\"a}iv{\"a}rinta, and S.~Siltanen.
\newblock Numerical computation of complex geometrical optics solutions to the
  conductivity equation.
\newblock {\em Applied and Computational Harmonic Analysis}, 29(1):391--403,
  2010.

\bibitem{Astala2006}
K.~Astala and L.~P{\"a}iv{\"a}rinta.
\newblock A boundary integral equation for {C}alder\'on's inverse conductivity
  problem.
\newblock In {\em Proc. 7th Internat. Conference on Harmonic Analysis,
  Collectanea Mathematica}, 2006.

\bibitem{Astala2006a}
K.~Astala and L.~P{\"a}iv{\"a}rinta.
\newblock {C}alder\'on's inverse conductivity problem in the plane.
\newblock {\em Annals of Mathematics}, 163(1):265--299, 2006.

\bibitem{Beals1981}
Richard Beals and Ronald~R. Coifman.
\newblock Scattering, transformations spectrales et \'equations d'\'evolution
  non lin\'eaires.
\newblock In {\em Goulaouic-Meyer-Schwartz Seminar, 1980--1981}, pages Exp. No.
  XXII,10. \'Ecole Polytech., Palaiseau, 1981.

\bibitem{Boiti1986}
M.~Boiti, J.~Jp. Leon, M.~Manna, and F.~Pempinelli.
\newblock On the spectral transform of a {K}orteweg-de {V}ries equation in two
  spatial dimensions.
\newblock {\em Inverse Problems}, 2(3):271--279, 1986.

\bibitem{Boiti1987}
M.~Boiti, J.~P. Leon, M.~Manna, and F.~Pempinelli.
\newblock On a spectral transform of a {KdV}-like equation related to the
  {S}chr{\"o}dinger operator in the plane.
\newblock {\em Inverse Problems}, 3:25--36, 1987.

\bibitem{Brown1997}
R.~M. Brown and G.~Uhlmann.
\newblock Uniqueness in the inverse conductivity problem for nonsmooth
  conductivities in two dimensions.
\newblock {\em Communications in Partial Differential Equations},
  22(5):1009--1027, 1997.

\bibitem{Calder'on1980}
A.-P. {C}alder{\'o}n.
\newblock On an inverse boundary value problem.
\newblock In {\em Seminar on {N}umerical {A}nalysis and its {A}pplications to
  {C}ontinuum {P}hysics ({R}io de {J}aneiro, 1980)}, pages 65--73. Soc. Brasil.
  Mat., Rio de Janeiro, 1980.

\bibitem{Cheney1999}
M.~Cheney, D.~Isaacson, and J.~C. Newell.
\newblock Electrical impedance tomography.
\newblock {\em SIAM Review}, 41(1):85--101, 1999.

\bibitem{Faddeev1966}
L.~D. Faddeev.
\newblock Increasing solutions of the {S}chr{\"o}dinger equation.
\newblock {\em Soviet Physics Doklady}, 10:1033--1035, 1966.

\bibitem{Gardner1967}
Clifford~S. Gardner, John~M. Greene, Martin~D. Kruskal, and Robert~M. Miura.
\newblock Method for solving the korteweg-devries equation.
\newblock {\em Phys. Rev. Lett.}, 19:1095--1097, Nov 1967.

\bibitem{Grinevich2000}
P.~G. Grinevich.
\newblock Scattering transformation at fixed non-zero energy for the
  two-dimensional {S}chr{\"o}dinger operator with potential decaying at
  infinity.
\newblock {\em Russian Mathematical Surveys}, 55(2):1015--1083, 2000.

\bibitem{Grinevich1986}
P.~G. Grinevich and S.~V. Manakov.
\newblock Inverse scattering problem for the two-dimensional {S}chr{\"o}dinger
  operator, the $\dbar$-method and nonlinear equations.
\newblock {\em Functional Analysis and its Applications}, 20:94--103, 1986.

\bibitem{Grinevich2012}
P.G. Grinevich and R.G. Novikov.
\newblock Faddeev eigenfunctions for point potentials in two dimensions.
\newblock {\em Physics Letters A}, 376(12–13):1102 -- 1106, 2012.

\bibitem{Hamilton2012}
S.J. Hamilton, C.N.L. Herrera, J.~L. Mueller, and A.~Von~Herrmann.
\newblock A direct {D}-bar reconstruction algorithm for recovering a complex
  conductivity in 2-{D}.
\newblock {\em Arxiv preprint}, 2012.

\bibitem{Isaacson2004}
D.~Isaacson, J.~L. Mueller, J.~C. Newell, and S.~Siltanen.
\newblock Reconstructions of chest phantoms by the {D}-bar method for
  electrical impedance tomography.
\newblock {\em IEEE Transactions on Medical Imaging}, 23:821--828, 2004.

\bibitem{Isaacson2006}
D.~Isaacson, J.L. Mueller, J.C. Newell, and S.~Siltanen.
\newblock Imaging cardiac activity by the {D}-bar method for electrical
  impedance tomography.
\newblock {\em Physiological Measurement}, 27:S43--S50, 2006.

\bibitem{Knudsen2003}
K.~Knudsen.
\newblock A new direct method for reconstructing isotropic conductivities in
  the plane.
\newblock {\em Physiological Measurement}, 24(2):391--403, 2003.

\bibitem{Knudsen2009}
K.~Knudsen, M.~Lassas, J.L. Mueller, and S.~Siltanen.
\newblock Regularized {D}-bar method for the inverse conductivity problem.
\newblock {\em Inverse Problems and Imaging}, 3(4):599--624, 2009.

\bibitem{Knudsen2004a}
K.~Knudsen and A.~Tamasan.
\newblock Reconstruction of less regular conductivities in the plane.
\newblock {\em Communications in Partial Differential Equations}, 29:361--381,
  2004.

\bibitem{Korteweg1895}
D.~J. Korteweg and G.~de~Vries.
\newblock Xli. on the change of form of long waves advancing in a rectangular
  canal, and on a new type of long stationary waves.
\newblock {\em Philosophical Magazine Series 5}, 39(240):422--443, 1895.

\bibitem{Lassas2007}
M.~Lassas, J.~L. Mueller, and S.~Siltanen.
\newblock Mapping properties of the nonlinear {F}ourier transform in dimension
  two.
\newblock {\em Communications in Partial Differential Equations},
  32(4):591--610, 2007.

\bibitem{Lassas2012b}
M.~Lassas, J.~L. Mueller, S.~Siltanen, and A.~Stahel.
\newblock The {N}ovikov-{V}eselov equation and the inverse scattering method:
  {II}. {C}omputation.
\newblock {\em Nonlinearity}, 25(6):1799--1818, 2012.

\bibitem{Lassas2012a}
M.~Lassas, J.~L. Mueller, S.~Siltanen, and A.~Stahel.
\newblock The {N}ovikov-{V}eselov equation and the inverse scattering method,
  {P}art {I}: {A}nalysis.
\newblock {\em Phys. D}, 241(16):1322--1335, 2012.

\bibitem{Mueller2003}
J.L. Mueller and S.~Siltanen.
\newblock Direct reconstructions of conductivities from boundary measurements.
\newblock {\em SIAM Journal on Scientific Computing}, 24(4):1232--1266, 2003.

\bibitem{Mueller2012}
J.L. Mueller and S.~Siltanen.
\newblock {\em Linear and Nonlinear Inverse Problems with Practical
  Applications}.
\newblock SIAM, 2012.

\bibitem{Murata1986}
Minoru Murata.
\newblock Structure of positive solutions to {$(-\Delta+V)u=0$} in {${\bf
  R}^n$}.
\newblock {\em Duke Math. J.}, 53(4):869--943, 1986.

\bibitem{Music2013}
M~Music, P~Perry, and S~Siltanen.
\newblock Exceptional circles of radial potentials.
\newblock {\em Inverse Problems}, 29(4):045004, 2013.

\bibitem{Nachman1988}
A.~Nachman, J.~Sylvester, and G.~Uhlmann.
\newblock An $n$-dimensional {B}org--{L}evinson theorem.
\newblock {\em Communications in Mathematical Physics}, 115:595--605, 1988.

\bibitem{Nachman1996}
A.~I. Nachman.
\newblock Global uniqueness for a two-dimensional inverse boundary value
  problem.
\newblock {\em Annals of Mathematics}, 143:71--96, 1996.

\bibitem{Novikov1992}
R.~G. Novikov.
\newblock The inverse scattering problem on a fixed energy level for the
  two-dimensional {S}chr{\"o}dinger operator.
\newblock {\em Journal of Functional Analysis}, 103(2):409--463, 1992.

\bibitem{Novikov1986}
S.~P. Novikov and A.~P. Veselov.
\newblock Two-dimensional {S}chr{\"o}dinger operator: inverse scattering
  transform and evolutional equations.
\newblock {\em Physica D}, 18:267--273, 1986.

\bibitem{Perry2012}
P.~A. {Perry}.
\newblock {Miura Maps and Inverse Scattering for the Novikov-Veselov Equation}.
\newblock {\em ArXiv e-prints}, January 2012.

\bibitem{Siltanen2000}
S.~Siltanen, J.~Mueller, and D.~Isaacson.
\newblock An implementation of the reconstruction algorithm of {A}. {N}achman
  for the 2-{D} inverse conductivity problem.
\newblock {\em Inverse Problems}, 16:681--699, 2000.

\bibitem{Siltanen1999}
Samuli Siltanen.
\newblock Electrical impedance tomography and {F}addeev {G}reen's functions.
\newblock {\em Annales Academiae Scientiarum Fennicae Mathematica
  Dissertationes}, 121(121):56, 1999.
\newblock Dissertation, Helsinki University of Technology, Espoo, 1999.

\bibitem{Sylvester1987}
J.~Sylvester and G.~Uhlmann.
\newblock A global uniqueness theorem for an inverse boundary value problem.
\newblock {\em Annals of Mathematics}, 125:153--169, 1987.

\bibitem{Tsai1989}
T.~Tsai.
\newblock {\em Nonlinear Evolutions of the {S}chr{\"o}dinger Operator in
  $\mathbb{R}^2$}.
\newblock PhD thesis, Yale University, 1989.

\bibitem{Tsai1993}
T.~Tsai.
\newblock The {S}chr{\"o}dinger operator in the plane.
\newblock {\em Inverse Problems}, 9:763--787, 1993.

\bibitem{Tsai1994}
T.~Tsai.
\newblock The associated evolution equations of the {S}chr{\"o}dinger operator
  in the plane.
\newblock {\em Inverse Problems}, 10:1419, 1994.

\bibitem{Vainikko2000}
Gennadi Vainikko.
\newblock Fast solvers of the {L}ippmann-{S}chwinger equation.
\newblock In {\em Direct and inverse problems of mathematical physics
  ({N}ewark, {DE}, 1997)}, volume~5 of {\em Int. Soc. Anal. Appl. Comput.},
  pages 423--440. Kluwer Acad. Publ., Dordrecht, 2000.

\bibitem{Veselov1984}
A.~P. Veselov and S.~P. Novikov.
\newblock Finite-zone, two-dimensional, potential {S}chr{\"o}dinger operators,
  explicit formulas and evolution equations.
\newblock {\em Sov. Math. Dokl}, 30:558--591, 1984.

\end{thebibliography}
