Helium Atom Numerical Solution TDSE Notes

             

Wave function expansion

　　 The angular basis is used as the eigenstate of the total orbital angular momentum of the two electrons and the component of $z$, which is the generalized spherical harmonic function $\mathcal{Y}_{l_1,l_2}^{L,M}( \hat{\boldsymbol{\mathbf{r}}} _1, \hat{\boldsymbol{\mathbf{r}}} _2)$.

　　 The space where the total wave function is located can be regarded as the tensor product space of angular space and radial space. Now that we have a base in angular space, the total wave function can be expanded on this base

Among them, $\lambda$ is the serial number after sorting all combinations of $(l_1,l_2,L,M)$. $ \left\lvert R_\lambda \right\rangle $ is the two-dimensional radial wave function That is, the total wave function is

　　 From the conclusion in the tensor product space , we can find the Hamiltonian $H$ about the "matrix element" $H_{\lambda, \lambda'}$ of the angular basis, and each matrix element is an operator in the radial space. List the "matrix form" of the Schrodinger equation, and you get a set of Schrodinger equations for the radial wave function of couple

　　 The total Hamiltonian can be expressed as

Where $H_i \ \ (i = 1, 2)$ is the Hamiltonian of a single electron, opposite to the angular matrix The second term is a diagonal matrix because the base of $\mathcal Y$ is the eigenfunction of $L_i^2$. $V_{12}$ is the Coulomb potential energy between two electrons, and $V_{Fi}$ is the potential energy of the electric field on each electron.

Symmetry of wave function

　　 We generally assume that the electron spin state is always singlet (antisymmetric), so the total wave function satisfies the particle exchange symmetry $\Psi( \boldsymbol{\mathbf{r}} _1, \boldsymbol{\mathbf{r}} _2) = \Psi( \boldsymbol{\mathbf{r}} _2, \boldsymbol{\mathbf{r}} _1)$. The symmetry of the partial wave can be derived from eq. 4

In theory, using this formula can reduce the repetitive partial wave in the program by half, and the coupling evolution of the partial wave can also be simplified, but it has not been realized due to the trouble of programming.

　　 In addition, it can be proved (actually discovered by numerical calculations) that $l_1 + l_2 - L$ is an odd partial wave (odd partial wave) and $l_1 + l_2 - L$ is an even partial wave (even partial wave) that evolve independently in the online polarization electric field, that is, $V_{12}$ and $z_1, z_2$ will not couple subwaves with different parities. So if the helium atom evolves from the ground state (including only the even partial wave: $l_1 = l_2, L = 0$) in the linear polarization electric field, then the odd partial wave term will always be zero and can be removed in the program.

Electronic Interaction Item

This is equivalent to decomposing the linear operator in the tensor product space into the sum of the tensor product of radial and angular operators. The second step uses eq. 9 , that is

　　 To calculate the matrix element $ \left\langle \mathcal Y_\lambda \middle| V_{12} \middle| \mathcal Y_{\lambda'} \right\rangle $, the linear combination of the products of the six spherical harmonic functions must be integrated twice. Using eq. 12 , the above formula can be expressed as a linear combination of six CG coefficients (double summation).

Now try to sum $m$ and use the choice theorem (the sum of the second line of the 3j symbol is 0), we get (note the range of $m$) at most only one item in the sum is not 0, so Among them, the sum requirement of $m_1, m_2$

　　 Since the CG coefficients are all real numbers, this must be a real number matrix. Selection rules are as follows

• According to the conservation of total angular momentum ($1/r_{12}$ operator and $L, L_z$ operator are exchangeable), when $L \ne L'$ or $M \ne M'$, the matrix element is zero.
• According to parity CG coefficients, when $l_1' + l + l_1$ or $l_2' + l + l_2$ is an odd number, the above formula is equal to zero.
• According to $\Pi \mathcal{Y} = (-1)^{l_1 + l_2} \mathcal{Y}$, considering that $1/r_{12}$ has even parity, the parity of $l_1 + l_2$ and $l_1' + l_2'$ must be different.

　　 In Aihua's paper, eq. 11 is represented by the 9j symbol, but it is limited to $M \ne M'$. However, after numerical verification, I found that it should be possible to extend Aihua’s formula to

Note that when $M = M'$, the result has nothing to do with $M$.

Poisson integral

　　 According to the idea of ​​the FEDVR basis, we would naturally want to use Gaussian Lobatto integral to calculate the $ \left\langle r_i r_j \middle| r_ < ^l/r_ > ^l \middle| r_{i'} r_{j'} \right\rangle $ after the expansion of the $V_{12}$ operator. However, this is not accurate, because the expansion of $r_ < ^l/r_ > ^l$ with a finite-order polynomial will result in error. Using the Poisson integral method can reduce the calculation error and still maintain the form of the diagonal matrix, that is

According to my guess (check it if you have time!), the result of Poisson integral integration should be the same as the result of direct numerical integration. Derivation refer to McCurdy, Solving the three-body Coulomb breakup problem using exterior complex scaling formula 68. Where $ \boldsymbol{\mathbf{T}} _l$ is the matrix of single electron kinetic energy operators Among them, $\phi_i$ is the orthogonal normalized base of the one-dimensional FEDVR grid point at $r_i$, and $r_{max}$ is the boundary of the last FE when calculating $T_l$. Although eq. 15 contains $r_{max}$, the result has nothing to do with $r_{max}$. The reason is that $T_l$ is also related to $r_{max}$, and the two items are offset by $r_{max}$.

Electric field action term

　　 For each electron,

Same as in the case of hydrogen atoms, there are For the linearly polarized electric field, we only need the last item, taking the first electron as an example, The last step uses the properties of $m_1 + m_2 + m_3 = 0$ in the 3j symbol and $m_1 + m_2 = M$ in the CG coefficient. The two delta functions verify that the $z$ directional electric field acting on the first electron will not affect the second electron’s angular momentum and the $z$ component of the total angular momentum. In the derivation, the conclusion of $m_1' = m_1, m_2' = m_2$ has been successively drawn, that is, the electric field will not affect the angular momentum $z$ component of each electron (so the sum of $m_1', m_2'$ is eliminated). In the final summation of $m_1,m_2$, $m_1+m_2=M$ is required, and $ \left\lvert m_1 \right\rvert < \min \left\{l_1, l'_1 \right\} $.

　　 In the same way, for the second electron

　　 In Aihua’s paper, $M = M' = 0$ is represented by the 9j symbol as (the numerical verification is the same)

Nonlinear polarization electric field

　　 For the nonlinear polarization electric field, the first two items of eq. 17 are not zero, but are only the same as the result of the $z$ component (only the last 3j coefficient becomes two items, and then divide by $\sqrt{2}$ to normalize)

Binding State

　　 The lower-energy bound state can be easily obtained by imaginary time evolution, that is, replacing $ \exp\left(- \mathrm{i} H \Delta t\right) $ with $ \exp\left(- H \Delta t\right) $, which is equivalent to changing $\Delta t$ into an imaginary number $- \mathrm{i} \Delta t$, hence the name. The ground state of helium atom only contains the partial wave of $L = 0$, which requires $l_1 = l_2 = 0, 1, 2, \dots$ (refer to fig. 2 ).

wave function evolution

　　 Evolution (propagation) also uses split operator

It's just that $H_{int} = V_{12} + V_F$ here is an item of non-diagonal (operator) matrix after all $ \left\langle \mathcal{Y} \middle| \ \middle| \mathcal{Y} \right\rangle $ acts, that is, the interaction of electrons and the effect of electric field on each electron.

　　 $H_0 = H_1 + H_2$ Since $H_1$ and $H_2$ are exchangeable, we can evolve them independently

That is, for the two-dimensional grid of each partial wave, use a method similar to hydrogen atoms to evolve each column and then evolve each row.

　　 According to the above derivation, it can be seen that after $H_{int}$ is expanded into a linear combination of the tensor product of the radial operator and the angular operator, the radial operators are all functions of $r$, which is all for the FEDVR base Diagonal matrix (even if Poisson integral is used). So when evolving $ \exp\left(- \mathrm{i} H_{int}\Delta t\right) $, each two-dimensional FEDVR grid can be performed independently, which provides the possibility of a large number of parallelism. In summary, the evolution in each direction of the three-dimensional wavefunction matrix can be performed independently for each column in the program.

　　 When evolving $ \exp\left(- \mathrm{i} H_{int}\Delta t\right) $, Aihua's approach is to further split into three operators

This is because $V_F = V_{F1} + V_{F2}$ is time dependent and proportional to the electric field component $E_z$.

　　 The reason for the further division of Aihua is that his code prestores the eigenvalues ​​and eigenvector matrices after the diagonalization of $V_F$ and $V_{12}$. If you directly diagonalize $H_{int}$, then you need to re-align each step. Keratinization once. The memory required for pre-storage is proportional to the square of the number of subwaves, but when we need a lot of subwaves, the memory is far from enough. In order to solve the memory problem, I use Crank-Nicolson evolution both radially and angularly, and regenerate the $V_F + V_{12}$ matrix for each grid point at each time step, and one evolution does not need to be divided into three times. For the grid points whose absolute value of the wave function is lower than a certain threshold, we can even ignore and not evolve.

　　 Another point, if you only care about single ionization, you can use a rectangular grid (see Ossiander 2017 Nat. Phys.). But pay attention to adding absorber on the short side to prevent reflection. )

Partial wave convergence analysis

　　 Using $l_{1,max}, l_{2,max}, L_{max}$ to generate the partial wave list is sometimes wasteful, and there will be many partial waves with extremely low probability. The best way is to plot the probability of each partial wave into a graph, from which you can easily see whether the partial wave has converged and whether there is waste.

Partial wave error analysis (original)

　　 When we intercept a finite number of partial waves, new errors will be introduced during the evolution of the partial waves (eq. 29 ). After interception, $V_F, V_{12}$ is still a Hermitian matrix, so their propagators are also unitary, so the modulus length of the wave function cannot be used to estimate the error. So another method can be used to estimate the error. The error can be estimated every few steps. Let $ \boldsymbol{\mathbf{V}} '_F, \boldsymbol{\mathbf{V}} '_{12}$ be a matrix with twice the number of bases, then you can calculate the number of sub-waves more than $ \boldsymbol{\mathbf{(}} \mathrm{i} V'_F\Delta t) \boldsymbol{\mathbf{v}} $ than the base, and sum the absolute value of each sub-wave.

Binding State

　　 Theoretically, the eigenvalues ​​of the total Hamiltonian matrix can be solved numerically, but in practice this matrix is ​​so large that it is almost impossible to solve it (in fact, we never store the total Hamiltonian matrix directly, but store the split operator). We use imaginary time propagation .

　　 For the initial state of virtual-time evolution, I use the multiplication of the ground state of two He+ (real wave function), which is slightly different from Aihua, but this has little effect. Since our split operators are all real-numbered matrices, the exponential function is also a real-numbered matrix after using imaginary time, so the ground state wave function obtained by the imaginary time evolution must also be a pure real number. We only need to use all the partial waves of $L = 0, M = 0$.

Double electron approximation

　　 Argon’s number of electrons is $N = 18$, ground state ionization potential energy is $I_p = 15.6 \,\mathrm{eV} $ ($0.5733 \,\mathrm{au} $), using equivalent potential energy

Satisfy $V_{eff}(r\to\infty) \to -2/r$, when $V_{eff}(r\to 0) \to -N/r$. The smaller the parameter $a$, the smaller the radius of the shielding potential energy, and the closer to the helium atom, the greater the ionization potential energy.