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!

  • $\begingroup$ As a personal opinion, I'm not a big fan of such transparent boundary conditions, as they have to be adjusted to the problem at hand. On the other hand, the alternative in the paper, absorbing boundaries, is conceptually simple and can be readily applied to a wide range of problems (also your ADI-stuff). Moreover, the performance gain they obtained in the paper by comparing transparent vs. absorbing was only 5/4 ... taken this factor for granted, it's nothing to really aim for, and it will go under in the effects of all the other implementation decisions. $\endgroup$
    – davidhigh
    Commented Feb 9, 2022 at 11:35
  • $\begingroup$ Thanks, I will have a look at them just now. Yes, I am not that much interested in performance gains, I just want something to work reasonably well ... This method in the post which I wrote, however, seems not to have something problem - dependent! $\endgroup$
    – velenos14
    Commented Feb 9, 2022 at 11:50
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Commented Feb 9, 2022 at 16:47
  • $\begingroup$ is this automatic? the bot thing? everytime I post a question, this appears. I highlighted with bold, it says my question is: ... $\endgroup$
    – velenos14
    Commented Feb 9, 2022 at 17:18
  • $\begingroup$ @davidhigh, I just want to clarify if I am allowed to use 1d boundary conditions on my 2d ADI scheme. at the first half-timestep propagation, from $t_n$ to $t_{n+1/2}$, if $B$ is an operator with derivatives only on $y$-coordinate, then can I apply the BCs for a fixed x-coordinate, for the $y = y_{min}$ (lower) and $y = y_{max}$ (upper) nodes? Remember, these BCs are not Dirichlet. Similarly, for the second half-timestep propagation, from $t_{n+1/2}$ to $t_n$, if $A$ is acting only on $x$-coordinate, then can I apply BC's for a fixed y-coordinate, for the $x=x_{min}$ and $x=x_{max}$ nodes? $\endgroup$
    – velenos14
    Commented Feb 10, 2022 at 22:21


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.