Lippmann–Schwinger Equation for 2B Bound States in 2D

The starting point for a 2B bound state with reduced mass \(\mu\) is the time-independent Schrödinger equation \[(H_0 + V)\,|\Psi\rangle = E\,|\Psi\rangle\]where\(H_0\) is the free Hamiltonian, \(V\)is the interaction, \(E<0\) is the binding energy, and \(|\Psi\rangle\)is the corresponding bound state wave function. Rearranging the equation above gives the homogeneous LS equation for a bound state\[|\Psi\rangle = G_0(E)\,V\,|\Psi\rangle, \qquad G_0(E) \equiv \frac{1}{E - H_0}\] which expresses the state in terms of \(G_0(E)\)and the interaction. In a momentum basis\(|\mathbf p\rangle\), adopting the completeness the free propagator\(\int \! d^2 p\, |\mathbf p\rangle\langle \mathbf p| = 1\), the LS equation projects to\[\psi(\mathbf p) = \frac1E - \dfrac{p ^ 2}{2\mu}\, \int_0^\infty \! dp'\, p' \int_0^2\pi \! d\phi'\, V(\mathbf p,\mathbf p')\, \psi(\mathbf p')\]where \(\psi(\mathbf p) \equiv \langle \mathbf p| \Psi\rangle\) and\(V(\mathbf p,\mathbf p') \equiv \langle \mathbf p| V | \mathbf p' \rangle\). For central interactions in 2D, it is convenient to separate the angular dependence by a partial-wave (PW) expansion, which reduces LS equation to a one-dimensional integral equation for each angular-momentum channel \(m\)\[\psi_m(p) = \frac1{E - \dfrac{p ^ 2}{2\mu}} \int_0^\infty \! dp'\, p'\, V_m(p,p')\, \psi_m(p').\]The PW-projected kernel can be obtained from the angular integration in momentum space, \[\begin{aligned}V_m(p, p') = \int_0^{2\pi} d\phi' \, V(p, p', \phi') \cos(m \phi') \equiv \int_0^\infty dr r \, J_m(pr) V(r) J_m(p'r).\end{aligned}\]with \(J_m\) the Bessel function of the first kind. The PW representation for each \(m\) reduces the 2D integral equation to a one-dimensional radial form. This framework is used throughout to compute 2D bound states, e.g., hydrogenic systems with Coulomb interactions and the deuteron with Yukawa-type interactions, and, in particular, excitons with the RK interaction, within a momentum-space LS scheme.

A convenient internal consistency check for a numerically obtained 2B binding energy and wave function is to compare the eigenvalue \(E\) from the LS equation with the Hamiltonian expectation value, \(\langle H \rangle = \langle H_0 \rangle + \langle V \rangle\), evaluated using the corresponding momentum-space wave function. The degree of agreement between \(E\) and \(\langle H \rangle\) provides a direct measure of the numerical accuracy of the method. For a given PW channel\(m\), with the wave function normalization \[\int \! d^2p\, \psi^2(\mathbf p) = 2\pi \int_0^\infty \! dp\, p\, \psi_m^2(p) = 1.\] The Hamiltonian expectation value, expressed in terms of the kinetic and potential contributions, reads\[\begin{aligned}\langle H \rangle &= \langle H_0 \rangle + \langle V \rangle, \\[3pt] \langle H_0 \rangle &= 2\pi \int_0^\infty \! dp\, p\;\frac{p ^ 2}{2\mu}\, |\psi_m(p)|^2, \\[3pt] \langle V \rangle &= 2\pi \int_0^\infty \! dp\, p \int_0^\infty \! dp'\, p'\; \psi_m(p)\, V_m(p,p')\, \psi_m(p') \,.\end{aligned}\]

Numerical Implementation

To numerically solve the LS integral equation for a 2B bound state in 2D, represented in the PW framework, the integral equation should be transformed into a discretized form that can be solved iteratively or through matrix diagonalization. To efficiently handle the integrals, variable transformations are introduced to map the infinite domains of momentum \(p\), spatial coordinate \(r\), and angle\(\phi\) onto finite intervals. The transformations are defined as: \[p,p', r = \frac{1 + x}{1 - x}, \quad \phi' = \pi (1+x), \quad x \in [-1,1],\]where \(x\) is a variable spanning the interval \([-1, 1]\), suitable for integration using Gaussian-Legendre quadrature. The discretized form of the integral equation for the wave function \(\psi_m\) becomes an eigenvalue problem:\[\lambda \ \psi_m = {\cal K}(E) \, \psi_m,\] where \({\cal K}(E)\)represents the integral kernel operator dependent on the binding energy \(E\). The eigenvalue \(\lambda\) is used as a criterion to verify the solution. The binding energy \(E\) is determined by solving for the eigenvalue \(\lambda = 1\). This involves iteratively searching over initial guesses for \(E\)until convergence is achieved. A tolerance of \(10^-6\) is set for the eigenvalue condition, ensuring high accuracy in the computed energy.
The iterative search procedure can be outlined as follows:

1. Select two initial guesses for the binding energy E.

2. Construct the kernel \({\cal K}(E)\) for each chosen energy.

3. Solve the eigenvalue problem \(\lambda \, \psi_m = {\cal K}(E) \, \psi_m\) to obtain \(\lambda\) and the corresponding wave function \(\psi_m\) for both initial guesses.

4. Using these two initial pairs \((E, \lambda)\), perform a linear fit to predict the next energy value for which \(\lambda = 1\).

5. Iterate the process with the updated energy values until \(\lambda\) converges to \(1\) within the specified tolerance.

The kernel \({\cal K}(E)\) is diagonalized at each step, which allows for the extraction of both the eigenvalue \(\lambda\) and the wave function \(\psi_m\). Efficient numerical techniques such as matrix diagonalization, implemented through libraries like LAPACK, or iterative solvers like the Lanczos algorithm, can be employed for this purpose.

Hydrogenic atom in 2D

We begin with the 2D hydrogenic system, a well-established case that provides an ideal benchmark for testing the accuracy of our numerical method. In contrast to its three-dimensional analogue, the 2D system exhibits distinctive features in its energy spectrum and wave function behavior arising from the reduced dimensionality. In this section, we present both analytical and numerical treatments of the problem, working in natural units where \(\hbar c = e = m = 1\). The 2B potential energy between a nucleus of charge \(Z\) and a single electron in configuration space is given by \[V(r) =-\dfrac{Z}{r},\] which can be expressed in momentum space as\[V({\mathbf p}, \mathbf p') =\dfrac{-Z}{2\pi \vert {\mathbf q}\vert}, \quad \mathbf q= \mathbf p-{\mathbf p}'.\]The analytical solutions of the LS integral equation for the 2B binding energy levels are given in Ref. \[E_{exact} = -\dfrac{Z ^ 2 m}{4\left(n + 1/2\right)^2}, \quad n = 0, 1, 2, \ldots\] The exact energy levels for 2D hydrogenic atoms offer a clear reference for comparison. Evaluating the agreement between our computed values and the known analytical results allows us to assess the reliability of the approach and identify potential issues at an early stage. Establishing this validation is essential before extending the method to more complex systems, such as excitons or deuterons. The results of this comparison are presented in Table 1, showing that our numerical values closely match the exact energy levels of the 2D hydrogenic atom. The relative errors are very small, less than 0.4% for the ground state and the first three excited states, demonstrating the accuracy and reliability of our method.

Table 1 Comparison of analytical (\(E_{num}\)) and numerical (\(E_{num}\)) 2B binding energies, obtained from Eqs. (14) and (4), respectively, for the 2D hydrogenic atom with partial wave \(m=0\) and nuclear charge \(Z=2\). Calculations employ \(500\) grid points for the magnitude of the relative momentum and \(40\) points for the angular variables.

n		\(E_exact\)		\(E_{Num}\)		\(\left(\dfrac{E_{exact} - E_{num}}{E_{exact}}\right) \times 100 \%\)
0		-4.00000		-3.99928		0.01800
1		-0.44444		-0.44441		0.00700
2		-0.16000		-0.16017		0.10625
3		-0.08163		-0.08193		0.36751

Fig. 1 complements the binding‐energy results in Table 1 by showing momentum-space hydrogenic wave functions for \(n=0-3\). The plots depict \(\psi^{(n)}_{m}(p)\) versus relative momentum \(p\), highlighting how the momentum distribution varies with the energy level. The ground state exhibits a broader high-\(p\)distribution, whereas excited states (\(n=1,2,3\)) develop nodes and shift weight toward lower \(p\), consistent with their larger spatial extent in configuration space.

Figure 1 Ground (\(n=0\)) and first three excited (\(n=1-3\)) momentum-space wave functions \(\psi^{(n)}_{m}(p)\) of the 2D hydrogenic atom for partial wave \(m=0\) and nuclear charge \(Z=2\) as functions of relative momentum\(p\). States correspond to the numerical binding energies in Table 1

Deuteron in 2D

This section presents the numerical simulation of the deuteron, the bound state of one proton and one neutron, in 2D using the MT potential, written as a sum of repulsive and attractive Yukawa terms, \[V(r) = V_r \frac{e ^ {-\mu_r r}}{r} + V_a\frac{e ^ {-\mu_a r}}{r},\] where \(V_r\) and \(V_a\) are the strengths of the repulsive and attractive components, respectively, and \(\mu_r\) and \(\mu_a\) are the corresponding range parameters. The potential transforms to momentum space as \[V(\mathbf p, \mathbf p') = \frac{1}{2\pi}\left( \frac{V_r}{\sqrt{q ^ 2 + \mu_r^2}} + \frac{V_a}{\sqrt{q ^ 2 + \mu_a^2}} \right), \quad q = | \mathbf p- \mathbf p' |.\]The parameter sets for the two MT models used in our calculations are listed in Table and the corresponding momentum-space potentials \(V(q)\) are plotted in Fig. as functions of the momentum transfer \(q\).

Table 2 Parameter sets for the two MT potential models used in our calculations.

MT-model	\(V_a \, \text{(MeV fm)}\)	\(\mu_a \, (\text{fm}^{-1})\)	\(V_r \, \text{(MeV fm)}\)	\(\mu_r \, (\text{fm}^{-1})\)
Model-1	\(-600.00\)	\(1.550\)	\(1438.7228\)	\(3.21\)
Model-2	\(-600.00\)	\(1.550\)	\(0\)	\(0\)

Figure 2 Comparison of the two MT potentials as functions of the momentum transfer \(q\).

In our calculations, we use the nucleon mass with \(\hbar^2/m = 41.47\ MeV \text·fm^2\). As shown in Table , the computed binding energies for the 2D deuteron differ markedly between the two MT models. This is expected: Model 1, which superposes repulsive and attractive Yukawa terms, supports only a single (ground) bound state, whereas Model 2, which is purely attractive, supports two deeper bound states (ground and first excited). The corresponding momentum-space wave functions are shown in Fig.The left panel (Model 1) displays a single ground-state wave function, \(\psi(p)\), versus relative momentum \(p\). The right panel (Model 2) presents two ground and excited bound states, labeled \(n=0\) and \(n=1\).

Table 3 Calculated 2D deuteron binding energies for the two MT potential models, using \(200\) grid points for the 2B relative-momentum magnitude and \(101\) grid points for the angular variables.

\(E_{d}\) (MeV)		\(E^{(0)}_{d}\) (MeV)	\(E^{(1)}_{d}\) (MeV)
Model-1		Model-2
\(-6.246\)		\(-7806\)	\(-312.4\)

Figure 3 2D momentum-space deuteron wave functions \(\psi(p)\), computed with the MT potential: Model 1 (left) and Model 2 (right), plotted versus relative momentum\(p\).

Excitons in 2D Materials

Having validated the framework on 2D hydrogenic atoms and the deuteron, we now apply it to excitons-bound electron-hole pairs formed by Coulomb attraction that play a central role in the optical response of 2D materials. In this section we solve the LS integral equation for excitons using the RK potential. This yields exciton binding energies and wave functions, providing insight needed to predict exciton behavior in practical devices and to inform the design of materials with targeted optical properties. The RK potential accounts for nonlocal, distance-dependent dielectric screening in 2D materials, and thus models electron-hole interactions more accurately than a bare Coulomb potential. The explicit forms of the RK potential in configuration and momentum space are\[\begin{aligned}V(r) &=& -\dfrac{\alpha \hbar c \pi}{2 r_0 \kappa}\left[H_0\left(\dfrac{r}{r_0}\right)-Y_0 \left(\dfrac{r}{r_0}\right)\right], \cr V(\mathbf p,\mathbf p') &=& - \dfrac{1}{4\pi^2}\left( \dfrac{1}{4\pi \epsilon_0 \kappa} \dfrac{2\pi e^2}{q(1 + r_0 q/\kappa )}\right), \quad q = | \mathbf p- \mathbf p' |\end{aligned}\]where \(H_0\)and \(Y_0\) denote, respectively, the order-zero Struve function and the Bessel function of the second kind. Here \(r\) denotes the electron-hole separation, \(r_0\) is the screening length, and \(\kappa=(\epsilon_1+\epsilon_2)/2\) is the average environmental dielectric constant. The factor \(\epsilon_0/e^2=\big(4\pi\,\alpha\,\hbar c\big)^-1\) converts between the fine-structure constant and energy units. We use \(\alpha^-1=137.035999084\), \(\hbar c=1973.269804\ \text{eV}\!\cdot\!\text{Å}\), and the electron rest mass \(m_0=0.510998950\ \text{MeV}\). The material-dependent RK parameters used in our calculations are summarized in Table 4.

Our exciton binding energies for monolayer MoS\(_2\), MoSe\(_2\), WS\(_2\), and WSe\(_2\), computed with the RK momentum-space potential using the parameters in Table 4 and obtained by solving the LS integral equation are reported in Table 5. As a numerical self-consistency check, we also compare each binding energy with the Hamiltonian expectation value \(\langle H\rangle=\langle H_0\rangle+\langle V\rangle\) evaluated from the corresponding wave function. In all four materials, the relative percentage differences are well below \(3 \cdot 10^-6\%\), indicating that the momentum/angle discretization, kernel construction, and diagonalization are converged and the solutions are numerically reliable.

In Table 6, we show the convergence of exciton binding energies for the TMDs as a function of the number of mesh points in the relative momentum magnitude \(N_p\), together with an extrapolation to\(1/N_p \to 0\) (last column). The extrapolated values are in excellent agreement with Refs.  across all materials, with relative percentage differences below \(0.3\%\).