Two- and Three-Particle Complexes with Logarithmic Interaction: Compact Wave Functions for Two-Dimensional Excitons and Trions

TL;DR

This paper develops compact wave functions for two- and three-particle complexes with logarithmic interactions in two dimensions, accurately describing excitons and trions in TMDC monolayers and enabling perturbative studies of magnetoexcitons.

Contribution

It introduces new wave functions for 2D excitons and trions with logarithmic interactions, achieving high accuracy and applicability for perturbation theory and experimental comparisons.

Findings

Wave functions yield 5-6 correct decimal digits for exciton energies.

Logarithmic potential predicts trion binding energies within 25% of experimental data.

Wave function structure at small distances is characterized for Rytova-Keldysh interaction.

Abstract

Assuming a logarithmic interaction between constituent particles, compact and locally accurate wave functions that describe bound states of the two-particle neutral and three-particle charged complexes in two dimensions are designed. Prime examples of these complexes are excitons and trions that appear in monolayers of Transition-Metal DichalCogenides (TMDCs). In the case of excitons, these wave functions led to 5-6 correct decimal digits in the energy and the diamagnetic shifts. In addition, it is demonstrated that they can be used as zero-order approximations to study magnetoexcitons via perturbation theory in powers of the magnetic field strength. For the trion, making a comparison with experimental data for concrete TMDCs, we established that the logarithmic potential leads to binding energies $\lesssim 25\%$ greater than experimental ones. Finally, the structure of the wave function at small distances is established for excitons whose carriers interact via the Rytova-Keldysh potential.

Tables

Table 1. Table 1: Dimensionless optimal parameters and variational energies ε n ρ , 0 ( v ​ a ​ r ) superscript subscript 𝜀 subscript 𝑛 𝜌 0 𝑣 𝑎 𝑟 \varepsilon_{n_{\rho},0}^{(var)} of the first six S 𝑆 S -states. Diamagnetic shifts are also presented. Displayed numbers are rounded.

$n_{\rho }$	$A$	$B$	$C$	$D$	$F$	${\epsilon }_{n_{\rho },0}^{\left(var\right)}$	$\frac{1}{8}⟨{\rho }^2⟩$
0	-0.2259	0.7684	0.1879	1.9102	0.9561	0.179 935	0.1363
1	-2.1430	0.9613	-0.0655	0.7327	1.5413	1.314 677	1.0363
2	-3.5415	1.0256	-0.3938	0.6052	1.8695	1.830 608	2.8459
3	-4.6212	1.0659	-0.6910	0.5576	2.1177	2.168 874	5.5630
4	-5.4135	1.0935	-0.9265	0.5393	2.3216	2.421 054	9.1873
5	-6.1751	1.1193	-1.1668	0.5214	2.4880	2.622 221	13.7150
6	-6.9278	1.1438	-1.4163	0.5044	2.6286	2.789 590	19.1070

Table 2. Table 2: Logarithmic Potential, the Ground State ( 0 , 0 ) 0 0 (0,0) : First 10 sums of energy corrections. Digits in bold are invariant with respect to the next order correction. The value of ε e ​ x ​ a ​ c ​ t subscript 𝜀 𝑒 𝑥 𝑎 𝑐 𝑡 \varepsilon_{exact} was calculated with the function ( 20 ).

Approximation	Value	Correction	Value
${\epsilon }_0+{\epsilon }_1$	0.179 935 4 52 016 501 408 229 152	-${\epsilon }_2$	5.1787914339321731$\times {10}^{-8}$
${\epsilon }_0+{\epsilon }_1+{\epsilon }_2$	0.179 935 400  228 587 068 907 420	${\epsilon }_3$	9.8170705093752$\times {10}^{-11}$
${\epsilon }_0+\dots +{\epsilon }_3$	0.179 935 400 32 6 757 774 001 172	-${\epsilon }_4$	8.61541999679$\times {10}^{-13}$
${\epsilon }_0+\dots +{\epsilon }_4$	0.179 935 400 325  896 232 001 493	${\epsilon }_5$	9.200491318$\times {10}^{-15}$
${\epsilon }_0+\dots +{\epsilon }_5$	0.179 935 400 325 905  432 492 811	-${\epsilon }_6$	1.16392216$\times {10}^{-16}$
${\epsilon }_0+\dots +{\epsilon }_6$	0.179 935 400 325 905 31 6 100 595	${\epsilon }_7$	1.681331$\times {10}^{-18}$
${\epsilon }_0+\dots +{\epsilon }_7$	0.179 935 400 325 905 317 7 81 927	-${\epsilon }_8$	2.7053$\times {10}^{-20}$
${\epsilon }_0+\dots +{\epsilon }_8$	0.179 935 400 325 905 317 75 4 874	${\epsilon }_9$	4.76$\times {10}^{-22}$
${\epsilon }_0+\dots +{\epsilon }_9$	0.179 935 400 325 905 317 755 3 50	-${\epsilon }_{10}$	9.0$\times {10}^{-24}$
${\epsilon }_0+\dots +{\epsilon }_{10}$	0.179 935 400 325 905 317 755 341	${\epsilon }_{11}$	1.8$\times {10}^{-25}$
${\epsilon }_{exact}$	0.179 935 400 325 905 317 755 341

Table 3. Table 3: Energies ( ε 𝜀 \varepsilon ) of the first Low-Lying States with quantum numbers n ρ ≤ 4 subscript 𝑛 𝜌 4 n_{\rho}\leq 4 and | m | ≤ 5 𝑚 5 |m|\leq 5 , see ( 6 ). Results were obtained with N = 12 𝑁 12 N=12 , see ( 20 ). Underlined digits correspond to those digits reproduced by the variational energies. All printed digits are exact: confirmed with N = 13 𝑁 13 N=13 .

	0	1
0	0.179 935 400 325 905 317 755 341	1.039 612 607 367 968 583 608 037
1	1.314 677 846 047 317 635 438 844	1.662 901 190 508 306 406 113 371
2	1.830 608 839 744 414 785 298 073	2.047 765 063 110 404 237 580 088
3	2.168 874 146 054 584 411 366 434	2.326 094 048 304 208 876 062 166
4	2.421 054 965 033 637 757 825 116	2.544 033 274 577 971 448 108 208
	2	3
0	1.497 798 460 867 032 070 612 310	1.811 273 253 112 008 598 564 854
1	1.929 287 879 273 176 751 640 227	2.141 542 186 466 426 635 213 589
2	2.233 478 680 963 778 182 582 480	2.392 481 446 049 754 044 603 719
3	2.467 896 772 583 752 523 081 303	2.594 392 795 360 084 586 946 462
4	2.658 389 492 811 470 794 957 093	2.763 108 473 682 027 259 596 914
	4	5
0	2.049 706 164 599 668 581 995 749	2.242 142 115 820 400 875 545 543
1	2.317 307 924 417 713 851 799 756	2.467 097 474 123 876 210 404 164
2	2.530 696 526 629 787 083 928 815	2.652 645 890 317 996 675 838 214
3	2.707 805 246 985 471 371 167 548	2.810 282 701 535 865 116 991 903
4	2.859 009 046 649 596 541 890 608	2.947 167 432 155 959 952 276 143

Table 4. Table 4: Ground state (0,0): Numerical coefficients ε n subscript 𝜀 𝑛 \varepsilon_{n} , n = 0 , 1 , … , 11 𝑛 0 1 … 11 n=0,1,...,11 of the series expansion ( 23 ) for ε 𝜀 \varepsilon calculated in non-linearization procedure.

$n$	${\epsilon }_n$	$n$	${\epsilon }_n$
0	0.179 935 400	6	$-\mathrm{0.018\hspace{0.17em}000\hspace{0.17em}868}$
1	0.136 337 679	7	0.030 034 032
2	$-\mathrm{0.023\hspace{0.17em}813\hspace{0.17em}717}$	8	$-\mathrm{0.057\hspace{0.17em}913\hspace{0.17em}120}$
3	0.0129 299 044	9	0.127 105 525
4	$-\mathrm{0.011\hspace{0.17em}145\hspace{0.17em}815}$	10	$-\mathrm{0.314\hspace{0.17em}135\hspace{0.17em}328}$
5	0.012 765 279	11	0.867 136 787

Table 5. Table 5: Trion dimensionless energy ε ​ ( σ ) 𝜀 𝜎 \varepsilon(\sigma) vs the mass ratio σ 𝜎 \sigma , see ( 36 ). Variational energies obtained through trial functions ( 37 ) and ( 40 ). For the latter, calculations correspond to N = 1 𝑁 1 N=1 , which were confirmed with N = 2 𝑁 2 N=2 . Results rounded to first display decimal digits.

	Wave function			Wave function
$\sigma$	Eq. (37)	Eq. (40)	$\sigma$	Eq. (37)	Eq. (40)
0	0.0871	0.0791	1	0.42265	0.4182
0.1	0.1365	0.1285	2	0.6091	0.6058
0.2	0.1809	0.1729	3	0.7450	0.7411
0.3	0.2202	0.2165	4	0.8506	0.8473
0.4	0.2563	0.2539	5	0.9378	0.9347
0.5	0.2891	0.2868	6	1.0117	1.0090
0.6	0.3194	0.3166	7	1.0763	1.0737
0.7	0.3476	0.3443	8	1.1333	1.1308
0.8	0.3740	0.3703	9	1.1845	1.1821
0.9	0.3990	0.3949	10	1.2309	1.2286

Table 6. Table 6: Negative charged trion binding energy in meV for different suspended TMDC monolayers. Results rounded to first displayed significant digits. Values of ρ 0 subscript 𝜌 0 \rho_{0} , σ 𝜎 \sigma , and experimental binding energies taken from reference [ 45 ] .

Material			$E_b$ (meV)
	${\rho }_0$ (Å)	$\sigma$	Present Results	Experiments
MoS2	39	0.81	41	34, 35
MoSe2	40	0.86	39	30
WS2	38	0.84	41	34, 36
WSe2	45	0.85	35	30

Equations

$$V_{ij}(\rho )=-\frac{q_i{q}_j}{{\rho }_0}\ln \left(\frac{\rho }{{\rho }_0}\right).$$

$${\rho }_0=\frac{d({\epsilon }_{||}-1)}{2},$$

$$-\frac{{\hslash }^2}{2\mu }\left({\partial }_{\rho }^2\psi +\frac{1}{\rho }{\partial }_{\rho }\psi \right)+\left(\frac{{\hslash }^2{m}^2}{2\mu {\rho }^2}+\frac{e^2}{{\rho }_0}\ln \left(\frac{\rho }{{\rho }_0}\right)\right)\psi =E\psi ,\rho \in [0,\infty ),$$

$$\rho \to {\left(\frac{{\hslash }^2{\rho }_0}{\mu e^2}\right)}^{-\frac{1}{2}}\rho ,$$

$$-\frac{1}{2}\left({\partial }_{\rho }^2\psi +\frac{1}{\rho }{\partial }_{\rho }\psi \right)+\left(\frac{m^2}{2{\rho }^2}+\ln \left(\rho \right)\right)\psi =\epsilon \psi .$$

$$E=\frac{e^2}{{\rho }_0}\epsilon -\frac{e^2}{2{\rho }_0}\ln \left(\frac{\mu e^2{\rho }_0}{{\hslash }^2}\right).$$

$${\int }_0^{\infty }|\psi (\rho ){|}^2\rho d\rho <\infty$$

$${\psi }_{n_{\rho },m}(\rho )={\rho }^{|m|}{P}_{n_{\rho },m}({\rho }^2)\exp (-{\Phi }_{n_{\rho },m}(\rho )).$$

$$P_{n_{\rho }}({\rho }^2)=\prod _{i=1}^{n_{\rho }}({\rho }^2-{\rho }_i^2),P_0(\rho )\equiv 1.$$

$${\psi }_{0,0}(\rho )=\exp (-{\Phi }_{0,0}(\rho )).$$

$${\Phi }_{0,0}(\rho )=\sum _{i=1}^{\infty }\sum _{j=0}^i{a}_{ij}{\rho }^{2i}{\ln }^j(\rho ),\rho \to 0,$$

$${\psi }_{0,0}(\rho )=\sum _{i=0}^{\infty }\sum _{j=0}^i{b}_{ij}{\rho }^{2i}{\ln }^j(\rho ),\rho \to 0,$$

$${\Phi }_{0,0}(\rho )=\rho \left(\sqrt{2\ln \rho }+\frac{1-\epsilon }{\sqrt{2\ln \rho }}+\mathcal{O}\left((\ln \rho {)}^{-3/2}\right)\right)+\frac{1}{2}\ln \rho +...,\rho \to \infty$$

$${\psi }_{n_{\rho },m}(\rho )={\rho }^{|m|}\sum _{i=0}^{\infty }\sum _{j=0}^i{b}_{ij}{\rho }^{2i}{\ln }^j(\rho ),\rho \to 0,$$

$${\Phi }_{n_{\rho },m}^{(approx)}(\rho )=\frac{A+(1-2B){\rho }^2\ln \rho +(C+B{\rho }^2)\ln (1+D{\rho }^2)}{\sqrt{F^2(1+D{\rho }^2)+\frac{1}{2}{\rho }^2\ln \rho }}-\frac{A}{F},$$

$${\psi }_{n_{\rho },m}^{(approx)}(\rho )={\rho }^{|m|}{P}_{n_{\rho }}({\rho }^2)\times \exp \left(-\frac{A+(1-2B){\rho }^2\ln \rho +(C+B{\rho }^2)\ln (1+D{\rho }^2)}{\sqrt{F^2(1+D{\rho }^2)+\frac{1}{2}{\rho }^2\ln \rho }}+\frac{A}{F}\right),$$

$$⟨{\psi }_{n_{\rho ,m}}^{(approx)}|{\psi }_{n_{\rho }^{\prime },m}^{(approx)}⟩=0,n_{\rho }^{\prime }=0,1,...,n_{\rho }.$$

$${\psi }_{n_{\rho },m}^{(exact)}={\psi }_{n_{\rho },m}^{(approx)}(1-{\varphi }_1-{\varphi }_2...).$$

$${\epsilon }_{n_{\rho },m}^{exact}={\epsilon }_{n_{\rho },m}^{(var)}+{\epsilon }_2+{\epsilon }_3+...,$$

$${\psi }_{n_{\rho },m}(\rho )\approx {\psi }_{0,m}^{(approx)}(\rho )\times \sum _{i=0}^N\sum _{j=0}^i{c}_{i,j}^{(n_{\rho },m)}{\rho }^{2i}{\ln }^j(\rho )$$

$$-\frac{1}{2}\left({\partial }_{\rho }^2\psi +\frac{1}{\rho }{\partial }_{\rho }\psi \right)+\left(\ln \left(\rho \right)+\frac{{\gamma }^2}{8}{\rho }^2\right)\psi =\epsilon \psi .$$

$$\gamma =\frac{B}{B_0},B_0=\frac{c{e}^3{\mu }^2}{{\hslash }^3}.$$

$$\epsilon (\gamma )={\epsilon }_0+\sum _{n=1}^{\infty }{\epsilon }_n{\gamma }^{2n}.$$

$$\varepsilon_{1}\ =\ \frac{1}{8}\left\langle\rho^{2}\right\rangle\$$

$$\Delta \epsilon (\gamma )=\sum _{n=1}^{\infty }{\epsilon }_n{\gamma }^{2n}$$

$${\hat{H}}_T\psi =E\psi ,\psi =\psi (u_1,u_2,u_3),$$

$$\displaystyle\hat{H}_{T}=\ -\frac{\hbar^{2}}{2}\Delta\ -\ \frac{e^{2}}{\rho_{0}}\ln\left(\frac{u_{1}}{\rho_{0}}\right)\ +\ \frac{e^{2}}{\rho_{0}}\ln\left(\frac{u_{2}}{\rho_{0}}\right)+\ \frac{e^{2}}{\rho_{0}}\ln\left(\frac{u_{3}}{\rho_{0}}\right)\$$

$$\Delta =$$

$$+\frac{1}{m_1}\frac{u_2^2+u_3^2-u_1^2}{u_2{u}_3}{\partial }_{u_3}{\partial }_{u_2}+\frac{1}{m_2}\frac{u_1^2+u_3^2-u_2^2}{u_1{u}_3}{\partial }_{u_1}{\partial }_{u_3}+\frac{1}{m_3}\frac{u_1^2+u_2^2-u_3^2}{u_1{u}_2}{\partial }_{u_1}{\partial }_{u_2}.$$

$${\mu }_{ij}=\frac{m_i{m}_j}{m_i+m_j},i,j=1,2,3.$$

$$dV\propto u_1{u}_2{u}_3{S}^{-1}d{u}_{1}d{u}_{2}d{u}_3,$$

Full text

Introduction

It was long-time ago established [1, 2] that a logarithmic interaction between free electrons $(e^{-})$ and holes $(h)$ may occur in two-dimensional semiconductors as a limiting case of the Rytova-Keldysh potential222See Section III for a detailed discussion based on excitons.. Let us consider a thin film/layer with a dielectric surrounding made of dielectric substrates. If the thickness is sufficiently small compared with the exciton Bohr radius, a logarithmic interaction between carriers confined to the film emerges as a result of the polarization of atomic orbitals [3]. Monolayers of Transition-Metal DichalCogenides (TMDCs) are relevant examples of this kind of film-substrate configurations. Such monolayers can be considered as two-dimensional since their thickness is a few Angstroms [4]. The low dimensionality (planar) and dielectric screening in such materials result in a strong electrostatic interaction. It allows the existence of stable bound states for complexes composed by electrons and holes, such as excitons ($e^{-},h$), and positively and negatively charged trions: $(e^{-},h,h)$ and $(e^{-},e^{-},h)$, respectively. Monolayers of TMDCs have recently received special attention for being ideal candidates for potential applications in optoelectronics [5, 6], valleytronics [7], enhanced photoluminescence [8, 9], and systems with pronounced many-body effects [5]. In particular, the optical response of these materials is explained in terms of complexes [3]. For a multilayer configuration of TMDCs, the pair-wise interaction between two carriers of charge $q_{i}$ and $q_{j}$ within the same layer remains logarithmic333For a detailed discussion and derivation of the interaction see [10], the supplemental material of [3]., having the same form as for monolayers:

[TABLE]

Here, $\rho$ denotes the relative distance between the carriers, whose electric charges are $q_{i}$ and $q_{j}$, and

[TABLE]

where $\varepsilon_{||}$ is the in-plane component of the dielectric permittivity tensor of the bulk layered material, and $d$ is the distance between layers.

In the effective-mass approximation, the quantum mechanical description of complexes is governed by the Schrödinger equation. In this context, the variational method has shown to be an adequate tool to study bound states of excitons and trions, see [11, 12, 13, 14, 15, 16] and references therein. However, as mentioned in [17], difficulties in constructing a reasonable wave function Ansatz in the case of larger complexes than excitons hinders the straightforward extension of the variational consideration. To overcome this drawback, some recent advances have been made in constructing adequate wave functions for intralayer trions [18, 14]. However, most of the variational functions favor the simplicity in calculations leading to reasonable results in terms of energy, instead of a correct description of the wave function. The main goal of this work is to show that, for two/three-particle complexes, both accurate energies and wave functions can be simultaneously achieved by using adequate compact trial functions.

In the present study, we construct compact Ansätze (trial functions) associated with the intralayer exciton of multi and monolayers made of TMDCs using the internal structure of the exact wave functions in the logarithmic regime of interaction. For unclear reasons to the authors, this approach has not been studied so far in a variational consideration. We focus on the construction of locally accurate approximations of the wave functions. They are valuable not only for finding energies and expectation values with high accuracy, but also for shedding some light on physical properties of the exact wave functions. In fact, the search for compact wave functions describing few-particle charged systems is an active field of research [19]. For example, they are widely used in atomic physics due to their usefulness to compute efficiently scattering cross sections [20] and matrix elements of singular operators [21]. In quasi-two-dimensional materials, the slow convergence of CI (configuration interaction) functions [22, 23, 24] has motivated the search for simple yet accurate wave functions for the description of hole-electron interaction.

As we show, our approximate solutions (taken as zero order approximation) lead to a convergent perturbation series to the exact solution. Two concrete physically relevant examples of application are discussed: (i) magneto-excitons in a weak field regime; and most importantly (ii) the construction of three-particle wave functions.

The present work is organized as follows. In Section I, we discuss the construction of compact parameter-dependent exciton wave functions. Concrete variational calculations for the spectrum of the first low-lying states are presented. We investigate the accuracy of the energy estimates using the non-linearization procedure [25] and an alternative variational trial function. Then, we study the effect of a weak uniform magnetic field to the energy spectrum using perturbation theory taking our compact functions as zero order approximation. In this line, we discuss the (re)summation of perturbation series using Padé approximants. In Section II, taking as building block the trial function constructed for the exciton, we propose a trial function for the ground state of the trion. We study the binding energy of the complex for concrete TMDCs and compare it with experimental results. A simple formula for the trion energy is provided. In Section III, we discuss extensions of our consideration to the Rytova-Keldysh potential. Finally, we summarized our results in Conclusions.

I Two-Particle System: Excitons

Consider the neutral system made of two charged particles (hole and electron) interacting through potential (1). After separating the motion of the center of mass and using polar coordinates, we arrive at the familiar two-dimensional radial Schrödinger equation for the relative motion,

[TABLE]

where $\mu$ is the reduced mass444By definition $\mu\ =\ \dfrac{m_{e}\,m_{h}}{m_{e}+m_{h}}$, where $m_{e}$ and $m_{h}$ are the electron and hole effective masses, respectively., $e$ denotes the charge of the hole, and $\hbar$ is the reduced Planck constant. Any energy and wave function can be labeled by ($n_{\rho},m$), but for simplicity we drop such labels for now. The radial quantum number takes the values $n_{\rho}=0,1,...$, meanwhile the magnetic quantum number $m=0,\pm 1,...$ . Using the transformation via the dimensionless coordinate

[TABLE]

we remove the explicit presence of the constants $\hbar$, $\mu$, $\rho_{0}$, and $e$ from the Schrödinger equation which is transformed into its dimensionless form

[TABLE]

In this equation, $\varepsilon$ plays the role of dimensionless energy, and it is related to $E$ through

[TABLE]

The second term in (6) only provides the reference point to measure energies, and it has no relevant role. From (6), it is clear that the energy difference between two arbitrary states does not depend on the reduced mass $\mu$. Since we are interested in bound states, we impose boundary conditions on (5) such that the normalizability requirement

[TABLE]

is fulfilled. Under these considerations, equation (5) does not admit an exact solution: energies and wave functions can only be found in approximate form. Interestingly, the same spectral problem defined by equation (5) can appear in another context. Indeed, the effect of a wiggly cosmic string for both mass-less and massive particle propagation along the string axis is governed by (5), see [26].

I.1 Ground State

Relevant information concerning the structure of the wave function can be revealed using asymptotic analysis. It is convenient to adopt the exponential representation of the wave function, namely

[TABLE]

The unknown function $\Phi_{n_{\rho},m}(\rho)$ is called phase, while $P_{n_{\rho}}$ is a polynomial of degree $n_{\rho}$. Explicitly,

[TABLE]

For a given state, this polynomial is determined by the position of the nodes $\rho_{i}$, $i=1,...,n_{\rho}$. Thus, representation (8) is unambiguous [25]. The asymptotic series of the ground state phase, $\Phi_{0,0}(\rho)$, shares properties with those for excited states. Hence, we focus on the quantum numbers ($n_{\rho}=0,m=0)$ from now on. In this case, (8) takes the form

[TABLE]

We construct the asymptotic series around two relevant points of the domain: $\rho=0$ and $\rho=\infty$; small and large relative distances, respectively. At $\rho=0$, it can be demonstrated that the asymptotic series of the phase has the following structure

[TABLE]

where $a_{ij}$ are coefficients with dependence on $\varepsilon$. As a consequence, the wave function has a similar asymptotic series,

[TABLE]

where $b_{ij}$ are $\varepsilon$-dependent coefficients. On the other hand, the first terms of the asymptotic series of the phase at $\rho=\infty$ are

[TABLE]

Note that the dominant term, $\mathcal{O}(\rho\sqrt{\ln\rho})$, does not depend on the energy. Therefore, the same leading term is expected for any bound state. Contrary to the series at $\rho=0$, see (11), we were unable to find the general structure of the phase at $\rho=\infty$. However, for the particular purpose of this work, this piece of information is not required. At $\rho\rightarrow 0$, the phase of any state has the same structure of series (11). On the other hand, the wave function of any state has the following asymptotic expansion,

[TABLE]

with coefficients $b_{ij}$ depending on $\varepsilon_{n_{\rho},m}$, and $|m|$. Using series (14), the real solutions of the equation $\psi_{n_{\rho},m}(L)=0$, with $L>0$ sufficiently large, define low-accuracy approximations to the low-lying energies and wave functions.

I.2 Compact Trial Function

Based on series (11), (12), and (13), we constructed a parameter-dependent approximation for the wave function of an arbitrary state $(n_{\rho},m)$. We followed the prescription described in [27], where it was applied successfully to the quartic anharmonic oscillator. Such prescription establishes the following: the approximate phase is the result of matching the series (11) and (13) in a minimal way, reproducing as many dominant growing terms in (13) as possible. This procedure is almost unambiguous, and it leads to

[TABLE]

where $\{A,B,C,D,F\}$ are five $(n_{\rho},m)$-dependent dimensionless parameters. By construction, the phase (15) reproduces functionally (same structure, but different coefficients) all the terms in the expansion at small distances (11), but only the leading one in the expansion at large distances (13). The approximate trial wave function of an arbitrary state ($n_{\rho},m$) is given ultimately by

[TABLE]

where the polynomial $P_{n_{\rho}}(\rho^{2})$ is of the form (9), and it carries the information about the nodes. To fix the value of free parameters, we use the variational method imposing the orthogonalization constraints

[TABLE]

These constraints define the position of the nodes. Thus, after fixing $m$, we move sequentially from the ground state $(n_{\rho}=0)$ to a higher excited state $n_{\rho}>0$. Under these constraints, we calculate the parameter-dependent expectation value of the Hamiltonian associated with (5), usually called variational energy and denoted by $\varepsilon^{(var)}_{n_{\rho},m}$. Then, using an optimization procedure, we can find the configuration of parameters that minimize the variational energy. Only for states with quantum numbers $(n_{\rho}=0,m)$, the variational principle guarantees that the variational energy is an upperbound of the exact energy. To ensure the square-integrability of (16), the constraint $D>0$ is imposed.

Concrete numerical calculations were carried out for the first low-lying states with $n_{\rho}\leq 6$ and $|m|\leq 4$. As a result, the variational energy was found with relative accuracy $\sim 10^{-6}$ or less. In Table 1, we present the optimal parameters and the variational energies for the first seven $S$-states. Meanwhile, nodes are obtained with 5 exact significant digits. It was confirmed by using the non-linearization procedure and making alternative accurate variational calculations, see below. The optimized variational energies reach and sometimes overcome the best results found in the literature [28, 29, 30, 3].

By construction, our wave functions $\psi_{n_{\rho,m}}^{(approx)}$ are locally accurate, mainly due to the correct asymptotic behavior of the wave function at small and large distances. To check this, we can estimate the local accuracy of $\psi_{n_{\rho},m}^{(approx)}(\rho)$ using the non-linearization procedure, via the perturbation series [25]:

[TABLE]

Furthermore, the non-linearization procedure dictates that

[TABLE]

if $\psi_{n_{\rho},m}^{(approx)}$ is chosen according to the above mentioned prescription. Numerical calculations for $S$-sates established that corrections $\phi_{n}$ in (18) decrease as $n$ grows, specifically $|\phi_{n+1}/\phi_{n}|\lesssim 10^{-1}$ for all $\rho$. Therefore, the local accuracy of wave functions is guaranteed. In turn, $|\varepsilon_{n+1}/\varepsilon_{n}|\lesssim 10^{-2}$, which indicates a fast rate of convergence of (19). In particular, the value of $\varepsilon_{2}$ suggests that our variational calculations lead to energies with 5 - 6 exact decimal digits. For the states considered, corrections in (18) and (19) were calculated using the Mathematica codes described in [31]. In Table 2, we present explicit values of the first 11 corrections $\varepsilon_{n}$ in (19) for the ground state (0,0). For this state, it is enough to consider the first three terms in (19) to reach the accuracy provided by the finite-element calculations, see [32].

I.3 Alternative Trial Function

The orthogonalization procedure used to determine the position of the nodes is accurate but impractical for highly excited states. For a given $m$, the $n_{\rho}$ excited state requires $n_{\rho}-1$ constraints in order to fulfill (17). To overcome this drawback, an alternative and efficient procedure is discussed in this Section. This approach only requires the knowledge of the ground state functions $(0,m)$. Using the wave functions $\psi_{0,m}^{(approx)}(\rho)$ with optimal parameters, we construct the expansion

[TABLE]

to describe any state. The factor $\psi^{(approx)}_{0,m}$ guarantees the correct the asymptotic dominant behavior at large $\rho$ in our trial function (20). Based on the series at small $\rho$, an additional factor is introduced in the form of a partial sum of (12) in order to improve the small-$\rho$ behavior. In this representation, any wave function (20) contains $(N+1)(N+2)/2$ terms. Using the linear variational principle, it is known that the energies and coefficients $c_{i,j}^{(m)}$ are determined by the secular equations (Rayleigh-Ritz method). Thus, this alternative procedure is a two-step variational consideration. To solve the secular equations, we use the Löwdin orthogonalization procedure [33] since the set of functions $\{\rho^{2i}\ln^{j}\rho\}$ is not orthogonal. In Table 3, we present the energies of the low-lying states with quantum numbers $n_{\rho}\leq 4$ and $|m|\leq 5$ using $N=12$. For all states considered, numerical results indicate that the rate of convergence is about 3-4 correct digits with an increment of $N$ to $N+1$. In fact, using $N=12$, we established 24 exact decimal digits for energies of all states considered. Comparing with our variational results, we established that our locally accurate approximations (16) lead to energies with 5-6 exact decimal digits, which is in agreement with our calculations via the non-linearization procedure.

I.4 Magnetoexcitons: S-States

Consider an exciton subjected to a time-independent magnetic field $\boldsymbol{B}=B\hat{\boldsymbol{z}}$.555Therefore, the magnetic field is transversal to the thin film. As long as the momentum of the center of mass is zero, it can be shown [34, 35, 36] that the Schrödinger equation666In the symmetric gauge. for the relative motion describing $S$-states is

[TABLE]

where

[TABLE]

Equation (21) is written in variable (4). Note that $B_{0}$ is the exciton unit of magnetic field while $\varepsilon$ is defined through (6). For weak magnetic fields, perturbation theory for $\varepsilon({\gamma})$ can be constructed in the form

[TABLE]

The first order correction ($\varepsilon_{1}$), given by

[TABLE]

is called the diamagnetic shift. For $S$-states, the experimental measurement of $\varepsilon_{1}$ reveals physical properties such as the exciton mass, size, and spin [37]. Taking our compact approximations, we calculated the diamagnetic shifts of the first $S$-states with $0\leq n_{\rho}\leq 6$ with accuracy of 5 significant digits. Results are presented in Table 1. To calculate higher order corrections $\varepsilon_{n>1}$, we used the non-linearization procedure as described in [38]. There, it was applied to construct the strong coupling expansion of the cubic anharmonic oscillator, see Subsection A. In this way, we calculated higher corrections $\varepsilon_{n>1}$. Using as zero-order approximation (16), we established corrections with 9 exact decimal digits. In Table 4, we present the first eleven corrections in (23) for the ground state. Following Dyson’s argument, it is expected to be a divergent series. However, numerical results suggest that the second term in (23), namely

[TABLE]

is an alternating series and Padé resummable. For example, using a Padé approximant $P_{5}^{6}(\gamma^{2})$ based on the first 11 corrections leads to accurate results in the domain $0\leq\gamma\leq 1$ with 4 exact significant digits. In Fig. 1, we presented the plots of $\Delta\varepsilon(\gamma)$ for the first seven $S$-states calculated through Padé approximants $P_{5}^{6}(\gamma^{2})$.

II Three-Particle System: Trions

In this Section, we consider the lowest bound state of the two-dimensional complex made of three charged particles of (effective) masses $(m_{1},m_{2},m_{3})$ and charges $(-e,e,e)$, respectively. The pairwise interaction between them is assumed to be logarithmic, and it is given by (1). In Fig. 2, the geometrical setting of the system is shown.

The ground state is a square-integrable eigenfunction with dependence only on the relative distances $u_{1}=r_{23}$, $u_{2}=r_{13}$, and $u_{3}=r_{12}$, and it is the lowest eigenfunction of the Hamiltonian777The ranges of the three variables $u_{i}$ are coupled and satisfy a triangle condition: their lengths must be such that they can form a triangle [39].

[TABLE]

[TABLE]

with [40]

[TABLE]

where

[TABLE]

are the reduced masses. The operator $\hat{H}_{T}$ is self-adjoint with respect to the volume element

[TABLE]

We set $m_{2}=m_{3}$, as it usually appears in trions888It has been established that for monolayers of TMDCs based of Tungsten, electrons in negative trions may have different effective masses, see [41].. Under this assumption,

[TABLE]

Next, we introduce the following change of variables

[TABLE]

and the parameter

[TABLE]

to remove in (26) the appearance of physical constants (except for the masses $m_{1}$ and $m_{2}=m_{3}$). Therefore, the trion ground state wave function corresponds to the nodeless solution of the dimensionless equation

[TABLE]

where

[TABLE]

In (34), $\varepsilon(\sigma)$ is related to the total energy ($E$), as follows

[TABLE]

In contrast to the energy levels of an exciton, the difference in energy of two arbitrary levels does depend on the masses through $\sigma$.

II.1 Trial Function

Based on the functions (16) constructed for excitons in Section I, we design compact trial functions for a trion in its ground state. Since we assumed that $m_{2}=m_{3}$, the positively charged trion $(e^{-},h,h)$ and the negative one ($h,e^{-},e^{-}$) are described by the same Schrödinger equation999When the constituent particles have different masses, they are distinguishable. Consequently, no symmetry requirement is imposed to the total wave function. In particular, in (37) the symmetrizer is no longer needed.. Their total wave functions, which include spin and valley quantum numbers, are antisymmetric with respect to the permutation of identical particles, see [18]. In particular, the spatial part of the ground state wave function must be symmetric101010According to our variational calculations, the anti-symmetric $S$-state is repulsive and does not correspond to a bound state. This feature was already established for the Rytova-Keldysh potential in the non-logarithmic regime, see [18, 42]. under this permutation. Therefore, we propose the following symmetric wave function

[TABLE]

where the symmetrizer operator is given by

[TABLE]

Here $\hat{P}_{23}$ is the permutation between the identical particles $2\rightarrow 3$ in $u$-variables. The form of the phase is given in (15). Note that $\Phi_{0,0}$ is given by formula (15) with optimized parameters found in Table 1. The construction of the previous wave function is motivated as follows. The factor $e^{-\Phi_{0,0}(\alpha^{2}\,u_{2})-\Phi_{0,0}(\beta^{2}\,u_{3})}$ describes the trion when the repulsive term of the interaction is absent. In turn, the factor $1+\gamma u_{1}^{2}\ln(u_{1}^{2})$ takes into account the polarization of the complex and describes the repulsion between equally charged carriers. This fact can be easily established using asymptotic analysis in (34) at small $u_{1}$, which leads to the expansion

[TABLE]

Parameters $\{\alpha,\beta\}$ were introduced in such a way that (i) they admit the meaning of screening charges if $m_{1}$ is kept fixed; and (ii) the square-integrability of the trial function is guaranteed. Altogether, the trial function (37) contains only three free real parameters: $\{\alpha,\beta,\gamma\}$. To fix their values, we use the variational method to find the optimal configuration that minimizes the variational energy. Numerical calculations were performed in perimetric coordinates using a modification of the FORTRAN code described in [43].

We carried out calculations for $\sigma\in[0,10]$, thus, covering representative TMDCs (see below). In Fig. 3, we present the optimal parameters as functions of $\sigma$ in the domain $\sigma\in[0,10]$. They have a smooth behavior, especially at large $\sigma$. We found that $\beta>\alpha$ in the domain considered. It hints that one of the two equally charged carriers is closer to the opposite-charged one. Compared to $\alpha$ and $\beta$, parameter $\gamma$ is smaller. It reaches it maximum value at $\sigma\sim 1-2$. In Table 5, optimized variational energy $\varepsilon(\sigma)$ is presented.

II.2 Alternative Trial Function

To investigate the accuracy of the variational calculations shown in the previous Section, we use an alternative trial function to estimate the energy of the trion in the lowest $S$-state with higher accuracy. In a similar way it was done for the exciton, see (20), we proposed

[TABLE]

where $c_{i,j,k}$ are coefficients to be determined. Without loss of generality, we set $c_{0,0,0}=1$ as normalization of the approximate wave function. Note that at $N=0$, trial functions (37) and (40) coincide. Meanwhile at $N>0$, they differ due to the insertion of the polynomial prefactor $\sim u_{1}^{2i}u_{2}^{2j}u_{3}^{2k}$. This insertion corrects the behavior of the wave function at small distances, leading to a higher accuracy for the ground state energy. To fix the value of the parameters $\{\alpha,\beta,\gamma\}$ and $c_{i,j,k}$’s, we followed a two-step variational consideration. First, we fix the value of parameters $\{\alpha,\beta,\gamma\}$ according to the variational calculations of the previous Section, in which the polynomial prefactor was not considered. The remaining free parameters, $c_{i,j,k}$’s, are defined by means of the secular equations. To solve them, we used the Löwdin orthogonalization procedure. In Table 5, we present the value of $\varepsilon(\sigma)$ calculated with $N=1$ in (40) for representative values of $\sigma$. Using $N=2$, we confirmed the results for $N=1$ up to the first four decimal digits. It indicates a fast rate of convergence with respect to $N$. In this way, we concluded that the wave function (37) provides an accuracy of two decimal digits for all $\sigma>0$.

The plot of $\varepsilon(\sigma)$, calculated through trial function (40) with $N=1$, is shown in Fig. 4. The curve described by $\varepsilon(\sigma)$ can be easily interpolated with an accuracy of 2 decimal digits by the simple function

[TABLE]

in the domain $\sigma\in[0,10]$.

II.3 Comparison with Experimental Data

In order to compare our results with experimental data, we calculated the binding energy ($E_{b}$). It is defined as the amount of energy needed to dissociate a trion into a neutral complex (exciton) and a free carrier. Therefore,

[TABLE]

where the energies correspond to ground states. According to (6), (36), and (41), we have a compact expression for $E_{b}$,

[TABLE]

Hence, the binding energy for trions is a function of $\sigma$.

In Table 6, we compare theoretical binding energies given by (43) with experimental ones for negatively charged trions111111Thus, we take $m_{1}=m_{h}$ and $m_{2}=m_{e}$. in molybdenum (Mo) and tungsten (W) dichalcogenide materials. As previously established [44], we found that logarithmic potential overbinds the trion, leading to larger binding energies than those reported in experiments. However, our results do not confirm the statement that the logarithmic potential leads to binding energies 50% larger, see [44]. The largest deviation occurs for MoSe2, reaching a deviation in energies of 30%. In turn, the smallest deviation is found for WS2, for which is 16%. Interestingly, either formula (43) or experimental results predicts equal binding energy for MoS2 and WS2.

III Toward the Rytova-Keldysh potential

If $\rho_{0}$ in (1) is comparable to the exciton Bohr radius, the logarithmic potential (1) is no longer suitable to describe the interaction between carriers inside monolayers of TMDCs. Furthermore, it overbinds three-particle complexes as we have seen. Under these circumstances, an adequate description of the interaction is given by the celebrated Rytova-Keldysh potential [1, 2], namely

[TABLE]

where

[TABLE]

Here $H_{0}$ and $Y_{0}$ are the Struve and Bessel function of second kind, respectively. At small $\rho$, the Rytova-Keldysh potential is reduced to the logarithmic potential as follows

[TABLE]

where $\Gamma=\gamma-\ln 2$; here $\gamma$ denotes the Euler-Mascheroni constant. The constant $\Gamma$ fixes the reference point for the energy. Therefore, it plays no relevant physical role when studying binding energies of complexes interacting within the logarithmic regime. On the other hand, at large distances

[TABLE]

Therefore, in this limit, the dominant interaction is given by the Coulomb one. The radial Schrödinger equation that describes excitons now reads

[TABLE]

To construct a compact trial function for excitons, we can follow the approach presented in Section I.B, in which the main ingredient is the asymptotic series. Using the transformation (4), we obtain the dimensionless Schrödinger equation

[TABLE]

where121212Here $a_{0}$ denotes the exciton Bohr radius defined as $a_{0}=\dfrac{\hbar^{2}}{\mu\,e^{2}}$. Meanwhile, $\mathcal{E}_{0}$ represents the exciton Rydberg constant.

[TABLE]

Note that $\lambda$ corresponds to the ratio of the two length scales of the system. According to parameters found in [45], $\lambda\lesssim 8\times 10^{-2}$ for all materials presented in Table 6. Therefore, it suggests taking the potential in equation (49) and expanding it in powers of $\lambda$,

[TABLE]

Therefore, the logarithmic potential is dominant131313The same conclusion for the $N$-particle complex can be established using similar arguments. at small $\lambda$. For arbitrary $\lambda$, an asymptotic analysis establishes that the phase of the exciton wave function has the asymptotic series

[TABLE]

where $a_{i,j}$ and $b_{ij}$ are coefficients, cf. (11). Consequently, the wave function behaves as

[TABLE]

with $A_{ij}$ and $B_{ij}$ coefficients. At large distances, the phase acquires the following form

[TABLE]

where $c_{k}$’s are coefficients. In contrast with (13), for the Keldysh-Rytova it is possible to find the structure exciton wave function at large distances. To the best of the authors’ knowledge, this piece of information is absent in the literature, and therefore it has not been exploited for the construction of trial functions, see for example [14]. Certainly, it will lead to an improvement in terms of local accuracy of wave function and variational energies compared with the current trial functions based on hydrogen-like orbitals, see [11]. In our present approach, the next step to construct a compact wave function is the interpolation between series (53) and (54). A minimal interpolation results in

[TABLE]

where $\{A,B,C,D,F,G\}$ are free parameters. Once introduced in (8), it leads to the approximate form of the exciton wave function. Its accuracy in the framework of the variational method will be studied elsewhere. Finally, the three-particle compact wave functions can be constructed following the same procedure shown in this work: going from two to three particles using as building block the exciton wave function.

IV Conclusions

In this article, locally accurate wave functions were constructed for the bound states of two-particle system (exciton) in two-dimensions whose constituent particles interact through a logarithmic potential. For states with quantum numbers $n_{\rho}\leq 4$ and $|m|\leq 5$, they allowed us to reproduce energies with 5-6 exact decimal digits. It was checked using the non-linearization procedure and an alternative two-step variational approach. For magnetoexcitons at rest, it was demonstrated that those functions (used as a zero-order approximation) lead to highly accurate coefficients of the weak coupling energy expansion in powers of the magnetic field strength. Numerical results for $S$-states, suggest that the weak coupling expansion is Padé (re)summable.

We constructed a compact approximation of the two-dimensional three-particle system (trion) ground state wave function using as a building block the exciton ground state one. It only depends on three non-linear free parameters that are fixed through the variational method. It contrasts with the $10^{6}$ linear variational parameters that were recently used in [42]. Our wave function led to binding energies in good agreement with experimental results for trions in concrete TMDCs made of Mo, W, S, and Se. It was shown that the logarithmic potential leads to binding energies $\lesssim 30\%$ larger than those coming from experiments. We find a simple formula for the binding energy as a function of the mass ratio of the constituent particles.

Finally, the structure of the exciton wave function at small distances whose carriers interact via the Rytova-Keldysh potential was established. This new piece of information may lead to the construction of improved variational wave functions used to study larger complexes in TMDCs like bi-excitons.

Acknowledgments

The authors thank Prof. A.V. Turbiner for drawing our attention to logarithmic interactions as well as valuable comments and suggestions. We thank J.C. López-Vieyra for the support with the variational calculations. J.C. del Valle thanks D.A Bonilla-Moreno for useful remarks and discussions. The authors thank M. Szyniszewski for his interest and for providing additional information on Ref. [45]. During the last stage of this work, J.C. del Valle was partially supported by the SONATABIS-10 grant no. 2019/34/E/ST1/00390. D.J.N acknowledges support by Fulbright COMEXUS Project NO. P000003405.