Transfer Matrix Formulation of Scattering Theory in Two and Three Dimensions

Scattering of waves has been one of the most important subjects of research in physics for over a century. Most students of physics begin their study of this phenomenon in undergraduate courses of quantum mechanics where they learn time-independent scattering theory and its application in solving simple one-dimensional problems. In more advanced courses they are exposed to the $\mathbf{S}$-matrix formulation of scattering theory which is of basic importance in quantum field theories and particle physics weinberg .

In one dimension, there is an alternative to the $\mathbf{S}$-matrix, called the transfer matrix, that similarly to the $\mathbf{S}$-matrix stores all the scattering information about the scattering potential. To give the definition of the transfer matrix, we recall that for a scattering potential $v$, every solution of the time-independent Schrödinger equation,

$$-\psi^{\prime\prime}(x)+v(x)\psi(x)=k^{2}\psi(x),$$

(1)

has the asymptotic form:

$$\psi(x)=A_{\pm}e^{ikx}+B_{\pm}e^{-ikx}~~{\rm for}~~x\to\pm\infty.$$

The transfer matrix of $v$ is the $2\times 2$ matrix $\mathbf{M}$ satisfying

$$\left[\begin{array}[]{c}A_{+}\\ B_{+}\end{array}\right]=\mathbf{M}\left[\begin{array}[]{c}A_{-}\\ B_{-}\end{array}\right].$$

(2)

It is easy to show that the reflection and transmission coefficients of $v$ are given by the entries of $\mathbf{M}$, razavy ; jpa-2009a ; sanchez .

The main advantage of the transfer matrix over the $\mathbf{S}$-matrix is that it satisfies an extremely useful composition rule which allows for the reduction of the scattering problem for $v$ to that of any set of constituent potentials $v_{1},v_{2},\cdots,v_{n}$ provided that $v=\sum_{j}v_{j}$ and the support $I_{j}$ of $v_{j}$ (i.e., the smallest interval outside which $v_{j}$ vanishes) do not intersect except possibly at a single point. If $\mathbf{M}_{j}$ is the transfer matrix of $v_{j}$ and $I_{j}$ lies to the left of $I_{j+1}$, the transfer matrix of $v$ is given by

$$\mathbf{M}=\mathbf{M}_{n}\mathbf{M}_{n-1}\cdots\mathbf{M}_{1}.$$

(3)

This relation makes the transfer matrix into an indispensable tool for dealing with a variety of scattering problems sanchez . Typical examples are the study of multilayer optical systems yeh-LPP-1 ; yeh-LPP-2 ; yeh-LPP-3 ; yeh-LPP-4 , acoustic waves accoustics-1 ; accoustics-2 ; accoustics-3 , transport phenomena in mesoscopic systems MK , photonic band structure calculations pendry-92 ; pendry-94 , optical potential engineering pra-2014bc-2015b-2 ; pra-2014bc-2015b-3 ; pra-2014bc-2015b-4 , etc. These applications are however limited to (effectively) one-dimensional physical systems. This is because the very definition of a multidimensional transfer matrix, that contains all the information about the scattering features of the potential and fulfils the composition property (3), has been lacking. The purpose of the present article is to introduce such a notion of transfer matrix and use it to develop an alternative to standard approach to scattering theory sakurai ; weinberg-qm that is capable of dealing with a variety of intricate scattering problems.

The basic conceptual framework for the present work is provided by the recent observation that the transfer matrix in one dimension can be expressed as

$$\mathbf{M}=\mathbf{U}(\infty,-\infty),$$

(4)

where $\mathbf{U}(x,x_{0})$ is the evolution operator for a fictitious non-unitary two-level quantum system with $x$ playing the role of the evolution parameter ap-2014 ; pra-2014a . This was originally motivated by the fact that the evolution operator of every quantum system has the same composition (group) property as the transfer matrix, i.e., (3). It turns that one can define a two-component state vector $\Psi(x)$ and an $x$-dependent effective matrix Hamiltonian ${\mathbf{H}}(x)$ of the form

	$\displaystyle\Psi(x)$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\left[\begin{array}[]{c}e^{-ikx}[\psi(x)-ik^{-1}\psi^{\prime}(x)]\\[3.0pt] e^{ikx}[\psi(x)+ik^{-1}\psi^{\prime}(x)]\end{array}\right],$		(7)
	$\displaystyle{\mathbf{H}}(\tau)$	$\displaystyle:=$	$\displaystyle\frac{v(x)}{2k}\left[\begin{array}[]{cc}1&e^{-2ikx}\\ -e^{2ikx}&-1\end{array}\right],$		(10)

such that the time-independent Schrödinger equation (1) becomes equivalent to the time-dependent Schrödinger equation

$$i\partial_{x}\Psi(x)=\mathbf{H}(x)\Psi(x),$$

(11)

with $x$ playing the role of “time,” ap-2014 .

Clearly, there is an infinity of choices for $\Psi(x)$ (and the corresponding $\mathbf{H}(x)$) that allow for expressing (1) in the form (11), FV-1958 ; jmp-1998 . Note however that among these, (7) is the only choice that fulfils the asymptotic boundary conditions:

$$\Psi(x)\to\left[\begin{array}[]{c}A_{\pm}\\ B_{\pm}\end{array}\right]~~{\rm for}~~x\to\pm\infty.$$

(12)

This relation is the key to the identification of the transfer matrix $\mathbf{M}$ with $\mathbf{U}(-\infty,\infty)$, where $\mathbf{U}(x,x_{0})$ is the evolution operator for the Hamiltonian (10), i.e., the operator satisfying

$$\Psi(x)=\mathbf{U}(x,x_{0})\Psi(x_{0}).$$

(13)

Equation (4) is a direct consequence of (2), (12), and (13) for $x_{0}\to-\infty$ and $x\to\infty$.

Another important consequence of (7) is that the corresponding effective Hamiltonian, namely (10), vanishes for the values of $x$ where $v(x)=0$. We can use this observation to give a derivation of the composition rule for transfer matrices. Let us denote by $\mathbf{H}_{j}(x)$ and $\mathbf{U}_{j}(x,x_{0})$ the effective Hamiltonian and the evolution operator for the constituent potentials $v_{j}$, respectively. Suppose that the support of $v_{j}$ is given by $I_{j}=[a_{j},a_{j+1}]$ for some real numbers $a_{j}$ and $a_{j+1}$. Then for all $x\notin I_{j}$, $\mathbf{H}_{j}(x)=0$. This in turn implies that whenever the interval $(x_{0},x)$ does not intersect $I_{j}$, we have $\mathbf{U}_{j}(x,x_{0})=\mathbf{I}$, where $\mathbf{I}$ is the $2\times 2$ identity matrix. In particular,

$$M_{j}=\mathbf{U}_{j}(\infty,-\infty)=\mathbf{U}_{j}(a_{j+1},a_{j})=\mathbf{U}(a_{j+1},a_{j}),$$

(14)

where the last equality follows from the fact that $v_{j}(x)$ is the restriction of $v(x)$ to $I_{j}$. Next, we recall that because $I_{j}$ lies to the left of $I_{j+1}$, $a_{1}\leq a_{2}\leq a_{3}\leq\cdots\leq a_{n+1}$. This together with (4), (13), and (14) imply

$$\mathbf{U}(\infty,\infty)=\mathbf{U}(a_{n+1},a_{n})\mathbf{U}(a_{n},a_{n-1})\cdots\mathbf{U}(a_{2},a_{1}).$$

In light of (4) and (14), this equation is equivalent to (3).

The idea of expressing the transfer matrix of one-dimensional scattering theory in terms of the evolution operator for a two-level quantum system has revealed an intriguing manifestation of geometric phases in semiclassical scattering jpa-2014a . It has also been used to devise a method of computing higher order corrections to semiclassical scattering calculations jpa-2014b and a generalization of the composition rule (3) to situations where the constituting potentials $v_{j}$ have overlapping support ap-2015 . In this article, we use this idea as a guiding principle for developing a transfer matrix formulation of the scattering theory in two and three dimensions. Since the basic difficulty that has hindered the progress towards solving this problem is already present in two dimensions, we give a more thorough treatment of scattering theory in two-dimensions. In order to focus on the conceptual developments and their concrete applications, we give the technical details of most of the calculations in the appendices.

Consider the Schrödinger equation,

$$\left[-\partial_{x}^{2}-\partial_{y}^{2}+v(x,y)\right]\psi(x,y)=k^{2}\psi(x,y),$$

(15)

where $v$ is a real- or complex-valued scattering potential and $(x,y)$ are cartesian coordinates in the plane $\mathbb{R}^{2}$. Let us identify the $x$-axis with the scattering axis, use $\mathbf{r}$ to denote the position vector associated with $(x,y)$, and take the Fourier transform of (15) with respect to $y$. This gives

$$\left[-\partial_{x}^{2}+v(x,i\partial_{p})\right]\tilde{\psi}(x,p)=\omega(p)^{2}\tilde{\psi}(x,p),$$

(16)

where $\tilde{\psi}(x,p):=\int_{-\infty}^{\infty}dye^{-ipy}\psi(x,y)$ is the Fourier transform of $\psi(x,y)$ over $y$ and $\omega(p):=\sqrt{k^{2}-p^{2}}$. We identify $v(x,i\partial_{p})$ with the integral operator defined by $v(x,i\partial_{p})\phi(p):=\frac{1}{2\pi}\int_{-\infty}^{\infty}dq\,\tilde{v}(x,p-q)\phi(q)$, where $\phi$ is a test function.

For the scattering problem, we are interested in the solutions of (16) that tend to a plane wave for $x\to\pm\infty$. These exist provided that $|p|\leq k$ and $|v(x,y)|$ has a sufficiently rapid asymptotic decay rate. We therefore demand that as $r:=\sqrt{x^{2}+y^{2}}\to\infty$, $|v(x,y)|$ tends to zero at such a rate that every solution of (15) has the following asymptotic form for $x\to\pm\infty$:

$$\frac{1}{2\pi}\int_{-k}^{k}dp\,e^{ipy}\left[A_{\pm}(p)e^{i\omega(p)x}+B_{\pm}(p)e^{-i\omega(p)x}\right].$$

Here $A_{\pm}$ and $B_{\pm}$ are coefficient functions vanishing outside the interval $[-k,k]$.

We identify the transfer matrix of the potential $v$ with the $2\times 2$ matrix operator $\mathbf{M}$ satisfying (2). As we show in the sequel, this notion of transfer matrix not only fulfills the composition rule (3), but also possesses complete information about the scattering properties of $v$.

We begin our analysis of the properties of $\mathbf{M}$ by introducing the following two-dimensional generalization of the two-component state vector (7) and effective Hamiltonian operator (10).

	$\displaystyle\Psi(x,p)$	$\displaystyle:=\frac{1}{2}\left[\begin{array}[]{c}e^{-i\omega(p)x}[\tilde{\psi}(x,p)-i\omega(p)^{-1}\tilde{\psi}^{\prime}(x,p)]\\[3.0pt] e^{i\omega(p)x}[\tilde{\psi}(x,p)+i\omega(p)^{-1}\tilde{\psi}^{\prime}(x,p)]\end{array}\right],$		(19)
	$\displaystyle\boldsymbol{\mathbf{H}}(x,p)$	$\displaystyle:=\frac{1}{2\omega(p)}e^{-i\omega(p)x\boldsymbol{\sigma}_{3}}v(x,i\partial_{p})\,\boldsymbol{\mathcal{K}}\,e^{i\omega(p)x\boldsymbol{\sigma}_{3}},$		(20)

where $\tilde{\psi}^{\prime}(x,p):=\partial_{x}\tilde{\psi}(x,p)$, $\boldsymbol{\sigma}_{i}$ are the Pauli matrices, and $\boldsymbol{\mathcal{K}}:=\boldsymbol{\sigma}_{3}+i\boldsymbol{\sigma}_{2}$. $\Psi(x,p)$ has the following two important properties:

1. 1.

   It satisfies $i\partial_{x}\Psi(x,p)=\boldsymbol{\mathbf{H}}(x,p)\Psi(x,p)$, which is equivalent to (16). In other words, $\Psi(x,p)=\boldsymbol{\mathbf{U}}(x,x_{0})\Psi(x_{0},p)$, where $\boldsymbol{\mathbf{U}}(x,x_{0})$ is the evolution operator for $\boldsymbol{\mathbf{H}}(x,p)$.

2. 2.

   It tends to $\left[\begin{array}[]{c}A_{\pm}(p)\\ B_{\pm}(p)\end{array}\right]$ for $x\to\pm\infty$.

By virtue of these properties and Eq. (2), we can express $\mathbf{M}$ in the form

	$\displaystyle\mathbf{M}$	$\displaystyle=$	$\displaystyle\boldsymbol{\mathbf{U}}(\infty,-\infty)=\mathscr{T}\exp\int_{-\infty}^{\infty}\!\!-i\boldsymbol{\mathbf{H}}(x,p)dx$
		$\displaystyle=$	$\displaystyle\mathbf{I}-i\int_{-\infty}^{\infty}\!\!\!dx_{1}\boldsymbol{\mathbf{H}}(x_{1},p)$
			$\displaystyle+(-i)^{2}\int_{-\infty}^{\infty}\!\!\!dx_{2}\int_{-\infty}^{x_{2}}\!\!\!dx_{1}\boldsymbol{\mathbf{H}}(x_{2},p)\boldsymbol{\mathbf{H}}(x_{1},p)+\cdots,$

where $\mathscr{T}$ is the “time-ordering” operation weinberg-qm with $x$ playing the role of “time.”

According to (20), for each $x$ satisfying $v(x,y)=0$, we have $\boldsymbol{\mathbf{H}}(x,p)=0$. Similarly to the case of one dimension, we can employ this property of $\boldsymbol{\mathbf{H}}(x,p)$ together with Eq. (II) to establish the composition rule (3). This rule holds true if $I_{j}$ do not have overlapping projections onto the scattering ($x$-) axis except possibly in a discrete set of points. For example, we can slice the support of $v$ along the $x$-axis, label the slices by $I_{j}$ (with $I_{j}$ to the left of $I_{j+1}$ as depicted in Fig. 1), and take

$$v_{j}(\mathbf{r})=\left\{\begin{array}[]{cc}v(\mathbf{r})&{\rm for}~\mathbf{r}\in I_{j},\\ 0&{\rm for}~\mathbf{r}\notin I_{j},\end{array}\right.$$

where $\mathbf{r}:=(x,y)$.

In order to see how $\mathbf{M}$ relates to the scattering amplitude of the potential $v$, we consider a left-incident wave and write the corresponding scattered wave in the form

$$\psi_{\rm scat}(\mathbf{r})=\Theta(-x)\psi_{-}(\mathbf{r})+\Theta(x)\psi_{+}(\mathbf{r})~~~~{\rm as}~r\to\infty,$$

(22)

where $\Theta$ is the Heaviside step function, i.e.,

$$\Theta(x):=\left\{\begin{array}[]{ccc}1&{\rm for}&x\geq 0,\\ 0&{\rm for}&x<0,\end{array}\right.$$

and $\psi_{-}(\mathbf{r})$ and $\psi_{+}(\mathbf{r})$ are respectively the reflected and transmitted waves given by

	$\displaystyle\psi_{\pm}(\mathbf{r}):=\frac{1}{2\pi}\int_{-k}^{k}dp\,T_{\pm}(p)\,e^{ipy}e^{\pm i\omega(p)x}~~{\rm for}~r\to\infty,$		(23)
	$\displaystyle T_{+}(p):=A_{+}(p)-A_{-}(p),~~~~~~~~~T_{-}(p):=B_{-}(p).$		(24)

The functions $T_{\pm}(p)$ yield the momentum representation of the reflected and transmitted waves $\psi_{\pm}(\mathbf{r})$.

For a left-incident plane wave with a definite momentum $k$ along the positive $x$-axis, we have

$\displaystyle A_{-}(p)=2\pi\delta(p),$

$\displaystyle B_{+}(p)=0,$

(25)

where $\delta$ stands for the Dirac delta-function. Substituting (25) in (2) and using (24), we can express $T_{\pm}(p)$ in terms of the entries $M_{ij}$ of $\mathbf{M}$ according to

		$\displaystyle T_{+}(p)=2\pi(M_{11}-1-M_{12}M_{22}^{-1}M_{21})\delta(p),$		(26)
		$\displaystyle T_{-}(p)=-2\pi M_{22}^{-1}M_{21}\delta(p).$

Next, we write the scattered wave in its standard form,

$$\psi_{{\rm scat}}(\mathbf{r})=\sqrt{\frac{i}{kr}}\,e^{ikr}f(\theta)~~~~{\rm as}~r\to\infty,$$

(27)

where $(r,\theta)$ are the polar coordinates of $\mathbf{r}$, and $f(\theta)$ is the scattering amplitude of the potential adhikari-86 . A careful analysis of the right-hand side of (23) gives a closed form asymptotic expression for $\psi_{\pm}(\mathbf{r})$ (See Appendix A.) Using this in (22) and comparing the result with (27), we obtain

$$f(\theta)=-\frac{i\omega(p)T_{\pm}(p)}{\sqrt{2\pi}}=-\frac{ik|\cos\theta|}{\sqrt{2\pi}}\,T_{\pm}(k\sin\theta),$$

(28)

where $\pm:={\rm sgn}(\cos\theta)$.

Equations (26) and (28) show that $\mathbf{M}$ includes all the information about the scattering properties of the potential $v$. In other words, the solution of the scattering problem can be achieved by evaluating $\mathbf{M}$ using (II) and computing $T_{\pm}(p)$ and $f(\theta)$ as given by (26) and (28).

Because the effective Hamiltonian $\mathbf{H}(x,p)$ scales with the potential $v$, Eq. (II) provides a convenient method of perturbative solution of the scattering problem (similarly to its one-dimensional analog pra-2014a ). We have checked that if we perform the first order perturbation theory using our approach, we obtain the known expression for the scattering amplitude in the first Born approximation Lapidus-1982 . Here we do not wish to discuss the utility of our method in perturbation theory. Instead, we demonstrate its application in providing an exact (non-perturbative) solution of the old and difficult problem of delta-function potential in two dimensions. The perturbation series solution of the Lippmann-Schwinger equation sakurai for the delta-function potential in two dimensions is plagued with the emergence of infinities. The same problem arises in the study of the corresponding spectral problem and has led to the development of intricate renormalization schemes that amount to the introduction of a length scale for the problem thorn ; exner ; henderson-1997 ; henderson-1998 ; cabo-1998 ; phillips ; rajeev-99 ; frederico ; Camblong-2001a ; camblong-2002 . The problem of divergences also arises in the standard treatment of the delta-function potential in 3D. As we show in the sequel, our approach avoids this problem in both two and three dimensions.

Consider

$$v(x,y)=\mathfrak{z}\,\delta(x)\delta(y),$$

(29)

where $\mathfrak{z}$ is a real or complex coupling constant. Then $\mathbf{H}=\mathfrak{z}\,\omega(p)^{-1}\delta(x)\delta(i\partial_{p})\boldsymbol{\mathcal{K}}/2$ and (II) gives

$\displaystyle\mathbf{M}=\mathbf{I}-\frac{i\mathfrak{z}}{2\omega(p)}\,\delta(i\partial_{p})\boldsymbol{\mathcal{K}}.$

(30)

As we show in Appendix B, the calculation of $T_{\pm}(p)$ is more conveniently carried out by substituting (30) in (2) and making use of (24) and (25). This gives $T_{\pm}(p)=-2i\mathfrak{z}/(4+i\mathfrak{z})\omega(p)$ which together with (28) yield the following exact expression for the scattering amplitude of (29).

$\displaystyle f(\theta)=-\sqrt{\frac{2}{\pi}}\frac{\mathfrak{z}}{4+i\mathfrak{z}}.$

(31)

The first Born approximation corresponds to the linear approximation of this equation, i.e., $f(\theta)\approx-\mathfrak{z}/2\sqrt{2\pi}$, which agrees with the results of Ref. Lapidus-1982 ; Camblong-2001a .

According to (31), $f(\theta)$ diverges for $\mathfrak{z}=4i$. This marks a zero-width resonance that is associated with what mathematicians call a spectral singularity prl-2009 . In optical applications pra-2011a , $v$ is related to the complex relative permittivity $\varepsilon$ of the medium according to $v=k^{2}(1-\varepsilon)$. In particular, (29) corresponds to a thin wire of nonmagnetic material that is aligned along the $z$-axis and described by $\varepsilon=1+i\zeta\delta(x)\delta(y)$ for some real parameter $\zeta$. Denoting the unit vector along the positive $z$-axis by $\mathbf{e}_{z}$, we can easily check that

$$\mathbf{E}(x,y,z,t)=e^{-ickt}\psi(x,y)\mathbf{e}_{z}$$

(32)

solves Maxwell’s equations for this system provided that $\psi$ satisfies (15) for the delta-function potential (29) with $\mathfrak{z}=-i\zeta k^{2}$. In terms of $\zeta$, the condition for the emergence of a spectral singularity takes the form $k^{2}=-4/\zeta$. This shows that the wire begins lasing pra-2011a for $k=2/\sqrt{-\zeta}$ whenever $\zeta<0$ (i.e., it is made of gain material), and functions as a coherent perfect absorber CPA ; CPA-Science ; longhi-1 ; longhi-3 for $k=2/\sqrt{\zeta}$ whenever $\zeta>0$ (it is made of lossy material).

As a simple application of the composition property of our transfer matrix, we address the scattering problem for electromagnetic waves interacting with an infinite planar slab of optically active material with a delta-function line defect on one of its faces as depicted in Fig. 2.

Identifying the location of the slab and the line defect respectively with $0\leq x\leq L$ and $x=y=0$ ($z$-axis) and considering the linearly polarized waves of the form (32), we can reduce the problem to the two-dimensional scattering problem for the potential $v_{\delta}+\tilde{v}$, where $v_{\delta}$ is the delta-function potential (31) which represents the defect, and $\tilde{v}$ is a barrier potential of the form

$$\tilde{v}(x,y):=\left\{\begin{array}[]{cc}\tilde{\mathfrak{z}}&{\rm for}~0\leq x\leq L,\\ 0&{\rm otherwise},\end{array}\right.$$

(33)

which represents the slab pra-2011a . Here $\tilde{\mathfrak{z}}:=k^{2}(1-\varepsilon)$, $\varepsilon$ is the complex relative permittivity of the slab, and $L$ is its thickness.

Because, along the $x$-axis, the support of $v_{\delta}$ is to the left of that of $\tilde{v}$, the transfer matrix $\mathbf{M}$ of $v_{\delta}+\tilde{v}$ is given by $\mathbf{M}=\tilde{\mathbf{M}}\mathbf{M}_{\delta}$, where $\tilde{\mathbf{M}}$ and $\mathbf{M}_{\delta}$ are respectively the transfer matrices of $\tilde{v}$ and $v_{\delta}$.

We have already computed $\mathbf{M}_{\delta}$. It has the form (30). We give the details of the calculation of $\tilde{\mathbf{M}}$ in Appendix C. Its entries have the following expression.

		$\displaystyle\tilde{M}_{11}(\omega)=\tilde{M}_{22}(-\omega)=[\cos({\mathfrak{n}}L\omega)+i{\mathfrak{n}}_{+}\sin({\mathfrak{n}}L\omega)]e^{-i\omega L},$		(34)
		$\displaystyle\tilde{M}_{12}(\omega)=\tilde{M}_{21}(-\omega)=i{\mathfrak{n}}_{-}\sin({\mathfrak{n}}L\omega)e^{-i\omega L},$

where ${\mathfrak{n}}:=\sqrt{1-\tilde{\mathfrak{z}}/\omega^{2}}=\sqrt{(\varepsilon-1)(k/\omega)^{2}+1}$ and ${\mathfrak{n}}_{\pm}:=({\mathfrak{n}}\pm{\mathfrak{n}}^{-1})/2$. For $p=0$, $\omega=k$, and ${\mathfrak{n}}$ coincides with the complex refractive index of the slab, namely $\sqrt{\varepsilon}$. In this case, (34) gives the known result obtained in one dimension pra-2011a ; pra-2015c .

Having determined $\mathbf{M}_{\delta}$ and $\tilde{\mathbf{M}}$, we can use the expression $\mathbf{M}=\tilde{\mathbf{M}}\mathbf{M}_{\delta}$ together with (2), (24) and (25) to compute $T_{\pm}$. As we show in Appendix D, this gives

	$\displaystyle\omega\,T_{-}$	$\displaystyle=2\pi k[\mathcal{X}(k)-1]\delta(p)-\frac{i\mathfrak{z}\,\mathcal{X}(k)\mathcal{X}(\omega)}{\mathcal{Y}(k)},$		(35)
	$\displaystyle\omega\,T_{+}$	$\displaystyle=2\pi k\left[\frac{1}{\tilde{M}_{22}(k)}-1\right]\delta(p)-\frac{i\mathfrak{z}\,\mathcal{X}(k)}{\mathcal{Y}(k)\tilde{M}_{22}(\omega)},$

where

	$\displaystyle\mathcal{X}(\omega)$	$\displaystyle:=1-\frac{\tilde{M}_{21}(\omega)}{\tilde{M}_{22}(\omega)}=\frac{2\left(e^{-2i{\mathfrak{n}}L\omega}+\frac{{\mathfrak{n}}-1}{{\mathfrak{n}}+1}\right)}{({\mathfrak{n}}+1)\mathcal{Z}(\omega)},$
	$\displaystyle\mathcal{Y}(k)$	$\displaystyle:=2+\frac{i\mathfrak{z}}{\pi}\int_{0}^{k}\frac{\mathcal{X}(\omega)d\omega}{\sqrt{k^{2}-\omega^{2}}},$
	$\displaystyle\mathcal{Z}(\omega)$	$\displaystyle:=e^{-2i{\mathfrak{n}}L\omega}-\left(\frac{{\mathfrak{n}}-1}{{\mathfrak{n}}+1}\right)^{2}.$

Substituting (35) in (28), we find the scattering amplitude for the system.

In the absence of the line defect, $\mathfrak{z}=0$, and the condition for the emergence of a spectral singularity takes the form $\tilde{M}_{22}(k)=0$ which is equivalent to $\mathcal{Z}(k)=0$. This equation turns out to be equivalent to the laser threshold and phase conditions for the slab pra-2011a . It is easy to see that the inclusion of the line defect yields other means of achieving a spectral singularity. In this case a spectral singularity occurs also for $\tilde{M}_{22}(\omega)=0$ which is equivalent to $\mathcal{Z}(\omega)=0$.

Recalling that the gain coefficient $g$ for the slab satisfies $g=-2k\kappa$, where $\kappa$ is the imaginary part of the refractive index $\sqrt{\varepsilon}$, we can follow the approach of pra-2011a to determine the threshold gain from $\mathcal{Z}(\omega)=0$. Using the fact that for generic active material, $|\kappa|\ll\eta-1$, where $\eta$ is the real part of $\sqrt{\varepsilon}$, we then find the following expression for the threshold gain coefficient

$$g=\frac{4\sqrt{\eta^{2}-\sin^{2}\theta}}{\eta L}\ln\left[\frac{\sqrt{\eta^{2}-\sin^{2}\theta}+|\cos\theta|}{\sqrt{\eta^{2}-1}}\right].$$

(36)

This is precisely the laser threshold condition for the oblique TE modes of the slab in the absence of the defect, if we identify $\theta$ with the incidence angle pra-2015a . Notice however that here we have a normally incident TE wave and $\theta$ corresponds to the direction of the scattered wave.

Fig. 3 shows a plot of the threshold gain coefficient (36) as a function of $\theta$ for a generic value of $\eta>1$.

As seen from this graph, $g$ attains its maximum value at $\theta=0,180^{\circ}$. If $\theta$ deviates from these values, $g$ decreases monotonically until it reaches its minimum value, which is zero, at $\theta=\pm 90^{\circ}$. This shows that if the slab consists of a gain medium, irrespective of how small its gain coefficient is, the defect will make it radiate laser light with a wavevector in some direction. This is actually easy to justify. The presence of the defect makes the waves scatter in different directions. Those that are scattered through angles closer to $\theta=\pm 90^{\circ}$ follow a longer optical path within the slab and lase more easily.

First we introduce some useful conventions and notions. We take the scattering axis of our system to be aligned along the $z$-axis of a cartesian coordinate system, label the position vector $\mathbf{r}$ with its cartesian coordinates $(x,y,z)$, and set $\vec{\rho}:=(x,y)$ and $(\vec{\rho},z):=(x,y,z)$. In general overarrows identify two-dimensional vectors lying in the $x$-$y$ plane.

Now, consider a potential $v:\mathbb{R}^{3}\to\mathbb{C}$ such that for $z\to\pm\infty$ every solution of the Schrödinger equation,

$$\left[-\partial_{x}^{2}-\partial_{y}^{2}-\partial_{z}^{2}+v(x,y,z)\right]\psi(x,y,z)=k^{2}\psi(x,y,z),$$

(37)

has the asymptotic form

$\displaystyle\frac{1}{(2\pi)^{2}}\int_{\mathscr{D}_{k}}d^{2}p\,e^{i\vec{p}\cdot\vec{\rho}}\left[A_{\pm}(\vec{p})e^{i\omega(\vec{p})z}+B_{\pm}(\vec{p})e^{-i\omega(\vec{p})z}\right],$

where $\mathscr{D}_{k}$ stands for the disc $\big{\{}\vec{p}\in\mathbb{R}^{2}\big{|}|\vec{p}|\leq k\big{\}}$, $A_{\pm}$ and $B_{\pm}$ are coefficient functions vanishing outside $\mathscr{D}_{k}$, and $\omega(\vec{p}):=\sqrt{k^{2}-\vec{p}^{\>2}}$. We can use (2) to define a transfer matrix $\mathbf{M}$ for the potential $v$, and check that it satisfies (II) if we replace $(x,p)$ by $(z,\vec{p})$, and define $\mathbf{H}(z,\vec{p})$ according to

$$\boldsymbol{\mathbf{H}}(z,\vec{p}):=\frac{1}{2\omega(\vec{p})}e^{-i\omega(\vec{p})z\boldsymbol{\sigma}_{3}}v(i\vec{\partial}_{p},z)\,\boldsymbol{\mathcal{K}}\,e^{i\omega(\vec{p})z\boldsymbol{\sigma}_{3}}.$$

(38)

Here $\vec{\partial}_{p}:=(\partial_{p_{x}},\partial_{p_{y}})$, $(p_{x},p_{y})$ are cartesian coordinates of $\vec{p}$, and $v(i\vec{\partial}_{p},z):=v(i\partial_{p_{x}},i\partial_{p_{y}},z)$. For details of the derivation of (38), see Appendix E.

According to (38), $\boldsymbol{\mathbf{H}}(z,\vec{p})$ vanishes for the values of $z$ for which $v(x,y,z)=0$. This together with Eq. (II) ensure that $\mathbf{M}$ satisfies the composition property (3).

In order to relate $\mathbf{M}$ to the scattering amplitude for the potential $v$, we introduce $T_{\pm}(\vec{p})$ by changing $p$ to $\vec{p}$ in (24). These have the same properties as their two-dimensional analog. We can express them in terms of $\mathbf{M}$ by enforcing (2) and setting

$\displaystyle A_{-}({\vec{p}})=4\pi^{2}\delta(p_{x})\delta(p_{y}),$

$\displaystyle B_{+}({\vec{p}})=0.$

In terms of $T_{\pm}(\vec{p})$ the reflected ($\psi_{-}$) and transmitted ($\psi_{+}$) waves take the form

$$\psi_{\pm}(\mathbf{r})=\frac{1}{4\pi^{2}}\int_{\mathscr{D}_{k}}d^{2}p\,T_{\pm}(\vec{p})e^{i\vec{p}\cdot\vec{\rho}}e^{\pm i\omega(\vec{p})z}~~~~{\rm as}~r\to\infty.$$

(39)

These in turn give the scattered wave $\psi_{\rm scat}$ according to (22) provided that we replace $x$ by $z$.

Next, we recall that $\psi_{\rm scat}(\mathbf{r})=f(\vartheta,\varphi)\,e^{ikr}/r$, where $f$ is the scattering amplitude and $(r,\vartheta,\varphi)$ are the spherical coordinates of $\mathbf{r}$ with $\vartheta$ and $\varphi$ respectively denoting the polar and azimuthal angles. Comparing this expression with the one given by (22) and (39) and using an argument given in Appendix F, we find

$$f(\vartheta,\varphi)=-\frac{ik|\cos\vartheta|}{2\pi}\,T_{\pm}(k\sin\vartheta\cos\varphi,k\sin\vartheta\sin\varphi),$$

(40)

where $\pm:={\rm sgn}(\cos\vartheta)$.

For the delta-function potential

$$v(x,y,z)=\mathfrak{z}\,\delta(x)\delta(y)\delta(z),$$

(41)

Eq. (38) gives $\boldsymbol{\mathbf{H}}(z,{\vec{p}})=\mathfrak{z}\,\omega(\vec{p})^{-1}\delta(z)\delta(i\partial_{p_{x}})\delta(i\partial_{p_{y}})\boldsymbol{\mathcal{K}}/2$. Substituting this relation in (II) and making use of (2), (25), and (40), we find

	$\displaystyle\mathbf{M}=\mathbf{I}-\frac{i\mathfrak{z}}{2\omega({\vec{p}})}\,\boldsymbol{\mathcal{K}}\delta(i\partial_{p_{x}})\delta(i\partial_{p_{y}}),$		(42)
	$\displaystyle f(\vartheta,\varphi)=\frac{-\mathfrak{z}}{4\pi+ik\mathfrak{z}}.$		(43)

In view of the latter relation the scattering length, i.e., $\xi:=-\lim_{k\to 0}f$, is given by $\xi=\mathfrak{z}/4\pi$. Therefore, we can express (43) as $f=-(\xi^{-1}+ik)^{-1}$. This coincides with Eq. (9) of Ref. rajeev-99 , which is obtained for general unitary interactions under the assumption that the range of the interaction is much smaller than the wavelength of the incident wave.

The standard physical application of delta-function potential is in the study of point defects in condensed matter systems. For a charged particle moving in a conducting medium, the presence of a point defect yields a non-vanishing resistivity, which is related to the scattering cross section of the defect $|f(\vartheta,\varphi)|^{2}$. Eq. (43) suggests that it is an isotropic function of space that is inversely proportional to $k^{2}+\mu^{2}$ where $\mu:=4\pi/|\mathfrak{z}|$.

Scattering of waves in one dimension is most conveniently studied in the framework of the transfer matrix formulation of scattering theory. This is mainly because of the practical advantages stemming from the composition property of the transfer matrix. Despite the enormous interest in the scattering of waves in two- and three-dimensions, there has been no successful attempt at developing a general transfer matrix formulation of scattering theory in these dimensions. In this article, we achieve this goal by introducing a genuine multi-dimensional notion of transfer matrix which enjoys the niceties of its one-dimensional analog.

In order to elucidate the practical aspects of our approach, we have explored its application in addressing the scattering problem for delta-function potentials in two and three dimensions. It is known that a proper treatment of the spectral problem for delta-function potential in two dimensions requires a coupling constant renormalization which makes it run with $k$ and introduces a length scale for the problem rajeev-99 ; camblong-2002 . As far as its physical applications are concerned, the delta-function potential is an abstraction for a physical potential whose range is much smaller than the relevant wavelengths. This provides an explanation for the length scale arising from the renormalization procedure and applies for the delta-function potential in three dimensions. It is remarkable that the direct application of our approach to delta function potentials yield finite results that do not rely on any approximation or renormalization schemes.

The composition property of the transfer matrix allows for extending our results to potentials consisting of finitely many delta-functions. In two dimensions, these describe finite arrays of parallel thin wires whose scattering properties are of direct interest in the study of metamaterials meta1 ; meta2 and nanophotonic devices nano1 ; nano2 ; nano3 . In three dimensions, delta-function potentials model point-like defects. This marks the significance of our approach for dealing with scattering problems in condensed matter physics. It is also not difficult to envisage a generalization of this approach that is capable of treating scattering of waves inside a waveguide. This will be of importance in the study of transport phenomena in mesoscopic systems MK .

To demonstrate the effectiveness of our approach in dealing with a realist physical system, we have used it to describe the scattering of electromagnetic ways from an optically active slab that has a surface line defect. Our results allow for the determination of the laser threshold condition for such a slab and reveal the fact that for any positive value of the gain coefficient the presence of the defect makes the slab begin emitting laser light in some direction. It is remarkable that this phenomenon does not depend on the strength of the delta-function interaction modeling the defect.

The main impact of the multidimensional transfer matrix formulation of scattering theory that we outline in this article will be in its numerical implementations. This is not only because of the composition property of the transfer matrix, which allows for the reduction of the scattering problem for a given potential to those of its truncations, but also because any discretization of this approach would yield matrix representations for the entries of the transfer matrix operator that could be substituted in (26) to determine $T_{\pm}$ and consequently the scattering amplitude.

In the following we give the details of various calculations that we refer to in the text of this article.