1 Introduction

Analytical results obtained for interacting fermions are of significant importance in various areas of physics, particularly in condensed matter physics, nuclear physics, and high-energy physics. Fermions represent a class of particles that obey Fermi–Dirac statistics and are subject to the Pauli exclusion principle. Understanding the behavior of interacting fermions is essential for a deeper understanding of many physical systems. Interacting fermions are often found in many-body quantum systems, in which the behavior of each particle is influenced by the presence of other particles, and these systems can be highly complex. Therefore, obtained analytical results can provide valuable insights into their behavior. These systems are often found in solid-state physics, where electrons in a crystal lattice interact with each other. The study of interacting fermions has played a pivotal role in understanding superconductivity and further such systems are central to quantum field theory, which is the theoretical framework that underlies the Standard Model of particle physics. While analytical results are valuable, it is important to note that, in many cases, analytical solutions are not possible due to the complexity of composite systems formed by interacting fermions. In such cases, some approximations, numerical simulations and computational methods play an important role in obtaining results. Analytical results obtained for interacting fermions are essential for understanding a wide range of physical phenomena, from the behavior of electrons in condensed matter systems to the structure of atomic nuclei and the fundamental forces of the universe. They provide the theoretical foundation for many areas of physics and often guide experimental research.

In non-relativistic quantum mechanics [1], determining the bound-states, scattering states, and resonance states of two-body systems typically involves solving the Schrödinger equation for the given potential. Bound states are characterized by having negative energy eigenvalues, which means the total energy of the system is less than the potential energy, and the two particles are bound together [1]. Resonance states are temporary states with a finite lifetime. They typically have complex energy eigenvalues with a negative real part and an imaginary part, representing the decay of the state. As resonance states have complex energies, the wave function will exhibit oscillatory behavior, representing the decaying and propagating character of the state. To distinguish between bound and resonance states, one must analyze the energy spectrum and the behavior of the associated wave function [1]. In relativistic quantum mechanics, particularly when dealing with interacting fermions, one typically works with the one-body Dirac equation in general. When dealing with a two-body system, one needs to consider the interaction between these two particles. The general structure of well-constructed two-body equations in relativistic quantum mechanics involves the extension of the Dirac equation to describe the fermion–fermion or fermion–antifermion systems. This is typically done using the concept of a multi-particle wave function. In such systems, the interaction terms are considered as the electromagnetic interaction, the strong nuclear force, or any other forces relevant to the specific system studied [1]. It is important to note that the form of the two-body equation may vary depending on the specific particles and interactions considered. In this framework, additionally, in practice, the two-body equations can become quite complex, and analytical solutions may not always be possible [1,2,3]. On the other hand, in general, phenomenologically established many-body equations are studied to explore the relativistic quantum motion of composite systems composed by interacting particles and the interaction potentials are chosen phenomenologically. That is one photon/boson exchange potential is considered in general. But, when constructing such two-body equations, time coordinates are given a priority. This essentially means ignoring the relative time difference. The first attempt to determine the relativistic dynamics of composite systems consisting of two interacting fermions using the relativistic two-body equation (accepted) took place shortly after the introduction of the Dirac equation. This equation is constructed by changing the Dirac matrices with the velocity of the particles and considering the modified Darwin potential as an inter-particle interaction potential [4]. But, this equation holds only in weak coupling regime since it does not hold when the interaction is long-range or velocities of the particles are high. To overcome this problem, a new formalism was introduced by starting from quantum field theory [5]. However, in this new formalism, the problem of relative time difference has emerged. So, this new equation can only be solved under certain approximations such as the instantaneous interaction. A two-body equation that is necessary to determine the dynamics of systems consisting of two interacting fermions must be a fully covariant one-time two-body equation that takes into account the kinematics of each particle, includes the most general electric and magnetic potentials, and includes accurate spin algebra. Such an equation was derived from quantum electrodynamics through the action principle [6], and then, it was shown that this equation can also be derived as an excited state of zitterbewegung [7]. This fully covariant two-body equation cannot be solved exactly for fermion–fermion or fermion–antifermion systems in \((3+1)\)-dimensions. This is because it gives 16-dimensional matrix equation resulting in a set of coupled equations consisting of two second-order differential equations even for the “well-known” systems such as one-electron atoms or quark–antiquark pairs [8, 9]. In recent years, this equation has been revisited, and it was shown that this one-time equation can be soluble for interacting fermions in \((2+1)\)-dimensions [10,11,12,13]. In this manuscript, we try to obtain results for a fermion–antifermion pair interacting via Coulomb-type inter-particle interaction potential by using the generators of the \(s\ell _{2}\) algebra, and then, we try to adapt the results to an exciton to determine real and damped oscillations for such an unstable system [14,15,16] in a monolayer medium in the vicinity of substrate [17].

An exciton in a monolayer material typically refers to an excited state of the material where an electron and a hole are bound together, and their energy levels allow them to interact with light, resulting in the emission or absorption of photons. Monolayer materials are two-dimensional materials that consist of a single layer of atoms or molecules, such as graphene or transition metal dichalcogenides [14, 15, 17]. In monolayer materials, excitons (bright) are those that efficiently couple with photons, which means they can easily emit or absorb light, making them observable in optical experiments. The energy levels of such excitons in monolayer materials depend on the band structure of the material. In transition metal dichalcogenides, for example, excitons are typically associated with the transition between the valence band and the conduction band. Hence, excitons play a crucial role in the optical properties of monolayer materials [18]. They are responsible for the features in the material’s absorption and photoluminescence spectra [14,15,16,17,18]. When a photon is absorbed, it can promote an electron from the valence band to the conduction band, creating an exciton. Conversely, when an exciton recombines, it emits a photon. Excitons in monolayer materials can have distinct spin and valley properties. The spin of an exciton refers to the orientation of the electron and hole spins, while the valley degree of freedom is related to the momentum of the electron and hole in the crystal’s momentum space. These properties can be manipulated for various applications, including valleytronics and spintronics [18]. Excitons in monolayer materials have led to the development of a range of optical devices, including light-emitting diodes and photodetectors. Their unique optical and electronic properties make them promising candidates for future technologies. Hence, determining the properties and behavior of excitons in monolayer materials is of great interest for both fundamental research and potential technological applications [19]. In this paper, we will try to explore the evolution of a fermion–antifermion pair interacting through Coulomb-type inter-particle interaction potential in \((2+1)\)-dimensions by searching for analytical solutions of the fully covariant two-body Dirac equation (see also [20]).

This paper is arranged as the following: In Sect. 2, we introduce the covariant two-body Dirac equation in \((2+1)\)-dimensional globally and locally flat spacetime background and derive a matrix equation for a static and spinless composite system formed by a fermion–antifermion pair bound together through an attractive Coulomb-type inter-particle interaction potential. In Sect. 3, we show non-perturbative second-order wave equation and determine its analytical solution for such a pair by using the generators of \(s\ell _{2}\) algebra. Then, we analyze the results in detail. In Sect. 4, we give a summary and discuss the effects of the surrounding medium on the real and damped modes of excitons.

2 Matrix equation

Let us start this section by introducing the fully covariant two-body Dirac equation. Then, we will derive a matrix equation for the relative motion of an interacting fermion–antifermion pair through an attractive Coulomb-type inter-particle interaction potential. In \((2+1)\)-dimensional flat spacetime background that can be represented through the metric: \({\text {d}}s^{2}=c^2{\text {d}}t^2-{\text {d}}x^2-{\text {d}}y^2\), the fully covariant two-body Dirac equation can be written in the following form [10,11,12,13]

$$\begin{aligned}&\left\{ \mathcal {H}_{f} \otimes \gamma ^{t^{\overline{f}}}+\gamma ^{t^{f}}\otimes \mathcal {H}_{\overline{f}} \right\} \Theta \left( {\textbf {x}}_{f},{\textbf {x}}_{\overline{f}}\right) =0,\nonumber \\&\mathcal {H}_{f}=\left[ \gamma ^{\mu ^{f}}\mathcal {D}_{\mu {f}}+i\mu _{f}{} {\textbf {I}}_2 \right] , \mathcal {H}_{\overline{f}}=\left[ \gamma ^{\mu ^{\overline{f}}}\mathcal {D}_{\mu _{\overline{f}}}+i\mu _{\overline{f}}{} {\textbf {I}}_2 \right] ,\nonumber \\&\mathcal {D}_{f}=\partial _{\mu }^{f}+i\frac{e^{f}\mathcal {A}_{\mu }^{\overline{f}}}{\hbar c}-\Gamma _{\mu }^{f},\quad \mathcal {D}_{\overline{f}}=\partial _{\mu }^{\overline{f}}+i\frac{e^{\overline{f}}\mathcal {A}_{\mu }^{f}}{\hbar c}-\Gamma _{\mu }^{\overline{f}},\nonumber \\&\mu _{f}=\frac{m_{f}c}{\hbar },\quad \mu _{\overline{f}}=\frac{m_{\overline{f}} c}{\hbar }. \end{aligned}$$
(2.1)

Here, the Greek indices indicate the coordinates of the background spacetime, the f (\(\overline{f}\)) refers to fermion (antifermion), \(\gamma ^{\mu }\) are the generalized Dirac matrices in \((2+1)\)-dimensions, \({\textbf {I}}_2\) denotes identity matrix, \(\mathcal {A}_{\mu }\) is the electromagnetic 3-vector potential, the letter e is the electric charge, \(\Gamma _{\mu }\) is the spinorial connection, m represents the rest mass, c is the light speed, \(\hbar\) is usual Planck constant, \(\Theta\) is the bi-spinor field depending on the spacetime position vectors (\({\textbf {x}}_{f}, {\textbf {x}}_{\overline{f}}\)) of the particles, and the \(\otimes\) symbol means the Kronecker product. The considered background spacetime is locally and globally flat, and hence, the spinorial affine connections do not change the dynamics of the system since they vanish [10]. The space-dependent Dirac matrices can be determined by using the following relation: \(\gamma ^{\mu }=e^{\mu }_{(k)}\overline{\gamma }^{(k)}\) in which \(e^{\mu }_{(k)}\) indicates the inverse tetrad fields and \(\overline{\gamma }^{(k)}\) are space-independent Dirac matrices. The \(\overline{\gamma }^{(k)}\) matrices are chosen in terms of the well-known Pauli spin matrices (\(\sigma _x, \sigma _y, \sigma _z\)) in \((2+1)\)-dimensions with respect to the signature (\(+,-,-\)) of the given metric, as \(\overline{\gamma }^{0}=\sigma _z\), \(\overline{\gamma }^{1}=i\sigma _x\) and \(\overline{\gamma }^{2}=i\sigma _y\) where \(i=\sqrt{-1}\) [10,11,12,13]. The tetrads can be obtained through \(g_{\mu \nu }=\)diag\((c^2,-1,-1)=e_{\mu }^{(k)}e_{\nu }^{(l)}\eta _{(k)(l)}\) in which \(g_{\mu \nu }\) is the contravariant metric tensor, \(e_{\mu }^{(k)}\) stand for the tetrads besides the flat Minkowski metric tensor \(\eta _{(k)(l)}=\)diag\((1,-1,-1)\). Accordingly, one can choose the tetrads as \(e^{(0)}_{0}=c\), \(e^{(1)}_{1}= 1\), \(e^{(2)}_{2}=1\) and hence one can obtain the inverse tetrads as the following \(e^{t}_{(0)}=1/c\), \(e^{x}_{(1)}=1\) and \(e^{y}_{(2)}=1\) through the orthogonality and orthonormality conditions: \(e^{\mu }_{(k)}e^{(k)}_{\nu }=\delta ^{\mu }_{\nu }\) and \(e^{(k)}_{\mu }e^{\mu }_{(l)}=\delta ^{(k)}_{(l)}\) where \(\mu ,\nu =t,x,y\) and \(k,l=0,1,2\). Here, we are interested in a fermion–antifermion pair interacting through an attractive Coulomb-type inter-particle interaction potential, namely \(\mathcal {A}_{t}=v({\textbf {x}}_{f}-{\textbf {x}}_{\overline{f}})\). Hence, we get that \(\mathcal {A}_{x}^{f,\overline{f}}=0\) and \(\mathcal {A}_{y}^{f,\overline{f}}=0\). Here, we try to determine real and damped oscillations for such a pair. To acquire this, we need to introduce the center of mass (\(\mathcal {R}\)) and relative motion coordinates (r) to arrive at a matrix equation for the relative motion of the considered fermion–antifermion pair for which \(e_{f}=-e_{\overline{f}}=e\) and \(m_{f}=m_{\overline{f}}=m\). As is usual with two-body problems, one can use the following equations [10,11,12,13] to derive a matrix equation for the relative motion of the considered pair

$$\begin{aligned}&r_{\mu }=x_{\mu }^{f}-x_{\mu }^{\overline{f}},\quad \mathcal {R}_{\mu }=\frac{x_{\mu }^{f}}{2}+\frac{x_{\mu }^{\overline{f}}}{2}, \quad x_{\mu }^{f}=\frac{r_{\mu }}{2}+\mathcal {R}_{\mu },\nonumber \\&x_{\mu }^{\overline{f}}=-\frac{r_{\mu }}{2}+\mathcal {R}_{\mu },\quad \partial _{x_{\mu }}^{f}=\partial _{r_{\mu }}+\frac{\partial _{\mathcal {R}_{\mu }}}{2},\nonumber \\&\partial _{x_{\mu }}^{\overline{f}}=-\partial _{r_{\mu }}+\frac{\partial _{\mathcal {R}_{\mu }}}{2}, \end{aligned}$$

which, of course, leads to \(\partial _{x_{\mu }}^{f}+\partial _{x_{\mu }}^{\overline{f}}=\partial _{\mathcal {R}_{\mu }}\). Through the obtained inverse tetrads, one finds the space-dependent Dirac matrices in Eq. (2.1) as the following \(\gamma ^{t^{f,\overline{f}}}=\sigma _z/c\), \(\gamma ^{x^{f,\overline{f}}}=i\sigma _x\) and \(\gamma ^{y^{f,\overline{f}}}=i\sigma _y\) even though the temporal parts depend on the speed of light c (constant). To observe any pairing effect, we can consider a static composite system whose center of mass is at rest at the center of the spacetime background. Here, we also consider that the interaction potential is time-independent. This fact allows us to practically factorize the resulting spinor as followsFootnote 1\(\Theta =\text {e}^{-i\omega t}\tilde{\Theta }({\textbf {r}})\). Here, \(\omega\) is the relativistic frequency and \(\tilde{\Theta }({\textbf {r}})\) is the spatial part of the resulting spinor field. At that rate, for the system under analysis, one can derive a matrix equation in the form of \(\hat{M}\tilde{\Theta }=0\) in which

$$\begin{aligned}\hat{M}=&\left( \begin{array}{cccc} \zeta (r)-\tilde{\mu }&{} \hat{\mathcal {D}}_{-} &{} -\hat{\mathcal {D}}_{-} &{} 0\\ -\hat{\mathcal {D}}_{+} &{} \zeta (r) &{} 0 &{} -\hat{\mathcal {D}}_{-}\\ \hat{\mathcal {D}}_{+} &{} 0 &{} \zeta (r) &{} \hat{\mathcal {D}}_{-}\\ 0 &{} \hat{\mathcal {D}}_{+} &{} -\hat{\mathcal {D}}_{+} &{} \zeta (r)+\tilde{\mu } \end{array} \right) ,\\ \tilde{\Theta }=&\left( \psi _{1}({\textbf {r}}) \ \psi _{2}({\textbf {r}}) \ \psi _{3}({\textbf {r}}) \ \psi _{4}({\textbf {r}})\right) ^{{\textbf {T}}}, \hat{\mathcal {D}}_{\pm }=\partial _{x}\pm i\partial _{y}, \end{aligned}$$
(2.2)

where \(\zeta (r)=\omega /c-v(r)\), \(\tilde{\mu }=2mc/\hbar\) and \({\textbf {T}}\) means transpose of the resulting spinor. Here, it is clear that the spinor components depend explicitly on the relative motion coordinates, as \(\psi _{j}(x,y), (j=1,2,3,4.)\). To exploit angular symmetry, we can consider polar coordinates by using the spin raising (+) and spin lowering (−) operators given by \(\hat{\mathcal {D}}_{\pm }=\text {e}^{\pm i\phi }\left( \pm \frac{i}{r}\partial _{\phi }+\partial _{r}\right)\) only for the transformed spinor [10]

$$\begin{aligned} \tilde{\Theta }=\left( \psi _{1}(r)\text {e}^{i(s-1)\phi } \ \psi _{2}(r)\text {e}^{is\phi } \ \psi _{3}(r)\text {e}^{is\phi } \ \psi _{4}(r)\text {e}^{i(s+1)\phi }\right) ^{\textbf {T}}. \end{aligned}$$

3 Real and damped oscillations

Fig. 1
figure 1

a Illustrates the impact of the effective dielectric constant on the exciton binding energy \(\mathcal {E}^{b}_{n}(\varepsilon )\) [eV]. b Depicts the variation in exciton decay time \(\tau _{n}(\varepsilon )\) [s] with respect to the effective dielectric constant (\(\varepsilon\)) of the surrounding medium

Full size image

In this section, we derive first a set of coupled equations, and then, we show a non-perturbative second-order wave equation for a static and spinless (\(s=0\)) system formed by a fermion–antifermion pair interacting through a Coulomb type inter-particle interaction potential v(r) [17]. Then, we use quasi-exactly solvable methods and adapt the results to an electron–hole pair so that we can observe the evolution of an exciton [17] in monolayer medium in the vicinity of the substrate. For the considered system, the matrix Eq. (2.2) leads to a set of coupled equations consisting of four equations. After some algebra, this set of equations can be clarified as the following

$$\begin{aligned}&\zeta (r)\psi _{+}-\tilde{\mu }\psi _{-}+4\psi _{0_{,r}}=0,\nonumber \\&\zeta (r)\psi _{-}-\tilde{\mu }\psi _{+}=0,\nonumber \\&\zeta (r)\psi _{0}+\psi _{+_{,r}}=0, \end{aligned}$$

where \(_{,r}\) means derivative with respect to r, \(\psi _{\pm }(r)=\psi _{1}\pm \psi _{4}\) and \(\psi _{0}=\psi _{2}=-\psi _{3}\). These equations can be solved in favor of \(\psi _{+}\). Accordingly, one can derive the following non-perturbative second-order wave equation

$$\begin{aligned} \psi _{+_{,rr}}-\frac{\zeta _{,r}}{\zeta }\psi _{+_{,r}}+\frac{\zeta ^2-\tilde{\mu }^2}{4}\psi _{+}=0. \end{aligned}$$
(3.1)

Here, we consider a fermion–antifermion pair interacting via Coulomb-type inter-particle potential in \((2+1)\)-dimensions. By choosing the interaction term as \(v(r)=-\frac{\alpha }{\varepsilon r}\), in which \(\alpha\) is the well-known fine structure constant (\(\sim 1/137\)) and \(\varepsilon\) is the effective dielectric constant of the surrounding medium [17], we will try to explore the real and damped modes of an exciton in monolayer medium in the vicinity of substrate. For such a system, through a new change of variable, \(\chi =\frac{\varepsilon \omega }{\alpha c}\ r\),Footnote 2 one obtains a wave equation which can be clarified, by considering an ansatz function \(\psi _{+}(\chi )=\chi ^{\frac{\alpha }{2\epsilon }}\text {e}^{\frac{\tilde{\alpha }}{2}\chi }\ \psi (\chi )\), as the following

$$\begin{aligned} \psi _{,\chi \chi }+\left[ \tilde{\alpha }+\frac{\beta +1}{\chi }+\frac{\gamma +1}{\chi -1}\right] \psi _{,\chi }+\left[ \frac{\bar{\mu }}{\chi }+\frac{\bar{\nu }}{\chi -1}\right] \psi =0. \end{aligned}$$
(3.2)

This wave equation is the simplest uniform shape of the confluent Heun differential equation [13, 21] and its regular solution around the regular singular point \(\chi =0\) is given as \(\psi (\chi )=\mathcal {N}\mathcal {H}_{\mathcal {C}}[\tilde{\alpha }, \beta , \gamma , \delta ,\tilde{\eta }; z]\) [13], in which \(\mathcal {N}\) is an arbitrary constant and

$$\begin{aligned} \tilde{\alpha }= & {} -\frac{\alpha }{\varepsilon \omega } \sqrt{\tilde{\mu }^2c^2-\omega ^2},\quad \beta =i\frac{\alpha }{\varepsilon },\quad \gamma =-2,\nonumber \\ \delta= & {} -\frac{\alpha ^2}{2\varepsilon ^2},\quad \tilde{\eta }=1+\frac{\alpha ^2}{2\varepsilon ^2},\nonumber \\ \bar{\mu }+\bar{\nu }=\, & {} \delta +\tilde{\alpha }\left[ \frac{\beta }{2}+\frac{\gamma }{2}+1\right] ,\nonumber \\ \tilde{\eta }= & {} \frac{1}{2}\left( \tilde{\alpha }+\tilde{\alpha }\beta -\beta -\beta \gamma -\gamma -\bar{\mu }\right) . \end{aligned}$$
(3.3)

According to quasi-exact solution method [22], one implements that

$$\begin{aligned} \mathcal {C}_{+}=\, & {} \tilde{\alpha },\quad -n\mathcal {C}_{+}=\bar{\mu }+\bar{\nu },\quad \mathcal {C}_{+-}=1,\quad \mathcal {C}_{0-}=-1,\nonumber \\ \mathcal {C}_{0}=\, & {} n-\tilde{\alpha }+2+\gamma +\beta ,\quad \mathcal {C}_{-}=-1-\frac{n}{2}-\beta ,\nonumber \\ \mathcal {C}= & {} -\delta +\frac{1}{2}[n-\tilde{\alpha }+2+\gamma +\beta ]n. \end{aligned}$$
(3.4)

According to Eqs. (3.3) and (3.4), the Lie-algebraic differential operator \(s\ell _{2}\) can be written as the following [22]

$$\begin{aligned} \mathcal {H}= &\, {} \mathcal {J}^{+}\mathcal {J}^{-}-\mathcal {J}^{0}\mathcal {J}^{-}+\tilde{\alpha }\mathcal {J}^{+}-\left( 1+\beta +\frac{n}{2} \right) \mathcal {J}^{-} \\{} & \, +\left( n+2-\tilde{\alpha }+\beta +\gamma \right) \mathcal {J}^{0}-\delta \\{} &\, + \frac{n}{2}\left( n+2-\tilde{\alpha }+\beta +\gamma \right) , \end{aligned}$$
(3.5)

where \(\mathcal {J}^{+}\), \(\mathcal {J}^{-}\) and \(\mathcal {J}^{0}\) are obtained through \(s\ell _{2}\) Lie algebra, by the generators [22]

$$\begin{aligned} \mathcal {J}^{+}=z^{2}\partial _{z}-n z,\,\, \mathcal {J}^{0}=z\partial _{z}-\frac{n}{2},\,\,\mathcal {J}^{-}=\partial _{z}. \end{aligned}$$
(3.6)

Hence, the spinor component \(\psi _{+}(\chi )\) can be determined as follows

$$\begin{aligned} \psi _{+}(\chi )=\chi ^{\frac{\alpha }{2\epsilon }}\text {e}^{\frac{\tilde{\alpha }}{2}\chi } \sum _{p=0}^{n}a_{p} \chi ^{p}, \end{aligned}$$
(3.7)

and three-term recurrence relation for the coefficients in Eq. (3.7) becomes

$$\begin{aligned} a_{p+1}= \,& {} \left[ \frac{p(p+1-\tilde{\alpha }+\beta +\gamma )-\delta }{(p+1) (p+1+\beta ) } \right] a_{p} -\left[ \frac{2 p \tilde{\alpha }}{(p+1 ) (p+1+\beta )} \right] a_{p-1}, \end{aligned}$$

if \(a_{-1}=0\) [13]. Via Eq. (3.3) and Eq. (3.4) one obtains the following expression for the relativistic frequency of the system under analysis

$$\begin{aligned} \omega _{n}= \frac{2mc^2}{\hbar }\sqrt{\frac{4\varepsilon ^2\left( n-\frac{i\alpha }{2\varepsilon }\right) ^2}{\alpha ^2+4\varepsilon ^2\left( n-\frac{i\alpha }{2\varepsilon }\right) ^2}}, \end{aligned}$$

where n is the radial quantum number (\(n=1,2,3,\ldots\)) [22]. Through power series expansion method, as is usual with a typical Coulomb problem, one can obtain an energy expression (\(\hbar \omega _{n}\)) in the form of

$$\begin{aligned}&\mathcal {E}_{n}=\mathcal {E}_{\mathcal {R}_n}+\mathcal {E}_{\mathcal {I}_n}, \end{aligned}$$
(3.8)

where

$$\begin{aligned} \mathcal {E}_{\mathcal {R}_n}\approx 2mc^2\left\{ 1-\frac{\alpha ^2/\varepsilon ^2}{8n^2}+\frac{15\alpha ^4/\varepsilon ^4}{128n^4}-\frac{105\alpha ^6/\varepsilon ^6}{1028n^6} \right\} , \end{aligned}$$

and

$$\begin{aligned} \mathcal {E}_{\mathcal {I}_n}\approx -im c^2\frac{\alpha ^3/\varepsilon ^3}{4n^3}+imc^2\frac{7\alpha ^5/\varepsilon ^5}{32n^5}. \end{aligned}$$

Accordingly, we can obtain the binding energy (\(\mathcal {E}^{b}_{n}\))

$$\begin{aligned} \mathcal {E}^{b}_{n}(\varepsilon )=-\frac{mc^2\alpha ^2/\varepsilon ^2}{4n^2}\left\{ 1-\frac{15\alpha ^2/\varepsilon ^2}{16n^2}+\frac{105\alpha ^4/\varepsilon ^4}{128n^4}\right\} , \end{aligned}$$
(3.9)

and decay time (\(\tau _{n}=\hbar /|\mathcal {E}_{\mathcal {I}_n}|\)) of the damped modes as the following

$$\begin{aligned} \tau _{n}(\varepsilon )=\frac{4\hbar n^3}{m c^2\alpha ^3/\varepsilon ^3\left\{ 1-\frac{7\alpha ^2/\varepsilon ^2}{8n^2}\right\} }, \end{aligned}$$
(3.10)

since \(\Theta \propto \text {e}^{-i\omega t}\). It can be seen that Eq. (3.10) gives decay time in the order of \(\sim 10^{-12}\) [s] for the ground state of an exciton in monolayer medium in the vicinity of substrate (see also [14,15,16,17,18,19]) if m is the usual electron mass (see also [20]). Here, we determine the discrete energy levels and decay time of the damped modes as a function of the effective dielectric constant of the surrounding medium. Here, we observe that the substrate material causes enhanced screening effects. That is, the substrate material can be responsible for lower binding energy and longer lifetime for excitonic states. Moreover, one can observe that the total lifetime of an exciton is dominated by the excited states because the decay time of the excited states can be quite longer than the ground state, especially for large values of the \(\varepsilon\). In the research [23, 24], exciton binding energies were measured differently for monolayer \(WS_{2}\) in the vicinity of different substrates (see also [25]). Our model clearly shows that the discrepancy between the previously announced results [23, 24] is due to the preferred substrates (Fig. 1).

4 Summary and discussions

Here, we have studied the quantum electrodynamics of a many-body system formed by an interacting fermion–antifermion pair in \((2+1)\)-dimensional spacetime background. By considering the inter-particle interaction potential between the particles as an attractive Coulomb-type potential, we search for solutions of a fully covariant two-body Dirac equation derived from quantum electrodynamics via the action principle (see also [20]). This equation has provided a non-perturbative second-order wave equation for such a static and spinless composite system, and we obtained its regular solution through the generators of \(s\ell _{2}\) algebra. This solution has allowed us to derive the quantization condition for the formation of such a pair, and accordingly, we determine the relativistic frequency modes (complex). Then, we have adapted the findings to an electron–hole pair in a monolayer medium in the vicinity of the substrate and show how an effective dielectric environment (\(\propto \varepsilon\)) modifies the real and damped modes of such an exciton. We have observed that the ground state binding energy of the exciton in question can vary from \(\sim 1\) eV to the order of a few meV with respect to varying values of the \(\varepsilon\) especially for \(\varepsilon >2\) (see also [14,15,16,17,18,19, 23,24,25]). Moreover, our model has shown that these states cannot be steady states and decay over time (exponentially). This provides a good ground to determine the decay time for each quantum state. We have observed that the decay time for the ground state of such an exciton is on the order of \(\sim 10^{-12}\) s. However, we have also shown that the decay process of an exciton is dominated by its excited states because the decay time of the excited states can be considerably longer than that of the ground state, especially for increasing \(\varepsilon\) values. Our results can help to explain the difference between the previously announced observations [23, 24]. Moreover, our findings impose that evolution of such an exciton in a monolayer medium can be actively controlled by arranging the effective dielectric constant (\(\varepsilon\)) of the surrounding medium.

Controlling exciton evolution in monolayer materials may lead to various benefits. By manipulating exciton properties, such as their binding energies and lifetimes, it may be possible to increase the efficiency of light absorption and emission in monolayer materials. This can lead to more efficient photodetectors, light-emitting devices, and photovoltaic cells. Also, monolayer materials are promising platforms for quantum technologies like quantum dots and quantum information processing. Precise control over exciton evolution is essential for the creation and manipulation of quantum states, making these materials potential candidates for quantum computing and communication. The ability to control exciton properties allows for the creation of tunable optoelectronic devices. This means that the optical and electronic characteristics of the materials can be adjusted to meet specific requirements.