I want to put (non-dirichlet) boundary conditions inside the code I wrote to solve the 2dim TDSE using the alternating direction implicit Peaceman - Rachford method. $$ (1 + iB\Delta t/2 ) \psi^{n+1/2} = (1 - iA\Delta t/2) \psi^n $$ $$ (1 + iA\Delta t/2 ) \psi^{n+1} = (1 - iB\Delta t/2) \psi^{n+1/2} $$, where $A$ contains operators only on the coordinate $x$ and $B$ contains operators only on the coordinate $y$ and $\psi = \psi(x,y,t)$.
This PR method, as applied to my particular case, is detailed in https://www.sciencedirect.com/science/article/pii/S0021999196955886 and in the Appendix of http://www.iap.tu-darmstadt.de/tqp/papers/BauMul99.pdf.
I have come across http://rrp.infim.ro/2003_55_4/carjan.pdf, from now on called Carjan 2003 which details a particular approach to implement Transparent Boundary Conditions for the 1 dim TDSE. Because the ADI basically solves a one-dim problem at each half timestep (from $t_n$ to $t_{n+1/2}$ and from $t_{n+1/2}$ to $t_n$) via Crank Nicolson in 1D, I believe I can implement BCs from 1dim TDSE studies. Am I correct?
If I am correct and I can implement them to my 2D problem, I am having trouble understanding how to treat the left boundary of a 1dimensional line (x-coordinate). Carjan 2003 only explains how to treat the right boundary (page 560, or 6/25 from the pdf). I have also went to the source of these Transparent Boundary Conditions described in Carjan 2003: the original (quoted as [11, 12] in the above link) is a treatment of the propagation of a laser beam in space [nothing related to quantum mechanics, but there is a nice linkage between the two, made in the publication from link above]. Still, nothing related to the left boundary, only the $x = x_{max}$ boundary (RHS one).
In the paper, the right boundary is treated by firstly assuming that the waves near it are of the form $\psi_0 \exp(i k_x x)$, where $i^2=-1$, and $k_x$ and $\psi_0$ complex constants to be determined. Then the continuity equation for the wavefunction $\psi$ is applied (1-dim problem) and integrated across the simulation domain: $$ \rho = \int_{a}^{b} \vert \psi \vert^2 dx $$ $$ \dfrac{\partial \rho}{\partial t} = \dfrac{i \hbar}{2 m} \left[ \psi^* \dfrac{\partial \psi}{\partial x} - \psi \dfrac{\partial \psi^*}{\partial x} \right] \vert_{a}^{b} \equiv -F_b + F_a $$
They say that $F_b$ signifies the "energy" flux leaving through the right boundary and $F_a$ is the "energy" flux entering through the left boundary. It is then wanted to obtain the $k_x$ such that $-F_b$ is positive. They say that this implies that the contribution to the overall change in energy from this RHS boundary will be always negative (radiative energy can only flow out of the problem region). They obtain a $k_x$ s.t. respect the above, in particular the condition is: real part of $k_x$ shall be positive to ensure that radiation only flows outward at the RHS boundary.
I have looked here https://www.youtube.com/watch?v=Ls5HS2MLXpg and I understand that the continuity equation arises due to a fluid flowing from the LHS into a volume (my line domain) and going out on the right (for no creation/destruction of the fluid). I can also understand that the derivation can be applied the other way around: the fluid coming into the volume from the right and leaving through the left boundary.
My 2D TDSE problem, my aim, for which I want to implement BC's, has a wavefunction $\psi(x,y,t)$ which spreads from its initial condition (localized at a central point) across the whole simulation domain as time passes. It delocalizes over time. I need to remove (absorb, not reflect, as much as possible) the wavefunction which going out of the boundaries of the simulation box. So the $x_{min}$ and $x_{max}$ boundaries shall both have such transparent boundary conditions. They both shall remove any incoming wavefunction.
My question is: in the derivation for the left boundary, how shall I proceed to do it? I am particularly concerned about the starting assumption: in the paper, for the RHS boundary, they said: waves near the boundary are of the form $\psi_0 \exp(i k_x x)$. Now, when treating the LSH boundary, shouldn't the waves be of the form: $\psi_0 \exp(- i k_x x)$. [Tangent question: where has the time dependence disappeared in their ansatz?] Then the continuity equation, after it's been integrated as (see the boundaries, changed): $$ \rho = \int_{b}^{a} \vert \psi \vert^2 dx $$ shouldn't it read: $$ \dfrac{\partial \rho}{\partial t} = \dfrac{i \hbar}{2 m} \left[ \psi^* \dfrac{\partial \psi}{\partial x} - \psi \dfrac{\partial \psi^*}{\partial x} \right] \vert_{b}^{a} \equiv -F_a + F_b $$. But if I assumed $exp(-i k_x x)$ for the waves near the LHS boundary, the derivation, from now on, is not similar to the derivation for the RHS boundary presented in the papers. Again, I need to make $F_a$ > 0, but this time i have that -1 in the exponential. What I want to emphasize is that I believe the derivation shall be more or less identical w.r.t. which boundary we consider, but I don't get this. In the end, eq. (3.7) from Carjan 2003 shows the result of putting a transparent boundary condition at the RHS boundary: it changes the last equation from the Crank-Nicolson scheme (i.e. the one for the node prior to the very last one, i.e. the node before the boundary one). When I work with the LHS boundary, the first equation from the Crank-Nicolson scheme will have to change. Thus my concerns because a mistake in the derivation will ruin everything.
Thank you!
