# How to solve the electric potential due to a charge distribution using fast Fourier transforms

I'm writing a toy program to simulate 2D wave functions. I'm using a split-operator method to solve the Schrödinger equation and have no problems with arbitrary potentials. However, I'd now like to have the wave functions react to the electric potential due to a 2D charge distribution.

The Poisson equation relates the potential with a charge distribution (I'm ignoring possible constants here since I'm not interested in actual units)

$$\nabla^2 V = -\rho$$.

If we take the Fourier transform we get

$$-k^2 \hat{V} = -\hat{\rho} \Longleftrightarrow \hat{V} = \hat{\rho}/k^2$$.

I'm well aware that at the origin $$k^2 = 0$$ and that the zero mode corresponds to the average value of $$V$$. Normally I'd just set $$\hat{\rho}/k^2$$ to zero there since the absolute value of the potential doesn't matter in most cases - it's its derivatives that appear in the equations.

Now, however, as far as I understand, in the Schrödinger equation we have

$$\left(-\nabla^2 + V\right)\psi$$,

meaning that it is actually the potential, not its derivatives, that shows up in the equation. How do I work out the proper value for $$\hat{\rho}/k^2$$ at the origin? Is it somehow related to the total charge in the system? My plan B is to solve $$V$$ using Jacobi over-relaxation, but I'd really like to stick with FFTs that are very efficient and noniterative.

I tried searching for previous questions and their answers here and elsewhere but couldn't find the answer I was looking for (or possibly identify one as such). Thanks a lot in advance!

$$\newcommand\ket[1]{|#1\rangle}$$The absolute part of $$V$$ is still irrelevant in the Schrödinger equation. If you have an initial state $$\ket\psi,$$ then its time evolution is $$\ket{\psi(t)}=\exp(-iHt/\hbar)\ket\psi=\exp(-i(-\nabla^2+V)t/\hbar)\ket\psi.$$ If you add a constant to $$V,$$ so $$V\mapsto V+E$$ (equivalently, add a constant to $$H$$ for $$H\mapsto H+E$$), then the new time evolution is just (remember $$\exp(A+B)=\exp A\exp B$$ when $$A,B$$ commute) $$\ket{\psi'(t)}=\exp(-iEt/\hbar)\ket{\psi(t)}.$$ Now remember that quantum states are only physical up to global phases: $$\ket\psi$$ and $$e^{i\theta}\ket\psi$$ always represent the same physical situation. In particular, all observables (Hermitian operators) have the same probability distributions. So at every time $$t,$$ the evolved state with the shifted potential $$\ket{\psi'(t)}$$ is the same physical state as the state resulting from the unshifted potential $$\ket{\psi(t)}$$: the absolute values of $$V$$ do not matter, even for the Schrödinger equation.

So the reason you didn't find many answers to your question is because there's no question to answer. You can do whatever is most convenient for handling the 0-frequency component of $$V$$ and you will still get all the right physics.

• Thanks a lot for the reply! But isn't $\theta$ a function of $t$ then? Won't that mean that the phase varies faster and the wave function has more momentum? What if I were to couple multiple wave functions together via their electric potentials (the potential for each wave function is a sum of those of all the others)? The charge distribution of each wave function could be $c \left| \psi \right|^2$ where $c$ is a constant specific to each wave function and integrating the norm over the entire system yields 1. I'd imagine I'd need at least a consistent way of handling the zero modes? Dec 1, 2022 at 16:51
• The momentum ($\hat p=-i\hbar\nabla$) no longer involves a time-derivative in QM. The distribution of any observable at a time $t$ depends only on the wavefunction at time $t.$ I'm not sure what you're doing with the coupled system. Without approximations, for $n$ coupled particles there is one wavefunction in $n$ variables. E.g. for two coupled particles you have a wavefunction $\psi(t,x_1,x_2)$ satisfying $i\hbar\frac{d[\psi(t,x_1,x_2)]}{dt}=(-k_1{\nabla_{x_1}}^2-k_2{\nabla_{x_2}}^2+kq_1q_2/|x_1-x_2|+q_1V(x_1)+q_2V(x_2))\psi(t,x_1,x_2).$ Shift $V$ by a constant and nothing happens.
– HTNW
Dec 1, 2022 at 20:55
• To be honest I'm not sure what I'm doing either but thanks again for the explanation. I assume the $k$'s are constants and the Laplacians are only with respect to the variables in their subscripts? Are the $V$'s particle-specific external potentials and the middle term gives the electric potential between the particles? Dec 2, 2022 at 22:15
• I was thinking of simulating $n$ wave functions for $n$ particles (again omitting constants for brevity/laziness): $i\partial\psi_i/\partial t=\left(-\nabla^2+\sum_{j\neq i}q_i q_j V_j \right)\psi_i$ where $V_j$ is given by the charge distribution $\rho_j$, which again is given by $c \left|\psi_j\right|^2$. This seems to give me something like particles bouncing off each other so that's quite promising. This might be completely wrong of course. I'm a bit confused how to interpret an $n$-dimensional wavefunction for $n$ particles when my system is 2-dimensional for any number of particles. Dec 2, 2022 at 22:21
• It only makes physical sense to have one external potential $V,$ shared between all particles! Otherwise you've interpreted me correctly. For $d$-dimensional space, the wavefunction of $n$ particles is a function of $n$ $d$-dimensional positions, for a total of $nd$ real inputs. Rarely, the particles will not be entangled and the wavefunction may separate as $\psi(x_1,\ldots,x_n)=\psi_1(x_1)\ldots\psi_n(x_n),$ but the presence of interactions means this kind of state should not last. Your method of simulation may get some qualitative details right but I don't think it can be correct.
– HTNW
Dec 3, 2022 at 1:06