Rigorous asymptotic study of the screened electrostatic potential in a thin dielectric slab

TL;DR

This paper rigorously analyzes the asymptotic behavior of the screened electrostatic potential in ultra-thin dielectric slabs, deriving exact effective 2D potentials to aid exciton energy calculations in 2D materials.

Contribution

It provides a rigorous asymptotic analysis and explicit expressions for the screened Coulomb potential in thin dielectric slabs, advancing theoretical understanding for 2D material exciton studies.

Findings

Derived regularized expressions for the potential.

Obtained exact effective 2D potential formulas.

Facilitates accurate exciton binding energy computations.

Abstract

The screened Coulomb potential plays a crucial role in the binding energies of excitons in a thin dielectric slab. The asymptotic behavior of this potential is studied when the thickness of the slab is very small as compared to the exciton Bohr radius. A regularized expression is given and the exact effective 2D potential is derived. These expressions may be useful for the computation of the exciton binding energy in 2D or quasi-2D materials.

Figures

Equations

$${\eta }_1=\frac{1}{2}\log \left(\frac{{\epsilon }_f+{\epsilon }_1}{{\epsilon }_f-{\epsilon }_1}\right),{\eta }_2=\frac{1}{2}\log \left(\frac{{\epsilon }_f+{\epsilon }_2}{{\epsilon }_f-{\epsilon }_2}\right),$$

$$Y_{\eta }=[-\eta /2,\eta /2{]}^{2}andY=[-1/2,1/2{]}^2.$$

$$V(\rho ,z_0,z_0^{\prime })=\frac{e{e}^{\prime }}{2\pi {\epsilon }_{f}d}I(\frac{\rho }{d},\frac{z_0}{d},\frac{z_0^{\prime }}{d}),$$

$$I(r,x,y)={\int }_0^{+\infty }W(u,x,y)J_0(ru)du,$$

$$W(u,x,y)=\frac{\cosh \left[u(\frac{1}{2}+y)+{\eta }_1\right]\cosh \left[u(\frac{1}{2}-x)+{\eta }_2\right]}{\sinh \left(u+{\eta }_1+{\eta }_2\right)},(x,y)\in Y.$$

$$-\frac{{\hslash }^2}{2m_1}{\Delta }_0\psi -\frac{{\hslash }^2}{2m_2}{\Delta }_0^{\prime }\psi +V(|{\mathbf{\rho }}_0-{\mathbf{\rho }}_0^{\prime }|,z_0,z_0^{\prime })\psi =E\psi$$

$$\left\{-{\Delta }_{||}-{\kappa }_1{\partial }_{z_0}^2-{\kappa }_2{\partial }_{z_0^{\prime }}^2+\frac{2\mu }{{\hslash }^2}V(\rho ,z_0,z_0^{\prime })\right\}\psi =\frac{2\mu }{{\hslash }^2}E\psi$$

$$\left\{-{\Delta }_{||}-{\kappa }_1{\partial }_{\tilde{z}}^2-{\kappa }_2{\partial }_{{\tilde{z}}^{\prime }}^2-\frac{2\mu d^2}{{\hslash }^2}V(r,z,z^{\prime })\right\}=\frac{2\mu d^2}{{\hslash }^2}E\psi .$$

$$a_0=4\pi {\hslash }^2{ϵ}_f/\mu e^2,$$

$$\left\{-{\Delta }_{||}-{\eta }^2{\kappa }_1{\partial }_{\zeta }^2-{\eta }^2{\kappa }_2{\partial }_{{\zeta }^{\prime }}^2-4\eta I(r,\zeta /\eta ,{\zeta }^{\prime }/\eta )\right\}\psi =E_{\eta }\psi ,$$

$$W(u,x,y)\sim \frac{1}{2}{e}^{-u|x-y|}asu\to +\infty .$$

$${\int }_0^{+\infty }{e}^{-\alpha x}{J}_0(\beta x)dx=\frac{{\beta }^{-\nu }[\sqrt{{\alpha }^2+{\beta }^2}-\alpha {]}^{\nu }}{\sqrt{{\alpha }^2+{\beta }^2}}.$$

$$I(r,x,y)={\int }_0^{+\infty }\left(W(u,x,y)-\frac{1}{2}{e}^{-u|x-y|}\right)J_0(ru)du+\frac{1}{2\sqrt{r^2+|x-y{|}^2}}.$$

$$I(r,x,y)={\int }_0^{+\infty }\left(W(u,x,y)-\frac{1}{2}{e}^{-u|x-y|}\right)J_0(ru)du+\frac{1}{2}{\int }_0^{+\infty }{J}_0(ru)e^{-u|x-y|}du$$

$${\int }_0^{+\infty }{e}^{-u|x-y|}{J}_0(ru)du=\frac{1}{\sqrt{r^2+|x-y{|}^2}}$$

$$-{\Delta }_{||}-4\eta I_0(r),$$

$$I_0(r)={\int }_{Y}I(r,x,y)dxdy.$$

$$\left\{-{\Delta }_{||}-{\eta }^2\mu /m_1{\partial }_{\zeta }^2-{\eta }^2\mu /m_2{\partial }_{{\zeta }^{\prime }}^2+\eta V(r,\zeta /\eta ,{\zeta }^{\prime }/\eta )\right\}{\psi }_{\eta }=E_{\eta }{\psi }_{\eta }.$$

$${\int }_{Y_{\eta }\times \mathbb{R}}{\nabla }_{||}{\psi }_{\eta }\nabla \varphi +{\eta }^2{\kappa }_1{\int }_{Y_{\eta }\times \mathbb{R}}{\partial }_{\zeta }{\psi }_{\eta }{\partial }_{\zeta }\varphi +$$

$${\eta }^2{\kappa }_2{\int }_{Y_{\eta }\times \mathbb{R}}{\partial }_{{\zeta }^{\prime }}{\psi }_{\eta }{\partial }_{\zeta }\varphi +\eta {\int }_{Y_{\eta }\times \mathbb{R}}V(r,\zeta /\eta ,{\zeta }^{\prime }/\eta )\psi \varphi =E_{\eta }{\int }_{Y_{\eta }\times \mathbb{R}}{\psi }_{\eta }\varphi .$$

$${\int }_{Y_{\eta }}{\partial }_{\zeta }{\psi }_{\eta }{\partial }_{\zeta }\varphi =O({\eta }^2),{\int }_{Y_{\eta }}{\partial }_{{\zeta }^{\prime }}{\psi }_{\eta }{\partial }_{{\zeta }^{\prime }}\varphi =O({\eta }^2),$$

$$\frac{1}{{\eta }^2}{\int }_{Y_{\eta }\times \mathbb{R}}V(r,\zeta /\eta ,{\zeta }^{\prime }/\eta ){\psi }_{\eta }\varphi d\zeta d{\zeta }^{\prime }=\int dr\left[{\int }_{Y}V(r,x,y)dxdy\right]{\psi }_{\eta }^0(r){\varphi }^0(r)+o(1),$$

$${\psi }_{\eta }^0(r)={\psi }_{\eta }(r,0,0),{\varphi }^0(r)=\varphi (r,0,0).$$

$$\int {\nabla }_{||}{\psi }_{\eta }^0{\nabla }_{||}{\varphi }^0+\eta \int dr\left[{\int }_{Y}V(r,x,y)dxdy\right]{\psi }_{\eta }^0(r){\varphi }^0(r)=E_{\eta }\int {\psi }_{\eta }^0{\varphi }^0.$$

$$V_{2D}(r)={\int }_Y\frac{1}{\sqrt{r^2+|z^{\prime }-z{|}^2}}dzd{z}^{\prime }$$

$$V_A=\frac{\lambda }{r},$$

$$V_K(\rho )=\frac{e{e}^{\prime }}{2\pi ϵd}{I}_K(\rho ),$$

$$I_K(\rho )={\int }_0^{+\infty }\frac{J_0(t)dt}{t+\frac{{ϵ}_1+{ϵ}_2}{{ϵ}_f}r}.$$

$$-\nabla \cdot (\epsilon (z)\nabla \phi )=\delta (z-z_0^{\prime })\otimes \delta (\mathbf{\rho })$$

$$\hat{\phi }(\mathbf{k},z;z_0^{\prime })=\int \phi (\mathbf{\rho },z;z_0^{\prime })e^{-i\mathbf{k}\cdot \mathbf{\rho }}{d}^2\mathbf{\rho },$$

$$\phi (\mathbf{\rho },z;z_0^{\prime })=\frac{1}{(2\pi {)}^2}\int \hat{\phi }(\mathbf{k},z;z_0^{\prime })e^{i\mathbf{k}\cdot \mathbf{\rho }}{d}^2\mathbf{k}$$

Full text

1. Introduction

The raise of 2D materials has generated a renewed interest in computing excitons binding energy in thin structures, such as a slab of dielectric material. In that situation, it has been well appreciated in the literature that the Coulomb potential is screened and that it plays a crucial role on the value of the binding energies [2, 3, 4, 5]. Indeed, within the effective mass approximation, the enveloppe of the exciton wave function satisfies a Schr dinger equation involving the screened Coulomb potential. Several approaches have been used to derive an expression for this potential. One of the oldest work on this is an article by Keldysh [1]. In order to apply this potential to the situation of 2D materials, several approaching have been put forward in order to obtain a strictly 2D potential. In [1], the space variables transverse to the material are abruptly put to 0. In [4, 5] the starting point is a 2D polarisability, whereas in [3] the charge distributions are described as lines of charge, in order to obtain an effective 2D potential. Therefore, it seems that a direct limit analysis of the 3D potential when the width of the slab is very small is lacking.

In this note, we propose a multiple scale approach to the study of the screened potential when the thickness of the slab is very small with respect to the Bohr radius. We obtain a regularized potential by exhibiting explicitly the singular part, i.e. the bare Coulomb potential. The 2D effective potential compares very well with the approximation given in [1].

2. Expression of the screened potential and regularization

2.1. The screened potential

Consider the geometry described in fig.1. Cylindrical coordinates $(\boldsymbol{\rho},z)$ are used. We consider the electrostatic interaction between two charges $e$ located respectively at $(\boldsymbol{\rho}_{0},z_{0})$ and $(0,z^{\prime}_{0})$. We assume that $z_{0}\geq z_{0}^{\prime}$. The particles are situated in a dielectric slab of permittivity $\varepsilon_{f}$ and width $d$ situated in the interval $z\in[-d/2,d/2]$. The slab is surrounded by two semi-infinite medium of permittivities $\varepsilon_{1}$ for $z<-d/2$ and $\varepsilon_{2}$ for $z>d/2$. We denote

[TABLE]

and

[TABLE]

The electrostatic potential energy between the particles is defined by: $V(\boldsymbol{\rho},z_{0},z^{\prime}_{0})=e^{\prime}\hat{\varphi}_{e}(0,z^{\prime}_{0})=e\hat{\varphi}_{e^{\prime}}(\boldsymbol{\rho},z_{0})$, where $\hat{\varphi}_{q}$ is the potentiel created by particle $q$.

In the appendix, we give a full derivation of the potential since in [1] the result is stated without details. After all calculations are performed (see the appendix) the final result is

[TABLE]

where

[TABLE]

and the kernel $W$ is given by

[TABLE]

Since the expression was obtained for $z_{0}>z^{\prime}_{0}$ (see appendix), these variables should be switched when $z^{\prime}_{0}>z_{0}$.

From now on, we assume $e=-e^{\prime}$, that is, the second particle is a hole.

2.2. Solving the Schrödinger equation-Defining new coordinates

In a generic case, the situation considered in fig.1 is a two-body problem. The particles have masses $m_{1}$ and $m_{2}$ and are parametrized in cylindrical coordinates by $(\boldsymbol{\rho}_{0},z_{0})$ and $(\boldsymbol{\rho}^{\prime}_{0},z^{\prime}_{0})$ . We look for a time-independent wave function $\psi(\boldsymbol{\rho}_{0},z,\boldsymbol{\rho}^{\prime}_{0},z^{\prime})$ satisfying:

[TABLE]

As usual in this kind of problem, we focus on the relative motion. Here however, we cannot do this for the height variables $z$ since the potential does depend separately on $z_{0}$ and $z^{\prime}_{0}$. Therefore, we denote $\boldsymbol{\rho}=\boldsymbol{\rho}_{0}-\boldsymbol{\rho}^{\prime}_{0}$ and we obtain after eliminating the movement of the center of gravity in the $xOy$ plane:

[TABLE]

where $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ is the reduced mass, $\kappa_{1}=\mu/m_{1},\,\kappa_{2}=\mu/m_{2}$ and $\rho=|\boldsymbol{\rho}|$.

Let us now normalize the variables relatively to the slab width $d$. We define the new variables $z_{0}=\tilde{z}d,\,z^{\prime}_{0}=\tilde{z}^{\prime}d$. These variables belong to the interval $[-1/2,1/2]$. We also normalize the in-plane variable and define: $\rho=rd$. Plugging these into the equation leads to

[TABLE]

These basic transformations allow to introduce the exciton Bohr radius

[TABLE]

quite naturally.

From now on, we denote $\eta=d/a_{0}$ which is the small parameter of our problem. The final variables are now $\zeta$ and $\zeta^{\prime}$, satisfying $\zeta=\tilde{z}\eta,\,\zeta^{\prime}=\tilde{z}^{\prime}\eta$. These new variables belong to $[-\eta/2,\eta/2]$.

The final spectral equation with a small parameter is

[TABLE]

where: $E_{\eta}=(2\mu d^{2}/\hbar^{2})E$.

2.3. Regularization of the potential

It is clear that when the slab has an infinite width, one should recover the usual Coulomb potential $V_{C}(r)=1/4\pi\varepsilon_{f}r$. Therefore we expect that $V$ be a perturbation of $V_{C}$. In this section, we exhibit the singular part of $V$. First, we establish two lemma.

Lemma 1**.**

The kernel of screened Coulomb potential has the following asymptotic behavior

[TABLE]

Proof.

This is quite elementary using the expression of $\sinh$ and $\cosh$ in terms of the exponential. ∎

In order to evaluate the singular part, we need the following result

Lemma 2**.**

The following result holds for $\Re\alpha>|\Im\beta|$

[TABLE]

Proof.

See [7, p.702] formula (6.23 (3)). ∎

We are now in a position to exhibit a regularized potential, that is, which is not singular at the origin

Theorem 1**.**

The following decomposition holds

[TABLE]

The kernel of the first integral tends exponentially fast towards 0 and it defines a function that is regular near the origin $x=y=0$.

Proof.

We substract the asymptotic behavior

[TABLE]

Using lemma (2) we get

[TABLE]

and the result follows. ∎

We have exhibited the screened electrostatic in a dielectric slab as the usual Coulomb potential plus a correcting term. This term is exponentially small as the slab width is large compare to the relative distance between the two electric charges. Furthermore our expression (5) is regularized and do not present any divergence as the height $z$ approaches zero. This expression is suitable for further numerical calculations aiming to compute the binding energy of an exciton in a 2D materials. In the following we show how a multiscale approach allows to obtain the 2D potential when the width of the slab is very small with respect to the exciton Bohr radius.

3. Asymptotics of the spectral problem

In the previous analysis two scales enter into the problem, namely the slab width $d$ and the height of the electric charges $z_{0},z^{\prime}_{0}$. Working with quantities normalized by the slab width $d$ helped us exhibiting the regular part of the 2D electrostatic potential. Finding the exciton binding energy introduces another scale: the exciton Bohr radius $a_{0}$. We are interested in problems where the Bohr radius is large as compared to the slab width: $a_{0}\gg d$, that is, we focus on two-dimensional problems. Since now two length scales contribute to the problem, we need to introduce a new parameter $\eta=d/a_{0}$ in order to vary one length-scale independently of the other. Letting $\eta$ approach zero allows to deal with a 2D problem, by considering the electrostatic problem for a 2D sheet.

3.1. The multiple scale approach

We are now in a position to obtain the limit behavior of the Hamiltonian when $\eta$ tends to [math].

Theorem 2**.**

As $\eta\rightarrow 0$, the asymptotic 2D expansionn to first order in $\eta$, of the Hamiltonian describing the exciton is given by

[TABLE]

where the effective 2D potential $I_{0}$ is given by

[TABLE]

Proof.

Consider the spectral problem

[TABLE]

In order to obtain the asymptotic behavior when $\eta\rightarrow 0$, we put this expression in variational form. To do so we use a test function $\phi(r,\zeta,\zeta^{\prime})$ such that, for each $r\in[0,+\infty[$, the function $(\zeta,\zeta^{\prime})\rightarrow\phi(r,\zeta,\zeta^{\prime})$ belongs to $D(Y)$, i.e. the Schwartz space of $C^{\infty}$ functions with compact support in $Y$.

It holds

[TABLE]

Let us now divide this equality by $\eta^{2}$. Using Lebesgue theorem, we know that, for a continuous summable function $f(x,y)$ defined on $Y_{\eta}$, one has: $\lim_{\eta\rightarrow 0}\eta^{-2}\int_{Y_{\eta}}f(x,y)=f(0,0)$. Therefore, we obtain

[TABLE]

and

[TABLE]

where

[TABLE]

We conclude that, up to order $\eta$, the variational relation is

[TABLE]

The result follows. ∎

This result shows that contrarily to what could be intuitively believed, the effective 2D potential is not obtained from the 3D one simply by putting $z=z^{\prime}=0$. Let us consider more specifically what happens with the strictly Coulomb part of the potential. The effective potential in that case is

[TABLE]

whereas by putting $z=z^{\prime}=0$ one obtains simply $V_{c}=1/r$. In fig. (8), we have plotted both $V_{2D}$ and $V_{c}$. There it can be seen that, as $r$ tends to $+\infty$, the effective potential behaves as the Coulomb potential. However, as $r$ tends to [math] it behaves as

[TABLE]

where $\lambda=\frac{8}{3}\,10^{-4}$, an expression that was obtained by fitting.

In [1], the following 2D approximation is introduced (the expression is adapted in order to take into account the change of unit system)

[TABLE]

where

[TABLE]

It can be expressed using Struve and Neumann functions [1, 6] but this representation is of little interest since the integral in (10) can be computed numerically very easily. This form of the 2D potential was used in several articles, see e.g. [6], in order to compute the binding energy of excitons. Several authors have derived this expression using various approaches [3, 5]. Let us compare our result (7) and this potential. Using $\epsilon_{1}=\epsilon_{2}=1$ and $\epsilon_{f}=10$. The ratio $I_{K}/I_{0}$ is plotted in fig. 3. An excellent agreement is seen on a large range of values. It is only near the singularity that the potentials differ largely from each other. Our expression for an effective 2D potential is exact, as it comes from a rigorous limit analysis. Its validity is not limited to specific values of the variable $r$. As such it can take into account precisely the interaction between particles at small distances.

4. Conclusion

We have derived the asymptotic behavior of the Hamiltonian for the exciton wave function in a thin dielectric slab pf width $d$. Our approach is based on a limit analysis of the Hamiltonian using a suitable small parameter $\eta=d/a_{0}$, where $a_{0}$ is the Bohr radius. Our result compares well with the known results of the literature, although it is not limited to asymptotic values of the permittivity contrast. The method used here could be extended to obtain an expansion of the Hamiltonian with respect to $\eta$ [8]. We have also provided a regularized expression for the full 3D potential which makes explicit to what extend this potential departs from the Coulomb ones. This could be especially useful for implementing a perturbation analysis.

5. Appendix: Derivation of the screened potential.

5.1. Expression of the Green function

5.1.1. The equation in the Fourier domain

Let us compute then the potential $\varphi(\boldsymbol{\rho},z;z^{\prime}_{0})$ created by a unit charge situated at $(0,z^{\prime}_{0})$. It satisfies the Poisson equation in the Schwartz distributions meaning:

[TABLE]

and the conditions at infinity: $\lim_{z\rightarrow\pm\infty}\hat{\varphi}=0$.

By invariance of the medium in the $\boldsymbol{\rho}$ directions, a partial Fourier transform is performed. The Fourier transform is defined by:

[TABLE]

and the inverse transform is:

[TABLE]

The Fourier transform of $\hat{\varphi}$ satisfies, in the distributional meaning, the following differential equation:

[TABLE]

where: $k=|\mathbf{k}|$. For simplicity, the dependence of $\hat{\varphi}$ with respect to $\mathbf{k}$ and $z^{\prime}_{0}$ is implicit: we denote $\hat{\varphi}(z)$ instead of $\hat{\varphi}(\mathbf{k},z;z^{\prime}_{0})$.

5.2. The boundary conditions

There are four regions to be considered:

1. (1)

$z>d/2$, $-\hat{\varphi}^{\prime\prime}+k^{2}\hat{\varphi}=0$, hence, taking into account the condition at infinity: $\hat{\varphi}(z)=A_{2}\,e^{-k(z-d/2)}$, 2. (2)

$z^{\prime}_{0}<z<d/2$, $\hat{\varphi}(z)=A^{+}_{f}\,e^{-k(z-d/2)}+B^{+}_{f}\,e^{k(z-d/2)}$, 3. (3)

$-d/2<z<z^{\prime}_{0}$, $\hat{\varphi}(z)=A^{-}_{f}\,e^{-k(z-d/2)}+B^{-}_{f}\,e^{k(z-d/2)}$, 4. (4)

$z<-d/2$, $\hat{\varphi}(z)=B_{1}\,e^{k(z+d/2)}$

The boundary conditions at the interfaces of each domain are implied by eq.(12): the function $\hat{\varphi}$ is continuous everywhere and the function $\varepsilon\partial_{x}\hat{\varphi}$ is continuous everywhere except at $z^{\prime}_{0}$ where it as a jump: $\partial_{x}\hat{\varphi}(z_{0}^{+})-\partial_{x}\hat{\varphi}(z_{0}^{-})=-1/\varepsilon$.

For $\hat{\varphi}$ this gives the relations:

• •

$z=d/2$, $\hat{\varphi}(\left.\frac{d}{2}\right|^{+})=\hat{\varphi}(\left.\frac{d}{2}\right|^{-})$:

[TABLE]

• •

$z=z^{\prime}_{0}$, $\hat{\varphi}(\left.z^{\prime}_{0}\right|^{+})=\hat{\varphi}(\left.z^{\prime}_{0}\right|^{-})$:

[TABLE]

• •

$z=-d/2$, $\hat{\varphi}(\frac{d}{2}^{+})=\hat{\varphi}(\frac{d}{2}^{-})$:

[TABLE]

For the derivative, we obtain:

• •

$z=d/2$,

[TABLE]

• •

$z=z^{\prime}_{0}$,

[TABLE]

• •

$z=-d/2$,

[TABLE]

From [(13),(16)] and [(15),(18)],we get:

[TABLE]

In order to clarify the derivation, we denote:

[TABLE]

From (13) and (19), we get:

[TABLE]

where

[TABLE]

This last expression can be simplified. Let us denote:

[TABLE]

From (21), we get:

[TABLE]

We are looking first for the expression of $B_{f}^{+}$. It is obtained from (17), by using [(19)(20)(22)]:

[TABLE]

this gives, upon using the expression (22) for $\tau$

[TABLE]

where:

[TABLE]

5.3. Expression of the energy of interaction

The electrostatic energy between both charges is given by:

[TABLE]

For $z_{0}\geq z^{\prime}_{0}$, it holds:

[TABLE]

this gives:

[TABLE]

Using (23), it comes:

[TABLE]

and finally:

[TABLE]

There is a typo in the expression given in the paper by Keldysh: in the integral defining the energy of interaction, the term $e^{2\mathbf{k}\cdot\boldsymbol{\rho}}$ should be replaced by $e^{i\mathbf{k}\cdot\boldsymbol{\rho}}$.

Consider the double integral in polar coordinates: $(\rho,\theta)$: $\boldsymbol{\rho}=\rho(\cos\theta,\sin\theta)$ and $k(\cos\psi,\sin\psi)$: $\boldsymbol{\rho}\cdot\mathbf{k}=\rho\,k\cos(\theta-\psi)$. It is possible without loss of generality to take $\theta=\pi/2$, from invariance of the problem under a rotation around axis $Oz$.

[TABLE]

Let us recall the generating series for Bessel functions $J_{n}(\rho)$:

[TABLE]

we obtain: $J_{0}(\rho)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{i\rho\sin\psi}\,d\psi$, therefore it holds:

[TABLE]

Finally, we change to the new variable: $u=k/d$ to get:

[TABLE]