Consider
$$ (-\frac{\partial^2}{{\partial x}^2}-\frac{\partial^2}{{\partial y}^2}+2c\delta (x - y))\psi(x,y) = E \psi(x,y)$$
with boundary conditions $\psi(0,y) = \psi(L,y) = \psi(x,0) = \psi(x,L) = 0$
The equation simplifies to not include the delta term in the regions $x < y$ and $x > y$. $\psi(x,y) = (Ae^{ik_x x}+Be^{-ik_x x})(Ce^{ik_y y} + De^{-ik_y y})$
Using the boundary conditions on the solution of the $x < y$ region $$\psi(x,y) = A(e^{ik_x x} - e^{-ik_x x})(e^{ik_y y}-e^{ik_y(L-y)})$$ Similarly in the region $y > x$ $$\psi(x,y) = B(e^{ik_x x} - e^{ik_x (L -x)})(e^{ik_y y}-e^{-ik_y y})$$
However I am struggling with imposing the boundary condition that occurs along the delta function discontinuity $y = x$.
Consider an integral along the y - direction passing the discontinuity at $x$, with limits $x - \epsilon$ and $x + \epsilon$, where $\epsilon$ is small.
$$\int_{x - \epsilon}^{x + \epsilon} (-\frac{\partial^2}{{\partial x}^2}-\frac{\partial^2}{{\partial y}^2}+2c\delta (x - y))\psi(x,y)dy = \int_{x - \epsilon}^{x + \epsilon} E\psi(x,y)dy$$ simplifies to $$\int_{x - \epsilon}^{x + \epsilon}-\frac{\partial^2}{{\partial x}^2}\psi(x,y)dy- \frac{\partial \psi(x,y)}{\partial y}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon}+2c\psi(x,y) =0$$ The first term on the LHS cannot be simplified further without use of the Leibnitz rule. The derivative with respect to $x$ cannot be pulled out of the integral naively since the limits of integration include dependence on $x$.
$$\frac{d}{dx}\int_{x - \epsilon}^{x + \epsilon} \psi(x,y) dy = \psi(x, x + \epsilon) - \psi(x, x - \epsilon) + \int_{x - \epsilon}^{x + \epsilon} \frac{\partial}{\partial x} \psi(x,y)dy$$
Differentiating with respect to $x$ again, $$\frac{d^2}{dx^2}\int_{x - \epsilon}^{x + \epsilon} \psi(x,y) dy = 2\frac{\partial \psi(x, x + \epsilon)}{\partial x}+\frac{\partial \psi(x, x + \epsilon)}{\partial y} - 2\frac{\partial \psi(x, x - \epsilon)}{\partial x} - \frac{\partial \psi(x, x - \epsilon)}{\partial y} + \int_{x - \epsilon}^{x + \epsilon} \frac{\partial^2}{{\partial x}^2} \psi(x,y)dy = 2\frac{\partial \psi(x,y)}{\partial x}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon} + \frac{\partial \psi(x,y)}{\partial y}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon}+\int_{x - \epsilon}^{x + \epsilon} \frac{\partial^2}{{\partial x}^2} \psi(x,y)dy$$ In the limit of vanishing $\epsilon$, $$\int_{x - \epsilon}^{x + \epsilon} \psi(x,y) dy \rightarrow 0$$ Hence $$-\int_{x - \epsilon}^{x + \epsilon} \frac{\partial^2}{{\partial x}^2} \psi(x,y)dy = 2\frac{\partial \psi(x,y)}{\partial x}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon} + \frac{\partial \psi(x,y)}{\partial y}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon}$$ Using the previously obtained expression of the integral over the discontinuity $$2\frac{\partial \psi(x,y)}{\partial x}\Bigg|_{y = x - \epsilon}^{y = x + \epsilon} + 2c\psi(x,y)=0$$ Is this condition correct? If yes, how do I proceed with imposing the boundary conditions on the obtained standing wave solutions?
If I attempt to impose continuity $$\psi(x,x) = A(e^{ik_x x} - e^{-ik_x x})(e^{ik_y x}-e^{ik_y(L-x)}) = B(e^{ik_x x} - e^{ik_x (L -x)})(e^{ik_y x}-e^{-ik_y x})$$ it seems that this will only work if $k_x = k_y$, and $A = B$. These conditions seem to be too strong to allow the discontinuity condition to be applied. Hints please?
