# Microscopic theory of the polarizability of transition metal dichalcogenides
excitons: Application to WSe2

J. C. G. Henriques1,2, M. F. C. Martins Quintela1, N. M. R. Peres1,2∗
1Department and Centre of Physics, and QuantaLab, 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 develop a fully microscopic theory of the polarizability of
excitons in transition metal dichalcogenides. We apply our method to the
description of the excitation $2$p dark states. These states are not
observable in absorption experiments but can be excited in a pump-probe
experiment. As an example we consider $2$p dark states in WSe2. We find a good
agreement between recent experimental measurements and our theoretical
calculations.

##  I Introduction

The optical properties of monolayer transition metal dichalcogenides (TMDs)
are of considerable interest _per se_ , but also from the point of view of
applications Lv  _et al._ (2015); Ma  _et al._ (2020); Mueller and Malic
(2018); Wang  _et al._ (2018); Mak  _et al._ (2018); Schneider  _et al._
(2018). It is by now well known that the optically bright exciton absorption
peaks correspond to the excitation of states in the $n$s series Schneider  _et
al._ (2018); Hsu  _et al._ (2019), with the 1s being energy split due to the
strong spin-orbit effect present in these systems. The excitons in the $n$p
series are optically dark in TMDs and a very faint d-exciton line has also
been predicted Chaves  _et al._ (2017) but never observed to date. The dark
excitons can, however, be controlled magnetically Zhang  _et al._ (2017). The
Berry phase-induced splitting of the 2p states in MoSe2 was studied in Ref.
Yong  _et al._ (2019).

The dielectric response of TMDs has two different regimes. The first, named
the _interband_ , occurs when an electron in the valence band is promoted to
the conduction band. Due to the attractive electrostatic interaction between
the hole left in the valence band and the electron promoted to the conduction
band, this excitation leads to the formation of excitonic states Steinleitner
_et al._ (2017); Hsu  _et al._ (2019). From those, the $ns$ are optically
active and can be seen in absorption measurements Koperski  _et al._ (2017).
The other possible regime, which we call the _intra-exciton transitions_ ,
consists of transitions between the occupied $1s$ state and the empty $np$
states of the excitonic energy levels Miyajima  _et al._ (2016); Berghäuser
_et al._ (2016). Each of the $1s\rightarrow np$ transitions Berghäuser  _et
al._ (2018); Yong  _et al._ (2019); Hsu  _et al._ (2019); Merkl  _et al._
(2019) is characterized by an in-plane polarizability (a quantity shown
relevant in other semiconductors Tian  _et al._ (2020); Wang  _et al._
(2006)), which in turn determines the dielectric response of the system in a
pump-probe experiment, using lasers of different frequencies. In this work, we
consider the case of the TMD WSe2 with which experiments similar to the ones
just described have been performed Poellmann  _et al._ (2015). In Ref. Cha
_et al._ (2016) MoS2 was studied. Also, valley-selective physics has been
measured in this system Sie  _et al._ (2015); Mak  _et al._ (2018). In Ref.
Poellmann  _et al._ (2015), a 90-fs laser pulse centered at $\lambda_{c}=742$
nm was shined on a WSe2monolayer leading to the population of the $1s$
excitonic state. At a certain variable time delay, the low-energy dielectric
response is probed by a phase-locked mid-infrared pulse.

The first experimental evidence of $1s\rightarrow np$ transitions between
exciton levels of the associated Rydberg series Fröhlich  _et al._ (1985) was
found in the solid state system Cu2O, at the temperature of 1.8 K Fröhlich
_et al._ (1985). More recently, the experiment of Poellmann _et al._ Poellmann
_et al._ (2015) showed that the same physics can be observed at room
temperature in WSe2, a work that sprung interest by other groups Zhang  _et
al._ (2015); Merkl  _et al._ (2019). The results of Poellmann _et al._ are
special due to the two-dimensional nature of this TMD, which entails a poorer
electrostatic screening between the electron and the hole brought about by the
reduced dimension of WSe2 when compared to the three dimensional case of Cu2O.

The field-resolved time-domain data allowed the retrieving of the full
dielectric response of the excited sample Poellmann  _et al._ (2015). The
dielectric response is characterized by by the real and imaginary parts of
dielectric function
$\Delta\epsilon(\omega)=\epsilon_{1}(\omega)+i\epsilon_{2}(\omega)$. The real
part follows a dispersive shape with a zero crossing at a photon energy of 180
meV, for a given exciton density (laser fluence). As the authors stress, this
spectrum is in striking contrast to the Drude-like response of free
electron–hole pairs observed after non-resonant above-band gap excitation, and
is, therefore, a clear signature of the resonant excitation of intra-exciton
transition $1$s$\rightarrow 2$p referred above. The authors have also mapped
the real part of the optical conductivity, $\Delta\sigma_{1}$, which directly
connects to the imaginary part of the dielectric function. The former quantity
shows pronounced resonance at the frequency value of 180 meV. As we shall see
further ahead, these two features can be retrieved from the polarizability
Wang  _et al._ (2006) of the excitonic states.

Figure 1: Pump-probe experiment for observing dark excitonic states in a
transition metal dichalcogenide. The high intensity pump populates the $1$s
ground state of the excitonic energy series, whereas the mid-IR low intensity
probe induces transitions between the excitonic ground state ($1$s) and dark
$n$p states, thus allowing to access this usually not accessible piece the of
the dielectric response of the excitonic gas. This procedure is represented by
the energy levels on the right part of the figure. In the left part of the
figure, a sample of WSe2 is depicted on which both the pump and the probe
laser impinge.

The position of excitonic energy levels are known to depend on whether the
system is confined. For example, it is well known that the Rydberg series of
the two-dimensional Hydrogen atom, which in the infinite systems is given by
$E_{n}=-1/[2(n-1/2)^{2}]$ with $n=1,2,3\ldots$, in atomic units, is strongly
affected by confining the atom in a disk of radius $R$ Chao-Cador and Ley-Koo
(2005). Recently, fabrication of transition metal dichalcogenide metamaterials
with atomic precision Munkhbat  _et al._ (2020) has opened the possibility of
studying the dependence of the excitonic energy series of TMDs upon geometric
confinement of different geometries. It is expected that in addition to
negative energy bound states, positive ones will also appear. As we show
ahead, our procedure may be applied to the study of excitons in circular boxes
(which can actually be obtained, approximately, even in a hexagonal lattice).

In this paper we establish the connection between the polarizability of an
exciton (a microscopic quantity) and the dielectric function, which is a
macroscopic quantity characterizing many particles, when many excitons,
considered non-interacting (and therefore, low density) form an excitonic gas.
We show that our microscopic theory can account well for the experimental
results of Poellmann _et al._ Poellmann  _et al._ (2015). The paper is
organized has follows: in section II we present the Fowler’s and Karplus’
method which will allow us to compute the optical polarizability without
evaluating a sum over states (whose direct approach is doomed to fail). This
method requires the solution of a differential equation, which is then used to
compute the matrix elements that define the polarizability. Afterwards, we
apply this technique to the simple case of the cylindrical potential well. The
goal of this section is to benchmark the method. In section III we explore the
problem of excitons in WSe2, presenting a way of solving the Fowler’s and
Karplus’ differential equation for excitonic problems, computing the
polarizability and, from it, the dielectric function and conductivity of a
macroscopic collection of excitons. This section ends with a comparison of our
theoretical results with experimental data and a good agreement is found. In
section IV we give our final remarks. An appendix showing how to solve the
Wannier equation in log-grid closes the paper.

##  II Karplus’ and Fowler’s method

In this section we will briefly present the Karplus’ and Fowler’s Podolsky
(1928); Karplus (1962); Karplus and Kolker (1963); Fowler (1984) method to
compute the dynamical polarizability of different systems. After the essential
features of this approach are laid out, we will explore the problem of the
infinite cylindrical well. This serves two purposes: on the one hand it allows
us to apply the formalism to a simple yet instructive case; and on the other,
it sets the stage for the problem of two-dimensional excitons that will be
treated in the subsequent section and involves the solution of the Wannier
equation.

###  II.1 Theory

Let us begin our discussion considering the following Hamiltonian written in
atomic units (a.u.; this system of units is used throughout the paper, except
when stated otherwise) in the dipole approximation:

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

where $\mu$ is a mass term, $\nabla^{2}$ is the Laplacian, $V(\mathbf{r})$ is
a potential energy term and $\mathbf{F}(t)$ is an external time dependent
harmonic electric field:

$\mathbf{F}(t)=F\left[e^{-i\omega t}+e^{i\omega t}\right],$ (3)

with $\omega$ the field’s frequency and $F$ its amplitude. This general
Hamiltonian can be used to describe the interaction of different systems with
an external electric field. Two examples which will be treated in this work
are the interaction of an electric field with a particle trapped inside a
cylindrical well, and the interaction of light with excitons in 2D materials.
The first will be treated ahead in the present section as an introductory
problem, while the latter is the main subject of this work and will be studied
in the following section.

Now, for simplicity, we consider the electric field to be applied along the
$-x$ direction. With this assumption, the time dependent SchrÃ¶dinger equation
reads:

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

where $|\psi(t)\rangle$ is the state vector describing the wave function of
the system in the presence of the electric field. Following Fowler’s and
Karplus’ approach, we expand the state vector in powers of $F$ as:

$\displaystyle|\psi(t)\rangle$ $\displaystyle=e^{-iE_{0}t}|\psi_{0}\rangle$
$\displaystyle+Fe^{-i(E_{0}-\omega)t}|\psi_{1}^{+}\rangle+Fe^{-i(E_{0}+\omega)t}|\psi_{1}^{-}\rangle+\mathcal{O}\left(F^{2}\right)$
(5)

where $E_{0}$ and $|\psi_{0}\rangle$ are the energy and state vector,
respectively, of the unperturbed system and $|\psi_{1}^{\pm}\rangle$ are the
first order corrections to the state vector. Inserting this expansion in Eq.
(4), and keeping only terms up to first order in $F$, we find:

$\displaystyle H_{0}|\psi_{0}\rangle$ $\displaystyle=E_{0}|\psi_{0}\rangle$
(6) $\displaystyle H_{0}|\psi_{1}^{\pm}\rangle+x|\psi_{0}\rangle$
$\displaystyle=\left(E_{0}\mp\omega\right)|\psi_{1}^{\pm}\rangle,$ (7)

where the first equation is nothing more than the SchrÃ¶dinger equation of the
unperturbed system, and the second equation defines the first order correction
$|\psi_{1}^{\pm}\rangle.$ The dynamical polarizability $\alpha(\omega)$ is
defined as Karplus and Kolker (1963):

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

such that:

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

We note that no term proportional to $F$ appears in systems with inversion
symmetry (like the ones we will be considering). Note that we have departed
here from the traditional Rayleigh-SchrÃ¶dinger perturbation theory, which
goes a step further and expands $|\psi_{1}^{\pm}\rangle$ in the basis of the
unperturbed Hamiltonian.

As just said, if Rayleigh-SchrÃ¶dinger time-dependent perturbation theory had
been used to compute the dynamical polarizability, the final result would
consist in a sum over states which, in principle, could not be exactly
computed, since it would require the computation of an arbitrary large number
of matrix elements whose form is unknown analytically. Also, a numerical
approach to the sum over states would fail, because of the difficulty of
describing the wave functions of the continuum of states. Usually, these sums
are truncated and only the first terms are considered. The quality of this
approximation depends on how fast the sum over states converges to a given
value Henriques and Peres (2020). In comparison, the presented approach
bypasses the sum over states, and requires the computation of only two matrix
elements in order to obtain the dynamical polarizability, as described by Eq.
(8), and, at the same time, accomplishes the summation over all states.

###  II.2 The infinite cylindrical well

In order to test Eq. (8), we will now explore a simple example: a particle
trapped inside a 2D cylindrical potential well. For simplicity let us consider
the particle to be an electron, such that $\mu=1$ (in atomic units) in Eq.
(2). The potential term reads:

$V(r)=\begin{cases}0,&r<R\\\ \infty,&r>R\end{cases},$ (10)

where $R$ is the radius of the cylindrical potential well.

Before the dynamical polarizability can be computed, the unperturbed system
has to be studied. Although this is a simple system with a well established
solution, we will give, for completeness, a brief description of the necessary
steps to find the wave functions, and the respective energies, of the
cylindrical potential well. First, the SchrÃ¶dinger equation is written in
cylindrical coordinates $(r,\theta)$, and a solution by separation of
variables is proposed. After this is done, two equations arise. The angular
equation is satisfied by a complex exponential $e^{im\theta},$ where $m$ is
the angular quantum number. The radial equation yields Bessel functions of the
first kind. After these two solutions are combined, the boundary condition is
imposed, that is, the requirement that the wave function must vanish when
$r=R$. This condition, similarly to what happens in the 1D infinite well,
defines the energy spectrum. In summary, the wave functions are:

$\psi_{nm}(r,\theta)=\frac{\mathcal{C}_{nm}}{\sqrt{2\pi}}J_{m}\left(\frac{q_{mn}r}{R}\right)e^{im\theta},$
(11)

where $J_{m}$ is a Bessel function of the first kind, $q_{mn}$ is the $n-$th
zero of $J_{m}(x)$ and
$\mathcal{C}_{nm}=\sqrt{2}\left[\left|J_{m+1}(q_{mn})\right|R\right]^{-1}$ is
a normalization constant. The energy spectrum reads:

$E_{nm}=\frac{1}{2}\left(\frac{q_{mn}}{R}\right)^{2}.$ (12)

In Table 1 we present the numerical value of a few energy levels for a system
with $R=2$ a.u..

$m=0$ | $m=1$
---|---
$n=1$ | $n=2$ | $n=3$ | $n=1$ | $n=2$ | $n=3$
0.723 | 3.809 | 9.361 | 1.835 | 6.152 | 12.937
Table 1: Numerical values of the energies (in atomic units) of the infinite
cylindrical potential well, for the first three states with $m=0$ and $m=1$.
Transitions from the ground state $(n=1$, $m=0)$ to the states with $m=1$
require the energies: 1.112, 5.429 and 12.214.

Now that we are in possession of the unperturbed wave functions and energy
spectrum, the computation of the dynamical polarizability due to the presence
of a time dependent external field can be performed. In agreement with the
previous section, we will consider a harmonic electric field aligned along the
$x-$direction. Furthermore, we will be concerned with the polarizability of
the ground state, which in this case corresponds to the state with $m=0$ and
$n=1$. To compute $\alpha(\omega)$ we will have to compute the matrix elements
$\langle\psi_{10}|x|\psi_{1}^{\pm}\rangle$. To do so we must first solve Eq.
(7) and find $|\psi_{1}^{\pm}\rangle$. For $0<r<R$ we have:

$\left(\frac{\mathbf{p}^{2}}{2}-E_{10}\pm\omega\right)\psi_{1}^{\pm}(\mathbf{r})=-r\cos\theta\psi_{10}(\mathbf{r}).$
(13)

We now set $\psi_{1}^{\pm}(\mathbf{r})$ as
$\left[\mathcal{P_{\pm}}(r)/\sqrt{r}\right]\cos\theta$. This transformation
has two advantages, as it takes care of the angular part of the differential
equation, while simultaneously eliminating any term proportional to a first
derivative in $r$. With this assumption we find:

$8r^{7/2}\psi_{10}(r)+\left[3+8r^{2}(-E_{0}\pm\omega)\right]\mathcal{P}_{\pm}(r)=4r^{2}\mathcal{P_{\pm}}^{\prime\prime}(r),$
(14)

subject to the boundary conditions $\mathcal{P}(0)=0$ and $\mathcal{P}(R)=0$.
A numerical solution to this differential equation can be easily obtained, and
once this is done, computing the polarizability is a trivial step, amounting
to a numerical quadrature. In Fig. 2 we plot the polarizability obtained after
the differential equation was solved numerically. In order to obtain the real
and imaginary parts of the polarizability, a small imaginary term was
subtracted to $\omega$, that is $\omega\rightarrow\omega-i\delta$, with
$\delta\ll\omega$; this parameter controls the linewidth of the resonances,
and should be chosen to match experimental results. Since our goal is not to
compare this simple example with experimental data we simply set $\delta=0.05$
a.u.. Inspecting the results in Fig. 2, we find that, as expected, clear
resonances appear at energies corresponding to transitions from the ground
state to states with $m=1$, as can be confirmed by recalling the values
presented in Table 1. The first resonance, corresponding to a transition from
the ground state to the state with $n=m=1$ is by far the dominant one, being
orders of magnitude more intense than the second and third peaks. Although not
depicted, at higher energies other peaks with even smaller oscillator
strengths appear.

Figure 2: Plot of the real (top) and imaginary (bottom) parts of the dynamical
polarizability for an infinite cylindrical well, obtained using Eq. (8), and
with $\psi_{1}^{\pm}$ computed numerically. Both plots present three
resonances whose positions exactly coincide with the energies associated with
transitions from the state with $m=0$ and $n=1$, to the three states with
$m=1$ presented in Table 1. In order to obtain the real and imaginary parts of
the polarizability, a small imaginary term, $\delta=0.05$ a.u., was subtracted
to $\omega$. We considered the particle to be an electron, confined to a disk
of radius $R=2$ a.u..

##  III Application to two-dimensional excitons

In the present section we will apply the same theoretical ideas developed so
far to compute the optical polarizability of excitons in 2D materials.
Contrary to the case of the infinite cylindrical potential well, the existence
of a Coulomb-like potential combined with an infinite geometry transforms this
problem into a rather involved one. To cope with the inherent difficulties, we
will consider the problem to be confined on a finite disk, and solve the
Fowler-Karplus’ differential equation using a variational approach. In the end
we will consider the radius of the disk to be sufficiently large so that the
obtained results are numerically identical to the ones we would obtain with an
infinite geometry, allowing us to compare our predictions with the
experimental results of Poellmann _et al_ Poellmann  _et al._ (2015).

###  III.1 Relation between polarizability, dielectric function, and optical
conductivity

Consider a gas with $N_{X}$ excitons per unit area in a two dimensional
material. The concentration $N_{X}$ is controlled by the interband excitation
process (fluence of the laser), the exciton recombination, and disorder. If
the electric field creating the exciton gas is weak, the excitons’
concentration is small and the gas can be considered non interacting. The
polarization of a 2D material due to _intra-exciton transitions_ is defined as

$\mathbf{P}_{{\rm ex}}^{{\rm intr.}}=N_{X}\mathbf{p}=N_{X}\alpha\mathbf{E},$
(15)

where $\mathbf{p}$ is the dipole moment of the exciton, in units of charge
times length and $\alpha$ is the polarizability, in unites of $\epsilon_{0}$
times volume, with $\epsilon_{0}$ the dielectric vacuum permittivity. The
susceptibility of the exciton gas associated with intra-excitons transitions
is defined has

$\chi_{{\rm ex}}^{{\rm intra.}}=\frac{\mathbf{P}_{{\rm ex}}^{{\rm
intr.}}}{\epsilon_{0}\mathbf{E}}=N_{X}\alpha/\epsilon_{0},$ (16)

and has units of length. We note that, in addition to $\chi_{{\rm ex}}^{{\rm
intra.}}$, there are other terms contributing to the total susceptibility of
the excitons gas: transition from the exciton vacuum to the $n$s levels and to
other even parity states, and the background susceptibility, accounting for
transitions to higher-energy bands, as well as other processes contributing to
the total dielectric function of a TMD Gu  _et al._ (2019). These, however,
are not the focus of this paper and have been addressed elsewhere Epstein  _et
al._ (2020). Thus, we introduce the dielectric function associated with intra-
exciton transition in a diluted medium (that is, low exciton density) as,

$\Delta\epsilon(\omega)=\chi_{{\rm ex}}^{{\rm
intra.}}/d=\alpha\frac{N_{X}}{\epsilon_{0}d}\frac{1}{g_{s}g_{v}},$ (17)

where $d$ is the effective thickness of the TMD layer and $g_{s}$ and $g_{v}$
are the spin and valley degeneracy factors.The real and imaginary parts of the
dielectric function due to intra-exciton transitions follow from this
equation. The relation between the dielectric function and the conductivity,
$\sigma$, follows from Kaindl  _et al._ (2009)
$\Delta\epsilon=i\sigma/(\epsilon_{0}\omega)$, and we obtain

$\sigma=-i\epsilon_{0}\omega\Delta\epsilon.$ (18)

The computation of these two quantities, $\Delta\epsilon(\omega)$ and
$\sigma$, for the $1$s$\rightarrow n$p excitonic transitions is our main goal.
(It would be interesting to derive a Clausius-Mossotti Ryazanov and Tishchenko
(2006), relation for two-dimensional excitons, although this is not the focus
of our work.)

###  III.2 The variational wave function to the ground state of the exciton

In what follows we will compute the optical polarizability of the excitonic
ground state ($1$s state). In order to do so, we need to first determine its
wave function. Finding an analytical expression for the $1s$ exciton state is
not a simple task Henriques  _et al._ (2020a, b), and because of that we will
consider the following double exponential variational ansatz to describe the
ground state Pedersen (2016) (see also Quintela and Peres (2020)):

$\psi_{0}=\mathcal{N}\left(e^{-ar}-be^{-\gamma ar}\right),$ (19)

where $\mathcal{N}$ is a normalization constant and $a,$ $b$, and $\gamma$ are
variational parameters obtained from energy minimization. This variational
ansatz is considerably more accurate than if a single exponential is used. In
the context of the excitonic problem, the Hamiltonian given in Eq. (2)
becomes:

$H=-\frac{1}{2\mu}\nabla^{2}+V_{{\rm RK}}(r)-\mathbf{r}\cdot\mathbf{F}(t),$
(20)

where $\mu$ is the reduced mass of the electron-hole system, $\mathbf{F}(t)$
is the external field previously defined in Eq. (3), and $V_{{\rm RK}}(r)$ is
the Rytova-Keldysh potential Rytova (1967); Keldysh (1979):

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

with $\kappa$ the mean dielectric constant of the media above and below the
TMD, $r_{0}$ is an intrinsic parameter of the 2D material which can be related
to the effective layer thickness $d$ Poellmann  _et al._ (2015), and
$\mathbf{H}_{0}$ and $Y_{0}$ are the Struve $H$-function of zero-th order and
the Bessel function of the second kind of the same order, respectively. This
potential is the solution of the Poisson equation for a charge embedded in a
thin film. For large distances the Rytova-Keldysh potential presents a Coulomb
tail, and diverges logarithmically near the origin.

Having determined a way to compute the excitonic ground state wave function
and binding energy, we turn our attention to solving the Fowler’s and Karplus’
Fowler (1984); Karplus and Kolker (1963) differential equation introduced in
Eq. (7). Comparing the excitonic problem with the simple example treated in
the previous section, we stress two key differences. While in the case of the
cylindrical well our problem had a bounded geometry, that is, the domain of
our problem was necessarily restricted to $r<R$, in the present scenario we
find an unbounded problem, whose domain extends up to infinity. Moreover, in
the previous section no potential term appeared in our calculations, while in
the present case we find the Rytova-Keldysh potential which, as previously
noted, presents a Coulomb$-1/r$ tail at large distance and a $\ln(r)$
divergence at short distances. This slow decaying behavior is known to be
problematic in numerical calculations. Considering these two remarks, it is
clear that finding the solution to the differential equation will not be as
straight forward a process as it was before.

The first step we take in order to proceed with the calculations is to confine
our system on a finite disk with large radius $R$. In practice, when numerical
values are introduced, we will consider $R$ to be sufficiently large in order
to obtain results which are numerically identical to the ones that one would
find for an unbounded problem, and thus be able to compare our results with
those from Ref. Poellmann  _et al._ (2015). This approach would also prove its
usefulness if a genuinely finite system was studied. The confinement is
reflected on the wave function of the ground state by the introduction of an
additional multiplicative term $(r-R)$, that is:

$\psi_{0}^{{\rm fin.}}(r)=\mathcal{N}^{{\rm fin}}\left(e^{-ar}-be^{-\gamma
ar}\right)\left(r-R\right).$ (22)

For a sufficiently large $R$ one finds:

$\lim_{R\rightarrow\infty}\psi_{0}^{{\rm fin.}}(r)=\psi_{0}(r).$ (23)

As before, the parameters $a,$ $b$ and $\gamma$ are determined from energy
minimization. This transformation does not allow us to immediately solve the
differential equation. However, it unlocks a new way of approaching the
problem.

###  III.3 The dynamic variational method

Since we are now effectively working on a finite disk, and we saw in the
previous section that Bessel functions are an appropriate complete set of
functions to describe a problem in such a geometry, we propose that
$\psi_{1}^{\pm}(\mathbf{r})$ can be written as:

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

where $J_{1}(z)$ is a Bessel function of the first kind, $z_{1n}$ is the
$n-$th zero of $J_{1}(z)$, $R$ is the radius of the disk, $N$ is the number of
Bessel function we choose to use, and $\left\\{c_{n}^{\pm}\right\\}$ are a set
of coefficients yet to be determined. As proposed in Ref. Karplus and Kolker
(1963); Montgomery and Rubenstein (1978); Yaris (1963, 1964) the values of
$\left\\{c_{n}^{\pm}\right\\}$ are determined from the minimization of the
following functional:

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

We now recall the following relation regarding the orthogonality of Bessel
functions on a finite disk of radius $R$ Jeffrey and Zwillinger (2014):

$\int_{0}^{R}J_{\nu}\left(\frac{\alpha_{\nu
m}r}{R}\right)J_{\nu}\left(\frac{\alpha_{\nu
n}r}{R}\right)rdr=\frac{R^{2}}{2}\delta_{nm}\left[J_{\nu+1}(\alpha_{\nu
m})\right]^{2},$ (26)

where $\alpha_{\nu m}$ is the $m-$th zero of $J_{\nu}(z).$ Using this
relation, one easily shows that the functional can be rewritten as:

$\displaystyle\mathcal{J}_{\pm}$ $\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}$
$\displaystyle+\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},$
(27)

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

$\displaystyle\mathcal{I}_{kn}$ $\displaystyle=\int
J_{1}\left(\frac{z_{1k}r}{R}\right)V_{{\rm
RK}}(r)J_{1}\left(\frac{z_{1n}r}{R}\right)rdr$ (28)
$\displaystyle\mathcal{S}_{n}$ $\displaystyle=\int
J_{1}\left(\frac{z_{1n}r}{R}\right)\psi_{0}^{{\rm fin.}}(r)r^{2}dr.$ (29)

Remembering our goal of minimizing the functional $\mathcal{J}_{\pm}$, we must
now differentiate it with respect to the different $c_{n}^{\pm}$ with
$n\in\left\\{1,2,...,N\right\\}$. Doing so, and noting in passing that
$\mathcal{I}_{kn}=\mathcal{I}_{nk}$, one finds:

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

with $j\in\left\\{1,2,...,N\right\\}$. It is now clear that this defines a
linear system of equations whose solution determines the values of the
coefficients $c_{n}^{\pm}$. Moreover, we can write this concisely using matrix
notation, as follows:

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

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)$ (32)
$\displaystyle\mathbf{S}^{{\rm T}}$
$\displaystyle=\left(\mathcal{S}_{1},\mathcal{S}_{2},\ldots,\mathcal{S}_{N}\right),$
(33)

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

$\left(\mathbb{M}\right)_{ij}=\begin{cases}g_{i}^{\pm}(\omega)+\mathcal{I}_{ii}&i=j\\\
\mathcal{I}_{ij}&i\neq j\end{cases},$ (34)

where $g_{i}^{\pm}(\omega)=R^{2}\left[z_{1i}^{2}/(2\mu
R^{2})-E_{0}\pm\omega\right]/2$. After $\mathbb{M}$ and $\mathbf{S}$ are
computed using Eqs. (28) and (29), respectively, the coefficients that
determine $\psi_{1}^{\pm}(\mathbf{r})$ are readily obtained:

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

and the solution of the differential equation is found. As was previously
mentioned, when specific values are introduced in the problem, the radius $R$
must be chosen large enough in order to produce results which are identical to
the ones of an unbounded system. As the value of $R$ increases, the value of
$N$ must also increase in order to accurately describe the desired physical
system. At first, this may seem problematic, since for each value of $\omega$
we would have to compute $\mathbb{M}$ (with $N^{2}$ entries) and solve the
linear system associated with it. However, two aspects prevent this approach
of becoming inadequately inefficient as $N$ increases: i) the matrix
$\mathbb{M}$ is symmetric, reducing the number of independent entries from
$N^{2}$ to $N(N+1)/2$; ii) more importantly, the most time consuming part of
the calculation, computing all the relevant $\mathcal{I}_{ij}$, is independent
of $\omega$, and does not need to be reevaluated every time a new energy is
used.

$\mu/m_{0}$ |  $r_{0}$ /$a_{0}$ | $\kappa$ |  $N_{X}$ (cm-2) | $d/a_{0}$ |  $R$ /$a_{0}$ | $N$
---|---|---|---|---|---|---
0.167 | $54$ | 3.32 | 1.8$\times 10^{12}$ | 6.04 | 800 | 35
Table 2: Parameters used to describe excitons in WSe2. The value of $R$ was
chosen in such a way that the wave functions of the first excited states are
effectively zero even before $r=R$. Inspecting Figure 4, we observe that the
states with $m=1$ and $n=2,3$ have wave functions which essentially vanish for
$r\sim 400a_{0}$, where $a_{0}\sim 0.53$ $\AA$ is the Bohr radius. The number
of Bessel functions ($N$) was chosen as the minimum number which upon addition
of more functions leaves the final result unchanged. The parameters related
with WSe2were taken from Ref. Poellmann  _et al._ (2015). The reduced mass
$\mu$ is given in units of the bare electron mass $m_{0}$.

At last, to compute the dynamical polarizability we look back at Eqs. (8) and
(24) and write:

$\alpha(\omega)=-2\pi\left(\mathbf{c}^{+}(\omega)+\mathbf{c}^{-}(\omega)\right)\cdot\mathbf{S},$
(36)

Now, to put together everything we developed so far to the test, let us
consider the specific case of excitons in WSe2. The different parameters that
characterize this problem are shown in Table 2. Note that the value of the
radius was chose as $R=800a_{0},$ with $a_{0}\sim 0.53$ $\AA$, the Bohr
radius. The choice of this value is understood from the inspection of Fig 4,
in the Appendix, where we present numerically obtained wave functions for the
first excited states of WSe2 (considering an unbounded system). There we
observe that the wave functions are essentially zero for $r\sim 400a_{0}$, and
thus the choice of $R=800a_{0}$ renders our problem effectively infinite.
Moreover, in Table 2 we also observe that $N=35$ Bessel functions were used.
This value was chosen since it corresponds to the minimum value of functions
upon which addition of more functions leaves the final result invariant.
Inspection of the values of $\mathbf{c}^{\pm}$ allowed us to conclude that the
method converged.

###  III.4 Results

In Fig. 3 we depict the real part of the dielectric function and the optical
conductivity, obtained from the real and imaginary parts of the the
polarizability using Eqs. (17) and (18), for two different linewidths, 10 meV
and 50 meV, and compare our theoretical prediction with experimental results
for the latter case. The linewidths are introduced in the problem by
subtracting a small imaginary part (10 or 50 meV in the present case) to the
photon energy. This is also what allows us to extract both the real and
imaginary parts from Eq. (17). To model the problem the values of Table 2 are
used.

Let us start by looking at the cases where a broadening of 10 meV was used.
Similarly to what was found for the cylindrical well, we observe two clearly
visible resonances at the energies corresponding to transitions from the
excitonic ground state to the first two excited states with $m=1$, that is the
$2$p and $3$p states. As expected, the resonance associated with the
$1s\rightarrow 2p$ transition is significantly more intense than the others.
More resonances exist after the second one, but are not resolved due to the
value of the broadening parameter and their smaller oscillator strength. No
resonances appear above the excitonic ground state energy. If the value of $R$
was chosen smaller this may not have been the case due to the appearance of
states with positive energy as a consequence of the confinement Aquino  _et
al._ (2005), however, these extra resonances would be barely visible due to
the smaller oscillator strengths. The values of the binding energies of the
different excitonic states were computed numerically using the Numerov-
shooting method described in the Appendix.

When $\Delta=50$ meV is used (a value fixed by the phenomenological fit
performed in Ref. Poellmann  _et al._ (2015)) we observe a different scenario
from the previous one. Due to the increased value of the linewidth, and the
way we introduced this parameter in the calculations, the $1s\rightarrow 2p$
resonance is shifted away from its theoretical value (computed numerically
using the shooting method described in the Appendix). This artifact needs to
be manually corrected. Contrary to the previous case, this time only the main
resonance is resolved. This is a consequence of the high value of the
broadening parameter considered, as well as the proximity of the different
resonances and their small oscillator strength when compared to the
1s$\rightarrow$2p transition. Moreover, for this specific scenario, we also
depict the experimental results measured in Ref. Poellmann  _et al._ (2015),
where the authors found that the resonance is well described by
phenomenologically modeling the $1s\rightarrow 2p$ transition with a
Lorentzian oscillator. A good agreement between our theoretical prediction and
the experimental result is visible. Our microscopic theory captures both the
position and magnitude of the measured resonances.

Figure 3: Real part of the dielectric function (top) and optical conductivity
(bottom) for two different linewidths. 10meV (left) and 50meV (right). When a
broadening of 10 meV is considered two resonances are clearly visible,
appearing at energies corresponding to transitions from the ground state
(1$s$) to the states $2p$ and $3p$. When a broadening of 50 meV is used, only
the main resonance is resolved. For this case we also depict the experimental
values measured in Ref. Poellmann  _et al._ (2015) (large black dots). The
agreement between our theoretical prediction and the experimental points is
remarkable. The parameters of Table 2 were used. The energies corresponding to
the transitions above mentioned are given in Table 3. The vertical dashed
lines represent the energy difference $E=E_{n{\rm p}}-E_{1{\rm s}}$ as
determined from the solution of the Wannier equation. As shown in the
experiments Poellmann  _et al._ (2015), the there is a relation between the
level broadening and the exciton’s concentration. Therefore, the panels on the
left-hand side of this figure do not represent a realistic experimental
situation but rather show the role that reducing the broadening has in
revealing higher energy dark states.

##  IV Conclusions

In this paper we have used Fowler’s and Karplus’ methods to access the
polarizability of 2D excitons when transitions from bright ($n$s) states to
dark ($n$p) states take place.

For benchmarking this method we first applied it to a single particle in a two
dimensional disk of radius $R$ with vanishing Dirichlet boundary conditions
(we note that the optical properties of condensed matter systems are sensitive
to the boundary conditions Montgomery and Pupyshev (2016)). We showed that the
peaks in the polarizability are in good agreement with the energy difference
between the ground state and the excited states of the particle in the disk.

Afterwards, we applied Fowler’s and Karplus’ method to the central problem of
our work, that of finding the dielectric function of the exciton gas. To do
so, we looked for the transitions from bright ($n$s) to dark ($n$p) states,
which can be accessed in pump-probe experiments. We chose the transition metal
dichalcogenide WSe2 as our test subject, since measurements of the dielectric
function due to the exciton gas have been made recently Poellmann  _et al._
(2015). We found a good agreement between the results obtained from our
microscopic theory and the experimental data. Indeed the Fowler-Karplus
approach allows us to access the excitonic response of the system without much
work, only requiring the knowledge of the ground state of the excitonic
problem. This is a major advantage over the overwhelming fully computational
calculations based on _ab-initio_ methods. For making our work mostly
analytical, we use an analytical variational wave function to describe the
ground state of the exciton problem. The numerical part of the calculation
consists on computing and solving a well behaved linear system of equations.
The role of disorder on the visibility of $1$s$\rightarrow n$p transitions was
also discussed. We showed that for small disorder, $\Delta=10$ meV, several
transitions to different $n$p states are visible in the polarizability
spectrum, whereas for large disorder, $\Delta=50$ meV, the first transition
$1$s$\rightarrow 2$p is rather broad, as seen in the experiments, and higher
order transition are masked by the linewidth of the peak.

Finally, we speculate that a similar experiment made on WSe2 encapsulated in
boron-nitride will allow the revealing of the transition to higher $n$p rather
than just the $1$s$\rightarrow 2$p transition. The nonlinear response of the
exciton gas Mossman  _et al._ (2016); Kocherzhenko  _et al._ (2019) will be
the subject of a future work.

## Acknowledgements

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. 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 and PTDC/NAN-
OPT/29265/2017.

##  Appendix A Numerov shooting method and the SchrÃ¶dinger equation in a log-
grid

In this appendix we give a description of the method used to compute the
binding energies of excitons in WSe2. To compute these energies one needs to
solve the Wannier equation (in atomic units):

$-\frac{1}{2\mu}\nabla^{2}\psi(\mathbf{r})+V_{{\rm
RK}}(r)\psi(\mathbf{r})=E\psi(\mathbf{r}),$ (37)

where, as in the main text, $\mu$ is the reduced mass of the electron-hole
system, $V_{{\rm RK}}$ is the Rytova-Keldysh potential defined in Eq. (21),
$\psi(\mathbf{r})$ is the exciton wave function, and $E$ its binding energies.

In order to solve this equation, we propose the usual solution by separation
of variables, such that $\psi(\mathbf{r})=R(r)\Theta(\theta)$. The angular
contribution trivially yields:

$\Theta(\theta)=\frac{e^{im\theta}}{\sqrt{2\pi}},$ (38)

where $m=0,\pm 1,\pm 2,...$ is the angular quantum number, and the
$\sqrt{2\pi}$ is a normalization factor. Considering the definition
$u(r)=R(r)\sqrt{r}$, the radial equation becomes:

$-\frac{1}{2\mu}\frac{d^{2}u}{dr^{2}}+\left[\frac{m^{2}}{2\mu
r^{2}}-\frac{1}{8\mu r^{2}}+V_{{\rm RK}}(r)\right]u(r)=Eu(r).$ (39)

Now, we found it useful to introduce the following change of variable:

$x=\ln r,$ (40)

which effectively transforms our problem from a linear to a logarithmic one.
This transformation proved to be useful in the stabilization of the numerical
calculations due to the divergent behavior of the Rytova-Keldysh potential at
the origin. Moreover, we introduce an auxiliary function, $f,$ defined as:

$u(x)=f(x)e^{x/2}.$ (41)

Doing so, the radial differential equation is transformed into:

$-\frac{1}{2\mu}f^{\prime\prime}(x)+\left[\frac{m^{2}}{2\mu}+e^{2x}V_{{\rm
RK}}(e^{x})\right]f(x)=e^{2x}Ef(x),$ (42)

with boundary conditions $f(\pm\infty)=0$. This is the equation that, when
solved, defines the binding energies for the excitons in WSe2. This equation
performs outstandingly for the usually problematic $m=0$ energy levels, where
the centrifugal barrier is absent.

To solve this equation we use the shooting method coupled to a Numerov
algorithm Izaac and Wang (2018); Berghe  _et al._ (1989); Johnson (1978,
1977); Pillai  _et al._ (2012). Given an initial energy guess, we integrate
Eq. (39) from left to right, and from right to left, and match the logarithmic
derivatives of two solutions somewhere sufficiently away from the edges. The
matched wave function will only represent a true bound state, with a correct
binding energy, when both the wave function itself as well as its derivative
are continuous across the whole domain. If this condition is not satisfied, a
new energy guess must be used and the process repeated. The choice of the
starting points of integration is of crucial importance to obtain faithful
results. The staring point on the left should be chosen as negative as
possible in order to accurately capture the behavior near $r=0$ (or
$x\rightarrow-\infty$). In practice, values such as $x_{{\rm min}}=-10$
already produce excellent results. The starting point on the right should be
chosen such that the wave function has effectively reached zero significantly
before $x_{{\rm max}}.$ Using semi-analytical methods, such as the one studied
in Ref. Henriques  _et al._ (2020a), allows one to compute a first guess of
the wave function, and thus determine what value of $x_{{\rm max}}$ should be
chosen in the application of the shooting-Numerov method. The method of Ref.
Henriques  _et al._ (2020a) is also a good starting point for the initial
guesses of the binding energies. In Ref. Izaac and Wang (2018) an efficient
algorithm to find the binding energies of Coulomb-like problems is described.
The method is numerically stable and fast.

Figure 4: Probability radial density of the the ground state ($n=1$ and $m=0$)
and the first excited states with $m=1$, of excitons in WSe2, obtained using
the shooting method. The parameters of Table 2 were used. The radius is given
in units of Bohr’s radius $a_{0}\approx 0.53$ $\AA$.

Using the previously described method, and the parameters presented in Table 2
of the main text, we were able to compute the wave functions of the ground
state, and the first two excited states with $m=1$, of excitons in WSe2
depicted in Figure 4. As expected, these wave functions resemble those one
would obtain with the Coulomb potential. However, since the Rytova-Keldysh
potential originates states with smaller binding energies than the Coulomb
potential, these wave functions are more spread out in space. Although only
three wave functions are presented, this method also allows the computation of
higher excited states.

$m=0$ | $m=1$
---|---
$n=1$ | $n=2$ | $n=3$ | $n=4$ | $n=5$
$-255$ | $-75$ | $-29$ | $-15$ | $-9$
Table 3: Binding energies (in meV) of a few excitonic states in WSe2, obtained
with the shooting method. The parameters of Table 2 were used. The energy
difference $\Delta E=-75+255=180$ meV coincides very well with the value
observed in the experiment at moderate densities of the exciton gas.

The binding energies of the ground state, as well as the first states with
$m=1$, are presented in Table 3. These energies, and in particular those
corresponding to the transitions $1s\rightarrow 2p$ and $1s\rightarrow 3p$,
are in agreement with Figure 3 of the main text.

## References

  * Lv  _et al._ (2015) R. Lv, J. A. Robinson, R. E. Schaak, D. Sun, Y. Sun, T. E. Mallouk, and M. Terrones, Accounts of Chemical Research 48, 56 (2015).
  * Ma  _et al._ (2020) Q. Ma, G. Ren, K. Xu, and J. Z. Ou, Advanced Optical Materials , 2001313 (2020).
  * Mueller and Malic (2018) T. Mueller and E. Malic, npj 2D Materials and Applications 2, 1 (2018).
  * Wang  _et al._ (2018) G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Reviews of Modern Physics 90, 021001 (2018).
  * Mak  _et al._ (2018) K. F. Mak, D. Xiao, and J. Shan, Nature Photonics 12, 451 (2018).
  * Schneider  _et al._ (2018) C. Schneider, M. M. Glazov, T. Korn, S. Höfling, and B. Urbaszek, Nature Communications 9, 2695 (2018).
  * Hsu  _et al._ (2019) 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, 2D Materials 6, 025028 (2019).
  * Chaves  _et al._ (2017) A. J. Chaves, R. M. Ribeiro, T. Frederico, and N. M. R. Peres, 2D Materials 4, 025086 (2017).
  * Zhang  _et al._ (2017) 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, Nature Nanotechnology 12, 883 (2017).
  * Yong  _et al._ (2019) C.-K. Yong, M. I. B. Utama, C. S. Ong, T. Cao, E. C. Regan, J. Horng, Y. Shen, H. Cai, K. Watanabe, T. Taniguchi, S. Tongay, H. Deng, A. Zettl, S. G. Louie, and F. Wang, Nature Materials 18, 1065 (2019).
  * Steinleitner  _et al._ (2017) P. Steinleitner, P. Merkl, P. Nagler, J. Mornhinweg, C. Schüller, T. Korn, A. Chernikov, and R. Huber, Nano Letters 17, 1455 (2017).
  * Koperski  _et al._ (2017) M. Koperski, M. R. Molas, A. Arora, K. Nogajewski, A. O. Slobodeniuk, C. Faugeras, and M. Potemski, Nanophotonics 6, 1289 (01 Nov. 2017).
  * Miyajima  _et al._ (2016) K. Miyajima, K. Sakaniwa, and M. Sugawara, Physical Review B 94, 195209 (2016), publisher: American Physical Society.
  * Berghäuser  _et al._ (2016) G. Berghäuser, A. Knorr, and E. Malic, 2D Materials 4, 015029 (2016).
  * Berghäuser  _et al._ (2018) G. Berghäuser, P. Steinleitner, P. Merkl, R. Huber, A. Knorr, and E. Malic, Physical Review B 98, 020301 (2018).
  * Merkl  _et al._ (2019) 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, Nature Materials 18, 691 (2019).
  * Tian  _et al._ (2020) T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos, Nano Letters 20, 841 (2020).
  * Wang  _et al._ (2006) F. Wang, J. Shan, M. A. Islam, I. P. Herman, M. Bonn, and T. F. Heinz, Nature Materials 5, 861 (2006).
  * Poellmann  _et al._ (2015) C. Poellmann, P. Steinleitner, U. Leierseder, P. Nagler, G. Plechinger, M. Porer, R. Bratschitsch, C. Schüller, T. Korn, and R. Huber, Nature Materials 14, 889 (2015).
  * Cha  _et al._ (2016) S. Cha, J. H. Sung, S. Sim, J. Park, H. Heo, M.-H. Jo, and H. Choi, Nature Communications 7, 10768 (2016).
  * Sie  _et al._ (2015) E. J. Sie, J. W. McIver, Y.-H. Lee, L. Fu, J. Kong, and N. Gedik, Nature Materials 14, 290 (2015).
  * Fröhlich  _et al._ (1985) D. Fröhlich, A. Nöthe, and K. Reimann, Physical Review Letters 55, 1335 (1985).
  * Zhang  _et al._ (2015) X.-X. Zhang, Y. You, S. Y. F. Zhao, and T. F. Heinz, Physical Review Letters 115, 257403 (2015).
  * Chao-Cador and Ley-Koo (2005) L. Chao-Cador and E. Ley-Koo, International Journal of Quantum Chemistry 103, 369 (2005).
  * Munkhbat  _et al._ (2020) B. Munkhbat, A. B. Yankovich, D. G. Baranov, R. Verre, E. Olsson, and T. O. Shegai, Nature Communications 11, 4604 (2020).
  * Podolsky (1928) B. Podolsky, Proceedings of the National Academy of Sciences 14, 253 (1928).
  * Karplus (1962) M. Karplus, The Journal of Chemical Physics 37, 2723 (1962).
  * Karplus and Kolker (1963) M. Karplus and H. J. Kolker, The Journal of Chemical Physics 39, 1493 (1963).
  * Fowler (1984) P. W. Fowler, Molecular Physics 53, 865 (1984).
  * Henriques and Peres (2020) J. C. G. Henriques and N. M. R. Peres, Phys. Rev. B 101, 035406 (2020).
  * Gu  _et al._ (2019) H. Gu, B. Song, M. Fang, Y. Hong, X. Chen, H. Jiang, W. Ren, and S. Liu, Nanoscale 11, 22762 (2019).
  * Epstein  _et al._ (2020) I. Epstein, B. Terres, A. J. Chaves, V.-V. Pusapati, D. A. Rhodes, B. Frank, V. Zimmermann, Y. Qin, K. Watanabe, T. Taniguchi, H. Giessen, S. Tongay, J. C. Hone, N. M. R. Peres, and F. H. L. Koppens, Nano Letters 20, 3545 (2020).
  * Kaindl  _et al._ (2009) R. A. Kaindl, D. Hägele, M. A. Carnahan, and D. S. Chemla, Phys. Rev. B 79, 045320 (2009).
  * Ryazanov and Tishchenko (2006) M. I. Ryazanov and A. A. Tishchenko, Journal of Experimental and Theoretical Physics 103, 539 (2006).
  * Henriques  _et al._ (2020a) J. C. G. Henriques, G. B. Ventura, C. D. M. Fernandes, and N. M. R. Peres, Journal of Physics: Condensed Matter 32 (2020a), 10.1088/1361-648X/ab47b3.
  * Henriques  _et al._ (2020b) J. C. G. Henriques, G. Catarina, A. T. Costa, J. Fernández-Rossier, and N. M. R. Peres, Phys. Rev. B 101, 045408 (2020b).
  * Pedersen (2016) T. G. Pedersen, Physical Review B 94, 125424 (2016).
  * Quintela and Peres (2020) M. F. C. M. Quintela and N. M. R. Peres, The European Physical Journal B 93, 222 (2020).
  * Rytova (1967) S. Rytova, Moscow University Physics Bulletin 22 (1967).
  * Keldysh (1979) L. Keldysh, Sov. J. Exp. and Theor. Phys. Lett. 29, 658 (1979).
  * Montgomery and Rubenstein (1978) H. E. Montgomery and T. G. Rubenstein, Chemical Physics Letters 58, 295 (1978).
  * Yaris (1963) R. Yaris, The Journal of Chemical Physics 39, 2474 (1963).
  * Yaris (1964) R. Yaris, The Journal of Chemical Physics 40, 667 (1964).
  * Jeffrey and Zwillinger (2014) A. Jeffrey and D. Zwillinger, eds.,  _Table of Integrals, Series, and Products_ , 8th ed. (Academic Press, Amsterdam ; Boston, 2014).
  * Aquino  _et al._ (2005) N. Aquino, G. Campoy, and A. Flores-Riveros, International Journal of Quantum Chemistry 103, 267 (2005).
  * Montgomery and Pupyshev (2016) H. E. Montgomery and V. I. Pupyshev, Physica Scripta 92, 015401 (2016).
  * Mossman  _et al._ (2016) S. Mossman, R. Lytel, and M. G. Kuzyk, JOSA B 33, E31 (2016).
  * Kocherzhenko  _et al._ (2019) A. A. Kocherzhenko, S. V. Shedge, X. Sosa Vazquez, J. Maat, J. Wilmer, A. F. Tillack, L. E. Johnson, and C. M. Isborn, The Journal of Physical Chemistry C 123, 13818 (2019).
  * Izaac and Wang (2018) J. Izaac and J. Wang,  _Computational Quantum Mechanics_, Undergraduate Lecture Notes in Physics (Springer International Publishing, 2018).
  * Berghe  _et al._ (1989) G. V. Berghe, V. Fack, and H. E. D. Meyer, Journal of Computational and Applied Mathematics 28, 391 (1989).
  * Johnson (1978) B. R. Johnson, The Journal of Chemical Physics 69, 4678 (1978).
  * Johnson (1977) B. R. Johnson, The Journal of Chemical Physics 67, 4086 (1977).
  * Pillai  _et al._ (2012) M. Pillai, J. Goglio, and T. G. Walker, American Journal of Physics 80, 1017 (2012).


# Cross-Domain Generalization Through Memorization:
A Study of Nearest Neighbors in Neural Duplicate Question Detection

Yadollah Yaghoobzadeh*, Alexandre Rochette*11footnotemark: 1, Timothy J. Hazen
Microsoft Turing
<EMAIL_ADDRESS>*Equal contribution.

###### Abstract

Duplicate question detection (DQD) is important to increase efficiency of
community and automatic question answering systems. Unfortunately, gathering
supervised data in a domain is time-consuming and expensive, and our ability
to leverage annotations across domains is minimal. In this work, we leverage
neural representations and study nearest neighbors for cross-domain
generalization in DQD. We first encode question pairs of the source and target
domain in a rich representation space and then using a k-nearest neighbour
retrieval-based method, we aggregate the neighbors’ labels and distances to
rank pairs. We observe robust performance of this method in different cross-
domain scenarios of StackExchange, Spring and Quora datasets, outperforming
cross-entropy classification in multiple cases. We will release our codes as
part of the publication.

##  1 Introduction

Figure 1: An illustration of our $k-$NN. Question pairs from both Source and
Target domains are encoded by in a common representation space, which is the
result of unsupervised adaptation of BERT on Target data. The top k similar
items from Source are then aggregated based on their label and distance to
provide a final confidence.

Duplicate question detection (DQD) is an important application in information
retrieval and NLP (Burke et al., 1997; Jeon et al., 2005; Lei et al., 2016;
Nakov et al., 2016; Rücklé et al., 2019). It allows systems to recognize when
two questions share an answer. This is significant for community forums, such
as StackExchange111https://stackexchange.com/ (SE) to increase their
effectiveness in avoiding redundant questions and displaying relevant answers
to search questions. It is also important for FAQ retrieval question answering
systems (Sakata et al., 2019).

To learn DQD models for SE, question pairs are usually annotated with
duplication information that is extracted from community-provided meta-data.
Such annotations are sparse for most domains, e.g., a new SE forum providing
support for a new product. Therefore, leveraging other training signals either
from unsupervised data or supervised data from other domains is important
(Shah et al., 2018; Poerner and Schütze, 2019).

Pre-trained language models (PLMs) like BERT (Devlin et al., 2018) and RoBERTA
(Liu et al., 2019) are great unsupervised textual representations. Several
recent efforts adapt PLMs for the domains of interest by self-supervised fine-
tuning on unsupervised domain data, which has shown to be promising in several
scenarios (Lee et al., 2019; Beltagy et al., 2019; Han and Eisenstein, 2019;
Gururangan et al., 2020). We follow that and tune BERT on SE domains to obtain
richer representations for the task of DQD.

Recently, $k$-nearest neighbors ($K-$NNs) is applied on the PLM
representations for language modeling (Khandelwal et al., 2019) and dialogue
(Fan et al., 2020). We extend this line of study and apply $k-$NN for cross-
domain generalization in DQD, where the models are trained on data from a
_source_ domain, and applied on data from a _target_ domain. To do so, we
represent pairs from source and target in a common representation space and
then score target pairs using nearest neighbors in the source pairs. Figure 1
shows an illustration of this procedure.

Our study on AskUbuntu as target and source datasets of (Shah et al., 2018),
which include several domains of SE and also Quora and Sprint, reveals that
$k-$NN is more effective compared to cross-entropy classification if (i) the
pair representation space from PLMs is rich for the target domain, i.e.,
adapted on the unsupervised data from target or similar domains; or (ii)
source and target domains have large distributional shifts.

We make the following contributions: (i) We present the first study of
combining strengths of $k-$NN and neural representations for cross-domain
generalization in a sentence matching task, i.e., DQD. (ii) Our experimental
results on cross-domain DQD demonstrate that $k-$NN on rich question-pair
representations advances the results of cross-entropy classification,
especially when shifts in source to target domains is substantial.

##  2 Related Work

Sparsity in DQD labeled examples in the SE domains is tackled by leveraging
the unsupervised data (Poerner and Schütze, 2019), the supervised data from
other domains, or both (Shah et al., 2018; Uva et al., 2018; Rücklé et al.,
2019; Rochette et al., 2019). We follow these approaches and learn
representations from unsupervised data and apply them for better
generalization when external supervised data in other domains is used.

A combination of $k-$NN with neural representations is the subject of several
earlier work, mostly in image classification (Papernot and McDaniel, 2018;
Cohen et al., 2018) Papernot and McDaniel (2018) show that $k-$NN is more
robust to out-of-distribution examples. More related to our work, Khandelwal
et al. (2019) apply $k-$NN on neural representations computed from and applied
on language modeling task and interpolates its scores with Softmax. They
validate that this is effective in different scenarios, including domain
adaptation. Here we do not compute representations for $k-$NN on the same task
as the one we apply; the representations are computed by language modeling and
applied on DQD.

##  3 Cross-Domain DQD

DQD is to identify pairs of questions answered by the same information. We
address DQD for community forums like SE. In SE, several domains are built to
address diverse user needs. We address DQD in a scenario where we perform on a
_Target_ domain with no task supervised data, by leveraging its unsupervised
text if exists, and labeled examples from a _Source_ domain.

###  3.1 Question Pair Representations

Our work leverages BERT (Devlin et al., 2018), a transformer-based language
model, pre-trained on general text from Wikipedia and BookCorpus. We fine-tune
BERT using its self-supervised objectives on unlabeled questions in SE
formatted as (Title, Body). Specifically, we concatenate the Title and Body of
the questions of a domain to form documents for BERT. We adopt the terminology
of Gururangan et al. (2020) and call this process _domain-adaptive pre-
training ( DAPT)_ and the resulting model BERT${}_{\textrm{DAPT}}$. DAPT
tailors the representation towards the particular vocabulary, syntax, and
semantics of the problem space (DQD). It also adapts the representation to a
specific domain (e.g. AskUbuntu). Therefore, we expect
BERT${}_{\textrm{DAPT}}$ to produce richer initial representations for our
data.

To obtain a representation for a pair of questions $(q_{1},q_{2})$, we
concatenate the text of both items and feed that as input to
BERT${}_{\textrm{DAPT}}$. The two items are separated by the [SEP] token. We
regard the embedding of the first token ([CLS]) at the final layer as our pair
representation. This representation is then utilized in two ways to produce
DQD predictions, as described in the following.

###  3.2 Classification (CLF)

For CLF, we follow the standard BERT training for sequence classification
starting from BERT${}_{\textrm{DAPT}}$ as a better initial PLM for cross-
domain DQD. In this setting, the CLS embedding is the input to a
classification layer with a cross-entropy loss. The gradients are back-
propagated to the BERT${}_{\textrm{DAPT}}$ parameters through CLS embedding
and tuned for the DQD task.

###  3.3 $k-$Nearest Neighbors ($k-$NN)

We leverage the self-supervised representation $f_{t}$, corresponding to the
CLS embedding in BERT${}_{\textrm{DAPT}}$ and encode each pair $q_{s_{i}}$ in
the Source training set $D_{s}$ using $f_{t}$ and preserve its associated
label $y_{i}\in\\{duplicate,\neg duplicate\\}$, as illustrated in Figure 1. A
distance function $d$ between two vectors $u$ and $v$ is selected in order to
establish the nearest neighbors. We use the cosine distance for this purpose:
$d(u,v)=1-cosine(u,v)$

One score for each potential label $y$ of a test pair $q_{t}$ is then computed
using $s$, representing the fraction of the mass of $1-d_{i}$ in K of each.

$s(y)=\frac{\sum_{i}^{k}{\mathds{1}_{y_{i}=y}(1-d(f_{t}(q_{s_{i}}),f_{t}(q_{t}))}}{\sum_{i}^{k}{1-d(f_{t}(q_{s_{i}}),f_{t}(q_{t}))}}$
(1)

##  4 Datasets

Dataset | Questions | Train | Dev | Test
---|---|---|---|---
AskUbuntu | 305,769 | 9,106 | 1,000 | 1,000
SuperUser | 390,378 | 9,106 | - | -
Sprint | 31,768 | 9,106 | - | -
Quora | 537,211 | 9,100 | - | -

Table 1:  The number of duplicates taken from Shah et al. (2018). For
AskUbuntu, as our Target, Train, Dev and Test, and for others Train numbers
are shown. Note that all AskUbuntu questions are considered for
BERT${}_{\textrm{DAPT}}$.

We experiment on cross-domain DQD datasets of Shah et al.
(2018)222github.com/darsh10/qra_code (See Table 1). For AskUbuntu and
SuperUser, the positive examples are taken from the duplicate marks in SE. For
Sprint, three paraphrases are generated by annotators for each question in a
set of FAQ. In these three datasets, 100 negatives are sampled randomly per
each positive. The annotation of Quora comes from the released Quora question
pairs dataset (Kaggle, 2017).

We further extract all questions for SE domains from the dump
files.333archive.org/details/stackexchange. These are integrated in our
unsupervised adaptations. For our analysis in §5.2, we create two additional
unsupervised corpora from SE. The first is from Academia consisting of around
27K questions. The second is from 33 different SE domains (See Table 5 in
Appendix), composed of around 1.5M questions.

###  4.1 Lexical Similarity Statistics

| dup |  $\neg$dup |  $S$(AskUbuntu)
---|---|---|---
AskUbuntu | 0.16 | 0.03 | 1.00
SuperUser | 0.19 | 0.03 | 0.22
Sprint | 0.37 | 0.04 | 0.03
Quora | 0\. 47 | 0.30 | 0.12
Table 2: Lexical similarity (measured by Jaccard index) between question pairs
within domains grouped by classes (first two columns) and between vocabulary
of Sources and AskUbuntu.

We select AskUbuntu as our only Target. In Table 2 (last column), we show the
lexical similarity between each Source and the Target. Accordingly, Sprint and
Quora hold low similarity with AskUbuntu, and SuperUser is the most similar
domain.

We also present the similarity between paired questions in each class of
duplicate and non-duplicate in Table 2. We observe that a duplicate pair in
Sprint and Quora has higher word-overlap on average compared to SE datasets.
Quora has another significant difference: its negative pairs are selected to
have a high lexical overlap. This means that the labeling function in Quora is
different from others.

##  5 Experiments

To obtain BERT${}_{\textrm{DAPT}}$, we fine-tune bert-base-cased
(BERT${}_{\textrm{BASE}}$) using language modeling scripts in Transformers
(Wolf et al., 2019) for $3$ epochs using default hyperparameters. To clarify
the effects of DAPT, we experiment with BERT${}_{\textrm{BASE}}$ as well. For
CLF, we fine-tune BERT on task training data for $10$ epochs with a learning
rate of 5e-5, early stopping on the Target dev set. All our $k-$NN experiments
are done with with the $k=100$ using Faiss library Johnson et al. (2017).

#### Evaluation metric

Since the annotations are incomplete in SE, Shah et al. (2018) propose to use
AUC as the metric for DQD performance. They report the normalized AUC(.05),
which is the area under the curve of the true positive rate as function of the
false positive rate ($fpr$), from $fpr=0$ to $fpr=.05$. We follow the same
protocol and use AUC(.05) metric.

###  5.1 Results

|  | Target = AskUbuntu
---|---|---
| Model | AskUbuntu | SuperUser | Sprint | Quora
| _BERT ${}_{\textrm{DAPT}}$_
1 | CLF | .923 | .870 | .749 | .609
2 |  $k-$NN | .936 | .908 | .753 | .800
| _BERT ${}_{\textrm{BASE}}$_
3 | CLF | .899 | .779 | .562 | .515
4 |  $k-$NN | .871 | .755 | .649 | .621
| _From Shah et al. (2018)_ |
5 | BiLSTM | .858 | .796 | .615 | .446
Table 3: Comparing AUC(.05) results of our models and the baseline. Four
Sources are evaluated for the Target, i.e., AskUbuntu.
BERT${}_{\textrm{DAPT}}$ is the adaptation of BERT${}_{\textrm{BASE}}$ on
AskUbuntu unsupervised data.

In Table 3, we present the performance of our models for AskUbuntu evaluation
set given Source data from AskUbuntu, SuperUser, Sprint, or Quora. We add in-
domain (AskUbuntu as Source) results for comparison. We include the results of
Shah et al. (2018) obtained on the same data by learning domain-adversarial
BiLSTM models in line 5.

The first block corresponds to BERT${}_{\textrm{DAPT}}$: the adapted BERT on
AskUbuntu unsupervised data. We see that $k-$NN outperforms CLF in all cases
(line 2 vs. 1), confirming that $k-$NN is more robust if the pair
representation is rich. The most obvious improvement belongs to Quora as
Source (.609 to .800), where the labeling function shifts significantly (See
§4.1).

In the second block, BERT${}_{\textrm{BASE}}$ results in consistently worse
models (line 3-4) compared to BERT${}_{\textrm{DAPT}}$. Here for Sprint and
Quora, $k-$NN again outperforms CLF, giving more evidence about robustness of
$k-$NN in the case of domain shifts. However, for SuperUser, a closely related
domain to AskUbuntu and also AskUbuntu itself, $k-$NN underperforms. Given
that the input representations to $k-$NN are not tuned on SE in the case of
BERT${}_{\textrm{BASE}}$, this is not surprising. CLF fine-tunes the
representations as part of its task training on the Source data, which in this
case is a related or same domain as the Target.

| Unsupervised Data | CLF |  $k-$NN
---|---|---|---
1 | None | .779 | .755
2 | Source | .855 | .886
3 | Target | .870 | .908
4 | Unrelated | .778 | .782
5 | 33 SE domains | .891 | .917
Table 4: AUC(.05) of training on Source = SuperUser and evaluating on Target =
AskUbuntu, as a function of the input corpus for BERT${}_{\textrm{DAPT}}$.
None corresponds to BERT${}_{\textrm{BASE}}$ and Unrelated to Academia.
(Similar results but with Source = AskUbuntu is in Appendix, Table 6.)

###  5.2 Domain of Unsupervised Data in BERT${}_{\textrm{DAPT}}$

Here we aim to understand more about the impact of the domain of unsupervised
data on the quality of BERT${}_{\textrm{DAPT}}$. In Table 4, we report results
for SuperUser as Source and AskUbuntu as Target, and vary the unsupervised
corpus, starting from no data (i.e., BERT${}_{\textrm{BASE}}$).

We choose these domains (lines 2-5): Source (SuperUser), Target (AskUbuntu),
Unrelated (Academia) as a lexically distant domain to Target, and a set of 33
SE domains including Source and Target (See Table 5 for full list).

Table 4 demostrates that adaptation on Target data (line 3) is better than
either of Source (line 2) or the unrelated domain (line 4). Adaptation on a
large number of domains (line 5) is the best; the information across a diverse
set of domains is complementary for the task. Notably, we observe that the
representation is more critical for $k-$NN compared to CLF: the difference
between the best representation (33 domains) and worst
(BERT${}_{\textrm{BASE}}$) is much greater in $k-$NN compared to CLF. This
behavior is understandable as CLF further updates the representation during
its task training, while for $k-$NN, it remains fixed.

##  6 Conclusion

In this work, we studied applying $k-$NN in DQD cross-domain generalization.
We compared $k-$NN and a cross-entropy classifier when different question-pair
representations are available. Our results showed that domain-adaptive pre-
training on target data gives rich representations, and $k-$NN is more robust
against distributional shifts compared to classification if question pairs are
encoded by these rich representations.

We plan to extend our study to other tasks and understand better the strengths
of memorization in learning robust models where rich PLM embeddings are
utilized to represent examples. We believe concurrently that the promising
results and findings of this presented study could benefit other NLP research
to explore this direction more.

## References

  * Beltagy et al. (2019) Iz Beltagy, Arman Cohan, and Kyle Lo. 2019.  Scibert: Pretrained contextualized embeddings for scientific text.  _CoRR_ , abs/1903.10676.
  * Burke et al. (1997) Robin D Burke, Kristian J Hammond, Vladimir Kulyukin, Steven L Lytinen, Noriko Tomuro, and Scott Schoenberg. 1997.  Question answering from frequently asked question files: Experiences with the faq finder system.  _AI magazine_ , 18(2):57–57.
  * Cohen et al. (2018) Gilad Cohen, Guillermo Sapiro, and Raja Giryes. 2018.  DNN or k-nn: That is the generalize vs. memorize question.  _CoRR_ , abs/1805.06822.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  BERT: pre-training of deep bidirectional transformers for language understanding.  _CoRR_ , abs/1810.04805.
  * Fan et al. (2020) Angela Fan, Claire Gardent, Chloe Braud, and Antoine Bordes. 2020.  Augmenting transformers with knn-based composite memory for dialogue.  _arXiv preprint arXiv:2004.12744_.
  * Gururangan et al. (2020) Suchin Gururangan, Ana Marasović, Swabha Swayamdipta, Kyle Lo, Iz Beltagy, Doug Downey, and Noah A Smith. 2020.  Don’t stop pretraining: Adapt language models to domains and tasks.  _arXiv preprint arXiv:2004.10964_.
  * Han and Eisenstein (2019) Xiaochuang Han and Jacob Eisenstein. 2019.  Unsupervised domain adaptation of contextualized embeddings: A case study in early modern english.  _CoRR_ , abs/1904.02817.
  * Jeon et al. (2005) Jiwoon Jeon, W Bruce Croft, and Joon Ho Lee. 2005.  Finding similar questions in large question and answer archives.  In _Proceedings of the 14th ACM international conference on Information and knowledge management_ , pages 84–90.
  * Johnson et al. (2017) Jeff Johnson, Matthijs Douze, and Hervé Jégou. 2017.  Billion-scale similarity search with gpus.  _arXiv preprint arXiv:1702.08734_.
  * Kaggle (2017) Kaggle. 2017.  Quora question pairs.
  * Khandelwal et al. (2019) Urvashi Khandelwal, Omer Levy, Dan Jurafsky, Luke Zettlemoyer, and Mike Lewis. 2019\.  Generalization through memorization: Nearest neighbor language models.  _arXiv preprint arXiv:1911.00172_.
  * Lee et al. (2019) Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang. 2019.  Biobert: a pre-trained biomedical language representation model for biomedical text mining.  _CoRR_ , abs/1901.08746.
  * Lei et al. (2016) Tao Lei, Hrishikesh Joshi, Regina Barzilay, Tommi S. Jaakkola, Kateryna Tymoshenko, Alessandro Moschitti, and Lluís Màrquez. 2016.  Semi-supervised question retrieval with gated convolutions.  In _NAACL HLT 2016, The 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, San Diego California, USA, June 12-17, 2016_ , pages 1279–1289.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Nakov et al. (2016) Preslav Nakov, Lluís Màrquez, Alessandro Moschitti, Walid Magdy, Hamdy Mubarak, Abed Alhakim Freihat, Jim Glass, and Bilal Randeree. 2016.  SemEval-2016 task 3: Community question answering.  In _Proceedings of the 10th International Workshop on Semantic Evaluation_ , pages 525–545.
  * Papernot and McDaniel (2018) Nicolas Papernot and Patrick McDaniel. 2018.  Deep k-nearest neighbors: Towards confident, interpretable and robust deep learning.  _arXiv preprint arXiv:1803.04765_.
  * Poerner and Schütze (2019) Nina Poerner and Hinrich Schütze. 2019.  Multi-view domain adapted sentence embeddings for low-resource unsupervised duplicate question detection.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 1630–1641.
  * Rochette et al. (2019) Alexandre Rochette, Yadollah Yaghoobzadeh, and Timothy J. Hazen. 2019.  Unsupervised domain adaptation of contextual embeddings for low-resource duplicate question detection.
  * Rücklé et al. (2019) Andreas Rücklé, Nafise Sadat Moosavi, and Iryna Gurevych. 2019.  Neural duplicate question detection without labeled training data.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 1607–1617.
  * Sakata et al. (2019) Wataru Sakata, Tomohide Shibata, Ribeka Tanaka, and Sadao Kurohashi. 2019.  FAQ retrieval using query-question similarity and bert-based query-answer relevance.  In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR 2019, Paris, France, July 21-25, 2019_ , pages 1113–1116. ACM.
  * Shah et al. (2018) Darsh Shah, Tao Lei, Alessandro Moschitti, Salvatore Romeo, and Preslav Nakov. 2018\.  Adversarial domain adaptation for duplicate question detection.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 1056–1063.
  * Uva et al. (2018) Antonio Uva, Daniele Bonadiman, and Alessandro Moschitti. 2018.  Injecting relational structural representation in neural networks for question similarity.  _arXiv preprint arXiv:1806.08009_.
  * Wolf et al. (2019) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, R’emi Louf, Morgan Funtowicz, and Jamie Brew. 2019.  Huggingface’s transformers: State-of-the-art natural language processing.  _ArXiv_ , abs/1910.03771.

academia android anime apple askubuntu astronomy aviation bicycles biology
bitcoin boardgames buddhism chemistry christianity cogsci cooking crypto cs
gaming hinduism islam judaism linguistics mechanics meta.stackexchange
meta.superuser philosophy politics skeptics sports superuser unix workplace
---
Table 5: The 33 StackExchange domains used in our unsupervised BERT adaptation analysis. | unsupervised Data | CLF |  $k-$NN
---|---|---|---
1 | None | .899 | .871
3 | Target | .923 | .936
4 | Unrelated | .899 | .890
5 | 33 SE domains | .942 | .942
Table 6: AUC(.05) of training on Source = AskUbuntu and evaluating on Target =
AskUbuntu, as a function of the input corpus for BERT${}_{\textrm{DAPT}}$.
None corresponds to BERT${}_{\textrm{BASE}}$ and Unrelated to Academia.