# Theory of third-order polarizability of interlayer excitons due to intra-
excitonic energy level transitions

M. F. C. Martins Quintela1,2, J. C. G. Henriques1,2, N. M. R. Peres1,2
1Department and Centre of Physics, University of Minho, Campus of Gualtar,
4710-057, Braga, Portugal  2International Iberian Nanotechnology Laboratory
(INL), Av. Mestre José Veiga, 4715-330, Braga, Portugal

###### Abstract

In this paper, we employ a fully microscopic approach to the study of
interlayer excitons in layered materials. We discuss the utilization of
Fowler’s and Karplus’ method to access the dynamical polarizability of
non–interacting interlayer excitons in a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$–based van der Waals heterostructure.
Following from the calculation of the linear polarizability, we consider
Svendsen’s variational method to the dynamic third–order polarizability. With
this variational method, we study both two–photon absorption and
third–harmonic generation processes for interlayer excitons in a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ hetero–bilayer, discussing the various
intra–excitonic energy level transitions observed.

##  I Introduction

The advent of the study of two–dimensional layers with atomic thickness has
renewed the interest in transition–metal dichalcogenides (TMDs) and their
optical/optoelectronic properties [1]. These materials, which have been
studied in their bulk form since the 1960’s [2, 3], have been shown to be good
candidates for various optical and optoelectronic applications [4, 5, 6, 7].

The optical response of these materials is dominated by excitons [8], and
optically bright exciton absorption peaks have been shown to correspond to the
excitation of states in the $ns$ series [9, 10]. As for $np$–series excitons,
these can be controlled magnetically despite being optically dark in TMDs
[11]. The selection rules of excitons in TMDs for absorption experiments have
been recently thoroughly studied[12].

The linear dielectric response of TMDs consists of two distinct regimes: the
interband regime, where electrons from the valence band are excited to the
conduction band, leaving being a hole [13, 6]. The attractive electrostatic
interaction between this newly–formed electron–hole pair leads to the
formation of a bound state (exciton), and its $ns$ states are optically active
and can be observed in absorption measurements [14]. The second regime,
characterized by intra-exciton transitions, consists on the transition between
the excitonic ground state ($1s$) and the empty $np$ states [15, 16]. Each of
the $1s\rightarrow np$ transitions [17, 18, 19] is characterized by a peak in
the dynamical polarizability which, in turn, determines the dielectric
response of the system in a pump–probe experiment [20, 21].

Interlayer excitons are formed when two monolayers are brought together with
type–II band alignment, where the conduction band minimum and the valence band
maximum are located in different layers [22, 23]. This allows the intralayer
electron–hole pairs to tunnel into interlayer excitons (Fig. 1–left), which
have a significantly longer life–time due to the small overlap of the
individual electron and hole wave functions [24, 25]. Additionally, interlayer
excitons exhibit luminescence at lower energies than their intralayer
counterpart [26] which, together with their lower binding energies, allows for
an easier identification of the specific species in question in polarizability
measurements [21]. Recently, Merkl _et al._[21] were able to observe the
transition between interlayer and intralayer excitonic phases in van der Waals
heterostructures. This observation was performed via the measurement of the
linear dielectric function of the excitons in a pump–probe experiment (Fig.
1–right), observing a significant shift in maximum of the optical
conductivity.

Figure 1: Left: Interlayer exciton in a van der Waals heterostructure. Right:
Diagram of pump-probe experiments for observing $np$ (dark) excitonic states.

The excitonic energy levels are known to depend on whether the system is
confined or not (for example, the Rydberg series for the two–dimensional
Hydrogen atom is strongly affected by radial confinement [27]). Following
recent advancements in the fabrication of TMD–based meta–materials with atomic
precision [28, 29], we focus our attention on radially confined interlayer
excitons in a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ bilayer.

From studying the dynamic polarizability of the excitonic states, the full
dielectric response of the excited sample can be extracted. This transition
from polarizability to dielectric response is one from the microscopic
(individual exciton) to the macroscopic (when many excitons form a low–density
exciton gas). The susceptibility of this exciton gas would be given, in a
simplified form, as
$\chi_{\mathrm{ex}}^{\mathrm{inter}}=\frac{N_{X}\alpha}{\epsilon_{0}}$, with
$N_{X}$ the exciton density in the two–dimensional material, $\alpha$ the
linear polarizability, and $\epsilon_{0}$ the vacuum dielectric constant.
Separating the real and imaginary part of the susceptibility, these can be
experimentally measured [30, 20, 31, 21] and, therefore, would provide a good
comparison point for the results we obtain for the linear polarizability.

The non–linear polarizabilities of the hydrogen atom have been thoroughly
studied throughout the years [32, 33], providing a good starting point for the
study of these same non–linear processes in excitons. Recently, the non–linear
optical effects from intra–excitonic transitions have been the focus of many
theoretical and experimental works. This has been performed for several TMDs,
including but not limited to $\mathrm{WSe}_{2}$ [34, 35] and $\mathrm{WS}_{2}$
[36, 37], as well as $\mathrm{hBN}$ homostructures [38]. Through doping
processes, the optical non–linear properties of TMDs can be tuned [39] which,
along with their large non–linear optical coefficients [12], increases their
feasibility for practical aplications [40] (e.g., as optical modulators).

We begin this paper by discussing the necessary modifications to the
Rytova–Keldysh potential [41, 42] when dealing with interlayer excitons,
followed by analyzing an harmonic oscillator approximation to this potential
which will motivate the chosen variational wave–function. After defining this
_ansatz_ , we outline the complete set of basis functions we will use to
approximate the excited states of the exciton. The necessary
material–dependent set of parameters characteristic of the specific van der
Waals heterostructure in question will also be discussed, obtaining the
desired energy levels.

In Section III we outline Fowler’s and Karplus’ method [43, 44] to compute the
dynamical polarizability of non–interacting excitons. We then apply this same
method to interlayer excitons in a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$–based
van der Waals heterostructure, comparing the obtained resonances against the
transitions calculated both numerically and variationally via a finite–basis
approach.

In Section IV, we turn to the non–linear response of the excitons [45, 46]. We
focus our discussion on the third–order polarizability, first outlining a
variational procedure based on the ideas of both Karplus and Svendsen [43, 47,
48] for the calculation of the intra–excitonic third–order polarizability.
This procedure was recently studied by Henriques _et al._ , applied to for the
study of the two–photon absorption for excitons in $\mathrm{WSe}_{2}$ [49].
After outlining this variational method, we look at the two–photon absorption
for excitons in $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$, discussing the various
intra–excitonic transitions observed. Finally, we consider the third–harmonic
generation, first with an extremely small broadening as to clearly
differentiate each individual resonance, and afterwards with two much larger
broadenings as to ascertain the feasibility of experimental detection of each
peak. A diagrammatic representation of the various transitions is also
presented for both the two–photon absorption and third–harmonic generation
processes.

We finish the paper with our closing remarks in Section V.

##  II Interlayer Rytova–Keldysh Potential

We begin this section by discussing the Hamiltonian in which we focus our
attention. This will lead to the discussion of the Rytova–Keldysh potential,
and of the necessary modifications when dealing with interlayer excitons.
Following from this discussion, we will approximate the interlayer potential
as an harmonic oscillator, which will allow us to obtain an analytical
solution that will serve as a starting point for the variational _ansatz_ of
the ground state we will use throughout this article. Considering the boundary
conditions of the circular dot in question, we will then define a complete set
of functions that will serve as the basis for our variational approach to the
excitonic states (starting by the $np$–series states for Section III and then
the $ns$– and $nd$–series for Section IV). Finally, we finish this section by
analyzing the material dependent parameters required to study this system,
comparing the various variational energies for the excitonic ground state
against the result from numerical integration of the Schrödinger equation in a
$\log$–grid [50, 51, 52] via the Numerov shooting method.

###  II.1 Model Hamiltonian in the Dipole Approximation

Let us now begin by considering the following Hamiltonian (in atomic units, as
will be used throughout this paper) in the dipole approximation

$\displaystyle H$
$\displaystyle=H_{0}-\mathbf{r}\cdot\mathbf{F}\left(t\right)$
$\displaystyle=-\frac{1}{2\mu}\nabla^{2}+V\left(r\right)-\mathbf{r}\cdot\mathbf{F}\left(t\right),$
(1)

where $\mu$ is the reduced mass of the electron–hole system, $\nabla^{2}$ is
the Laplacian operator (taken in polar coordinates throughout this paper),
$V\left(r\right)$ is a potential energy term and $\mathbf{F}\left(t\right)$ is
an external time–dependent field. This external field will be initially
ignored as to first outline the variational methods and wave functions that
will be used, and its action will be studied in Sections III and IV, when the
methods to obtain both the linear and the third–order dynamical polarizability
are described, respectively.

Considering this same electron–hole system in a van der Waals (vdW)
heterostructure, their potential energy is accurately modeled by the
Rytova–Keldysh potential [41, 42]

$V_{RK}\left(r\right)=-\frac{\pi}{2r_{0}}\left[\mathbf{H}_{0}\left(\kappa\frac{r}{r_{0}}\right)-\mathbf{Y}_{0}\left(\kappa\frac{r}{r_{0}}\right)\right],$
(2)

with $\kappa$ the mean dielectric constant of the media, $r_{0}$ an intrinsic
parameter of the 2D material (interpretable as an in–plane screening length),
and $\mathbf{H}_{0},\mathbf{Y}_{0}$ the Struve–H and Bessel–Y (second kind)
special functions of zero–th order, respectively.

When considering interlayer excitons, a minimum separation originating from
the physical distance between the two layers appears in the Rytova–Keldysh
potential. The interlayer modified Rytova–Keldysh potential then reads [53]

$\displaystyle V_{\mathrm{i-RK}}\left(r\right)$
$\displaystyle=-\frac{\pi}{2r_{0}}\left[\mathbf{H}_{0}\left(\kappa\frac{\sqrt{r^{2}+d^{2}}}{r_{0}}\right)-\right.$
$\displaystyle\left.\qquad-\mathbf{Y}_{0}\left(\kappa\frac{\sqrt{r^{2}+d^{2}}}{r_{0}}\right)\right],$
(3)

with $d$ this interlayer separation distance.

Presenting the same Coulomb tail behaviour as the Rytova–Keldysh potential
[Eq. (2)] at large distances, the finite interlayer separation $d$ eliminates
the logarithmic divergence at the origin. The absence of this divergence makes
it so that an inverse–exponential _ansatz_ for the excitonic wave function
[54] is not the most adequate. This, together with the parabolic nature of the
interlayer modified Rytova–Keldysh potential near $r=0$, motivates the
_ansatz_ that will be obtained in the following section.

###  II.2 Harmonic Oscillator Approximation: Gaussian–Based _Ansatz_

For a first approach to the interlayer modified Rytova–Keldysh potential,
given by Eq. (3), one can expand this potential as a power–series near $r=0$.
This will allow us to obtain an analytical solution in this regime which, in
turn, will motivate the considered variational _ansatz_.

Assuming $d\gg\sqrt{\left\langle r^{2}\right\rangle}$ (where the average of
the radius squared is taken with the exciton wave function), Eq. (3) can be
expanded up to second order in $r$ as

$V\left(r\right)\approx-V_{0}+\gamma r^{2},$ (4)

where

$\displaystyle V_{0}$
$\displaystyle=\frac{\pi}{2r_{0}}\left[\mathbf{H}_{0}\left(\kappa\frac{d}{r_{0}}\right)-\mathbf{Y}_{0}\left(\kappa\frac{d}{r_{0}}\right)\right],$
$\displaystyle\gamma$
$\displaystyle=-\frac{\pi\kappa}{4dr_{0}^{2}}\left[\mathbf{H}_{1}\left(\kappa\frac{d}{r_{0}}\right)-\mathbf{Y}_{-1}\left(\kappa\frac{d}{r_{0}}\right)\right].$
(5)

Considering this potential, the Hamiltonian will be

$H=-\frac{\hbar^{2}}{2\mu}\nabla^{2}-V_{0}+\gamma r^{2}$

and its solution is given by [55, 56, 57]

$\Psi\left(\eta,\varphi\right)=A_{n,l}\left(\frac{\hbar}{\mu\omega}\right)^{\left|l\right|/2}\eta^{\left|l\right|/2}e^{-\eta/2}L_{n}^{\left|l\right|}\left(\eta\right)\frac{e^{il\varphi}}{\sqrt{2\pi}},$
(6)

where $\eta=\mu\omega r^{2}/\hbar$, $\omega=\sqrt{2\gamma/\mu}$,
$n=\left\\{0,1,2,3,\ldots\right\\}$ and $l=\left\\{0,\pm 1,\pm
2,\ldots\right\\}$, $A_{n,l}$ is a normalization constant, and
$L_{a}^{b}\left(x\right)$ is the associate Laguerre polynomial of degree $b$
and order $a$. The energy levels of this wave function are given by

$E_{n}=\hbar\left(\frac{2\gamma}{\mu}\right)^{1/2}\left(2n+\left|l\right|+1\right)-V_{0}.$

For the ground-state $\left(n=0,\right.$ $\left.l=0\right)$, Eq. (6) reads
(discarding numerical factors)

$\Psi\left(\eta\right)\propto e^{-\eta/2}.$

As such, an _ansatz_ based on this solution is given by

$\psi\left(r\right)=\mathcal{C}e^{-\frac{r^{2}}{2\beta}},$ (7)

where $\mathcal{C}$ is a normalization constant and $\beta$ is a variational
parameter. As we are interested in studying interlayer excitons confined to a
finite radius, a multiplicative factor is added to impose boundary conditions,

$\psi_{0}\left(r\right)=\mathcal{C}e^{-\frac{r^{2}}{2\beta}}\left(R-r\right),$
(8)

where $R$ is the radius of the enclosure, considered $R=1200$ throughout this
paper (in atomic units). This radius is considered sufficiently large such
that the variational wave functions become zero significantly before the
boundary is reached (clearly seen in Fig. 2, where the wave functions are
already zero at around $\sim R/3$).

Similarly to the _ansatz_ defined by Pedersen in [54], we modify Eq. (8) into
a sum of two gaussians, given by

$\psi_{0}\left(r\right)=\mathcal{C}\left(e^{-\frac{r^{2}}{2a}}+be^{-\frac{r^{2}}{2c}}\right)\left(R-r\right),$
(9)

which will the be _ansatz_ of the excitonic ground state we will consider
throughout this article. Additionally, a visual comparison of Eq. (8) against
Eq. (9), as well as the approximate wave function obtained via the method
described in Sec. II.3, is present in Fig. 2. Clearly the _ansatz_ (9)
performs much better than the _ansatz_ (8).

Figure 2: Radial profile of the Double–Gaussian Ansatz [Eq. (9)], the Gaussian
Ansatz [Eq. (8)] and the Finite–Basis Approximation (Sec. II.3) of the
ground–state excitonic wave function.

###  II.3 Variational Approach from Boundary Conditions

In certain conditions, namely for a small interlayer separation $d$, the
obtained eigenvalues from minimization of the Hamiltonian

$H_{0}=-\frac{1}{2\mu}\nabla^{2}+V_{\mathrm{i-RK}}\left(r\right)$ (10)

with the _ansatz_ from Eq. (9) are not significantly different from those
obtained from numerical integration (a small difference is present for the
present case in Table 2). Additionally, obtaining higher energy states
involves orthogonalization against Eqs. (8) or (9), followed by normalization,
a process whose complexity increases substantially when increasing the number
of excited states.

A different approach to obtaining variational solutions to the Schrödinger
equation with the interlayer modified Rytova–Keldysh is by integrating the
Hamiltonian matrix elements for a complete basis whose elements obey the
necessary boundary conditions. Truncating this basis and diagonalizing the
Hamiltonian matrix, we obtain a set of variational wave functions that
converge towards the real wave functions of the system as the basis size
increases.

Following [58], a basis of functions can be constructed from the solutions of
the circular infinite well, given by

$\psi_{n,l}\left(r,\varphi\right)=\frac{C_{n,l}}{\sqrt{2\pi}}e^{il\varphi}J_{l}\left(z_{n,l}\frac{r}{R}\right),$
(11)

where $J_{l}\left(x\right)$ is the Bessel function of the first kind of order
$l$, and $z_{n,l}$ is the $n$-th zero of $J_{l}\left(z\right)$. $C_{n,l}$ is a
normalization constant given by

$C_{n,l}=\sqrt{\frac{2}{R^{2}J_{\left|l\right|+1}^{2}\left(z_{n,l}\right)}}$

and the eigenvalues of the infinite circular well are

$E_{n,l}=\frac{\hbar}{2\mu}\left(\frac{z_{n,\left|l\right|}}{R}\right)^{2}.$

As both $V_{RK}\left(r\right)$ and $V_{i-RK}\left(r\right)$ are invariant
under rotations, the quantum number $l$ is well-defined. As such, we can
integrate the Hamiltonian and diagonalize it for a finite number of functions
with fixed $l$. Knowing this, we will start the discussion by outlining the
necessary material dependent parameters for a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ van der Waals heterostructure.

###  II.4 Material Dependent Parameters in $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$

To perform the necessary calculations and obtain the wave functions and energy
eigenvalues, we must substitute the material–specific parameters
characteristic of a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ van der Waals (vdW)
heterostructure [21].

The reduced mass of interlayer excitons in this vdW heterostructure is
$\mu=0.15$ and, as such, the effective Bohr radius will be $a_{0}\approx
0.581\,\text{\AA}$. The thickness of each layer is considered
$d_{\mathrm{mat}}=5.7\,\text{\AA}$ and the gap distance is
$d_{\mathrm{gap}}=1\,\text{\AA}$. As such, the effective $d$ parameter will be

$d_{\mathrm{eff}}=\frac{1}{2}d_{\mathrm{mat}}+d_{\mathrm{gap}}+\frac{1}{2}d_{\mathrm{mat}}=6.7\,\text{\AA}.$

The permittivities for each material are given by

$\epsilon_{\left(n\right)}=\sqrt{\epsilon_{\left(n\right)}^{\bot}\epsilon_{\left(n\right)}^{\|}},$
(12)

where $\epsilon_{\left(n\right)}^{\|}$ and $\epsilon_{\left(n\right)}^{\bot}$
are the in- and out-of-plane relative permittivities of the material $n$. The
numerical values of each component of the relativie permittivities in each
material is present in Table 1.

| $\mathrm{WSe_{2}}$ | $\mathrm{WS_{2}}$
---|---|---
$\epsilon_{\left(n\right)}^{\bot}$ | $7.5$ | $6.3$
$\epsilon_{\left(n\right)}^{\|}$ | $13.36$ | $11.75$
$\epsilon_{\left(n\right)}$ | $10.01$ | $8.604$
Table 1: In–plane, out–of–plane, and total relative permittivities of
$\mathrm{WSe}_{2}$ and $\mathrm{WS}_{2}$.

Knowing the permittivities in each material, the average permittivity of the
vdW heterostructure is

$\kappa_{\mathrm{mat}}=\frac{\epsilon_{\mathrm{WSe_{2}}}+\epsilon_{\mathrm{WS_{2}}}}{2}=9.31.$

Regarding the screening length $r_{0}$, this parameter can be separated as a
sum for each layer (as described in [53])

$r_{0}\rightarrow r_{0}^{\left(1\right)}+r_{0}^{\left(2\right)},$

where $r_{0}^{\left(n\right)}$ denotes the screening length for the layer $n$.
The screening length for each individual layer can be obtained as [20]

$r_{0}^{\left(n\right)}=\frac{\epsilon_{n}^{2}-1}{2\epsilon_{n}}d_{n}\sqrt{\frac{\epsilon_{n}^{\|}}{\epsilon_{n}^{\bot}}},$
(13)

where $d$ is the thickness of the layer. As the in- and out-of-plane relative
permittivities of the two materials are known, the total screening length is

$\displaystyle r_{0}$ $\displaystyle=70.73\,\text{\AA}.$

Choosing a basis size of $120$ and the same enclosure radius as in Sec. II.2
($R=1200$), the ground–state energy eigenvalue is (in both atomic units and
meV)

$E_{1,0}=-0.00143378=-39.0152\,\mathrm{meV}.$

A comparison against both variational _ansatze_ [Eqs. (8) and (9)] and the
numerical results from considering the shooting method in a log–grid is given
in Table 2.

|  $1s$ State
---|---
Finite Basis | $-0.00143378$
Shooting Method | $-0.00144375$
Double–Gaussian [Eq. (9)] | $-0.00142879$
Gaussian [Eq. (8)] | $-0.00136314$
Table 2: Comparison of the different estimates for the excitonic ground–state
energy eigenvalue for a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ heterostructure.

##  III Dynamical Variational Method for Linear Polarizability

In this section we will briefly outline Fowler’s and Karplus’ [43, 44] method
to compute the dynamical polarizability of various systems. This discussion
will be along the same lines as in [59], serving as a quick outline of the
procedure as some results will be necessary further ahead. We apply this
method to calculate the linear polarizability of interlayer excitons in a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ circular dot, comparing the obtained peaks
with the calculated $1s\rightarrow np$ transitions.

Looking back to Eq. (1), we consider (without loss of generality) the external
field pointing along the $-x$ direction. Closely following the approach
presented in [59], first outlined by Fowler and Karplus [43, 44], the
time–dependent Schrödinger equation reads

$\left[H_{0}+xF\left(t\right)\right]\left|\psi\left(t\right)\right\rangle=i\frac{\partial}{\partial
t}\left|\psi\left(t\right)\right\rangle.$ (14)

where $\left|\psi\left(t\right)\right\rangle$ describes the wave function of
the system in the presence of the electric field $\mathbf{F}$. The dynamical
polarizability is defined as

$\alpha\left(\omega\right)=-\left\langle\psi_{0}\right|x\left|\psi_{1}^{+}\right\rangle-\left\langle\psi_{0}\right|x\left|\psi_{1}^{-}\right\rangle$
(15)

such that

$E\left(\omega\right)=E_{0}-\frac{1}{2}\alpha\left(\omega\right)F^{2}.$ (16)

As we are working on a finite disk of radius $R$, Bessel functions of the
first kind are an appropriate complete set of functions to describe a problem
in such a geometry. Using these functions as a basis, we write
$\psi_{1}^{\pm}\left(\mathbf{r}\right)$ as

$\psi_{1}^{\pm}\left(\mathbf{r}\right)=\cos\theta\sum_{n=1}^{N}c_{n}^{\pm}J_{1}\left(\frac{z_{1,n}r}{R}\right),$
(17)

where $J_{1}\left(z\right)$ is the Bessel function of the first kind of order
$1$, $z_{1,n}$ is the $n$–th zero of $J_{1}\left(z\right)$, $N$ is the number
of Bessel functions we choose to use, and $\left\\{c_{n}^{\pm}\right\\}$ are a
set of coefficients yet to be determined. As proposed in Refs. [43, 60, 61,
62], the values of $\left\\{c_{n}^{\pm}\right\\}$ are determined from the
minimization of the functional

$\displaystyle\mathcal{J}_{\pm}$ $\displaystyle=\int
d\mathbf{r}\,\psi_{1,n}^{\pm}\left(\mathbf{r}\right)\left[H_{0}-E_{0}\pm\hbar\omega\right]\psi_{1,n}^{\pm}\left(\mathbf{r}\right)+$
$\displaystyle+2\int
d\mathbf{r}\,\psi_{1,n}^{\pm}\left(\mathbf{r}\right)r\cos\theta\psi_{0}\left(\mathbf{r}\right).$
(18)

Recalling the orthogonality relation of Bessel functions on a finite disk of
radius $R$ [63],

$\displaystyle\int_{0}^{R}J_{\nu}\left(\frac{z_{vm}r}{R}\right)J_{\nu}\left(\frac{z_{vn}r}{R}\right)rdr=$
$\displaystyle=\frac{R^{2}}{2}\delta_{nm}\left[J_{\left|\nu\right|+1}\left(z_{vm}\right)\right]^{2},$
(19)

with $z_{vm}$ the $m-$th zero of $J_{\nu}\left(z\right)$, one easily shows
that the functional can be rewritten as

$\displaystyle\mathcal{J}_{\pm}=\pi\sum_{n=1}\sum_{k=1}c_{n}^{\pm}c_{k}^{\pm}\mathcal{I}_{kn}+2\pi\sum_{n=1}c_{n}^{\pm}\mathcal{S}_{n}+$
$\displaystyle+\frac{\pi
R^{2}}{2}\sum_{n=1}\left(c_{n}^{\pm}\right)^{2}\left[\frac{\hbar^{2}}{2\mu}\frac{z_{1n}^{2}}{R^{2}}-E_{0}\pm\hbar\omega\right]\left[J_{2}\left(z_{1n}\right)\right]^{2},$
(20)

where $\mathcal{I}_{kn}$ and $\mathcal{S}_{n}$ refer to the following
integrals involving one and two Bessel functions, respectively

$\mathcal{I}_{kn}=\int_{0}^{R}J_{1}\left(\frac{z_{1k}r}{R}\right)V_{{\rm
RK}}\left(r\right)J_{1}\left(\frac{z_{1n}r}{R}\right)r\,dr,$ (21)
$\mathcal{S}_{n}=\int_{0}^{R}J_{1}\left(\frac{z_{1n}r}{R}\right)\,r\,\psi_{0}^{{\rm
fin.}}\left(r\right)rdr.$ (22)

Differentiating $\mathcal{J}_{\pm}$ with respect to the different
$c_{n}^{\pm}$ coefficients, one finds

$\displaystyle
c_{j}^{\pm}\left\\{\frac{R^{2}\left[J_{2}\left(z_{1j}\right)\right]^{2}}{2}\left[\frac{z_{1j}^{2}}{2\mu
R^{2}}-E_{0}\pm\omega\right]+\mathcal{I}_{jj}\right\\}+$
$\displaystyle+\sum_{n\neq
j}^{N}c_{n}^{\pm}\mathcal{I}_{jn}=-\mathcal{S}_{j},$ (23)

with $j\in\left\\{1,2,...,N\right\\}$. This equation defines a linear system
of equations whose solution determines the values of the coefficients
$c_{n}^{\pm}$.

We can write Eq. (23) in a more concise manner, using matrix notation, as

$\mathbb{M}\cdot\mathbf{c}^{\pm}=-\mathbf{S},$ (24)

where $\mathbf{c}^{\pm}$ and $\mathbf{S}$ are column vectors defined as

$\displaystyle\left[\mathbf{c}^{\pm}\right]^{{\rm T}}$
$\displaystyle=\left(c_{1}^{\pm},c_{2}^{\pm},\ldots,c_{N}^{\pm}\right),$
$\displaystyle\mathbf{S}^{{\rm T}}$
$\displaystyle=\left(\mathcal{S}_{1},\mathcal{S}_{2},\ldots,\mathcal{S}_{N}\right),$
(25)

and $\mathbb{M}$ is a $N\times N$ matrix with:

$\left(\mathbb{M}\right)_{ij}=g_{i}^{\pm}\left(\omega\right)\delta_{ij}+\mathcal{I}_{ij},$
(26)

where $\delta_{ij}$ is the Kronecker delta and

$g_{i}^{\pm}\left(\omega\right)=\frac{R^{2}\left[J_{2}\left(z_{1i}\right)\right]^{2}}{2}\left[\frac{z_{1i}^{2}}{2\mu
R^{2}}-E_{0}\pm\omega\right].$ (27)

After $\mathbb{M}$ and $\mathbf{S}$ are computed using Eqs. (21) and (27), and
(22), respectively, the coefficients that determine
$\psi_{1}^{\pm}\left(\mathbf{r}\right)$ are readily obtained as

$\mathbf{c}^{\pm}=-\mathbb{M}^{-1}\cdot\mathbf{S},$ (28)

and the solution of the differential equation is found.

At last, to compute the dynamical polarizability we look back at Eq. (17) and
write

$\displaystyle\alpha\left(\omega\right)$
$\displaystyle=-g_{v}\pi\sum_{n=1}^{N}\left(c_{n}^{+}+c_{n}^{-}\right)\times$
$\displaystyle\qquad\times\int_{0}^{R}J_{1}\left(\frac{z_{1n}r}{R}\right)r\,\psi_{0}^{{\rm
fin.}}\left(r\right)r\,dr$
$\displaystyle=-g_{v}\pi\left(\mathbf{c}^{+}+\mathbf{c}^{-}\right)\cdot\mathbf{S}.$
(29)

where $g_{v}=2$ is the valley degeneracy. The $\omega$ dependence on the right
hand side is present in the coefficients $c_{n}^{\pm}$ [more explicitly, this
dependence is present in $\mathbb{M}$ from the
$g_{i}^{\pm}\left(\omega\right)$ term, as defined in Eq. (27)].

The polarizability for interlayer excitons in a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ heterostructure is plotted in Fig. 3,
normalized by its value at the maximum for the excitonic $1s\rightarrow 2p$
transition (with $E_{1s\rightarrow 2p}\sim 29.0576\,\mathrm{meV}$) for easier
visualization.

Figure 3: Real (blue) and imaginary (red) parts of the normalized
polarizability for interlayer excitons in a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$
heterostructure. The vertical black lines represent the excitonic
$1s\rightarrow np$ transitions.

Having outlined and obtained the linear polarizability in Fig. 3, we now turn
our attention to the third–order polarizability of interlayer excitons in this
same vdW heterostructure.

##  IV Svendsen’s Method for the Third–Order Polarizability

Having outlined Fowler’s and Karplus’ method to obtain the linear dynamical
polarizability in Section III, we will now consider a similar method for the
third–order polarizability. We will follow the procedure delineated by
Svendsen [47, 48] and recently applied by Henriques _et al_ [49] to intralayer
excitons in TMDs.

Maintaining the definition for the $\mathbf{c}^{\pm}$ coefficients given in
Eq. (28), and following closely [48, 43, 49], we write the functional

$\displaystyle\mathcal{K}$
$\displaystyle=\left\langle\xi_{\alpha\beta}\left(\omega_{a},\omega_{b}\right)\middle|H_{0}-E_{0}+\omega_{a}+\omega_{b}\middle|\xi_{\alpha\beta}\left(\omega_{a},\omega_{b}\right)\right\rangle+$
$\displaystyle\quad+\left\langle\xi_{\alpha\beta}\left(\omega_{a},\omega_{b}\right)\middle|\mathbf{d}_{\beta}\middle|\psi_{\alpha}\left(\omega_{a}\right)\right\rangle+\left\langle\psi_{\alpha}\left(\omega_{a}\right)\middle|\mathbf{d}_{\beta}\middle|\xi_{\alpha\beta}\left(\omega_{a},\omega_{b}\right)\right\rangle$
(30)

where $\mathbf{d}$ is the dipole operator defined previously, $E_{0}$ the
ground–state energy, $H_{0}$ the unperturbed Hamiltonian [defined in Eq. (10)]
and $\left|\psi_{\alpha}\right\rangle$ the basis obtained via the
$\mathbf{c}^{\pm}$ coefficients and the Bessel functions
$J_{1}\left(z_{1,n}\right)$ in Eq. (17). This functional is then minimized
with respect to $\left|\xi_{\alpha\beta}\right\rangle$, defined analogously to
Eq. (17) as

$\displaystyle\xi_{\alpha\beta}\left(\omega_{a},\omega_{b};r\right)=$
$\displaystyle=\sum_{n=1}^{+\infty}\sum_{\left|l\right|=0}^{+\infty}\zeta_{l,n}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)J_{l}\left(z_{l,n}\frac{r}{R}\right)\frac{e^{il\theta}}{\sqrt{2\pi}}.$

Focusing on the $xx$ component of $\xi_{\alpha\beta}$, and defining
$g_{n}^{\left(l\right)}\left(\omega\right)$ as a more general form of Eq. (27)

$g_{n}^{\left(l\right)}\left(\omega\right)=\frac{R^{2}\left[J_{\left|l\right|+1}\left(z_{l,n}\right)\right]^{2}}{2}\left[\frac{z_{l,n}^{2}}{2\mu
R^{2}}-E_{0}+\omega\right]$

and $\mathcal{V}_{nm}^{\left(l\right)}$ analogously to Eq. (21)

$\displaystyle\mathcal{V}_{nm}^{\left(l\right)}$ $\displaystyle=\left\langle
J_{l}\left(z_{l,n}\frac{r}{R}\right)\middle|V_{i-RK}\left(r\right)\middle|J_{l}\left(z_{l,m}\frac{r}{R}\right)\right\rangle$
$\displaystyle=\int_{0}^{R}J_{l}\left(z_{l,n}\frac{r}{R}\right)V_{i-RK}J_{l}\left(z_{l,m}\frac{r}{R}\right)rdr,$

the functional of Eq. (30) can then be rewritten using the both the
$c^{\pm}_{n}$ and the $\zeta_{l,n}$ coefficients as

$\displaystyle\mathcal{K}$
$\displaystyle=\sum_{n=1}^{+\infty}\sum_{\left|l\right|=0}^{+\infty}\left[\zeta_{l,n}\right]^{\dagger}g_{n}^{\left(l\right)}\left(\omega_{a}+\omega_{b}\right)\zeta_{l,n}+$
$\displaystyle\quad+\sum_{n,m=1}^{+\infty}\sum_{\left|l\right|=0}^{+\infty}\left[\zeta_{l,n}\right]^{\dagger}\mathcal{V}_{nm}^{\left(l\right)}\zeta_{l,m}+$
$\displaystyle\quad+\frac{1}{2}\sum_{n,m=1}^{+\infty}\left\\{\left[\zeta_{0,n}\right]^{\dagger}\mathcal{T}_{nm}^{l=0}\cdot\left(c_{m}^{+}+c_{m}^{-}\right)+\right.$
$\displaystyle\left.+\left[\zeta_{\left|l\right|=2,n}\right]^{\dagger}\left(\mathcal{T}_{nm}^{l=2}\cdot
c_{m}^{+}+\mathcal{T}_{nm}^{l=-2}\cdot c_{m}^{-}\right)+c.c.\right\\}.$ (31)

In Eq. (31), $\mathcal{T}_{nm}^{l}$ is the dipole transition amplitude between
the functions of the basis, denoted by
$\left\langle\psi_{\alpha}\left(\omega_{a}\right)\middle|\mathbf{d}_{\beta}\middle|\xi_{\alpha\beta}\left(\omega_{a},\omega_{b}\right)\right\rangle$
in Eq. (30) and given by

$\displaystyle\mathcal{T}_{nm}^{l}$ $\displaystyle=\left\langle
J_{l}\left(z_{l,n}\frac{r}{R}\right)\middle|\,r\,\middle|J_{1}\left(z_{1,m}\frac{r}{R}\right)\right\rangle$
$\displaystyle=\int_{0}^{R}J_{l}\left(z_{l,n}\frac{r}{R}\right)rJ_{1}\left(z_{1,m}\frac{r}{R}\right)rdr.$

The solution of this integral can be written in a closed form by applying the
more general expression

$\displaystyle\int_{0}^{1}J_{\nu}\left(\alpha r\right)rJ_{\nu+1}\left(\beta
r\right)rdr=\frac{\alpha
J_{\nu+1}\left(\alpha\right)}{\left(\alpha^{2}-\beta^{2}\right)^{2}}\times$
$\displaystyle\qquad\times\left[-2\beta
J_{\nu}\left(\beta\right)+\left(\alpha^{2}-\beta^{2}\right)J_{\nu+1}\left(\beta\right)\right],$
(32)

valid as long as $J_{\nu}\left(\alpha\right)=0$. Comparing with the definition
of $\mathcal{T}_{nm}^{l}$, Eq. (32) can be simplified further as
$J_{\nu+1}\left(\beta\right)=0$. Performing the necessary changes of variable,
and substituting $J_{\nu+1}\left(\beta\right)$ into Eq. (32),
$\mathcal{T}_{nm}^{l}$ can be written as

$\displaystyle\mathcal{T}_{nm}^{\left(l=0\right)}$
$\displaystyle=-2R^{3}\frac{z_{0,n}z_{1,m}J_{1}\left(z_{0,n}\right)J_{0}\left(z_{1,m}\right)}{\left(z_{0,n}^{2}-z_{1,m}^{2}\right)^{2}},$
$\displaystyle\mathcal{T}_{nm}^{\left(l=2\right)}$
$\displaystyle=-2R^{3}\frac{z_{1,n}z_{2,m}J_{2}\left(z_{1,n}\right)J_{1}\left(z_{2,m}\right)}{\left(z_{1,n}^{2}-z_{2,m}^{2}\right)^{2}}.$
(33)

As done in the appendix of [49], we separate the $l=0$ and the $l=2$ cases.
Looking first at the more complicated case of $l=0$, we must ensure that
$\left|\xi_{\alpha\beta}\right\rangle$ is orthogonal to the considered
ground–state [(in our case, the _ansatz_ of Eq. (9)]. As such, to ensure this,
we must have

$\displaystyle\left\langle\psi_{0}\middle|\xi_{\alpha\beta}\right\rangle$
$\displaystyle=\sum_{n=1}^{+\infty}\zeta_{0,n}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)\left\langle\psi_{0}\middle|J_{0}\left(z_{0,n}\frac{r}{R}\right)\right\rangle$
$\displaystyle=0,$ (34)

where non–zero angular momentum terms vanish upon angular integration due to
the isotropic nature of the ground–state _ansatz_. Separating the first term
of the sum, this condition can then be rewritten as

$\displaystyle\zeta_{0,1}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)\left\langle\psi_{0}\middle|J_{0}\left(z_{0,1}\frac{r}{R}\right)\right\rangle+$
$\displaystyle+\sum_{n=2}^{+\infty}\zeta_{0,n}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)\left\langle\psi_{0}\middle|J_{0}\left(z_{0,n}\frac{r}{R}\right)\right\rangle=0,$

which can then be rearranged as

$\displaystyle\zeta_{0,1}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)=$
$\displaystyle=-\sum_{n=2}^{+\infty}\zeta_{0,n}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)\frac{\left\langle\psi_{0}\middle|J_{0}\left(z_{0,n}\frac{r}{R}\right)\right\rangle}{\left\langle\psi_{0}\middle|J_{0}\left(z_{0,1}\frac{r}{R}\right)\right\rangle}.$
(35)

This means that $\zeta_{0,1}$ is no longer considered an independent variable
and, as such, we will focus our attention on the $n\geq 2$ terms. For
compactness, we define

$\displaystyle\mathbf{f}_{n}=\frac{\left\langle\psi_{0}\middle|J_{0}\left(z_{0,n}\frac{r}{R}\right)\right\rangle}{\left\langle\psi_{0}\middle|J_{0}\left(z_{0,1}\frac{r}{R}\right)\right\rangle},$
(36)

which lets us write Eq. (35) as

$\displaystyle\zeta_{0,1}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)=-\sum_{n=2}^{+\infty}\zeta_{0,n}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)\mathbf{f}_{n}.$
(37)

Substituting Eq. (37) into Eq. (31) and then differentiating the resulting
expression with respect to $\left[\zeta_{l,n}\right]^{\dagger}$, we obtain for
$n\geq 2$,

$\displaystyle\left[g_{1}^{\left(0\right)}\left(\omega_{a}+\omega_{b}\right)\mathbf{f}_{n}\mathbf{f}_{m}+\mathcal{V}_{1,1}^{\left(0\right)}\mathbf{f}_{n}\mathbf{f}_{m}\right]\zeta_{0,n}-\mathcal{V}_{1,m}^{\left(0\right)}\mathbf{f}_{n}\zeta_{0,n}-\mathcal{V}_{n,1}^{\left(0\right)}\mathbf{f}_{m}\zeta_{0,n}-\frac{1}{2}\mathbf{f}_{n}\sum_{m=1}^{+\infty}\mathcal{T}_{nm}^{l=0}\cdot\left(c_{m}^{+}+c_{m}^{-}\right)+$
$\displaystyle+\sum_{n=2}^{+\infty}g_{n}^{\left(0\right)}\left(\omega_{a}+\omega_{b}\right)\zeta_{0,n}+\sum_{n,m=2}^{+\infty}\mathcal{V}_{nm}^{\left(0\right)}\zeta_{0,m}+\frac{1}{2}\sum_{n,m=2}^{+\infty}\left\\{\mathcal{T}_{nm}^{l=0}\cdot\left(c_{m}^{+}+c_{m}^{-}\right)+c.c.\right\\}$
$\displaystyle=0.$ (38)

This means that the system can be easily solved algebraically as

$\displaystyle\mathbf{\zeta}_{0}^{x,x}\left(\omega_{a},\omega_{b}\right)=\left[\mathbb{F}+\mathbb{M}^{\left(0\right)}\left(\omega_{a}+\omega_{b}\right)\right]^{-1}\cdot\left[\mathbf{W}_{0}^{x,x}\left(\omega_{a}\right)+\mathbf{f}_{0}^{x,x}\left(\omega_{a}\right)\right],$
(39)

with $\mathbb{M}^{\left(0\right)}\left(\omega_{a}+\omega_{b}\right)$ defined
as in Eq. (25) and

$\mathbf{W}_{0}^{x,x}=-\mathbb{T}^{0}\cdot\mathbf{c},$
$\left(\mathbf{f}_{0}^{x,x}\right)_{n}=\mathbf{f}_{n}\sum_{m=1}^{N}c_{m}\mathcal{T}_{1,m}^{\left(0\right)}c_{m},$
(40) $\displaystyle\left(\mathbb{F}\right)_{i,j}$
$\displaystyle=\left[g_{1}^{\left(0\right)}\left(\omega_{a}+\omega_{b}\right)+\mathcal{V}_{1,1}^{\left(0\right)}\right]\mathbf{f}_{i}\mathbf{f}_{j}-$
$\displaystyle\qquad-\mathcal{V}_{1,j}^{\left(0\right)}\mathbf{f}_{i}-\mathcal{V}_{i,1}^{\left(0\right)}\mathbf{f}_{j}.$

The missing coefficient
$\zeta_{0,1}^{\alpha,\beta}\left(\omega_{a},\omega_{b}\right)$ is obtained by
calculating Eq. (37) with the solution of Eq. (39). Several expressions have
been already simplified by taking advantage of the parity of Bessel–J
functions,

$J_{-l}\left(x\right)=(-1)^{l}J_{l}\left(x\right).$

Looking now at the simpler $l=2$ case, and as no orthogonality–based
restrictions need to be applied to its coefficients, we can directly
differentiate Eq. (31) with respect to $\left[\zeta_{l,n}\right]^{\dagger}$
and then minimize the resulting expresion, obtaining

$\displaystyle\sum_{n=1}^{+\infty}g_{n}^{\left(2\right)}\left(\omega_{a}+\omega_{b}\right)\zeta_{2,n}+\sum_{n,m=1}^{+\infty}\mathcal{V}_{nm}^{\left(2\right)}\zeta_{2,m}+$
$\displaystyle+\frac{1}{2}\sum_{n,m=1}^{+\infty}\left\\{\mathcal{T}_{nm}^{l=2}\cdot
c_{m}^{+}+\mathcal{T}_{nm}^{l=-2}\cdot c_{m}^{-}\right\\}=$ $\displaystyle 0.$
(41)

This expression results in a similar system to the one in Eq. (28), given by

$\displaystyle\mathbf{\zeta}_{2}^{x,x}\left(\omega_{a},\omega_{b}\right)=\left[\mathbb{M}^{\left(2\right)}\left(\omega_{a}+\omega_{b}\right)\right]^{-1}\cdot\mathbf{W}_{2}^{x,x}\left(\omega_{a}\right),$
(42)

with $\mathbf{W}_{2}^{x,x}$ defined analogously to $\mathbf{W}_{0}^{x,x}$,
after simplifications from parity of Bessel functions.

Having defined both $\mathbf{\zeta}_{0}^{x,x}$ and $\mathbf{\zeta}_{2}^{x,x}$,
the third–order susceptibility follows from [49] (with
$\omega_{\sigma}=\omega_{1}+\omega_{2}+\omega_{3}$) as

$\displaystyle\chi_{x,x,x,x}^{\left(3\right)}\left(-\omega_{\sigma};\omega_{1},\omega_{2},\omega_{3}\right)$
$\displaystyle=\frac{1}{3!}\mathcal{P}\left\\{-\mathbf{c}^{\dagger}\left(\left(-\omega_{\sigma}\right)^{\dagger}\right)\cdot\left[\mathbb{T}_{0}\cdot\zeta_{0}^{x,x}\left(-\omega_{2},-\omega_{3}\right)+\mathbb{T}_{2}\cdot\zeta_{2}^{x,x}\left(-\omega_{2},-\omega_{3}\right)\right]\right.+$
$\displaystyle\left.+\left[\mathbf{S}\cdot\mathbf{c}\left(-\omega_{\sigma}\right)\right]\times\left[\sum_{n=1}^{N}c_{n}^{\dagger}\left(\left(-\omega_{2}\right)^{\dagger}\right)c_{n}\left(\omega_{1}\right)2\pi\frac{R^{2}}{2}\left[J_{2}\left(\mathrm{z}_{1,n}\right)\right]^{2}\right]\right\\},$
(43)

with $\mathcal{P}$ denoting the different permutations of the frequencies
$\left(-\omega_{\sigma};\omega_{1},\omega_{2},\omega_{3}\right)$, and
$\mathbf{S}$ as defined in Eq. (22). Having defined Eq. (43), we will now
utilize it to compute two different regimes for the third–order polarizability
of interlayer excitons in $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$.

###  IV.1 Two–Photon Absorption in $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$

Figure 4:
Left: Real (blue) and imaginary (red) parts of the two–photon absorption
third–order susceptibility (normalized by the peak value for the $1\rightarrow
3d$ transition at $\sim 17.5\,\mathrm{meV}$) as a function of the incident
photons’ energies for a broadening $\Delta=0.05\,\mathrm{meV}$ and basis size
$N=120$. The resonances correspond to $1s\rightarrow ns$ transitions (black
dashed lines) and to $1s\rightarrow nd$ transitions (dark–orange dashed
lines). The dominant $1s\rightarrow 2p$ transition is at $\sim
29.0576\,\mathrm{meV}$, outside of the frequency range plotted in this figure.
Right: Diagramatic representation of the $1s\rightarrow ns$ (blue arrows) and
$1s\rightarrow nd$ (red arrows) transitions. Topmost black line represents the
resonance when $\hbar\omega=E_{1s\rightarrow 2p}$.

We will now compute the $xxxx$ component of the two–photon absorption
third–order susceptibility, defined as
$\chi_{xxxx}^{\mathrm{TPA}}\left(\omega\right)=\chi_{xxxx}^{\left(3\right)}\left(-\omega;\omega,-\omega,\omega\right)$,
which will be given by

$\displaystyle\chi_{xxxx}^{\mathrm{TPA}}\left(\omega\right)$
$\displaystyle=\frac{1}{3!}\mathcal{P}\left\\{-\mathbf{c}^{\dagger}\left(\left(-\omega\right)^{\dagger}\right)\cdot\left[\mathbb{T}_{0}\cdot\zeta_{0}^{x,x}\left(\omega,-\omega\right)+\mathbb{T}_{2}\cdot\zeta_{2}^{x,x}\left(\omega,-\omega\right)\right]\right.+$
$\displaystyle\quad\left.+\left[\mathbf{S}\cdot\mathbf{c}\left(-\omega\right)\right]\times\left[\sum_{n=1}^{N}c_{n}^{\dagger}\left(\omega^{\dagger}\right)c_{n}\left(\omega\right)2\pi\frac{R^{2}}{2}\left[J_{2}\left(\mathrm{z}_{1,n}\right)\right]^{2}\right]\right\\}.$
(44)

As the two–photon absorption was considered in [49], we will not discuss it as
in–depth as we will with the third–harmonic generation. Similarly to what was
done in Fig. 3, we plot the real and imaginary parts of the two–photon
absorption third–order susceptibility for a broadening small enough
$\left(0.05\,\mathrm{meV}\right)$ such that each individual peak can be
clearly observed. This plot is visible in Fig. 4, normalized by the peak value
for the $1\rightarrow 3d$ transition at $\sim 16.6\,\mathrm{meV}$, along with
a diagram representing the $1s\rightarrow ns$ and $1s\rightarrow nd$
transitions.

###  IV.2 Third–Harmonic Generation

To finalize this paper, we will now compute the $xxxx$ component of the
third–order susceptibility third–harmonic generation of interlayer excitons in
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$, which is defined as
$\chi_{xxxx}^{\mathrm{THG}}\left(3\omega\right)=\chi_{xxxx}^{\left(3\right)}\left(-3\omega;\omega,\omega,\omega\right)$.
Following the same approach as in Sec. IV.1, we will begin by computing
$\chi_{xxxx}^{\mathrm{THG}}\left(3\omega\right)$ for a small enough broadening
such that each resonance is clear.

As previously, $\chi_{xxxx}^{\mathrm{THG}}\left(3\omega\right)$ can be quickly
obtained by modifying Eq. (43), reading

$\displaystyle\chi_{xxxx}^{\mathrm{THG}}\left(3\omega\right)$
$\displaystyle=\frac{1}{3!}\mathcal{P}\left\\{-\mathbf{c}^{\dagger}\left(\left(-3\omega\right)^{\dagger}\right)\cdot\left[\mathbb{T}_{0}\cdot\zeta_{0}^{x,x}\left(-\omega,-\omega\right)+\mathbb{T}_{2}\cdot\zeta_{2}^{x,x}\left(-\omega,-\omega\right)\right]\right.+$
$\displaystyle\quad\left.+\left[\mathbf{S}\cdot\mathbf{c}\left(-3\omega\right)\right]\times\left[\sum_{n=1}^{N}c_{n}^{\dagger}\left(\left(-\omega\right)^{\dagger}\right)c_{n}\left(\omega\right)2\pi\frac{R^{2}}{2}\left[J_{2}\left(\mathrm{z}_{1,n}\right)\right]^{2}\right]\right\\}.$
(45)

The left–hand side of Fig. 5 has been normalized by the peak value for the
resonance at $3\hbar\omega=E_{1s\rightarrow 2p}$, whilst the right–hand side
features a diagram representing the various transitions.

Figure 5:
Left: Real (blue) and imaginary (red) parts of the third–harmonic generation
third–order susceptibility (normalized by the peak value for the $1\rightarrow
2p$ resonance at $\sim 29.0576,\mathrm{meV}$) as a function of the incident
photons’ energies for a broadening $\Delta=0.05\,\mathrm{meV}$ and basis size
$N=120$. The resonances correspond to $1s\rightarrow np$ transitions (black
dashed lines), $1s\rightarrow ns$ transitions (green dashed lines) and
$1s\rightarrow nd$ transitions (dark–orange dashed lines). The $\sim
87.1727\,\mathrm{meV}$ resonance (black dot–dashed line) corresponds to the
response at frequency $3\hbar\omega=3E_{1s\rightarrow 2p}$.
Right: Diagramatic representation of the various transitions visible in the
plot on the left–hand side. The colors of the arrows represent the different
transitions: $1s\rightarrow np$ (black), $1s\rightarrow ns$ (green), and
$1s\rightarrow nd$ (dark–orange), while the number of arrows represents the
frequency of the transition: 3 arrows for $3\hbar\omega=E_{i\rightarrow f}$, 2
arrows for $2\hbar\omega=E_{i\rightarrow f}$, and 1 arrow for
$\hbar\omega=E_{i\rightarrow f}$.

Spectral broadening as low as $\Delta\sim 2\,\mathrm{meV}$ (i.e., the
broadening that was considered for Fig. 3) can be achieved for low temperature
encapsulated systems [64]. Considering much higher values of the broadening
than in Fig. 5, namely $\Delta=1\,\mathrm{meV}$ and $\Delta=2\,\mathrm{meV}$,
we recompute Eq. (45), obtaining Fig. 6.

Figure 6: Real (left) and imaginary (right) parts of the normalized
third–harmonic generation third–order susceptibility for different broadening
values: $\Delta=2\,\mathrm{meV}$ (solid lines) and $1\,\mathrm{meV}$ (dashed
lines). The basis size is fixed at $N=120$. The resonances correspond to
$1s\rightarrow np$ transitions (black dashed lines), $1s\rightarrow ns$
transitions (green dashed lines) and $1s\rightarrow nd$ transitions
(dark–orange dashed lines). The $\sim 87.1727\,\mathrm{meV}$ resonance (black
dot–dashed line) corresponds to the response at
$3\hbar\omega=3E_{1s\rightarrow 2p}$.

Although some $1s\rightarrow np$ resonances are already apparent at
$\Delta=2\,\mathrm{meV}$, namely at $\sim 40\,\mathrm{meV}$ ($1s\rightarrow
5p$) and $\sim 60\,\mathrm{meV}$ ($1s\rightarrow 8p$), most $1s\rightarrow ns$
and $1s\rightarrow nd$ transitions only become clearly visible in the sub–meV
regime. At those broadenings, the most clear peaks are at $3\hbar\omega\approx
50\,\mathrm{meV}$ ($1s\rightarrow 3d$) and at $3\hbar\omega\approx
56\,\mathrm{meV}$ ($1s\rightarrow 4s$). The two resonances at $\sim
72\,\mathrm{mev}$ and $\sim 75\,\mathrm{mev}$ are joined into one at
$\Delta=1\,\mathrm{meV}$, although their presence is noticeable from
comparison with the adjacent $1s\rightarrow np$ transition peaks.

The $1s\rightarrow np$ resonance peaks occur for
$3\hbar\omega=E_{1s\rightarrow np}$, while the $1s\rightarrow ns$ and
$1s\rightarrow nd$ peaks are situated at
$3\hbar\omega=\frac{3}{2}E_{1s\rightarrow ns/nd}$. Additionally, the peak at
$3\hbar\omega=87.1727\,\mathrm{meV}$ corresponds to
$3\hbar\omega=3E_{1s\rightarrow 2p}$. This is due to the response at both
$2\omega$ [65] and $\omega$ also being present in this process, described in
detail in Section 3.2 of [66], i.e. $2\hbar\omega=E_{1s\rightarrow ns/nd}$ and
$\hbar\omega=E_{1s\rightarrow 2p}$. No other peaks for
$\hbar\omega=E_{1s\rightarrow np}$ are visible, as they would already be
outside the considered frequency domain ($3E_{1s\rightarrow 3p}\sim
105.872\,\mathrm{meV}$, barely outside of the plot).

##  V Conclusions

In this paper, we discussed interlayer excitons in a
$\mathrm{WSe}_{2}/\mathrm{WS}_{2}$–based van der Waals heterostructure, for
which recent experimental measurements of the polarizability through
pump–probe experiments [21] have been performed. In these experiments, the
transition from bright to dark states has been accessed for both interlayer
and intralayer excitons. In this paper, we focused ourselves solely on the
transitions between intra–excitonic states for interlayer excitons in a
circular $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ dot.

Approximating the interlayer modified Rytova–Keldysh potential in a power
series, we motivate a Gaussian variational _ansatz_ for the ground state of
the system. After discussing the necessary material–dependent parameters, we
consider both this _ansatz_ and a complete set of basis functions the Bessel
functions of the first kind to estimate the excitonic wave functions through
variational methods.

We apply Fowler’s and Karplus’ method variational method to access the linear
polarizability of two–dimensional interlayer excitons, focusing on the
transition from bright to dark states (more specifically, $1s\rightarrow np$
transitions) in a circular dot of a $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$–based
van der Waals heterostructure with Dirichlet boundary conditions.

For the linear regime we consider a broadening of $\Delta=2\,\mathrm{meV}$,
allowing for the clear observation of the multiple peaks in the excitonic
polarizability. We observe an almost perfect agreement between the frequency
of each resonance and the energy differences between the ground state and the
excited states of the exciton. These energy differences were calculated both
numerically (via the Numerov shooting method) and variationally (via both a
Gaussian _ansatz_ for the ground–state and a finite basis of Bessel function
of the first kind), allowing for a greater confidence in their values.

Onto the third–order polarizability, we began by providing an in–depth
discussion of Svendsen’s variational method [47, 48]. After outlining this
method and analyzing the required modifications due to the necessity of
orthogonality between the various states, we work on the algebraization of the
method, arriving at a purely vectorial problem after the necessary integrals
are computed. The obtained expression is then tested first with the two–photon
absorption, and then with third–harmonic generation.

Starting by two–photon absorption, a process which was discussed in–depth by
Henriques _et al._ [49], we perform the necessary calculations with a
broadening small enough such that the resonance associated with each
individual intra–excitonic energy level transition can be clearly identifiable
$\left(\Delta=0.05\,\mathrm{meV}\right)$. The two–photon absorption
third–order susceptibility is then normalized by the value at the
$E_{1s\rightarrow 3d}$ resonance as to facilitate the comparison of each peak.
Of note is the amplitude of the resonance associated with the
$E_{1s\rightarrow 7d}$ transition, peaking at around $40\%$ of the maximum.

Finally turning to the third–harmonic generation in the third–order
susceptibility, we begin our discussion by again performing the calculation
with a very small broadening $\left(\Delta=0.05\,\mathrm{meV}\right)$. Each
individual transition was clearly identifiable, with resonances at
$\hbar\omega=E_{i\rightarrow f}$ and $3\hbar\omega=E_{i\rightarrow f}$ for
$1s\rightarrow np$ transitions, and at $2\hbar\omega=E_{i\rightarrow f}$ for
both $1s\rightarrow ns$ and $1s\rightarrow nd$ transitions. Considering
experimentally–obtainable values of the broadening, we compare the
third–harmonic generation third–order susceptibility for
$\Delta=1\,\mathrm{meV}$ and $\Delta=2\,\mathrm{meV}$, arguing for the
possibility of experimentally observing few well resolved resonance peaks,
namely those associated with $1s\rightarrow 3d$ and $1s\rightarrow 4s$
transitions.

## Acknowledgments

M.F.C.M.Q. acknowledges the International Nanotechnology Laboratory for the
Quantum Portugal Initiative grant. N.M.R.P acknowledges support by the
Portuguese Foundation for Science and Technology (FCT) in the framework of the
Strategic Funding UIDB/04650/2020. J.C.G.H. acknowledges the Center of Physics
for a grant funded by the UIDB/04650/2020 strategic project. N.M.R.P. also
acknowledges support from the European Commission through the project
“Graphene–Driven Revolutions in ICT and Beyond” (Ref. No. 881603, CORE 3),
COMPETE 2020, PORTUGAL 2020, FEDER and the FCT through projects
POCI-01-0145-FEDER- 028114, POCI-01-0145-FEDER-02888 and PTDC/NAN-
OPT/29265/2017.

## References

  * [1] N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi, and N. Marzari.  Two-dimensional materials from high-throughput computational exfoliation of experimentally known compounds.  Nat. Nano., 13(3):246–252, feb 2018.
  * [2] R. F. Frindt and A. D. Yoffe.  Phys. properties of layer structures : optical properties and photoconductivity of thin crystals of molybdenum disulphide.  Proc. R. Soc. Lond. A, 273(1352):69–83, apr 1963.
  * [3] E. Fortin and W.M. Sears.  Photovoltaic effect and optical absorption in MoS2.  J. of Phys. and Chem. of Sol., 43(9):881–884, jan 1982.
  * [4] R. Lv, J. A. Robinson, R. E. Schaak, D. Sun, Y. Sun, T. E. Mallouk, and M. Terrones.  Transition Metal Dichalcogenides and Beyond: Synthesis, Properties, and Applications of Single- and Few-Layer Nanosheets.  Acc. of Chem. Res., 48(1):56–64, 2015.
  * [5] T. Mueller and E. Malic.  Exciton physics and device application of two-dimensional transition metal dichalcogenide semiconductors.  npj 2D Mat. and App., 2(1):29, December 2018.
  * [6] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek.  Colloquium: Excitons in atomically thin transition metal dichalcogenides.  Rev. Mod. Phys., 90(2):021001, 2018.
  * [7] Q. Ma, G. Ren, K. Xu, and J. Z. Ou.  Tunable Optical Properties of 2D Materials and Their Applications.  Adv. Opt. Mat., page 2001313, 2020.
  * [8] H. Yu, X. Cui, X. Xu, and W. Yao.  Valley excitons in two-dimensional semiconductors.  Nat. Sci. Rev., 2(1):57–70, 2015.
  * [9] C. Schneider, M. M. Glazov, T. Korn, S. Höfling, and B. Urbaszek.  Two-dimensional semiconductors in the regime of strong light-matter coupling.  Nat. Comm., 9(1):2695, 2018.
  * [10] W.-T. Hsu, J. Quan, C.-Y Wang, L.-S. Lu, M Campbell, W-H Chang, L.-J. Li, X. Li, and C.-K. Shih.  Dielectric impact on exciton binding energy and quasiparticle bandgap in monolayer WS 2 and WSe 2.  2D Mat., 6:025028, 2019.
  * [11] X.-X. Zhang, T Cao, Z Lu, Y.-C. Lin, F. Zhang, Y Wang, Z. Li, J. C. Hone, J. A. Robinson, D. Smirnov, S. G. Louie, and T. F. Heinz.  Magnetic brightening and control of dark excitons in monolayer WSe 2\.  Nat. Nano., 12(9):883–888, 2017.
  * [12] A. Taghizadeh and T. G. Pedersen.  Nonlinear optical selection rules of excitons in monolayer transition metal dichalcogenides.  Phys. Rev. B, 99:235433, Jun 2019.
  * [13] M. V. Durnev and M. M. Glazov.  Excitons and trions in two-dimensional semiconductors based on transition metal dichalcogenides.  Physics-Uspekhi, 61(9):825–845, sep 2018.
  * [14] M. Koperski, M. R. Molas, A. Arora, K. Nogajewski, A. O. Slobodeniuk, C. Faugeras, and M. Potemski.  Optical properties of atomically thin transition metal dichalcogenides: observations and puzzles.  Nanophotonics, 6(6):1289 – 1308, 01 Nov. 2017.
  * [15] G. Berghäuser, A. Knorr, and E. Malic.  Optical fingerprint of dark 2p-states in transition metal dichalcogenides.  2D Mat., 4(1):015029, 2016.
  * [16] K. Miyajima, K. Sakaniwa, and M. Sugawara.  Optical transitions from the lowest to higher exciton and biexciton Rydberg states in CuCl.  Phys. Rev. B, 94(19):195209, 2016.  Publisher: American Physical Society.
  * [17] D. Fröhlich, A. Nöthe, and K. Reimann.  Observation of the Resonant Optical Stark Effect in a Semiconductor.  Phys. Rev. Let., 55(12):1335–1337, 1985.
  * [18] G. Berghäuser, P. Steinleitner, P. Merkl, R. Huber, A. Knorr, and E. Malic.  Mapping of the dark exciton landscape in transition metal dichalcogenides.  Phys. Rev. B, 98(2):020301, 2018.
  * [19] T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos.  Electronic Polarizability as the Fundamental Variable in the Dielectric Properties of Two-Dimensional Materials.  Nano Letters, 20(2):841–851, 2020.
  * [20] C. Poellmann, P. Steinleitner, U. Leierseder, P. Nagler, G. Plechinger, M. Porer, R. Bratschitsch, C. Schüller, T. Korn, and R. Huber.  Resonant internal quantum transitions and femtosecond radiative decay of excitons in monolayer WSe 2.  Nat. Mat., 14(9):889–893, 2015.
  * [21] P. Merkl, F. Mooshammer, P. Steinleitner, A. Girnghuber, K.-Q. Lin, P. Nagler, J. Holler, C. Schüller, J. M. Lupton, T. Korn, S. Ovesen, S. Brem, E. Malic, and R. Huber.  Ultrafast transition between exciton phases in van der waals heterostructures.  Nat. Mat., 18(7):691–696, April 2019.
  * [22] P. Rivera, K. L. Seyler, H. Yu, J. R. Schaibley, J. Yan, D. G. Mandrus, W. Yao, and X. Xu.  Valley-polarized exciton dynamics in a 2d semiconductor heterostructure.  Science, 351(6274):688–691, February 2016.
  * [23] C. Jin, J. Kim, M. I. B. Utama, E. C. Regan, H. Kleemann, H. Cai, Y. Shen, M. J. Shinner, A. Sengupta, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang.  Imaging of pure spin-valley diffusion current in WS2-WSe2heterostructures.  Science, 360(6391):893–896, May 2018.
  * [24] Y. Gong, J. Lin, X. Wang, G. Shi, S. Lei, Z. Lin, X. Zou, G. Ye, R. Vajtai, B. I. Yakobson, H. Terrones, M. Terrones, B.-K. Tay, J. Lou, S. T. Pantelides, Z. Liu, W. Zhou, and P. M. Ajayan.  Vertical and in-plane heterostructures from WS2/MoS2 monolayers.  Nat. Mat., 13(12):1135–1142, September 2014.
  * [25] B. Miller, A. Steinhoff, B. Pano, J. Klein, F. Jahnke, A. Holleitner, and U. Wurstbauer.  Long-lived direct and indirect interlayer excitons in van der Waals heterostructures.  Nano Lett., 17(9):5229–5237, aug 2017.
  * [26] F. Vialla, M. Danovich, D. A. Ruiz-Tijerina, M. Massicotte, P. Schmidt, T. Taniguchi, K. Watanabe, R. J. Hunt, M. Szyniszewski, N. D. Drummond, T. G. Pedersen, V. I. Fal’ko, and F. H. L. Koppens.  Tuning of impurity-bound interlayer complexes in a van der Waals heterobilayer.  2D Mat., 6(3):035032, may 2019.
  * [27] L. Chao-Cador and E. Ley-Koo.  Two-dimensional hydrogen atom confined in circles, angles, and circular sectors.  Int. J. of Q. Chem., 103(4):369–387, 2005.
  * [28] B. Munkhbat, A. B. Yankovich, D. G. Baranov, R. Verre, E. Olsson, and T. O. Shegai.  Transition metal dichalcogenide metamaterials with atomic precision.  Nat. Comm., 11(1):4604, 2020.
  * [29] F. Karimi, S. Soleimanikahnoj, and I. Knezevic.  Tunable plasmon-enhanced second-order optical nonlinearity in transition metal dichalcogenide nanotriangles.  Phys. Rev. B, 103:L161401, 2021.
  * [30] F. Wang, J. Shan, M. A. Islam, I. P. Herman, M. Bonn, and T. F. Heinz.  Exciton polarizability in semiconductor nanocrystals.  Nat. Mat., 5(11):861–864, 2006.
  * [31] S. Cha, J. H. Sung, S. Sim, J. Park, H. Heo, M.-H. Jo, and H. Choi.  1 s -intraexcitonic dynamics in monolayer MoS 2 probed by ultrafast mid-infrared spectroscopy.  Nat. Comm., 7(1):10768, 2016.
  * [32] D. P. Shelton.  Hyperpolarizability of the hydrogen atom.  Phys. Rev. A, 36(7):3032–3041, 1987.
  * [33] N. L. Manakov, S. I. Marmo, and E. A. Pronin.  Dynamic hyperpolarizabilities of excited states of hydrogen.  J. of Exp. and Theo. Phys., 98(2):254–271, February 2004.
  * [34] H. Zeng, G.-B. Liu, J. , Y. Yan, B. Zhu, R. He, L. Xie, S. Xu, X. Chen, W. Yao, and X. Cui.  Optical signature of symmetry variations and spin-valley coupling in atomically thin tungsten dichalcogenides.  Scientific Reports, 3(1), April 2013.
  * [35] H. G. Rosa, Y. W. Ho, I. Verzhbitskiy, M. J. F. L. Rodrigues, T. Taniguchi, K. Watanabe, G. Eda, V. M. Pereira, and J. C. V. Gomes.  Characterization of the second- and third-harmonic optical susceptibilities of atomically thin tungsten diselenide.  Scientific Reports, 8(1), July 2018.
  * [36] C. Janisch, Y. Wang, D. Ma, N. Mehta, A. L. Elías, N. Perea-López, M. Terrones, V. Crespi, and Z. Liu.  Extraordinary second harmonic generation in tungsten disulfide monolayers.  Scientific Reports, 4(1), July 2014.
  * [37] C. Torres-Torres, N. Perea-López, A. L. Elías, H. R. Gutiérrez, D. A. Cullen, A. Berkdemir, F. López-Urías, H. Terrones, and M. Terrones.  Third order nonlinear optical response exhibited by mono- and few-layers of WS2.  2D Mat., 3(2):021005, April 2016.
  * [38] K. Yao, N. R. Finney, J. Zhang, S. L. Moore, L. Xian, N. Tancogne-Dejean, F. Liu, J. Ardelean, X. Xu, D. Halbertal, K. Watanabe, T. Taniguchi, H. Ochoa, A. Asenjo-Garcia, X. Zhu, D. N. Basov, A. Rubio, C. R. Dean, J. Hone, and P. J. Schuck.  Enhanced tunable second harmonic generation from twistable interfaces and vertical superlattices in boron nitride homostructures.  Science Advances, 7(10):eabe8691, March 2021.
  * [39] V. Shahnazaryan, V. K. Kozin, I. A. Shelykh, I. V. Iorsh, and O. Kyriienko.  Tunable optical nonlinearity for transition metal dichalcogenide polaritons dressed by a fermi sea.  Phys. Rev. B, 102:115310, 2020.
  * [40] G. Wang, S. Zhang, X. Zhang, L. Zhang, Y. Cheng, D. Fox, H. Zhang, J. N. Coleman, W. J. Blau, and J. Wang.  Tunable nonlinear refractive index of two-dimensional MoS_2, WS_2, and MoSe_2 nanosheet dispersions [invited].  Phot. Res., 3(2):A51, March 2015.
  * [41] S. N. Rytova.  The screened potential of a point charge in a thin film.  Mosc. Un. Phys. Bul., 22(30), 1967.
  * [42] L. V. Keldysh.  Coulomb interaction in thin semiconductor and semimetal films.  Sov. J. Exp. and Theor. Phys. Lett., 29:658, 1979.
  * [43] M. Karplus and H. J. Kolker.  A Variation-Perturbation Approach to the Interaction of Radiation with Atoms and Molecules.  J. of Chem. Phys., 39(6):1493–1506, 1963.
  * [44] P. W. Fowler.  Energy, polarizability and size of confined one-electron systems.  Mol. Phys., 53(4):865–889, 1984.
  * [45] S. Mossman, R. Lytel, and M. G. Kuzyk.  Dalgarno–lewis perturbation theory for nonlinear optics.  J. of the Opt. Soc. of Am. B, 33(12):E31, September 2016.
  * [46] A. A. Kocherzhenko, S. V. Shedge, X. S. Vazquez, J. Maat, J. Wilmer, A. F. Tillack, L. E. Johnson, and C. M. Isborn.  Unraveling excitonic effects for the first hyperpolarizabilities of chromophore aggregates.  J. of Phys. Chem. C, 123(22):13818–13836, February 2019.
  * [47] E. N. Svendsen.  A variational method for calculation of dynamic third-order susceptibilities.  Int. J. of Qua. Chem., 34(S22):477–480, 1988.
  * [48] J.-R. Gabryl, C. Barbier, P. Lemaire, and E. N. Svendsen.  Variational calculation of the dynamic third-order susceptibility of water.  J. of Mol. Str. Theochem, 357(3):199–205, 1995.
  * [49] J. C. G. Henriques, Høgni C. Kamban, Thomas G. Pedersen, and N. M. R. Peres.  Calculation of the nonlinear response functions of intraexciton transitions in two-dimensional transition metal dichalcogenides.  Phys. Rev. B, 103:235412, Jun 2021.
  * [50] B. R. Johnson.  New numerical methods applied to solving the one-dimensional eigenvalue problem.  J. of Chem. Phys., 67(9):4086–4093, 1977.
  * [51] B. R. Johnson.  The renormalized Numerov method applied to calculating bound states of the coupled-channel Schrödinger equation.  J. of Chem. Phys., 69(10):4678–4688, 1978.
  * [52] G. V. Berghe, V. Fack, and H. E. de Meyer.  Numerical methods for solving radial Schrödinger equations.  J. of Comp. and Ap. Math., 28:391–401, 1989.
  * [53] H. C. Kamban and T. G. Pedersen.  Interlayer excitons in van der waals heterostructures: Binding energy, stark shift, and field-induced dissociation.  Scientific Reports, 10(1), March 2020.
  * [54] T. G. Pedersen.  Exciton stark shift and electroabsorption in monolayer transition-metal dichalcogenides.  Phys. Rev. B, 94:125424, Sep 2016.
  * [55] A. Edery and P. Laporte.  First and second-order relativistic corrections to the two and higher-dimensional isotropic harmonic oscillator obeying the spinless salpeter equation.  J. Phys. Comm., 2(2):025024, feb 2018.
  * [56] B. Cohen-Tannoudji, F. Laloe, and D. Claude.  Quantum Mechanics: volume 3. Wiley, 12 2019.
  * [57] M. Quintela and N. Peres.  A colloquium on the variational method applied to excitons in 2d materials.  Eur. Phys. J. B, 93(12):222, 2020.
  * [58] E. Castaño.  The confined two-dimensional hydrogen atom in the linear variational approach.  Rev. Mex. de Fís., 51:6, 01 2005.
  * [59] J. C. G. Henriques, M. F. C. Martins Quintela, and N. M. R. Peres.  Theoretical model of the polarizability due to transitions between exciton states in transition metal dichalcogenides: application to wse2.  J. Opt. Soc. Am. B, 38(7):2065–2074, Jul 2021.
  * [60] H. E. Montgomery and T. G. Rubenstein.  One-electron wavefunctions. Dynamic dipole polarizabilities.  Chem. Phys. Let., 58(2):295–297, 1978.
  * [61] R. Yaris.  Time-Dependent Perturbation—Variation Method.  J. of Chem. Phys., 39(10):2474–2477, 1963.
  * [62] R. Yaris.  Time-Dependent Perturbation-Variation Method. II.  J. of Chem. Phys., 40(3):667–668, 1964.
  * [63] A. Jeffrey and D. Zwillinger.  Table of Integrals, Series, and Products.  Academic Press, Amsterdam ; Boston, 8th edition edition, 2014.
  * [64] C. Robert, M. A. Semina, F. Cadiz, M. Manca, E. Courtade, T. Taniguchi, K. Watanabe, H. Cai, S. Tongay, B. Lassagne, P. Renucci, T. Amand, X. Marie, M. M. Glazov, and B. Urbaszek.  Optical spectroscopy of excited exciton states in MoS2 monolayers in van der waals heterostructures.  Phys. Rev. Mat., 2(1), January 2018.
  * [65] W. Kim, J. Y. Ahn, J. Oh, J. H. Shim, and S. Ryu.  Second-harmonic young’s interference in atom-thin heterocrystals.  Nano Letters, 20(12):8825–8831, 2020.
  * [66] R. Boyd.  Nonlinear Optics.  Academic Press, London, 2020.


# Hawking radiation from a Vaidya black hole: a semi-classical approach and
beyond

Haryanto M. Siahaana and Triyantab
Theoretical High Energy Physics and Instrumentation Division,
Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung,
Jalan Ganesha 10, Bandung 40132, Indonesia.
a
b<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract:

We derive the Hawking radiation for Vaidya black hole in the tunneling picture
from the corresponding single particle action by the use of the radial null
geodesic and the Hamilton-Jacobi method (beyond semi-classical approximation).
Both results are then analyzed and compared.

Hawking radiation, Vaidya spacetime, semi-classical method and beyond.

††preprint: arXiv:0811.1132 [gr-qc]

##  1 Introduction

In 1974, Hawking startled the physics community by proving that black hole
evaporates particles [1]. It contradicts with classical general relativistic
definition of a black hole, an object that nothing can escape from it [2],
[3]. Hawking’s derivation was a quantum field theoretically. The semi-
classical approach in modeling the black hole radiation as a tunneling effect
was developed in the past decade [4], [5], [6], [7] and attracted many
attentions [e.g. [16]]. This is mostly because of its bold physical picture
and it is easier, mathematically.
There are two ways to perform semi-classical analysis for a black hole
radiation. The first is by the use of radial null geodesic method developed by
Parikh and Wilczek [3], [4]. In this method, one has to get the expression
dr/dt from the radial null geodesic condition, $ds^{2}=d\Omega=0$, for a
metric that has the form
$ds^{2}=-a\left(r\right)dt^{2}+b\left(r\right)dr^{2}+r^{2}d\Omega^{2}$. Then,
the obtained expression is used to calculate the imaginary part of the action.
We then end up with the obtaining of the Hawking temperature
$T_{H}=\left({k\beta}\right)^{-1}$, after equating the transition rate for a
particle in a tunneling process $\Gamma\sim\exp\left[{2{\mathop{\rm
Im}\nolimits}S}\right]$ with the Boltzmann factor $\exp\left[{-\beta
E}\right]$.
The second method was developed by Padmanabhan et al. [5] and it attracted
some applications (e.g. in 2-d stringy black holes [17]). In the method, the
scalar wave function is determined by the ansatz
$\phi\left({r,t}\right)=\exp\left[{{{-iS\left({r,t}\right)}\mathord{\left/{\vphantom{{-iS\left({r,t}\right)}\hbar}}\right.\kern-1.2pt}\hbar}}\right]$
where $S\left({r,t}\right)$ is the action for a single scalar particle.
Inserting this ansatz into the Klein-Gordon equation in a gravitational
background, one yields an equation for the action $S\left({r,t}\right)$ which
can be solved by the Hamilton-Jacobi method. After obtaining the action, one
can get the probability for outgoing and ingoing particles,
$P_{out}=\left|{\phi_{{\rm{out}}}}\right|^{2}$ and
$P_{in}=\left|{\phi_{{\rm{in}}}}\right|^{2}$, respectively. The Hawking
temperature can be obtained by using the ’principle of detailed balance’ [5],
$P_{out}=\exp\left[{-\beta E}\right]P_{in}=\exp\left[{-\beta E}\right]$, since
all particles must be absorbed by the black hole.
An interesting improvement has recently been made by Banerjee and Majhi [8].
They expand the action for single particle in the power series of Planck
constant,
$S\left({r,t}\right)=\sum\nolimits_{n=0}{\alpha_{n}\hbar^{n}S_{n}\left({r,t}\right)}$.
By considering that, for all $n$, $S_{n}\left({r,t}\right)$ and
$S_{n+1}\left({r,t}\right)$ are proportional to each other with the same
proportionality, one can write
$S\left({r,t}\right)=\left({1+\sum\nolimits_{n=1}{\alpha_{n}\hbar^{n}}}\right)S_{0}\left({r,t}\right)$.
By using the method of Padmanabhan et al., the action gives the correction
value of the Hawking temperature. Even though the correction is in principle
negligible, due to the very small of the Planck constant, one could regard
this as an effect of quantum analysis in the semi-classical quantum gravity.
Several developments have been made, including fermionic consideration and the
relation of the method to trace anomaly [15].
In this paper, we consider a more general metric with the mass of black hole
is time ($t$) and radius ($r$) dependent. A semiclassical approach for a
dynamical black hole has also recently worked in [18]. As generally known,
this subject is well described by the Vaidya metric. To get an exact ($t-r$)
dependence of the Hawking temperature, we insist to work in using
Schwarzschild like metric for the Vaidya space time [9] and impose the
condition for a very slowly varying mass of the black hole, $\left|{\partial
m/\partial r}\right|\equiv\left|{m^{\prime}}\right|<<1$ and $\left|{\partial
m/\partial t}\right|\equiv\left|{\dot{m}}\right|<<1$. In this approximation,
the squared $m^{\prime}$ and $\dot{m}$ are well approximated to zero. This
consideration will be useful later in showing the proportionality between each
term of action’s expansion. Our aim is to see the ($t-r$) dependence of
Hawking temperature in the case of varying mass of black hole. We derive it
both in the radial null geodesic and Hamilton-Jacobi method. For Hamilton-
Jacobi method, we directly consider all of the expansion terms of action
(beyond the method of Padmanabhan et al.) as Banerjee and Majhi had done for
some metrics with constant black hole masses.
The organization of our paper is as follows. In the second section, we will
derive the Hawking temperature for a metric with a varying mass by the use of
radial null geodesic method. In the third section, the Hawking temperature is
obtained by the Hamilton Jacobi method, by considering all terms in action
expansion. In the last section, we give a conclusion for our work. For the
rest of this paper, we use the unit dimension: Newton constant, light velocity
in vacuum, and Boltzman constant, $G=c=k_{B}=1$.

##  2 Hawking Temperature from Radial Null Geodesic Method in Vaidya Space-
Time

We start with the metric derived by Farley and D’Eath [9] for Vaidya space-
time

$\displaystyle
ds^{2}=-\left({\frac{{\dot{m}}}{{x\left(m\right)}}}\right)^{2}\left({1-\frac{{2m}}{r}}\right)dt^{2}+\left({1-\frac{{2m}}{r}}\right)^{-1}dr^{2}+r^{2}d\Omega^{2}.$
(1)

In the above, the black hole mass $m$ varies with time $t$ and radius $r$,
$m\equiv m\left({r,t}\right)$. $x\left(m\right)$ is an arbitrary function of
mass and $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ is the metric
of 2-sphere. Later, rather than adopting the law of mass evolution by Hawking
[10], [11], $dm/dt=-C/m^{2}$, we choose a weaker condition, $\partial
m/\partial t\equiv\dot{m}=-C(m)/m^{2}$, which is a partial derivative, in
contrast to the Hawking law that is an exact differentiation. It will be shown
that this choice leads to a better form of Hawking temperature. The variable
$C\left(m\right)$ in this expression counts the number of particles emitted by
the black hole with mass $m$. This variable will increase for decreasing mass
[9] and for early stage of evaporation, we could consider that this variable
is a small valued-function because the evaporated mass is small compared to
the total mass. In this line of thought, putting $x\left(m\right)=-C(m)/m$ is
still valid, at least at the early stage of evaporation. Our purpose in
setting $x\left(m\right)$ in this form is to reduce the variable contained in
the metric and later on, for the case of fixed mass with respect to radius, to
obtain time dependent Hawking temperature [11],
$T_{H}\left(t\right)=\hbar/8\pi m(t)$.
The above considerations lead the metric (1) into the form

$\displaystyle
ds^{2}=-\left({\frac{1}{{m}}}\right)^{2}\left({1-\frac{{2m}}{r}}\right)dt^{2}+\left({1-\frac{{2m}}{r}}\right)^{-1}dr^{2}+r^{2}d\Omega^{2}.$
(2)

It has the general form

$\displaystyle
ds^{2}=-F\left({r,t}\right)dt^{2}+G\left({r,t}\right)^{-1}dr^{2}+r^{2}d\Omega^{2}.$
(3)

where in our case $F\left({r,t}\right)\equiv m^{-2}\left({1-2mr^{-1}}\right)$
and $G\left({r,t}\right)\equiv 1-2mr^{-1}$. The general form (3) will simplify
our next calculation. Unless for specific purposes, we will write
$F\left({r,t}\right)$ and $G\left({r,t}\right)$ as $F$ and $G$, respectively,
for the sake brevity.
It turns out that the metric (2) has a coordinate singularity at $r_{h}=2m$
which of course is time dependent. Painleve transformation which is used to
remove coordinate singularity for a metric with time-like Killing vector, is
also applicable in this analysis. Transforming

$\displaystyle dt\to dt-\sqrt{\frac{{1-G}}{{FG}}}dr$ (4)

the metric (3) changes into

$\displaystyle
ds^{2}=-F\left({r,t}\right)dt^{2}+2F\sqrt{\frac{{1-G}}{{FG}}}dtdr+dr^{2}+r^{2}d\Omega^{2},$
(5)

and therefore no coordinate singularity is found. Thus, by such a Painleve
transformation, it understandable that, in principle, a coordinate singularity
in general relativity can be removed only by changing coordinate, without
defining new physical condition or theorem. In the consideration that a
tunneled particle moves on the path with no singularity, (5) will be useful in
describing particle’s dynamics. For radial null geodesic,
$ds^{2}=d\Omega^{2}=0$, the differentiation of radius with respect to time can
be obtained from (5) as

$\displaystyle\frac{{dr}}{{dt}}=\sqrt{\frac{F}{G}}\left({\pm
1-\sqrt{1-G}}\right),$ (6)

where $+(-)$ signs denote outgoing(ingoing) radial null geodesics.
Near the horizon, we can expand the coefficient $F$ and $G$ by the use of
Taylor expansion. Since $F$ and $G$ are $\left({t-r}\right)$ dependent, and we
only need their approximation values for short distances from a point
(horizon), we could apply the Taylor expansion at a fixed time. So, we can
write

$\displaystyle\left.{F\left({r,t}\right)}\right|_{t}\simeq\left.{F^{\prime}\left({r,t}\right)}\right|_{t}\left({r-r_{h}}\right)+\left.{O\left({\left({r-r_{h}}\right)^{2}}\right)}\right|_{t},$
(7)

and

$\displaystyle\left.{G\left({r,t}\right)}\right|_{t}\simeq\left.{G^{\prime}\left({r_{h},t}\right)}\right|_{t}\left({r-r_{h}}\right)+\left.{O\left({\left({r-r_{h}}\right)^{2}}\right)}\right|_{t}.$
(8)

By the approximations (7) and (8) above, the dependence of radius to time in
(6) can be approached by

$\displaystyle\frac{{dr}}{{dt}}\simeq\frac{1}{2}\sqrt{F^{\prime}\left({r_{h},t}\right)G^{\prime}\left({r_{h},t}\right)}\left({r-r_{h}}\right).$
(9)

Now, we discuss the action of outgoing particle through the horizon. In the
original work by Parikh and Wilczek [3], the imaginary action is written as

$\displaystyle{\mathop{\rm Im}\nolimits}S={\mathop{\rm
Im}\nolimits}\int\limits_{r_{{\rm{in}}}}^{r_{{\rm{out}}}}{p_{r}dr}={\mathop{\rm
Im}\nolimits}\int\limits_{r_{{\rm{in}}}}^{r_{{\rm{out}}}}{\int\limits_{0}^{p_{r}}{dp_{r}^{\prime}dr}}={\mathop{\rm
Im}\nolimits}\int\limits_{r_{{\rm{in}}}}^{r_{{\rm{out}}}}{\int\limits_{0}^{H}{\frac{{dH^{\prime}}}{{{\textstyle{{dr}\over{dt}}}}}}}dr.$
(10)

The above expression is due to the Hamilton equation $dr/dt=dH/dp_{r}|_{r}$
where $r$ and $p_{r}$ are canonical variables (in this case, the radial
component of the radius and the momentum). As a reminder, the action of a
tunneled particle in a potential barrier higher than the energy of the
particle itself will be imaginary, $p_{r}=\sqrt{2m\left({E-V}\right)}$.
Different from discussions of several authors for a static black hole mass
[e.g. [3], [5], [6], [7], [8], [12], [13]], the outgoing particle’s energy
must be time dependent for black holes with varying mass. So, the
$dH^{\prime}$ integration at (10) is for all values of outgoing particle’s
energy, say from zero to $+E\left(t\right)$.
By using approximation (9), we can perform the integration (10). For $dr$
integration, we can perform a contour integration for upper half complex plane
to avoid the coordinate singularity $r_{h}$. The result is

$\displaystyle{\mathop{\rm Im}\nolimits}S=\frac{{2\pi
E\left(t\right)}}{{\sqrt{F^{\prime}\left({r_{h},t}\right)G^{\prime}\left({r_{h},t}\right)}}}.$
(11)

Since the tunneling probability is given by
$\Gamma\sim\exp\left[{-{\textstyle{2\over\hbar}}{\mathop{\rm
Im}\nolimits}S}\right]$, equalizing it with the Boltzmann factor
$\exp\left[{-\beta E\left(t\right)}\right]$ for a system with time dependent
energy we obtain

$\displaystyle
T_{H}=\frac{{\hbar\sqrt{F^{\prime}\left({r_{h},t}\right)G^{\prime}\left({r_{h},t}\right)}}}{{4\pi}}.$
(12)

As seen in expression (12), the Hawking temperature $T_{H}$ is time dependent.
We will see it later that it is also radius dependent. Inserting the values of
$F^{\prime}\left({r_{h},t}\right)$ and $G^{\prime}\left({r_{h},t}\right)$ with
the condition of $\left({m^{\prime}}\right)^{2}\simeq 0$, one has

$\displaystyle
T_{H}=\frac{{\hbar\sqrt{\left({m\left({r_{h},t}\right)}\right)^{2}-2m^{\prime}\left({r_{h},t}\right)\left({m\left({r_{h},t}\right)+2m\left({r_{h},t}\right)}\right)}}}{{8\pi\left({m\left({r_{h},t}\right)}\right)^{2}}}.$
(13)

In the above, the mass $m$ is written as $m\left({r_{h},t}\right)$ explicitly
to remind us that the mass and its derivative
$m^{\prime}\left({r_{h},t}\right)$ are evaluated at $r_{h}$ (the radius of
event horizon). Interestingly, for a black hole whose mass is only time
dependent, not radius dependent, $m^{\prime}=0$, we get

$\displaystyle T_{H}=\frac{\hbar}{{8\pi m\left(t\right)}}.$ (14)

and thus, we recover to the time-dependent Hawking temperature. Further
interesting investigation would be to understand the black hole model with
$x(m)=-C(m)/m$, however it will not be discussed in this paper.

##  3 Hawking Temperature from Hamilton-Jacobi Method in Vaidya Space-Time

In this section, we would work in scalar field theory with gravitational
background. We still work in the metric (1) and impose the condition
$\dot{m}={{-C\left(m\right)}\mathord{\left/{\vphantom{{-C\left(m\right)}{m^{2}}}}\right.\kern-1.2pt}{m^{2}}}$
later. In this section, the condition $\left({\dot{m}}\right)^{2}\simeq 0$ and
$\left({x\left(m\right)}\right)^{2}\simeq 0$ are very important in showing the
proportionality between each expansion terms of action. To avoid confusion
with $F\left({r,t}\right)$ and $G\left({r,t}\right)$ that have been used in
the previous section, we rewrite metric (1) as

$\displaystyle
ds^{2}=-f\left({r,t}\right)dt^{2}+g\left({r,t}\right)^{-1}dr^{2}+r^{2}d\Omega^{2},$
(15)

where $f\left({r,t}\right)\equiv\dot{m}(1-2mr^{-1})/x(m)$ and
$g\left({r,t}\right)\equiv\left({1-2mr^{-1}}\right)$.
Massless scalar particles under the gravitational background $g_{\mu\nu}$ obey
the Klein-Gordon equation

$\displaystyle\frac{{-\hbar^{2}}}{{\sqrt{-g}}}\partial_{\mu}\left[{g^{\mu\nu}\sqrt{-g}\partial_{\nu}}\right]\phi=0.$
(16)

For spherical symmetric black hole, we may reduce our attention only to
($r-t$) sector in the space-time, or in other words, we reduce to two
dimensional black hole problems. In this consideration, equation (16) under
the background metric (15) simplifies to

$\displaystyle\partial_{t}^{2}\phi-\frac{1}{{2fg}}\left({\dot{f}g+\dot{g}f}\right)\partial_{t}\phi-\frac{1}{2}\left({f^{\prime}g+fg^{\prime}}\right)\partial_{r}\phi-
fg\partial_{r}^{2}\phi=0.$ (17)

In the above, $f=f\left({r,t}\right)$ and $g=g\left({r,t}\right)$. By the
standard ansatz for scalar wave function
$\phi\left({r,t}\right)=\exp\left[{-{\textstyle{i\over\hbar}}S\left({r,t}\right)}\right]$,
equation (17) leads to the equation for the action $S\left({r,t}\right)$

$\displaystyle\begin{array}[]{l}\left({\frac{{-i}}{\hbar}\left({\frac{{\partial^{2}S}}{{\partial
t^{2}}}}\right)}\right)-\frac{1}{{\hbar^{2}}}\left({\frac{{\partial
S}}{{\partial
t}}}\right)^{2}-\frac{1}{{2fg}}\left({\dot{f}g+\dot{g}f}\right)\left({\frac{{-i}}{\hbar}}\right)\left({\frac{{\partial
S}}{{\partial t}}}\right)\\\
-\frac{1}{2}\left({f^{\prime}g+fg^{\prime}}\right)\left({\frac{{-i}}{\hbar}}\right)\left({\frac{{\partial
S}}{{\partial
r}}}\right)-fg\left({\frac{{-i}}{\hbar}\left({\frac{{\partial^{2}S}}{{\partial
r^{2}}}}\right)-\frac{1}{{\hbar^{2}}}\left({\frac{{\partial S}}{{\partial
r}}}\right)^{2}}\right)=0\\\ \end{array}.$ (20)

Now, our next step is to solve this equation. An approximation method can be
applied by expanding the action in the order of Planck constant power,

$\displaystyle
S\left({r,t}\right)=S_{0}\left({r,t}\right)+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}S_{n}\left({r,t}\right)},$
(21)

for $n=1,2,3,...$. The constant $\alpha_{n}$ is set to keep all the expansion
terms have the action’s dimension. Taking unit dimensions $G=c=k_{B}=1$,
$\alpha_{n}$ would have the dimension of $[m]^{-2n}$ which $m$ refers to the
mass. It is clear that this expansion would lead to a very long equation. Due
to the very small value of the Planck constant, many authors [5], [6], [12],
[13] neglect the terms for $n\geq 1$. This consideration is acceptable, and
including higher terms is just adding correction for semi-classical derivation
of Hawking temperature. By grouping all the terms into the same powers of
$\hbar$, we could write for some lowest rank as

$\displaystyle\hbar^{0}{\rm{}}:{\rm{}}fg\left({\partial_{r}S_{0}}\right)^{2}-\left({\partial_{t}S_{0}}\right)^{2}=0,$
(22)
$\begin{array}[]{l}\hbar^{1}{\rm{}}:{\rm{}}ifg\partial_{r}^{2}S_{0}-i\left({\partial_{t}^{2}S_{0}}\right)+2fg\left({\partial_{r}S_{0}}\right)\left({\partial_{r}S_{1}}\right)+\frac{i}{2}\left({f^{\prime}g+fg^{\prime}}\right)\partial_{t}S_{0}\\\
{\rm{}}+\frac{i}{2}\frac{{\left({\dot{f}g+\dot{g}f}\right)}}{{fg}}\partial_{t}S_{0}-2\left({\partial_{t}S_{0}}\right)\left({\partial_{t}S_{1}}\right)=0,\\\
\end{array}$
$\begin{array}[]{l}\hbar^{2}{\rm{}}:{\rm{}}ifg\partial_{r}^{2}S_{1}-i\partial_{t}^{2}S_{1}+2fg\left({\partial_{r}S_{0}}\right)\left({\partial_{r}S_{2}}\right)+fg\left({\partial_{r}S_{1}}\right)^{2}+\frac{i}{2}\left({f^{\prime}g+fg^{\prime}}\right)\partial_{t}S_{1}\\\
{\rm{}}+\frac{i}{2}\frac{{\left({\dot{f}g+\dot{g}f}\right)}}{{fg}}\partial_{t}S_{1}-2\left({\partial_{t}S_{0}}\right)\left({\partial_{t}S_{2}}\right)-\left({\partial_{t}S_{1}}\right)^{2}=0,\\\
\end{array}$
$\begin{array}[]{l}\hbar^{3}{\rm{}}:{\rm{}}ifg\partial_{r}^{2}S_{2}-i\partial_{t}^{2}S_{2}+2fg\left({\partial_{r}S_{0}}\right)\left({\partial_{r}S_{3}}\right)++2fg\left({\partial_{r}S_{1}}\right)\left({\partial_{r}S_{2}}\right)-2\left({\partial_{t}S_{0}}\right)\left({\partial_{t}S_{3}}\right)\\\
{\rm{}}-2\left({\partial_{t}S_{1}}\right)\left({\partial_{t}S_{2}}\right)+\frac{i}{2}\left({f^{\prime}g+fg^{\prime}}\right)\partial_{t}S_{2}+\frac{i}{2}\frac{{\left({\dot{f}g+\dot{g}f}\right)}}{{fg}}\partial_{t}S_{2}=0,\\\
\end{array}$ $...$ $...$

As obtained by Banerjee and Majhi [8] for the metric that has a time-like
Killing vector, the metric (15) also leads to such a relation:

$\hbar^{0}{\rm{}}:{\rm{}}\partial_{t}S_{0}=\pm\sqrt{fg}\partial_{r}S_{0},$
$\displaystyle\hbar^{1}{\rm{}}:{\rm{}}\partial_{t}S_{1}=\pm\sqrt{fg}\partial_{r}S_{1},$
(23)
$\hbar^{2}{\rm{}}:{\rm{}}\partial_{t}S_{2}=\pm\sqrt{fg}\partial_{r}S_{2},$
$...$ $...$

To get the benefits of the conditions for black holes with slowly varying
mass, $(m^{\prime})^{2}\simeq 0$ and $(\dot{m})^{2}\simeq 0$, the arbitrary
function $x\left(m\right)$ must be taken to be $-m^{\prime}$ [9]. This
mechanism means that we have used the dynamic equation for black hole mass as
$m^{\prime}=C\left(m\right)/m$. Unless these conditions, the set of equations
(23) can not be obtained. Then one can obtain a general pattern for arbitrary
$n$, thus for the terms of $\hbar^{n}$, that is
$\partial_{t}S_{n}=\pm\sqrt{fg}\partial_{r}S_{n}$. From the set of equations
(22), one can see that each $S_{n}\left({r,t}\right)$ is proportional to
$S_{0}\left({r,t}\right)$. By this evidence, we can write the expansion (21)
into

$\displaystyle
S\left({r,t}\right)=S_{0}\left({r,t}\right)+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}S_{0}\left({r,t}\right)}=\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)S_{0}\left({r,t}\right).$
(24)

The extra value
$\sum\nolimits_{n}{\alpha_{n}\hbar^{n}S_{0}\left({r,t}\right)}$ in (24) can be
regarded as the correction term of the semi-classical analysis. So, our next
step is to find the solution of $S_{0}\left({r,t}\right)$ satisfying
$\partial_{t}S_{0}=\pm\sqrt{fg}\partial_{r}S_{0}$.
In the standard Hamilton-Jacobi method, $S_{0}\left({r,t}\right)$ can be
written into two parts, the time part which has the form of $Et$ and the
radius part $\tilde{S}_{0}\left(r\right)$ which is in general a radius
dependent only. Since our metric coefficients are both radius and time
dependent, the standard method would not be applicable. We could generalized
the method by making an ansatz

$\displaystyle
S_{0}\left({r,t}\right)=\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}+\tilde{S}_{0}\left({r,t}\right).$
(25)

At the first sight, it seems that the ansatz is rather strange, that is
$S\left({r,t}\right)$ and $\tilde{S}_{0}\left({r,t}\right)$ are both $t$ and
$r$ dependent. The term $\int{E\left({t^{\prime}}\right)dt^{\prime}}$ is more
understandable, since the emitted particle’s energy is continuum and time
dependent. Let see how it works.
From (25), one can write that

$\displaystyle\partial_{t}S_{0}\left({r,t}\right)=E\left(t\right)+\partial_{t}\tilde{S}_{0}\left({r,t}\right)$
(26)

and

$\displaystyle\partial_{r}S_{0}\left({r,t}\right)=\partial_{r}\tilde{S}_{0}\left({r,t}\right).$
(27)

Since $\tilde{S}_{0}\left({r,t}\right)$ is $t$ and $r$ dependent, one can
write

$\displaystyle\frac{{d\tilde{S}_{0}\left({r,t}\right)}}{{dr}}=\partial_{r}\tilde{S}_{0}\left({r,t}\right)+\partial_{t}\tilde{S}_{0}\left({r,t}\right)\frac{{dt}}{{dr}}.$
(28)

Eliminating $dt/dr$ by the use of
${{dr}\mathord{\left/{\vphantom{{dr}{dt}}}\right.\kern-1.2pt}{dt}}=\pm\sqrt{fg}$,
equation (28) can be written as

$\displaystyle\frac{{d\tilde{S}_{0}\left({r,t}\right)}}{{dr}}\mp\frac{1}{{\sqrt{fg}}}\partial_{t}\tilde{S}_{0}\left({r,t}\right)=\partial_{r}\tilde{S}_{0}\left({r,t}\right).$
(29)

Combining the first equation of (23), with equations (26) and (29) where we
should note that the action equation
$-\left({fg}\right)^{-1/2}\partial_{t}S_{0}\left({r,t}\right)=\partial_{r}S_{0}\left({r,t}\right)$
is belong to outgoing particle and with the radial evolution is
$dr/dt=\sqrt{fg}$, then we could write

$\displaystyle\mp\left({fg}\right)^{-1/2}\left({E\left(t\right)+\partial_{t}\tilde{S}_{0}\left({r,t}\right)}\right)=\frac{{d\tilde{S}_{0}\left({r,t}\right)}}{{dr}}\mp\left({fg}\right)^{-1/2}\partial_{t}\tilde{S}_{0}\left({r,t}\right).$
(30)

From (30), we can get the exact differentiation of
$\tilde{S}_{0}\left({r,t}\right)$

$\displaystyle\frac{{d\tilde{S}_{0}\left({r,t}\right)}}{{dr}}=\mp\left({fg}\right)^{-1/2}E\left(t\right),$
(31)

and the solution of $\tilde{S}_{0}\left({r,t}\right)$ can be obtained by
integration

$\displaystyle\tilde{S}_{0}\left({r,t}\right)=\mp
E\left(t\right)\int{\frac{{dr}}{{\sqrt{fg}}}}.$ (32)

After inserting the mass evolution equation,
$\dot{m}=-{{C\left(m\right)}\mathord{\left/{\vphantom{{C\left(m\right)}{m^{2}}}}\right.\kern-1.2pt}{m^{2}}}$,
and the arbitrary function
$x\left(m\right)=-{{C\left(m\right)}\mathord{\left/{\vphantom{{C\left(m\right)}m}}\right.\kern-1.2pt}m}$,
one may identify that $f$ and $g$ are exactly equal to $F$ and $G$ stated in
the previous section. By this equality, the integration (32) can be evaluated
by adopting the value of $\int{\left({FG}\right)^{-1/2}dr}$ as in obtaining
expression (11) from (10) along with it’s approximation method (near horizon
Taylor expansion). The result for the integration (32) is

$\displaystyle\tilde{S}_{0}\left({r,t}\right)=\mp
E\left(t\right)\frac{{i\pi}}{{\sqrt{F^{\prime}G^{\prime}}}}.$ (33)

The equation gives the complete action

$\displaystyle
S\left({r,t}\right)=\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)\left({\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}\mp
E\left(t\right)\frac{{i\pi}}{{\sqrt{F^{\prime}G^{\prime}}}}}\right).$ (34)

The signs $+\left(-\right)$ in expression (34) refer to the action for ingoing
(outgoing) particle.
Back to our first ansatz for scalar wave function,
$\phi=\exp\left[{{\textstyle{{-i}\over\hbar}}S\left({r,t}\right)}\right]$, the
wave function for ingoing and outgoing massless scalar particle can be read of
as

$\displaystyle\phi_{in}\left({r,t}\right)=\exp\left[{-\frac{i}{\hbar}\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)\left({\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}+E\left(t\right)\frac{{i\pi}}{{\sqrt{F^{\prime}G^{\prime}}}}}\right)}\right]$
(35)

and

$\displaystyle\phi_{out}\left({r,t}\right)=\exp\left[{-\frac{i}{\hbar}\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)\left({\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}-E\left(t\right)\frac{{i\pi}}{{\sqrt{F^{\prime}G^{\prime}}}}}\right)}\right]$
(36)

respectively. Consequently, from (35) one can get the ingoing probability of
particle as below

$\displaystyle
P_{in}=\exp\left[{\frac{2}{\hbar}\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)\left({{\mathop{\rm
Im}\nolimits}\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}+\frac{{\pi
E\left(t\right)}}{{\sqrt{F^{\prime}G^{\prime}}}}}\right)}\right].$ (37)

This ingoing probability must be equal to unity since all particles including
the massless one are absorbed by the black hole. This consideration gives us
the relation

${\mathop{\rm
Im}\nolimits}\int\limits_{0}^{t}{E\left({t^{\prime}}\right)dt^{\prime}}=-\frac{{\pi
E\left(t\right)}}{{\sqrt{F^{\prime}G^{\prime}}}},$

which leads the outgoing probability

$\displaystyle P_{out}=\exp\left[{-\frac{{4\pi
E\left(t\right)}}{{\hbar\sqrt{F^{\prime}G^{\prime}}}}\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)}\right].$
(38)

Finally, to get the Hawking temperature from the outgoing probability (38), we
equate this probability expression with $\exp\left[{-\beta
E\left(t\right)}\right]$ which in [5] is called ’detailed balance’ principle.
It yields

$\displaystyle
T_{H^{\prime}}=\frac{{\hbar\sqrt{F^{\prime}\left({r_{h},t}\right)G^{\prime}\left({r_{h},t}\right)}}}{{4\pi\left({1+\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}}\right)}}.$
(39)

We use the index $H^{\prime}$ to distinguish the result (39) with that of
(12). Neglecting the correction term $\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}$
in the denominator of (39), one has recovered result (12) which has been
obtained by the radial null geodesics at the previous section.

##  4 Conclusion

We have worked out the Hawking temperature for the metric with no time-like
Killing vector and with the time and radius dependent coefficients. It has
been shown that by the use of mass evolution $\partial m/\partial
t=-C(m)/m^{2}$ and mass dependent function $x(m)=-C(m)/m$, the resulting
Hawking temperature coincides with that for the widely known time dependent
one in the case of radius independent mass. By neglecting the
$\sum\nolimits_{n}{\alpha_{n}\hbar^{n}}$ term in the denominator of the
resulting Hawking temperature for the case of beyond semi-classical approach
as given in Section 3, one recovers the result in the section 2. However, the
equal obtained results in (12) and (38) with neglected correction terms are
actually based on some subtleties. In section 2, we use the $\dot{r}$ that is
obtained from the coordinate after Painleve transformation which differs with
a factor half with untransformed one (directly obtained from Schwarzschild
like metric (3)). This is clear since the particle moves in the path with no
singularity111We thank B. Majhi and R. Banerjee for pointing this out.. Thus,
in section 3 we use the $\dot{r}$ in deriving the action’s solution (29) that
is obtained from Schwarzschild like metric (without Painleve transformation),
since the equation of action comes from it (17). It needs further
investigations to verify whether the results are still the same in the
suitable Painleve coordinate in deriving the wave equation (16) (parallel to
section 3.2 [8]).
The slowly varying mass with respect both to time and radius, considered in
Section 3, affects our result on the Hawking temperature expression and the
proportionality between each terms of action’s expansion. We have restored the
$x(m)$ as $-m^{\prime}$ in section 3 to get the proportionality by using our
slowly varying mass condition. This consideration of course yields an equation
for mass dynamics $m^{\prime}=C\left(m\right)/m$ which needs further
investigations. Of course, our analysis is not valid for the case of highly
radiated black holes. There are several models that have been proposed for
describing black hole radiation. But, the lack of experimental data does not
enable one to compare between the models.
Further investigation on the relationship between ingoing and outgoing
probabilities for the case of black hole with varying mass/energy by the use
of path integral method is intriguing. This work was first performed by Hartle
and Hawking [14] for a black hole with constant mass and constant energy of
outgoing particle. The relation proposed in the paper,
$P_{out}=\exp\left[{-\beta E\left(t\right)}\right]P_{in}$, is an intuitive
manner as it might not be that simple.

##  5 Acknowledgements

HMS thanks R. Kurniadi and A. Nugroho from the Faculty of Mathematics and
Natural Sciences, ITB for useful discussions. HMS also thanks B. Majhi and R.
Banerjee from S. N. Bose National Centre for Basic Sciences for useful
comments.

## References

  * [1] S. W. Hawking (1975), Particle Creation by Black Holes, Comm. Math. Phys. 43, 199.
R. M. Wald (1994), Quantum Field Theory in Curved Spacetime and Black Hole
Thermodynamics, The Univ. of Chicago Press.

  * [2] R. M. Wald (1984), General Relativity, The Univ. of Chicago Press.
  * [3] M. K. Parikh and F. Wilczek (2000), Hawking Radiation as Tunneling, Phys. Rev. Lett. 85, 5042 [arXiv:hep-th/9907001].
  * [4] M. K. Parikh (1998), Membrane Horizons : The Black Hole’s New Clothes, Ph.D. Thesis, Princeton University, Princeton, NJ [arXiv:hep-th/9907002].
  * [5] K. Srinivasan and T. Padmanabhan (1999), Particle Production and Complex Path Analysis, Phys. Rev. D 60, 024007 [arXiv:hep-th/9907002].
S. Shankaranarayanan, T. Padmanabhan, and K. Srinivasan (2002), Hawking
Radiation in Different Coordinate Settings : Complex Path Approach, Class.
Quant. Grav. 19, 2671, [arXiv:gr-qc/0010042].

  * [6] P. Mitra (2006), Hawking Temperature from Tunneling Formalism, Phys. Lett. B 648, 240 [arXiv:hep-th/0611265].
  * [7] E.T. Akhmedov, V. Akhmedova, and D. Singleton (2006), Hawking Temperature in the Tunneling Picture, Phys. Lett. B 642, 124 [arXiv:hep-th/0608098].
E.T. Akhmedov, V. Akhmedova, D. Singleton, and T. Pilling (2007), Thermal
Radiation of Various Gravitational Background, Int. J. Mod. Phys. A 22, 1705,
[arXiv:hep-th/0605137].
T. Pilling (2008), Black Hole Thermodynamics and the Factor 2 Problem, Phys.
Lett. B 660, 402, [arXiv:0709.1624].
B.D. Chowdhury (2008), Pramana 70, 593, [arXiv:hep-th/0605197].
A.J.M. Medved and Elias C. Vagenas (2005), On Hawking Radiation as Tunneling
with Back-Reaction, Mod. Phys. Lett. A 20 , 2449, [arXiv:gr-qc/0504113].
Ryan Kerner and Robert B. Mann (2006), Tunnelling, Temperature and Taub-NUT
Black Holes, Phys. Rev. D 73, 104010, [arXiv:gr-qc/0603019].
Michele Arzano, A.J.M. Medved and E. C. Vagenas (2005), Hawking Radiation as
Tunneling through the Quantum Horizon, JHEP 0509, 037, [arXiv:hep-th/0505266].
Rabin Banerjee and Bibhas R. Majhi (2008), Quantum Tunneling and Back
Reaction, Phys.Lett. B 662, 62, [arXiv:0801.0200].
Rabin Banerjee, Bibhas R. Majhi, and Saurav Samanta (2008), Noncommutative
Black Hole Thermodynamics, Phys.Rev. D 77, 124035, [arXiv:0801.3583].

  * [8] R. Banerjee and Bibhas R. Majhi (2008), Quantum Tunneling Beyond Semiclassical Approximation, JHEP 0806, 095, [arXiv:0805.2220].
  * [9] A. N. St. J. Farley and P. D. D’Eath (2006), Vaidya Space-Time in Black Hole Evaporation, Gen. Rel. Grav. 38, 425, [arXiv:gr-qc/0510040].
  * [10] S.W. Hawking (1974), Nature, 248, 30.
  * [11] S. Massar (1995), The Semi-Classical Back Reaction to Black Hole Evaporation, Phys. Rev. D 52, 5857, [arXiv:gr-qc/9411039].
  * [12] Marco Angheben, Mario Nadalini, Luciano Vanzo, and Sergio Zerbini (2005), Hawking Radiation as Tunneling for Extremal and Rotating Black Holes, JHEP 0505, 014, [arXiv:hep-th/0503081].
M. Nadalini, L. Vanzo and S. Zerbini (2005), Hawking radiation as tunneling:
The D dimensional rotating case, J. Phys. A 39, 6601 [arXiv:hep-th/0511250].

  * [13] E.T. Akhmedov, T. Pilling, and Douglas Singleton, Subtleties in the Quasi-Classical Calculation of Hawking Radiation, [arXiv:0805.2653].
  * [14] J.B. Hartle and S.W. Hawking (1976), Path Integral Derivation of Black Hole Radiance, Phys. Rev. D 13, 2188.
  * [15] Sujoy K. Modak (2008), Corrected Entropy of BTZ Black Hole in Tunneling Approach, [arXiv:0807.0959].
Bibhas R. Majhi (2008), Fermion Tunneling Beyond Semiclassical Approximation,
[arXiv:0809.1508].
Rabin Banerjee and Bibhas R. Majhi (2008), Quantum Tunneling and Trace
Anomaly, [arXiv:0808.3688].

  * [16] E.C. Vagenas (2000), Are extremal 2-D black holes really frozen?, Phys. Lett. B 503, 399 [arXiv:hep-th/0012134].
E.C. Vagenas (2001), Two-dimensional dilatonic black holes and Hawking
radiation, Mod. Phys. Lett. A 17, 609 [arXiv:hep-th/0108147].
E.C. Vagenas (2001), Semiclassical corrections to the Bekenstein-Hawking
entropy of the BTZ black hole via selfgravitation, Phys. Lett. B 533, 302
[arXiv:hep-th/0109108].
A. J. M. Medved and E. C. Vagenas (2005), On Hawking radiation as tunneling
with logarithmic corrections, Mod. Phys. Lett. A 20, 1723 [arXiv:gr-
qc/0505015].
Qing-Quan Jiang and Shuang-Qing Wu (2005), Hawking radiation of charged
particles as tunneling from Reissner-Nordstrom-de Sitter black holes with a
global monopole, Phys. Lett. B 635, 151 [arXiv:hep-th/0511123].
Qing-Quan Jiang, Shuang-Qing Wu and Xu Cai (2005), Hawking radiation as
tunneling from the Kerr and Kerr-Newman black holes, Phys. Rev. D 73 064003
[arXiv:hep-th/0512351].
R. Di Criscienzo, M. Nadalini, L. Vanzo, S. Zerbini and G. Zoccatelli (2007),
On the Hawking radiation as tunneling for a class of dynamical black holes,
Phys. Lett. B 657, 107 [arXiv:0707.4425].

  * [17] E.C. Vagenas (2001), Complex paths and covariance of Hawking radiation in 2-D stringy black holes, Nuovo Cim. B 117, 899 [hep-th/0111047].

  * [18] S.A. Hayward, R. Di Criscienzo, L. Vanzo, M. Nadalini and S. Zerbini (2008), Local Hawking temperature for dynamical black holes, [arXiv:0806.0014].