# PhysNLU: A Language Resource for Evaluating Natural Language Understanding
and Explanation Coherence in Physics

###### Abstract

In order for language models to aid physics research, they must first encode
representations of mathematical and natural language discourse which lead to
coherent explanations, with correct ordering and relevance of statements. We
present a collection of datasets developed to evaluate the performance of
language models in this regard, which measure capabilities with respect to
sentence ordering, position, section prediction, and discourse coherence.
Analysis of the data reveals equations and sub-disciplines which are most
common in physics discourse, as well as the sentence-level frequency of
equations and expressions. We present baselines that demonstrate how
contemporary language models are challenged by coherence related tasks in
physics, even when trained on mathematical natural language objectives.

Keywords: mathematical text, physics, natural language understanding,
discourse coherence

languageresourceLanguage Resources

PhysNLU: A Language Resource for Evaluating Natural Language Understanding and
Explanation Coherence in Physics

Jordan Meadows1,2, Zili Zhou1, André Freitas1,2
---
1Department of Computer Science, University of Manchester
2Idiap Research Institute, Switzerland
<EMAIL_ADDRESS>
{zili.zhou<EMAIL_ADDRESS>

Abstract content

## 1\. Introduction

Physics literature is a form of mathematical language which is unique beyond
simply domain vocabulary. How physicists use mathematics to reason and
explain, separates their field fundamentally from other disciplines, including
mathematics. Many of its sub-disciplines are situated between pure mathematics
and engineering, while others conjoin computer science and biology, with
physical methods acting as a well-travelled bridge between the formal and
natural sciences. It has not been proven, for example, that smooth solutions
[Pizzocchero, 2021, Gala et al., 2021, Miller, 2021] always exist for the
Navier-Stokes equations (a millenium problem) despite their widespread use in
simulating and engineering fluid dynamics, while biophysics demonstrates that
fundamental problems in ecology and evolution can be characterized by
computational complexity classes [Ibsen-Jensen et al., 2015]. Physics
discourse serves as a universal mechanism for generating empirically
falsifiable quantitative theory in the natural sciences and engineering [Smith
and Fleck, 2017, Coffey and Kalmykov, 2012], separate to both pure mathematics
and any downstream field.
Its core traits are reflected in unique literary devices and its mathematical
explanations will differ to those in the formal sciences as a result. A
concrete example is the physics derivation; a core explanatory or
argumentative device less rigorous and more informal than mathematical proofs
[Meadows and Freitas, 2021, Davis, 2019, Kaliszyk et al., 2015], which
generally results in predictive equations relating physical quantities, rather
than generating a truth value for a given conjecture (e.g., twin primes). Such
equations are central components of physics descriptions, with natural
language forming around them and their elements, and their relation to other
equations through derivations.
Mathematics as a whole, particularly logic, is less concerned with this
predictive modelling of real world systems, let alone when such systems are
quantum or relativistic, or both. Suggesting that mathematicians work at a
level of abstraction higher than that of physicists (i.e., proof frameworks
compared to specific derivations), Feynman famously states that
“Mathematicians are only dealing with the structure of reasoning…”.
Within the unique sphere of physics literature, we introduce a suite of
datasets which together gauge a model’s proficiency in recognising whether or
not a physics-related explanation is coherent. In parallel with tasks inspired
by DiscoEval [Chen et al., 2019], we aim to “evaluate the discourse-related
knowledge captured by pretrained sentence representations” in the physics
domain. From the proposed data we show that modern pretrained language models
are challenged by these tasks even after fine-tuning, in particular
demonstrating that a recent language model [Shen et al., 2021] trained on a
large corpus of mathematical text, is outperformed by even vanilla BERT-Base
and all popular non-mathematical language models considered in this work.
We contribute the following:

  1. 1.

We introduce PhysNLU; a collection of 4 core datasets related to sentence
classification, ordering, and coherence of physics explanations based on
related tasks [Chen et al., 2019]. Each dataset comprises explanations
extracted from Wikipedia including derivations and mathematical language. We
additionally present 2 parent datasets extracted from 6.3k articles related to
physics, in both raw Wikipedia data, and in a form that mimics WikiText-103
[Merity et al., 2016], which is a popular dataset used in related work [Iter
et al., 2020]. PhysNLU is avaliable
online111https://github.com/jmeadows17/PhysNLU.

  2. 2.

We provide analysis of linguistic features of physics text, including insights
such as sentence and example-level distribution of mathematical content across
the datasets, the frequency in which explained concepts relate to physics sub-
domains, and the most frequent equations in the discourse.

  3. 3.

We demonstrate how the state-of-the-art does not exhibit proficient inference
capabilities with respect to tasks concerning order, coherence, relative
position, and classification of sentence-level physics explanations, even when
approaches have been designed explicitly for mathematical language, through
baselines extracted from experiments involving a selection of pretrained
language models.

## 2\. Task Description

The tasks considered in this work probe model proficiency across 4 categories,
originally designed for general language, but here employed specifically for
physics discourse containing mathematics. Binary Sentence Ordering tests the
ability of a model to recognise order at the shortest possible scale, between
two sentences. Sentence Position tests this order and position recognition at
a larger scale, closer to that of full paragraphs. Discourse Coherence tests
whether a model can determine whether a sequence of statements in an
explanation are continuous and relevant. Sentence Section Prediction tests how
well a model can link individual sentences to a specific section of an
explanation. Together, in our context, they evaluate the discourse-related
knowledge captured by pretrained sentence representations, and physics
explanation coherence with respect to order and sentence relevance. We now
describe our method for data collection for each of the 6 datasets, including
the 4 directly used in the forthcoming experiments for each task as described
in Figure 2, while an overview of our contributions are displayed in Figure 1.

Figure 1: Workflow and overview. Figure 2: Each numbered box is a physics
statement. Explanation of how to obtain the differential form of Gauss’ law
via the divergence theorem, used to demonstrate how 4 evaluation tasks handle
physics explanations (extracted from Wikipedia, including errors). (SP)
Sentence position (top left) takes a random sentence from the description and
moves it to position 1, then the model predicts the true position of the new
first sentence, which in this case is 4. (DC) Discourse coherence (top right)
randomly replaces a sentence with another from a different physics
explanation, where the model predicts if the explanation is coherent. (BSO):
Binary sentence ordering (bottom left) swaps two consecutive sentences, and
the model predicts if the second entails the first. (SSP): Sentence section
prediction takes a random sentence, and the model predicts if it belongs to an
introduction or otherwise.

## 3\. Dataset Collection

### 3.1. PhysNLU-WikiRaw

Starting from an English Wikipedia XML dump, we select articles with a mention
of “physics” and contain at least one equation defined with a <math> tag.
After cleaning articles to contain mostly mathematical natural language, and
removing those which are predominantly tables, this results in a dataset
containing 6.3k articles. We include article titles and corresponding raw
unedited text, as well as wikipedia article categories.

### 3.2. PhysNLU-WikiText

This data mimics WikiText-103 [Merity et al., 2016] which is used for the
approach introducing the CONPONO objective [Iter et al., 2020] during a
preprocessing stage222https://github.com/google-
research/language/blob/master/language/conpono/create_pretrain_data/wiki_preproc_pipeline.py.
Among other similarities, we opt to nest section titles within equals signs
(e.g. “ = Title = ”, “ = = Section = = ”) and omit reference and “see more”
sections. The major linguistic differences between WikiText-103 and PhysNLU-
WikiText are the inclusion of mathematical content as well as structures which
may contain mathematical expressions such as tables, which are infrequent.
This core dataset is taken as the starting point from which to derive the
other datasets. We then extract 516k unique sentences for use in the following
datasets, where sentences are determined by splits on full stops which,
similarly to commas, are separated by a space from words (e.g. “end of
sentence .”). We correct for issues with names (e.g. J. J. Hopfield) and
abbreviations, and some instances where full stops should be present but are
omitted.

### 3.3. PhysNLU-BSO (Binary Sentence Ordering)

We take all pairs of consecutive sentences from sections, where selected pairs
overlap. Each pair has a 50% chance that the pair order is swapped and we
include a label to denote whether a swap (1) has occurred or not (0), suitable
to be framed as binary classification. The BSO dataset contains 459k examples.

### 3.4. PhysNLU-SP (Sentence Position)

We take the first 5 sentences from each applicable section, select a sentence
at random and move it to the first position (shifting the others down). The
number of the swapped sentence is the label corresponding to each set of 5
sentences, suitable for multiclass classification. The SP dataset contains 40k
examples.

### 3.5. PhysNLU-DC (Discourse Coherence)

The first 6 sentences from each applicable section are selected, then between
positions 2 and 5 inclusive a sentence is swapped with another article at
random, with 50% swap occurrence. Whether a swap has occurred or not is
included as a label for each example for binary classification, and the DC
dataset contains 35k examples.

### 3.6. PhysNLU-SSP (Sentence Section Prediction)

All sentences from the introduction sections of each article are selected and
an equal number of sentences are extracted from elsewhere at random from the
corpus. Introduction sentences are associated with a label (1) while non-
introductory sentences are associated with a separate label (0) for binary
classification. The SSP dataset contains 90k examples.

## 4\. Dataset Statistics

We now analyse our data with a focus on equations and mathematical natural
language. Table 1 shows an overview of notable features, such as the
proportion of examples in each dataset which contains mathematical
expressions, or specifically equations.
Figure 3 describes the proportion of sentences which contain at least $n$
mathematical elements for $n\in[1,6]$, where an element is identified via 3
separate tags: <math>, $\\{$math$|$, and $\\{$mvar$|$. The lighter bars
correspond to all sentences present in the evaluation data, while the darker
bars correspond to sentences in the SSP dataset which contain proportionally
less math. Introductory sentences make up half of the data for SSP and usually
they do not contain mathematical language or equations, which accounts for
this gap. The proportion of math in sentences from the BSO, SP, and DC
datasets are practically equivalent to the overall proportion. Figure 4 shows
the relative frequency and proportion of the top 8 Wikipedia categories
associated with each article. A single article can correspond to a large
number of categories, out of 12.5k categories in our case. Notably, fields
related to quantum mechanics are by far the most frequent, where 10% of the
data corresponds to either “Quantum mechanics”, “Quantum field theory”, or
“Condensed matter physics”.
Figure 5 displays how often specific equations are present in the corpora. One
might be tempted to claim that this demonstrates how physicists tend to argue
and explain using initial conditions with respect to time ($t=0$),
displacement ($x=0$, $r=0$, $z=0$), and angle ($\theta=0$), however this exact
string matching is biased towards simple equations. As the complexity of
equations increases to include multiple terms, and many terms are equivalent
in meaning but different in notation, there will be multiple equations in the
data which correspond to the same physics. A more accurate way to assess this
would involve classifying groups of equations with a good math retrieval model
[Peng et al., 2021] and counting group frequency. This analysis does offer
insight for simple equations however. For example, it reflects the convention
that people prefer to start counting from $n=1$ in physics, which occurs more
frequently than $n=0$, that the famous $E=mc^{2}$ is more prolific than the
similarly famous $F=ma$, and that the most frequently discussed Maxwell
equation is $\nabla\cdot\mathbf{B}=0$.
Figure 6 shows the proportion of examples from each evaluation dataset which
contain at least $n$ counts of either a <math> equation or non-equational
math, for $n\in[1,6]$. The darker bars represent solely equations. DC, SP, BSO
and SSP examples comprise 6, 5, 2, and 1 sentences respectively, so we expect
that the proportion of math included decreases in that order.

Dataset | Size | % with math | % with equations
---|---|---|---
DC | 35 k | 45 | 35
SP | 40 k | 36 | 29
BSO | 459 k | 24 | 17
SSP | 90 k | 12 | 7
Table 1: Metadata for each evaluation dataset. The percentage of examples with
math refers to the percentage of examples with LaTeX text between the XML tags
<math> and </math>. The percentage of equations refers to only math which
contains at least one equality. Figure 3: For 516k unique sentences used in
the evaluation tasks, the proportion of sentences which contain at least $n$
math elements is shown for $n\in[1,6]$. The two <math> variants correspond to
LaTeX text identified by the XML <math> tag which respectively do and do not
contain an equality. $\\{$math$|$ and $\\{$mvar$|$ correspond to any
mathematical text identified by each marker. The darker bars correspond to
only sentences extracted from the SSP dataset. Figure 4: Top 8 most frequent
article categories out of 12.5k, excluding categories containing the phrase
“Articles containing …”, where each article may correspond to multiple
categories. Figure 5: Top 8 most frequent <math> tagged equations by exact
string matching after accounting for spaces, commas, and full stops. Figure 6:
The proportion of examples which contain at least $n$ math elements is shown
for $n\in[1,6]$, where examples are sourced from the DC, SP, BSO, and SSP
datasets. The math is identified via the <math> tag, where lighter bars
correspond to at least $n$ math of any kind, while the darker bars correspond
to the inclusion of only equations.

## 5\. Results

We evaluate models on 4 tasks from the DiscoEval suite [Chen et al., 2019]. We
remove the PDTB and RST related tasks due to the lack of a linguistic
framework for describing discourse relations in the physics context. Table 1
gives additional information regarding the data used for each task.

### 5.1. Evaluation Tasks

DiscoEval is “designed to evaluate discourse-related knowledge in pretrained
sentence representations”. We briefly describe the 4 evaluation tasks from
DiscoEval considered in our work, with examples shown in Figure 2. We use the
same conventions [Chen et al., 2019] for representing concatenation of vectors
$[.,.,...]$.
Sentence position (SP) involves considering 5 consecutive sentences at a time,
moving a random sentence to the first position, then predicting the correct
position of the first sentence. We take the first 5 sentences from every
paragraph in our data. Classifiers are trained by encoding the 5 sentences to
vector representations $x_{i}$, then vectors $x_{1}-x_{i}$ are concatenated to
$x_{1}$ for $i\in[2,5]$ as input to the classifier as:
$[x_{1},\;x_{1}-x_{2},\;x_{1}-x_{3},\;x_{1}-x_{4},\;x_{1}-x_{5}]$.
Binary sentence ordering (BSO) involves taking pairs of contiguous sentences
from a paragraph, swapping the order 50% of the time and predicting if a swap
has occurred. A classifier is trained by concatenating $x_{1}$ and $x_{2}$
with their element-wise difference as: $[x_{1},x_{2},x_{1}-x_{2}]$.
Discourse coherence (DC) involves taking 6 consecutive sentences, replacing a
sentence from positions 2-5 inclusive with a sentence from a random article
with 50% frequency and predicting if a swap has occurred. We take the first 6
sentences from each paragraph for this task. Each vector $x_{i}$ is
concatenated for input to the classifier as:
$[x_{1},x_{2},x_{3},x_{4},x_{5},x_{6}]$.
Sentence section prediction (SSP) involves sampling a sentence from either the
abstract of a scientific article or elsewhere with equal probability, and
predicting if the sentence belongs to the abstract. In our case, we sample
from article introductions as we do not have abstracts. The original task
[Chen et al., 2019] involves the omission of equations which increases the
difficulty of the task, but due to the nature of our problem space we leave
them in. The classifier input is just the vector representation $x_{1}$.

### 5.2. Baselines

We include 7 baseline transformer-architecture models [Vaswani et al., 2017]
in this study, BERT-base-uncased, BERT-large-uncased [Devlin et al., 2019],
RoBERTa-base [Liu et al., 2019], MathBERT [Shen et al., 2021], MegatronBERT
[Shoeybi et al., 2019], CONPONO (K=2), CONPONO (K=4) [Iter et al., 2020].
BERT-base-uncased, BERT-large-uncased are two BERT models pretrained on large-
scale common domain text corpora. Based on the 12 or 24 encoder layers, the
models achieved state-of-the-art performances on several NLP tasks such as
sentence classification, next sentence prediction, token classification, etc.
BERT-large-uncased has a larger parameter size than BERT-base-uncased. The
smaller model outperforms the larger in our case.
RoBERTa-base has the same architecture as BERT but is pretrained with more
data, different hyperparameters, and only full-length sentences, and shows the
BERT models were undertrained.
The MathBERT model uses the BERT architecture and further trains on the masked
language modelling objective using a large corpus of mathematical text,
including arXiv paper abstracts and textbooks covering pre-k to graduate-
level. Their model shows notable improvement over BERT in an array of
mathematical tasks.
MegatronBERT model improves the architecture of the original BERT models to
enable model deployment across distributed GPU environments, while
simultaneously improving the accuracy of the model.
The CONPONO models use the encoder architecture (and same data) from BERT to
encode text segments, but are pretrained instead with the CONPONO objective
together with Masked Language Modelling (MLM). The purpose is to let the
models learn discourse relationships between sentences with respect to order
and distance, while using negatives to increase sentence representation
quality. We use 2 versions of CONPONO as baselines in this paper, K=2 means
considering a maximum of 2 sentences before or after the anchor segment during
pretraining, K=4 means considering a maximum of 4 sentences before or after.
70-30 split testing and 5-fold split testing are conducted for each baseline.
We list the results in Tables 2-5 and observe the following:

  * •

MathBERT is outperformed by all other models in each task, including base and
large vanilla BERT. Given the improved performance of the model over BERT on 3
tasks related to mathematical text [Shen et al., 2021], and MathBERT being
specialized for mathematical tasks, this is a surprising result.

  * •

BERT-base-uncased model outperforms BERT-large-uncased in each task.

  * •

Comparing CONPONO K=2 and CONPONO K=4, the performances are similar in each
task, with K=2 being marginally better.

  * •

MegatronBERT outperforms all models in all tasks except for discourse
coherence (DC), but its parameter size is larger than most of other baselines
including BERT-based-uncased and CONPONO models.

  * •

CONPONO K=2 model outperforms BERT-base-uncased, RoBERTa-base, and MathBERT on
SP and BSO tasks.

  * •

All baselines perform poorly on the DC task.

  * •

All baselines perform relatively well on the SSP task, but the accuracy
performances are close, with no massive margin even between all-around worst
performer MathBERT and generally good MegatronBERT.

| 70-30 Split | K-fold
---|---|---
| Acc | F1 | AP | ROC | Acc | F1 | AP | ROC
BERT-base-uncased | 0.332 | 0.320 | 0.346 | 0.676 | 0.340$\pm$0.010 | 0.314$\pm$0.007 | 0.349$\pm$0.004 | 0.678$\pm$0.004
BERT-large-uncased | 0.306 | 0.254 | 0.318 | 0.651 | 0.306$\pm$0.006 | 0.269$\pm$0.017 | 0.317$\pm$0.004 | 0.651$\pm$0.005
RoBERTa-base | 0.281 | 0.227 | 0.283 | 0.611 | 0.287$\pm$0.008 | 0.217$\pm$0.011 | 0.291$\pm$0.003 | 0.615$\pm$0.003
MathBERT | 0.258 | 0.248 | 0.256 | 0.571 | 0.259$\pm$0.003 | 0.245$\pm$0.005 | 0.257$\pm$0.003 | 0.574$\pm$0.004
MegatronBERT | 0.647 | 0.646 | 0.729 | 0.893 | 0.647$\pm$0.003 | 0.646$\pm$0.003 | 0.729$\pm$0.002 | 0.894$\pm$0.001
CONPONO K=2 | 0.536 | 0.536 | 0.577 | 0.820 | 0.532$\pm$0.005 | 0.528$\pm$0.005 | 0.574$\pm$0.006 | 0.816$\pm$0.003
CONPONO K=4 | 0.527 | 0.523 | 0.571 | 0.816 | 0.526$\pm$0.008 | 0.522$\pm$0.009 | 0.572$\pm$0.009 | 0.815$\pm$0.004

Table 2: Sentence Position

| 70-30 Split | K-fold
---|---|---
| Acc | F1 | AP | ROC | Acc | F1 | AP | ROC
BERT-base-uncased | 0.535 | 0.530 | 0.543 | 0.550 | 0.517$\pm$0.010 | 0.444$\pm$0.067 | 0.539$\pm$0.001 | 0.546$\pm$0.001
BERT-large-uncased | 0.521 | 0.477 | 0.529 | 0.532 | 0.517$\pm$0.010 | 0.467$\pm$0.069 | 0.528$\pm$0.001 | 0.532$\pm$0.001
RoBERTa-base | 0.531 | 0.529 | 0.537 | 0.546 | 0.531$\pm$0.002 | 0.525$\pm$0.007 | 0.538$\pm$0.001 | 0.546$\pm$0.001
MathBERT | 0.511 | 0.443 | 0.519 | 0.518 | 0.513$\pm$0.001 | 0.492$\pm$0.023 | 0.520$\pm$0.002 | 0.520$\pm$0.002
MegatronBERT | 0.762 | 0.762 | 0.855 | 0.853 | 0.762$\pm$0.001 | 0.762$\pm$0.001 | 0.855$\pm$0.001 | 0.853$\pm$0.001
CONPONO K=2 | 0.671 | 0.671 | 0.741 | 0.742 | 0.670$\pm$0.002 | 0.670$\pm$0.003 | 0.741$\pm$0.002 | 0.741$\pm$0.002
CONPONO K=4 | 0.665 | 0.664 | 0.734 | 0.735 | 0.666$\pm$0.001 | 0.665$\pm$0.001 | 0.734$\pm$0.001 | 0.735$\pm$0.001

Table 3: Binary Sentence Ordering

| 70-30 Split | K-fold
---|---|---
| Acc | F1 | AP | ROC | Acc | F1 | AP | ROC
BERT-base-uncased | 0.529 | 0.527 | 0.534 | 0.541 | 0.522$\pm$0.011 | 0.491$\pm$0.044 | 0.535$\pm$0.004 | 0.539$\pm$0.004
BERT-large-uncased | 0.505 | 0.402 | 0.519 | 0.522 | 0.514$\pm$0.006 | 0.435$\pm$0.056 | 0.520$\pm$0.007 | 0.524$\pm$0.008
RoBERTa-base | 0.517 | 0.498 | 0.529 | 0.539 | 0.515$\pm$0.012 | 0.447$\pm$0.044 | 0.533$\pm$0.004 | 0.541$\pm$0.005
MathBERT | 0.507 | 0.501 | 0.506 | 0.506 | 0.506$\pm$0.004 | 0.494$\pm$0.015 | 0.509$\pm$0.004 | 0.512$\pm$0.005
MegatronBERT | 0.529 | 0.515 | 0.539 | 0.546 | 0.524$\pm$0.003 | 0.514$\pm$0.012 | 0.534$\pm$0.002 | 0.541$\pm$0.003
CONPONO K=2 | 0.534 | 0.532 | 0.538 | 0.545 | 0.533$\pm$0.003 | 0.532$\pm$0.004 | 0.536$\pm$0.003 | 0.544$\pm$0.004
CONPONO K=4 | 0.531 | 0.531 | 0.539 | 0.545 | 0.529$\pm$0.004 | 0.524$\pm$0.005 | 0.539$\pm$0.005 | 0.545$\pm$0.005

Table 4: Discourse Coherence

| 70-30 Split | K-fold
---|---|---
| Acc | F1 | AP | ROC | Acc | F1 | AP | ROC
BERT-base-uncased | 0.695 | 0.694 | 0.756 | 0.763 | 0.693$\pm$0.003 | 0.693$\pm$0.003 | 0.754$\pm$0.002 | 0.762$\pm$0.002
BERT-large-uncased | 0.670 | 0.670 | 0.723 | 0.733 | 0.666$\pm$0.003 | 0.665$\pm$0.003 | 0.718$\pm$0.002 | 0.728$\pm$0.002
RoBERTa-base | 0.649 | 0.648 | 0.706 | 0.713 | 0.651$\pm$0.002 | 0.650$\pm$0.002 | 0.709$\pm$0.002 | 0.715$\pm$0.002
MathBERT | 0.618 | 0.618 | 0.679 | 0.677 | 0.619$\pm$0.006 | 0.619$\pm$0.006 | 0.677$\pm$0.003 | 0.676$\pm$0.004
MegatronBERT | 0.705 | 0.705 | 0.778 | 0.782 | 0.704$\pm$0.003 | 0.703$\pm$0.003 | 0.776$\pm$0.001 | 0.779$\pm$0.002
CONPONO K=2 | 0.687 | 0.686 | 0.761 | 0.756 | 0.685$\pm$0.004 | 0.684$\pm$0.004 | 0.746$\pm$0.002 | 0.754$\pm$0.003
CONPONO K=4 | 0.684 | 0.683 | 0.758 | 0.755 | 0.683$\pm$0.003 | 0.683$\pm$0.003 | 0.743$\pm$0.002 | 0.752$\pm$0.002

Table 5: Sentence Section Prediction

## 6\. Related Work

DiscoEval [Chen et al., 2019] is a suite of evaluation tasks with the purpose
of determining whether sentence representations include information about the
role of a sentence in its discourse context. They build sentence encoders
capable of modelling discourse information via training objectives that make
use of natural annotations from Wikipedia, such as nesting level, section and
article titles, among others.
Other core work [Iter et al., 2020] involves pretraining on both MLM and a
contrastive inter-sentence objective (CONPONO), where they achieve state-of-
the-art benchmarks for five of seven tasks in the DiscoEval suite,
outperforming BERT-Large despite equalling the size of BERT-Base and training
on the same amount of data. BERT-Base pretrained additionally on BSO in place
of CONPONO, and BERT-large, claim the remaining two benchmarks.
Using full encoder-decoder transformer [Vaswani et al., 2017] architecture and
an additional masked attention map which incorporates relationships between
nodes in operator trees (OPTs) of equations [Davila et al., 2016, Davila and
Zanibbi, 2017], MathBERT [Peng et al., 2021] approach pretrains with three
objectives on arXiv data each extracting a specific latent aspect of
information. MLM learns text representations, context correspondence
prediction learns the latent relationship between formula and context, and
masked substructure prediction learns semantic-level structure of formulas by
predicting parent and child nodes in OPTs. This model obtains state-of-the-art
results in retrieval-based math tasks, but their model is unpublished. We use
an alternative MathBERT [Shen et al., 2021] that is trained on a larger corpus
of mathematical text, and demonstrates mathematical proficiency over regular
BERT.
A data extraction pipeline [Ferreira and Freitas, 2020] collects 20k entries
related to mathematical proofs from the ProofWiki
website333https://proofwiki.org/wiki/MainPage, such as definitions, lemmas,
corollaries, and theorems. They evaluate BERT and SciBERT by fine-tuning on a
pairwise relevance classification task with their NL-PS dataset, where they
classify if one mathematical text is related to another. As we highlight in
the Introduction, physics and mathematics literature differ in their
overarching considerations, and more specifically, the unstructured informal
physics Wikipedia explanations that we present in our data naturally differ
from the structured proofs present in NL-PS.
Their work builds on previous efforts applying NLP to general mathematics. One
early approach [Zinn, 2003] proposes proof representation structures via
discourse representation theory, including a prototype for generating formal
proofs from informal mathematical discourse. Another approach [Cramer et al.,
2009] focuses on development of a controlled natural language for mathematical
texts which is compatible with existing proof verification software. Since
these early developments, natural language-based problem solving and theorem
proving have progressed significantly, but are still below human-level
performance. For example, efforts towards building datasets for evaluating
math word problem solvers [Huang et al., 2016] concludes the task as a
significant challenge, with more recent large-scale dataset construction and
evaluation work [Amini et al., 2019, Miao et al., 2021] confirming that model
performance is still well below the gold-standard. This difficulty extends to
approaches involving pre-university math problems and geometric quantities
[Matsuzaki et al., 2017, Lu et al., 2021]. For automated theorem proving and
mathematical reasoning, various datasets and accompanying approaches have been
proposed [Kaliszyk et al., 2017, Bansal et al., 2019] including more recent
work with equational logic [Piepenbrock et al., 2021] and language models
[Rabe et al., 2020, Han et al., 2021].
A dataset construction approach and accompanying heuristic search for
automating small physics derivations has been developed [Meadows and Freitas,
2021] which allows published results in modern physics to be converted into a
form interpretable by a computer algebra system [Meurer et al., 2017], which
then accommodates limited informal mathematical exploration. Detailed physics
derivation data is scarce, and others have tackled such issues via synthetic
data [Aygün et al., 2020], though not yet in physics. Reinforcement learning
has been employed [Luo and Liu, 2018] to solve differential equations in
nuclear physics with a template mapping method, and proof discovery and
verification has been explored in relativity [Govindarajalulu et al., 2015].
Theorem proving and derivation automation in physics remains elusive, with
detailed discussions available in the literature [Kaliszyk et al., 2015,
Davis, 2019].
With language models demonstrating logical capabilities with respect to type
inference, missing assumption suggestion, and completing equalities [Rabe et
al., 2020], as well as state-of-the-art performance in math retrieval and
tasks related to equation-context correspondence [Peng et al., 2021], we
believe our present work will contribute towards physics natural language /
equational reasoners capable of generating coherent mathematical explanations
and derivations.

## 7\. Conclusion

Within the domain of physics, we present 2 parent datasets for general use and
4 specific datasets corresponding to discourse evaluation tasks [Chen et al.,
2019] collectively referred as PhysNLU. The presented data frequently features
equations, formulae, and mathematical language. Our analysis reveals that
concepts related to quantum mechanics are most commonly discussed as
determined by Wikipedia article category, that equations related to initial
and boundary conditions are the most frequently considered when considering
near-exact string matching, and we report the proportion of sentences and
dataset examples which contain equations and mathematical terms identified by
annotation frameworks native to Wikipedia. Finally, we present baseline
results for popular non-mathematical language models and demonstrate that,
despite expensive pretraining efforts and specialised training objectives for
learning various aspects of mathematical text, such efforts do not improve the
performance of language models in tasks related to sentence ordering,
position, and recognising whether physics explanations are coherent. Future
work will involve developing objectives which aid performance in this regard.

## Acknowledgements

This work was partially funded by the SNSF project NeuMath (200021_204617).

## 8\. Bibliography

## References

  * Amini et al., 2019 Amini, A., Gabriel, S., Lin, P., Koncel-Kedziorski, R., Choi, Y., and Hajishirzi, H.  (2019).  Mathqa: Towards interpretable math word problem solving with operation-based formalisms.  arXiv preprint arXiv:1905.13319.
  * Aygün et al., 2020 Aygün, E., Ahmed, Z., Anand, A., Firoiu, V., Glorot, X., Orseau, L., Precup, D., and Mourad, S.  (2020).  Learning to prove from synthetic theorems.  arXiv preprint arXiv:2006.11259.
  * Bansal et al., 2019 Bansal, K., Loos, S., Rabe, M., Szegedy, C., and Wilcox, S.  (2019).  Holist: An environment for machine learning of higher order logic theorem proving.  In International Conference on Machine Learning, pages 454–463. PMLR.
  * Chen et al., 2019 Chen, M., Chu, Z., and Gimpel, K.  (2019).  Evaluation benchmarks and learning criteria for discourse-aware sentence representations.  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 649–662, Hong Kong, China, November. Association for Computational Linguistics.
  * Coffey and Kalmykov, 2012 Coffey, W. and Kalmykov, Y. P.  (2012).  The Langevin equation: with applications to stochastic problems in physics, chemistry and electrical engineering, volume 27.  World Scientific.
  * Cramer et al., 2009 Cramer, M., Fisseni, B., Koepke, P., Kühlwein, D., Schröder, B., and Veldman, J.  (2009).  The naproche project controlled natural language proof checking of mathematical texts.  In International Workshop on Controlled Natural Language, pages 170–186. Springer.
  * Davila and Zanibbi, 2017 Davila, K. and Zanibbi, R.  (2017).  Layout and semantics: Combining representations for mathematical formula search.  In Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval, pages 1165–1168.
  * Davila et al., 2016 Davila, K., Zanibbi, R., Kane, A., and Tompa, F. W.  (2016).  Tangent-3 at the ntcir-12 mathir task.  In NTCIR.
  * Davis, 2019 Davis, E.  (2019).  Proof verification technology and elementary physics.  In Algorithms and Complexity in Mathematics, Epistemology, and Science, pages 81–132. Springer.
  * Devlin et al., 2019 Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  (2019).  Bert: Pre-training of deep bidirectional transformers for language understanding.  In NAACL.
  * Ferreira and Freitas, 2020 Ferreira, D. and Freitas, A.  (2020).  Natural language premise selection: Finding supporting statements for mathematical text.
  * Gala et al., 2021 Gala, S., Galakhov, E., Ragusa, M. A., and Salieva, O.  (2021).  Beale–kato–majda regularity criterion of smooth solutions for the hall-mhd equations with zero viscosity.  Bulletin of the Brazilian Mathematical Society, New Series, pages 1–13.
  * Govindarajalulu et al., 2015 Govindarajalulu, N. S., Bringsjord, S., and Taylor, J.  (2015).  Proof verification and proof discovery for relativity.  Synthese, 192(7):2077–2094.
  * Han et al., 2021 Han, J. M., Rute, J., Wu, Y., Ayers, E. W., and Polu, S.  (2021).  Proof artifact co-training for theorem proving with language models.  arXiv preprint arXiv:2102.06203.
  * Huang et al., 2016 Huang, D., Shi, S., Lin, C.-Y., Yin, J., and Ma, W.-Y.  (2016).  How well do computers solve math word problems? large-scale dataset construction and evaluation.  In Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 887–896.
  * Ibsen-Jensen et al., 2015 Ibsen-Jensen, R., Chatterjee, K., and Nowak, M. A.  (2015).  Computational complexity of ecological and evolutionary spatial dynamics.  Proceedings of the National Academy of Sciences, 112(51):15636–15641.
  * Iter et al., 2020 Iter, D., Guu, K., Lansing, L., and Jurafsky, D.  (2020).  Pretraining with contrastive sentence objectives improves discourse performance of language models.
  * Kaliszyk et al., 2015 Kaliszyk, C., Urban, J., Siddique, U., Khan-Afshar, S., Dunchev, C., and Tahar, S.  (2015).  Formalizing physics: automation, presentation and foundation issues.  In International Conference on Intelligent Computer Mathematics, pages 288–295. Springer.
  * Kaliszyk et al., 2017 Kaliszyk, C., Chollet, F., and Szegedy, C.  (2017).  Holstep: A machine learning dataset for higher-order logic theorem proving.  arXiv preprint arXiv:1703.00426.
  * Liu et al., 2019 Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., and Stoyanov, V.  (2019).  Roberta: A robustly optimized bert pretraining approach.  ArXiv, abs/1907.11692.
  * Lu et al., 2021 Lu, P., Gong, R., Jiang, S., Qiu, L., Huang, S., Liang, X., and Zhu, S.-C.  (2021).  Theorem-aware geometry problem solving with symbolic reasoning and theorem prediction.
  * Luo and Liu, 2018 Luo, M. and Liu, L.  (2018).  Automatic derivation of formulas using reforcement learning.  arXiv preprint arXiv:1808.04946.
  * Matsuzaki et al., 2017 Matsuzaki, T., Ito, T., Iwane, H., Anai, H., and Arai, N. H.  (2017).  Semantic parsing of pre-university math problems.  In Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 2131–2141.
  * Meadows and Freitas, 2021 Meadows, J. and Freitas, A.  (2021).  Similarity-based equational inference in physics.  Physical Review Research, 3(4), Oct.
  * Merity et al., 2016 Merity, S., Xiong, C., Bradbury, J., and Socher, R.  (2016).  Pointer sentinel mixture models.  arXiv preprint arXiv:1609.07843.
  * Meurer et al., 2017 Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., et al.  (2017).  Sympy: symbolic computing in python.  PeerJ Computer Science, 3:e103.
  * Miao et al., 2021 Miao, S.-Y., Liang, C.-C., and Su, K.-Y.  (2021).  A diverse corpus for evaluating and developing english math word problem solvers.  arXiv preprint arXiv:2106.15772.
  * Miller, 2021 Miller, E.  (2021).  A survey of geometric constraints on the blowup of solutions of the navier–stokes equation.  Journal of Elliptic and Parabolic Equations, pages 1–11.
  * Peng et al., 2021 Peng, S., Yuan, K., Gao, L., and Tang, Z.  (2021).  Mathbert: A pre-trained model for mathematical formula understanding.
  * Piepenbrock et al., 2021 Piepenbrock, J., Heskes, T., Janota, M., and Urban, J.  (2021).  Learning equational theorem proving.  arXiv preprint arXiv:2102.05547.
  * Pizzocchero, 2021 Pizzocchero, L.  (2021).  On the global stability of smooth solutions of the navier–stokes equations.  Applied Mathematics Letters, 115:106970.
  * Rabe et al., 2020 Rabe, M. N., Lee, D., Bansal, K., and Szegedy, C.  (2020).  Mathematical reasoning via self-supervised skip-tree training.  arXiv preprint arXiv:2006.04757.
  * Shen et al., 2021 Shen, J. T., Yamashita, M., Prihar, E., Heffernan, N., Wu, X., Graff, B., and Lee, D.  (2021).  Mathbert: A pre-trained language model for general nlp tasks in mathematics education.  arXiv preprint arXiv:2106.07340.
  * Shoeybi et al., 2019 Shoeybi, M., Patwary, M. A., Puri, R., LeGresley, P., Casper, J., and Catanzaro, B.  (2019).  Megatron-lm: Training multi-billion parameter language models using model parallelism.  ArXiv, abs/1909.08053.
  * Smith and Fleck, 2017 Smith, R. W. and Fleck, C.  (2017).  Derivation and use of mathematical models in systems biology.  In Pollen Tip Growth, pages 339–367. Springer.
  * Vaswani et al., 2017 Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I.  (2017).  Attention is all you need.
  * Zinn, 2003 Zinn, C.  (2003).  A computational framework for understanding mathematical discoursexy.  Logic Journal of IGPL, 11(4):457–484.


# Controlled coherent-coupling and dynamics of exciton complexes in a MoSe2
monolayer

Aleksander Rodek1, Thilo Hahn2, James Howarth3, Takashi Taniguchi4, Kenji
Watanabe5, Marek Potemski1,6,7, Piotr Kossacki1, Daniel Wigger8, and Jacek
Kasprzak1,9,10 1Faculty of Physics, University of Warsaw, ul. Pasteura 5,
02-093 Warszawa, Poland  2Institute of Solid State Theory, University of
Münster, 48149 Münster, Germany  3National Graphene Institute, University of
Manchester, Booth St E, M13 9PL United Kingdom  4International Center for
Materials Nanoarchitectonics, National Institute for Materials Science, 1-1
Namiki, Tsukuba 305-0044, Japan  5Research Center for Functional Materials,
National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
6Laboratoire National des Champs Magnétiques Intenses, CNRS-UGA-UPS-INSA-EMFL,
25 Av. des Martyrs, 38042 Grenoble, France  7CENTERA Labs, Institute of High
Pressure Physics, PAS, 01- 142 Warsaw, Poland  8School of Physics, Trinity
College Dublin, Dublin 2, Ireland  9Univ. Grenoble Alpes, CNRS, Grenoble INP,
Institut Néel, 38000 Grenoble, France  10Walter Schottky Institut and TUM
School of Natural Sciences, Technische Universität München, 85748 Garching,
Germany<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantifying and controlling the coherent dynamics and couplings of optically
active excitations in solids is of paramount importance in fundamental
research in condensed matter optics and for their prospective optoelectronic
applications in quantum technologies. Here, we perform ultrafast coherent
nonlinear spectroscopy of a charge-tunable MoSe2 monolayer. The experiments
show that the homogeneous and inhomogeneous line width and the population
decay of exciton complexes hosted by this material can be directly tuned by an
applied gate bias, which governs the Fermi level and therefore the free
carrier density. By performing two-dimensional spectroscopy, we also show that
the same bias-tuning approach permits us to control the coherent coupling
strength between charged and neutral exciton complexes.

  * January 2023

##  1 Introduction

The ability to isolate monolayer flakes of semiconducting transition-metal
dichalcogenides (TMDs), such as MoSe2, and the discovery of their enhanced
light-matter interaction at the monolayer limit a decade ago [1, 2]
established a novel benchmark for semiconductor optics. The combination of
heavy effective masses and spin-valley locking – both inherent in the
monolayer’s band structure – with a large out-of-plane dielectric contrast
result in particularly strong Coulomb interactions among free carriers hosted
by TMD monolayers [3, 4]. For the same reason, the exciton transitions display
a non-hydrogenic excitation spectrum [5] and have large binding energies and
oscillator strengths [6, 7], such that they dominate the optical response of
TMD materials even at room temperature.

Further to that, the ability to stack different sorts of layered van der Waals
materials, such as TMDs, hexagonal boron nitride (hBN), graphene or graphite,
into van der Waals heterostructures [8, 9], facilitates the fabrication of
optoelectronic devices [10, 11, 12, 13, 14]. Sandwiching layered TMDs or other
materials between high quality hBN diminishes the electronic disorder and
protects them from degradation by ambient agents. In general, the hBN-
encapsulation procedure reduces the exciton’s spectral inhomogeneous
broadening $\hbar\sigma$, thus improving the optical performance of TMDs [15,
16, 17, 18]. It even enables reaching the homogeneous limit for the exciton’s
spectral line shape [19, 20]. This has been a milestone for monolayer
spectroscopy, as it enhances the visibility of phenomena previously blurred by
different types of inhomogeneity, with the appearance of a novel destructive
photon echo being an example [21].

We are employing a four-wave mixing (FWM) spectroscopy approach to study the
coherent optical properties of TMD exciton complexes. This method is recently
gaining more popularity for the investigation of layered semiconductors:
Besides the fundamental studies of the exciton coherence and population
dynamics in different TMD systems [22, 23, 19] it was also used to probe
intervalley scattering processes [24], valley coherence times [25], and
mechanisms of exciton broadening [26]. In particular it allowed for first
observations of biexcitons [27, 28] or coupling between different exciton
states [28] in various TMD materials.

Electron beam lithography and metal deposition can further be employed
complement van der Waals heterostructures with electronic contacts and gates
[29, 30]. In this way, one can inject carriers into TMD monolayers, permitting
us to control the free electron density $n_{\rm e}$ in our sample. This is
crucial, not only for optoelectronic applications, but also when exploring the
fundamental physics of correlated many-body systems in solids; landmark
examples being the recent revelation of Wigner crystals in a gate-tunable
MoSe2 monolayer [31] or optical sensing of the quantum Hall effect in graphene
[32]

The presence of free carriers modifies the excitonic optical response
substantially [33, 34], similarly as was shown and partially explained for
semiconductor quantum wells [35, 36, 37, 38]. Firstly, it is supposed that
free carriers affect the relative absorption (oscillator or effective dipole
strength) and the possible couplings between neutral excitons (X) and
(negatively) charged excitons (X-), also called trions. Secondly, the
exciton’s dephasing and therefore the spectral line shape should be sensitive
to the electron density $n_{\rm e}$ [39], in an analogous way as it depends on
the total exciton density, through the mechanisms called excitation induced
dephasing [40, 41, 20], sometimes described by a local field effect [42, 43].
Thirdly, with increasing $n_{\rm e}$ the exciton complexes are screened more
efficiently from intrinsic disorder, which in turn should reduce their
inhomogeneous broadening [3, 44], further impacting the excitons’ radiative
decay rates [45]. Here, we address the basic interplay between free carriers
$n_{\rm e}$, excitons X, and trions X- by performing coherent nonlinear
spectroscopy [46, 47] of a gated and hBN-encapsulated MoSe2 monolayer.

##  2 Sample & Experiment

A microscopy image of our device is presented in Fig. 1(a), where the MoSe2
monolayer flake (marked in red) is encapsulated between thin bottom (dashed
white) and top (solid white) hBN films. The bottom graphite flake (dashed
black) lies under the hBN spacer, while the top few-layer graphene flakes
(solid black) are adjacent to the MoSe2 as shown in the layer schematic in the
top right inset. The graphite and graphene operate as bottom and top gates,
respectively. Details regarding the fabrication of our device are provided in
the Supporting Information (SI). We have fabricated and investigated two
samples of the same design, which show very similar results throughout the
studied properties.

The device permits us to control the free electron density $n_{\rm e}$ by
applying a gate bias $U$ via the gold bands, as shown in yellow in Fig. 1(a),
connecting with both graphite and graphene layers acting as electrodes. When
increasing the bias and therefore $n_{\rm e}$ we observe modifications of the
optical response of the monolayer, both in the transition energy and the
spectral line shapes, as shown in the white-light reflection contrast spectra
presented in Fig. 1(b). We note that by varying $n_{\rm e}$ we can alter the
relative spectral intensities of X and the X-, as quantified in (c). We see
that already the linear optical response, as observed in reflectance, is tuned
when increasing $n_{\rm e}$, as the oscillator strength shifts from X to X- –
we will come back to this point in the last section.

Figure 1: Investigated MoSe2 device and optical characterization. (a) Optical
image of the heterostructure. The individual components are marked by colored
lines and the white dots mark the positions of the laser beams for the FWM
measurements and the inset shows the layer structure of the sample. (b)
Reflectance contrast spectrum of the sample as a function of the applied
voltage showing the neutral (X) and negatively-charged (X-) exciton
transitions. (c) Individual peak intensities X in blue and X- in red from (b)
and their ratio in green. (d) Imaging of the time-integrated FWM at
$\tau_{12}=1.2$ ps. The black area in the lower-left corner of the image is an
experimental artefact generated when the reference beam impinges on the highly
reflecting metallic contact $\sim 10$µm below the edge of the MoSe2 monolayer.
The outline of the MoSe2 flake is marked in red and the line where the
reference beam $\mathcal{E}_{\rm ref}$ crosses the bottom edge of the top hBN
flake in white.

To infer the coherent nonlinear response of excitons, we perform FWM
microscopy, similarly as in our recent works [22]. We use three, co-linearly
polarized pulses of 150 femtosecond duration, here labeled as ${\cal
E}_{1,2,3}$, that can be tuned into resonance with the excitonic transitions.
The beams are focused down to a diffraction-limited spot on the sample’s
surface, placed in a helium flow cryostat, setting the temperature at $T=8$ K
for all experiments. The complex-valued FWM signal, generated in a standard,
so-called photon echo, configuration is retrieved in reflectance by combining
optical heterodyning and spectral interferometry with a reference pulse [48].
The reference is focussed through the same objective (${\rm NA}=0.65$) below
${\cal E}_{1,2,3}$ and is reflected from the neighboring hBN without MoSe2, as
shown in Fig. 1(a) by the white dots.

To investigate the exciton’s coherence dynamics we measure the FWM amplitude
as a function of time delay $\tau_{12}$ between the first and the second
arriving pulses, while fixing $\tau_{23}$. Conversely, the exciton’s
population dynamics is probed in FWM when varying the delay $\tau_{23}$, while
fixing $\tau_{12}$. These pulse sequences are depicted in Fig. 2(a) and Fig.
3(a), respectively. We note that for all experiments we keep the same
excitation conditions of $0.3$ µW average power for each of the beams,
generating an exciton density of a few $10^{10}$cm-2 [22, 20]. To confirm the
expected linear scalings of the FWM amplitude with the pulse areas of ${\cal
E}_{1,2,3}$ for small powers, we checked the FWM power dependence (see SI Fig.
S1), also monitoring that the excitation induced dephasing can be neglected
for this range of exciton densities [20].

To characterize the overall device, we set $U=-0.5\,$V and perform a FWM
amplitude mapping across the MoSe2 flake for $\tau_{12}=1.2\,$ps, by scanning
the objective’s position. The resulting time-integrated FWM image is presented
in Fig. 1(d) (for detailed analysis of this imaging see SI Figs S2, S3, and
S4). The dominant signals reflect the third-order optical susceptibility of
the MoSe2 monolayer (red line). Note, that the generated FWM is weaker in the
upper part of the flake, which is due to the reduced reflectance of the
reference pulse $\mathcal{E}_{\rm ref}$ when impinging on the top hBN. This
crossover between low and high FWM amplitude regions (marked in white)
faithfully follows the shape of the top hBN’s edge when $\mathcal{E}_{\rm
ref}$ leaves the top hBN, while the other pulses $\mathcal{E}_{1,2,3}$ still
probe the fully encapsulated MoSe2 monolayer.

In time-resolved FWM, a typical case presented in Fig. 2(b), without
introducing free carriers we observe a more or less pronounced photon echo
[26, 22, 23], which is a consequence of inhomogeneous broadening
$\hbar\sigma$. The latter induces rephasing of the signal for $t=\tau_{12}$,
such that normally the echo is aligned with the diagonal in the
$(t,\tau_{12})$-space, as depicted as dashed diagonal in (b). However, for
short delays $\tau_{12}<0.5\,$ps we spot deviations from the typical echo. In
this range the signal is retarded in time, such that the transient is not
following the diagonal, and is instead curved, into a comet-like shape. The
echo distortion is readily reproduced, when taking into account the local
field effect [21], representing an effective exciton-exciton interaction, and
is described in leading order in the local field coupling by

$\displaystyle|p_{\rm FWM}(t,\,\tau_{12})|$ $\displaystyle\sim$
$\displaystyle\Theta(\tau_{12})\Theta(t)te^{-\beta(t+\tau_{12})-\frac{1}{2}\sigma^{2}(t-\tau_{12})^{2}}$
(1) $\displaystyle+$
$\displaystyle\Theta(-\tau_{12})\Theta(t+\tau_{12})(t+\tau_{12})e^{-\beta(t-\tau_{12})-\frac{1}{2}\sigma^{2}(t-\tau_{12})^{2}}\,.$

The fitted $(t,\,\tau_{12})$-dynamics is presented in Fig. 2(c), from which we
extract homogeneous $\beta$ and inhomogeneous $\sigma$ dephasing rates and
corresponding full-width at half-maxima (FWHM) spectral broadenings
$\hbar\beta=0.86\,$meV and
$2\sqrt{2\ln(2)}\hbar\sigma=\hbar\tilde{\sigma}=3.9\,$meV, respectively. By
inspecting spectrally-resolved FWM mappings, we reveal typical fluctuations of
the exciton’s central energy ($\approx\pm 5\,$meV) and line widths
($\approx\pm 1\,$meV), as shown in Fig. S2 in the SI, attributed to varying
strain generated during the sample assembly [49, 50, 51]. We therefore fix the
excitation spot for the entire experiment, as marked in Fig. 1(a), to probe an
area of relatively small $\sigma$ within the distribution. The same mapping is
performed for the time-resolved FWM signal, which allows to extract
homogeneous and inhomogeneous broadenings, summarized in SI Fig. S3. We find
an expected and pronounced correlation between $\sigma$ and the FWHM of the
FWM spectra, as shown in SI Fig. S4. There, we also find a correlation between
the exciton transition energy and the inhomogeneous broadening.

##  3 Coherence dynamics of neutral & charged excitons

We now proceed to the investigation of the exciton homogeneous linewidth
$\hbar\beta$ depending on the induced electron density $n_{e}$. For this
purpose, we center the excitation at the exciton energy, which for the
investigated area equals 1640 meV, and measure the photon echo as a function
of the gate bias $U$. Exemplary exciton coherence dynamics, i.e., time-
integrated FWM amplitudes as a function of $\tau_{12}$, are shown in Fig. 2(d)
(blue dots). Due to the local field, the measured dynamics deviate from simple
combined exponential and Gaussian decays according to Eq. (1). The
corresponding fits are shown as pale blue lines and reproduce the experiments
well. The data clearly indicate a growth of the dephasing rate when increasing
the electron density (top to bottom). All fitted homogeneous broadenings are
presented in (f, blue filled circles) as function of the gate bias, which
quantifies the impression of a strong rise of $\hbar\beta$ with increasing
$n_{e}$. The fitted dynamics also yield $\hbar\tilde{\sigma}$ as a function of
$n_{e}$ (f, blue open circles). Note, that for the fits we used the entire
signal dynamics in real and delay time from which we also extract the time-
integrated data in Fig. 2(d, e). Interestingly, the inhomogeneous broadening
is suppressed efficiently, when injecting more electrons. The FWM transient
thus evolves from the photon echo to free-induction decay like, as we
explicitly show in SI Fig. S5. This strikingly shows that the excitons become
less sensitive to the underlying disorder, which can be static due to the
strain or dynamic via fluctuating charges. With increasing $n_{e}$, the
disorder thus gets screened more efficiently and the applied voltage can
neutralize fluctuating charges [44], both mechanisms making the system less
inhomogeneous. We have previously found that the reduced inhomogeneity goes
hand-in-hand with a shortening of the excitons’ radiative lifetime [19]. The
small rise of the X’s inhomogeneous broadening at $U>4$ V is an artefact from
the weak optical response which results in an increased uncertainty of the
fit. This is also reflected by the significant increase of the error bars
(shaded area).

Figure 2: X and X- dephasing versus applied gate bias. (a) Scheme of the
three-pulse FWM sequence probing coherence dynamics via $\tau_{12}$. (b)
Measured time-resolved FWM amplitude showing a clear photon echo. (c)
Theoretical fit of (b) employing the local field model with the fitted
homogeneous and inhomogeneous dephasing rates $\beta$ and $\sigma$ as given in
the plot. (d, e) Exemplary time-integrated FWM amplitude dynamics as a
function of $\tau_{12}$ for different gates biases. (d) For X and (e) for X-.
(f) Homogeneous (filled circles) and inhomogeneous (open circles) FWHM line
widths extracted from (d, e) and time-resolved data as a function of gates
bias with X in blue and X- in red. The shaded areas show the uncertainties.

The second resonance, occurring at 1620 meV is the charged exciton transition
X-, which we selectively address by centering the laser pulse spectrum to this
energy. We carry out the same routine as before to determine the charged
exciton’s homogeneous and inhomogeneous broadening and their dependence on
$n_{\rm e}$. The results are presented in Fig. 2(e) and (f) in red colors.
When comparing X’s and X-’s homogeneous broadenings in (f), we observe that
for the same electron density the X broadening is significantly larger than
that of X-. This finding is similar to experiments performed on non-
intentionally doped samples [25]. Following an intermediate drop with
increasing $n_{\rm e}$, we eventually also observe an increase of the charged
exciton’s homogeneous broadening accompanied by a decrease in inhomogeneous
broadening (open red circles). The increased line width of both exciton
complexes is attributed to the dephasing due to interactions with the Fermi-
sea of electrons, similarly as in past studies on semiconductor quantum wells
[40, 39].

With this, we evidence the control the inhomogeneous broadening and excitonic
dephasings and thus the spectral line shape of the optical transitions, via a
tunable gate bias introducing free carriers into the MoSe2 monolayer. We
distinguish this line broadening mechanism from the ones previously
investigated in TMDs, i.e., phonon-induced [26, 22, 23, 52, 49] and
excitation-induced dephasing [19, 20].

At this point, we remark that the measured behavior of dephasings does not
match with a simple picture of the dipole (oscillator) strength transfer from
X to X- when increasing $n_{\rm e}$. Assuming that the dephasing was entirely
governed by the radiative lifetime, a reduction for the X (an increase for X-)
dipole strength should increase (reduce) the lifetime and consequently the
dephasing time. However, the opposite trend is observed in Fig. 2(f). Next to
the suppression of $\sigma$ via screening, the second possible cause for this
behavior is the presence of other decay channels that increase with a growing
$n_{\rm e}$, which would result in an accelerated dephasing while remaining
close to a lifetime-limited condition. To shed light on the dominating decay
mechanism, we thus measure FWM as a function of $\tau_{23}$, monitoring the
coherent population dynamics of X and X-.

##  4 Population dynamics of neutral & charged excitons

In Fig. 3(b) we present the time-integrated FWM signals of X as a function of
$\tau_{23}$ for selected gate biases (amplitude as red and phase as blue
dots). As the FWM response is measured in a coherent fashion, we can retrieve
its amplitude and phase, improving the insight into the involved decay
processes affecting the excitonic populations that occur. A natural choice to
describe this dependence is to consider a coherent superposition of several
exponential decays [53, 19] via

$\displaystyle S_{\rm FWM}(\tau_{23},t)$ $\displaystyle=$ $\displaystyle
A_{\rm off}\exp(i\varphi_{\rm off})+\exp(i\varphi_{\rm dr}t)\Bigg{\\{}A_{\rm
nr}\exp\left(i\varphi_{\rm nr}-\frac{\tau_{23}^{2}}{\tau_{0}^{2}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n}A_{n}\left[1+\mathrm{erf}\left(\frac{\tau_{23}}{\tau_{0}}-\frac{\tau_{0}}{2\tau_{n}}\right)\right]\exp\left(i\varphi_{n}+\frac{\tau_{0}^{2}}{2\tau_{n}^{2}}-\frac{\tau_{23}}{\tau_{n}}\right)\Bigg{\\}}$

where $t$ is real time of the FWM transient, ($A_{\rm nr},\varphi_{\rm nr}$)
are amplitude and phase of the two-photon absorption,
($A_{n},\tau_{n},\varphi_{n}$) are the amplitude, characteristic time, and
phase of a given decay processes, $\varphi_{\rm dr}$ is the phase drift during
the measurement, ($A_{\rm off}$, $\varphi_{\rm off}$) are the amplitude and
phase of the complex offset, and $\tau_{0}$ is the pulse duration of the
femto-second laser of around 150 fs.

As exemplarily represented by the fitted pale curves in Fig. 3(b, top), we
find that for gate bias values $U\leq 0$, which corresponds to the charge-
neutrality regime, i.e., no free electrons, the exciton’s population dynamics
can be fitted with a bi-exponential decay with characteristic timescales
$\tau_{1}<1\,$ps and $\tau_{2}\approx 7\,$ps. We plot the extracted timescales
of the relaxation processes in Fig. 3(c). This result is in agreement with
previous studies [22, 53], with the faster component being attributed to the
exciton’s decay, which can be due to several mechanisms: radiative
recombination, non-radiative scattering into various momentum-dark states and
localized states generated by the disorder, and finally X to X- conversion.
The longer decay stems from the relaxation of thermalized higher $k$-vector
exciton states, which scatter back into the light-cone (i.e. states with
momentum values $|k|<n\omega/c$) through non-radiative processes and
subsequently contribute to the FWM signal.

With increasing gate bias $U$ and thus $n_{\rm e}$, the description of the
observed population dynamics requires the introduction of a third, slower
relaxation process $\tau_{3}$. Earlier FWM experiments on samples naturally
doped with electrons also indicate the presence of this long-lived component
and attribute it to a decaying population of excitons with spin-forbidden
transitions [53]. As the disorder is screened more efficiently with increasing
$n_{\rm e}$, we tentatively suggest that the observed faster dynamics could be
also due to a weaker exciton localization.

The most pronounced change in the population dynamics, caused by the increase
of $n_{\rm e}$, can be observed during the first few picoseconds in Fig. 3(a).
For positive $U$ in Fig. 3(c), we observe a rapid increase of the fastest
relaxation rate. This behaviour, fortified by the shift of the FWM phase,
occurs in the same voltage range as the previously mentioned shortening of the
exciton dephasing time (see blue dots in Fig. 2(f)). Interestingly, when
comparing the retrieved $\tau_{1}$ values with $T^{\rm X}_{2}/2=1/(2\beta)$
(gray line in Fig. 3(c, bottom)), we see the same trend with increasing bias.
In particular, we find that for large gate biases the measurements approach
$\tau_{1}=T^{\rm X}_{2}/2$ for $U\approx 3\,$V and therefore the special
situation of lifetime-limited dephasing. This result is also in line with the
previous finding that the exciton looses nearly all inhomogeneous broadening
for $U>3\,$V and therefore the dominant dephasing mechanism for smaller $U$.

In Fig. 3(d, e) we present analogous data from the measurements where the
excitation energy is tuned to the charged exciton transition X-. In this case,
we only observe two distinct relaxation processes [53, 54] contributing to the
investigated range of $\tau_{23}$. X- is the lowest energy state and – in
contrast to X – cannot relax to other excitonic states. Consequently, the
initial decay of X- of around 6 ps, reflects its radiative and non-radiative
recombination. We note here also, that the simple extension of the light-cone
escape process as for the neutral excitons is now not valid for X-. The free
carrier leftover after the recombination allows for the additional momentum
transfer and recombination of X- from non-zero k-vectors [55, 56].

Figure 3: X and X- population dynamics versus gate bias. (a) Scheme of the
three-pulse FWM probing the population dynamics. (b) Exemplary FWM amplitude
(red dots) and phase (blue dots) dynamics of X as a function of $\tau_{23}$,
with fitted curves in pale colors, according to Eq. 1. (c) Extracted
characteristic decay times of the three identified relaxation channels as a
function of gate bias. The gray curve shows the blue data from Fig. 2(f) in
the form $T_{2}^{\rm X}=1/\beta$. (d, e) Same as (b, c) but for X-. Figure 4:
Quantum beat of the X-X- Raman coherence revealed in the population dynamics.
(a) Amplitude of the FWM spectra of X and X- in blue together with the
amplitude of reference spectrum in red. (b) Spectral dynamics of the FWM
amplitude as a function of $\tau_{23}$ showing a clear beating of X-. (c)
Theoretical simulation of (b). (d) Dynamics of the X- single after spectral
integration over the respective peak. Experiment as green dots and simulation
as pale line.

##  5 Controlled coherent coupling between neutral & charged excitons

So far, we have investigated the coherence dynamics of X and X- separately, by
selectively addressing the respective resonances. Conversely, when the two
complexes are excited simultaneously, both exciton species coexist and
interact with each other. To trigger the interplay between X and X- and reveal
their coupling, we excite them in tandem, as shown spectrally in Fig. 4(a) by
the laser spectrum in red and the FWM spectrum in blue. The measured
spectrally resolved population dynamics are presented in (b). Interestingly,
for the initial delays $\tau_{23}<0.5\,$ps, the FWM displays a beating
particularly pronounced on the X- resonances. The first driving beam induces
the Raman coherence between X and X- as depicted schematically in (b). Due to
the spectral splitting $\delta$ between X and X-, the density grating
generated by the second beam oscillates with the period
$2\pi\hbar/\delta\approx 0.14\,$ps. The FWM released by the probe therefore
shows the Raman quantum beat [57], as confirmed by the simulation in (c), and
directly shown in (d) as time traces after spectrally integrating over the X-
peak. This result indicates that X and X- are Raman-coupled, as they share a
common ground state. Note, that the simulation (pale green line) only
considers a single exciton decay channel, while in the experiment (green dots)
the interplay between different relaxation paths leads to the rising signal
for $\tau_{23}>0.5\,$ps (see Fig. 3).

A beating of a similar origin is observed in the coherence dynamics, probed by
the $\tau_{12}$-dependence, as shown in SI Fig. S6. When coherently coupled,
the first-order absorptions of X and X-, created by the leading pulse, are
mutually converted into the FWM of X- and X, respectively, by the following
two pulses, as schematically shown in Fig. 5(a). To pinpoint this phenomenon
we perform two-dimensional (2D) FWM spectroscopy [58, 59]. In this approach,
originating from nuclear magnetic resonance spectroscopy, Fourier
transformations are performed along two time axes: the direct time axis and
the indirect delay axis [60, 61, 62]. In our case, the transform along the
direct axis is automatically performed by the optical spectrometer, yielding
the FWM spectra. Conversely, the Fourier transform of the indirect axis has to
be recovered from the $\tau_{12}$-sequence. To track the FWM phase when
varying $\tau_{12}$, we apply the phase-referencing scheme [63], which
overcomes the need for an active phase-stabilization and permits us to
accurately perform the Fourier transform along $\tau_{12}$, yielding the
energy axis $\hbar\omega_{12}$.

Figure 5: Phase-referenced two-dimensional FWM spectroscopy of X and X-. (a)
The pulse sequence used in the 2D FWM experiment. $\tau_{12}$ is scanned,
while the other delay is set to $\tau_{23}=0.13\,$ps maximizing the coherent
coupling. (b) Examples of 2D FWM spectra for different gate biases as given in
the plots. (c) Integrated peak amplitudes of the four peak X, X, XX-, and X-X
(marked in (b)) as function of gate bias. (d) Characteristic peak rations
extracted from (c) as dots and from Fig. 1(c) as gray line. The shaded areas
in (c, d) mark the uncertainty ranges.

A conclusive display for the coherent coupling between X and X- is presented
in Fig. 5 (b), showing 2D FWM amplitude spectra for selected gate voltages as
labeled in the pots. The delay between pulses 2 and 3 was chosen to
$\tau_{23}=0.13\,$ps, which locates us in the coherent coupling range as
demonstrated in the SI Fig. S7. We have checked that for $\tau_{23}>1\,$ps the
X and X- coupling is present (see Fig. S8), although dominated by an
incoherent population transfer between the two complexes [26]. Conceptually,
the 2D FWM spectra of our V-shaped system consist of only four peaks. For the
diagonal pair, labeled X- and X, detection energy and $\hbar\omega_{12}$
energy are identical: FWM emission arises from the same absorption.
Decisively, we clearly detect off-diagonals, labeled XX- and X-X. This means
that the FWM of the charged exciton is also driven by the neutral exciton’s
first-order absorption, and vice-versa, respectively. Such coherent coupling
was previously reported in semiconductor nanostructures [48, 64, 59],
including TMDs [65, 66, 28]. Here, thanks to the tunability of $n_{\rm e}$ in
this gated MoSe2, we obviously find that the coherent nonlinear optical
response of the X - X- complex can be controlled simply by applying an
external bias.

The variation of the peak amplitudes in the 2D spectra with changing gate bias
has two potential reasons: (i) Change of the dipole strengths of X and X- due
to additional free charges, as already demonstrated in the reflectivity
measurement in Fig. 1(b), (ii) The coherence transfer between both exciton
species is affected by the free carriers. To disentangle these two effects we
need to quantify the strength of the coherent coupling depending on $n_{\rm
e}$. In panel (c) we plot the peaks’ integrated amplitudes versus the gate
bias, where we expectedly find that X (blue) clearly drops, while X- (red)
builds up when increasing the carrier density. At the same time the off-
diagonal peaks (green and yellow) – representing the coupling – show no clear
trend. As mentioned before, the variation of the reflectivity spectrum upon an
applied gate bias indicates a change of the dipole strengths of X and X- and
therefore of the pulse areas $\theta_{\rm X}$ and $\theta_{{\rm X}^{-}}$
applied to the two transitions. Note, that we do not consider any specific
origin for the change of the strength of optical response of X and X-. Now, in
order to extract the strength of the X - X- coupling found in the 2D spectra,
we can extract a quantity that does not depend on the pulse areas. To achieve
this, we note that the amplitudes of diagonal peaks scale as $\theta_{\rm
X}^{3}$ and $\theta_{{\rm X}-}^{3}$, whereas the off-diagonals follow the
scaling $\theta_{\rm X}^{2}\theta_{{\rm X}-}$ for XX- and $\theta_{\rm
X}\theta_{{\rm X}-}^{2}$ for X-X. Plotting now the expression,
XX${}^{-}\cdot$X-X/(X$\cdot$X-) in Fig. 5(d, top) the pulse areas and
therefore the dipole strengths cancel out. We see though, that this quantity
increases with the electron density, indicating the increase of the coupling
strength.

Another insightful quantity that can be extracted from the reflectivity
spectra and the 2D spectra is the sum of peak ratios. For the 2D spectra we
calculate XX-/X+X-X/X-, which translates into pulse areas as $\theta_{{\rm
X}^{-}}/\theta_{\rm X}+\theta_{\rm X}/\theta_{{\rm X}^{-}}$. To extract an
equivalent quantity from the reflectivity spectra we need to calculate
$\sqrt{R_{{\rm X}^{-}}}/\sqrt{R_{\rm X}}+\sqrt{R_{\rm X}}/\sqrt{R_{{\rm
X}^{-}}}$ from the peak intensities $R_{\rm X}$ and $R_{{\rm X}^{-}}$. The two
quantities are plotted (normalized to unity) in Fig. 5(d) as turquoise dots
for 2D and gray line for reflectivity. The behaviour of both curves is rather
unspecific for negative and small gate voltages, where X- is only weakly
addressed if present at all. Note, that in this range of bias values the X- is
very weak which is leading to the relatively large uncertainties (shaded
areas). Therefore, the peak ratios in this range are not particularly
significant and we do not expect to draw conclusions from this bias range.
However, the 2D spectra result shows a clearly growing trend for increasing
charge densities, while the reflectivity results shrink. This discrepancy is a
second hint that the variation of the dipole strengths with increasing gate
bias alone, cannot fully explain the variation of the coherent coupling
between X and X- observed in 2D FWM.

##  6 Conclusions & Outlook

FWM spectroscopy methods are powerful tools that have led in the recent years
to a significant improvement of the understanding of the rich exciton physics
in TMDs. This has established the method as a versatile technique for studies
of atomically thin materials and TMD-based heterostructures as detailed in the
Introduction. However, so far a particular focus on the important issue of the
free carrier influence on the studied phenomena has been missing. A factor
that fundamentally affects the optical response of TMDs, and governs even the
most basic characteristics like the number and respective intensities of
different exciton states.

Using FWM micro-spectroscopy, we have shown that the homogeneous linewidth and
population lifetime of excitonic complexes hosted by a MoSe2–monolayer-based
van der Waals heterostructure can be tuned via the free electron density,
which is injected by a gate bias applied to our device. With increasing gate
bias, the exciton’s inhomogeneous broadening decreases, reflecting the
screening of the disorder via free the electron gas, which also increases the
excitons’ radiative decay rates. Conversely, the homogeneous broadening
increases, which is attributed to the combined increase of the radiative decay
rate and the conversion rate of the neutral exciton towards the charged one.
By exciting the neutral (X) and charged exciton (X-) simultaneously and
probing the population dynamics we demonstrated that the coherence between X
and X- leads to a characteristic quantum beat in the FWM signal. By then
performing two-dimensional FWM spectroscopy for a variety of applied gate
biases, we have further demonstrated that the X - X- coherent coupling can be
controlled by the gate voltage, and hence by the free electron density.
Through considering specific peak ratios, we were able to demonstrate that the
change of coherent coupling can be disentangled from the variation of dipole
strengths arising from the injection of free carries. An increase of the
coupling strength with $n_{\rm e}$ could be linked with the screening of
disorder via the electron gas. This illustrates the utility and versatility of
ultrafast nonlinear spectroscopy in investigating optical responses of
excitonic systems, going beyond the capabilities of linear methods.

In future developments, by performing FWM with spatially separated driving
beams, while using devices hosting highly diffusive excitons, it will be
possible to achieve non-local coherent coupling in a two-dimensional
semiconductor. Combining this approach with two-color FWM spectroscopy would
also permit to selectively address the coherence transfer between neutral and
charged excitons. Our findings yield exciting prospects for forthcoming
investigations of coherent phenomena in the context of recent discoveries of
strongly-correlated exciton phases in solids [31], optically probed quantum
Hall states [32], moiré superlattices [67], and magnetic two-dimensional
materials [68]. In practice, coherent nonlinear spectroscopy could be used to
optically infer the dephasing processes of the Umklapp branches of TMD Wigner
crystals and fractional quantum Hall states in graphene.

This work was supported by the Polish National Science Centre under decisions
DEC-2020/39/B/ST3/03251. The Warsaw team (A.R., P.K. and M.P.) acknowledges
support from the ATOMOPTO project (TEAM program of the Foundation for Polish
Science, co-financed by the EU within the ERDFund), CNRS via IRP 2D Materials,
EU Graphene Flagship. M.P acknowledges support by the Foundation for Polish
Science (MAB/2018/9 Grant within the IRA Program financed by EU within SG OP
Program). The Polish participation in the European Magnetic Field Laboratory
(EMFL) is supported by the DIR/WK/2018/07 of MEiN of Poland. A. R.
acknowledges support of this work from the Diamentowy Grant under decision
DI2017008347 of MEiN of Poland. K. W. and T. T. acknowledge support from the
Elemental Strategy Initiative conducted by the MEXT, Japan, (grant no.
JPMXP0112101001), JSPS KAKENHI (grant no. JP20H00354), and the CREST
(JPMJCR15F3), JST. J. H. acknowledges support from EPSRC doctoral prize
fellowship. D. W. was supported by the Science Foundation Ireland (SFI) under
Grant 18/RP/6236. We thank Ronald Cox and Julien Renard for their constructive
comments on the manuscript, Rafał Bożek for the AFM measurement of the hBN
thicknesses, and Tilmann Kuhn for helpful discussions.

## References

## References

  * [1] Splendiani A, Sun L, Zhang Y, Li T, Kim J, Chim C Y, Galli G and Wang F 2010 Nano Lett. 10 1271–1275
  * [2] Mak K F, Lee C, Hone J, Shan J and Heinz T F 2010 Phys. Rev. Lett. 105(13) 136805
  * [3] Raja A, Chaves A, Yu J, Arefe G, Hill H M, Rigosi A F, Berkelbach T C, Nagler P, Schüller C, Korn T, Nuckolls C, Hone J, Brus L E, Heinz T F, Reichman D R and Chernikov A 2017 Nat. Commun. 8 15251
  * [4] Van Tuan D, Yang M and Dery H 2018 Phys. Rev. B 98 125308
  * [5] Chernikov A, Berkelbach T C, Hill H M, Rigosi A, Li Y, Aslan O B, Reichman D R, Hybertsen M S and Heinz T F 2014 Phys. Rev. Lett. 113(7) 076802
  * [6] Wang G, Chernikov A, Glazov M M, Heinz T F, Marie X, Amand T and Urbaszek B 2018 Rev. Mod. Phys. 90(2) 021001
  * [7] Stepanov P, Vashisht A, Klaas M, Lundt N, Tongay S, Blei M, Höfling S, Volz T, Minguzzi A, Renard J, Schneider C and Richard M 2021 Phys. Rev. Lett. 126(16) 167401
  * [8] Geim A K and Grigorieva I V 2013 Nature 499 419–425
  * [9] Cadiz F, Courtade E, Robert C, Wang G, Shen Y, Cai H, Taniguchi T, Watanabe K, Carrere H, Lagarde D, Manca M, Amand T, Renucci P, Tongay S, Marie X and Urbaszek B 2017 Phys. Rev. X 7(2) 021026
  * [10] Mak K F, He K, Shan J and Heinz T F 2012 Nat. Nanotechnol. 7 494–498
  * [11] Zeng H, Dai J, Yao W, Xiao D and Cui X 2012 Nat. Nanotechnol. 7 490–493
  * [12] Mak K F, McGill K L, Park J and McEuen P L 2014 Science 344 1489–1492
  * [13] Kioseoglou G, Hanbicki A T, Currie M, Friedman A L, Gunlycke D and Jonker B T 2012 Appl. Phys. Lett. 101 221907
  * [14] Cao T, Wang G, Han W, Ye H, Zhu C, Shi J, Niu Q, Tan P, Wang E, Liu B and Feng J 2012 Nat. Commun. 3 887
  * [15] Wierzbowski J, Klein J, Sigger F, Straubinger C, Kremser M, Taniguchi T, Watanabe K, Wurstbauer U, Holleitner A W, Kaniber M, Müller K and Finley J J 2017 Sci. Rep. 7 12383
  * [16] Ajayi O A, Ardelean J V, Shepard G D, Wang J, Antony A, Taniguchi T, Watanabe K, Heinz T F, Strauf S, Zhu X Y and Hone J C 2017 2D Mater. 4 031011
  * [17] Manca M, Glazov M M, Robert C, Cadiz F, Taniguchi T, Watanabe K, Courtade E, Amand T, Renucci P, Marie X, Wang G and Urbaszek B 2017 Nat. Commun. 8 14927
  * [18] Lindlau J, Selig M, Neumann A, Colombier L, Förste J, Funk V, Förg M, Kim J, Berghäuser G, Taniguchi T, Watanabe K, Wang F, Malic E and Högele A 2018 Nat. Commun. 9 2586
  * [19] Jakubczyk T, Nayak G, Scarpelli L, Liu W L, Dubey S, Bendiab N, Marty L, Taniguchi T, Watanabe K, Masia F, Nogues G, Coraux J, Langbein W, Renard J, Bouchiat V and Kasprzak J 2019 ACS Nano 13 3500–3511
  * [20] Boule C, Vaclavkova D, Bartos M, Nogajewski K, Zdražil L, Taniguchi T, Watanabe K, Potemski M and Kasprzak J 2020 Phys. Rev. Mater. 4(3) 034001
  * [21] Hahn T, Vaclavkova D, Bartos M, Nogajewski K, Potemski M, Watanabe K, Taniguchi T, Machnikowski P, Kuhn T, Kasprzak J and Wigger D 2022 Adv. Sci. 9 2103813
  * [22] Jakubczyk T, Delmonte V, Koperski M, Nogajewski K, Faugeras C, Langbein W, Potemski M and Kasprzak J 2016 Nano Lett. 16 5333–5339
  * [23] Jakubczyk T, Nogajewski K, Molas M R, Bartos M, Langbein W, Potemski M and Kasprzak J 2018 2D Mater. 5 031007
  * [24] Hao K, Moody G, Wu F, Dass C K, Xu L, Chen C H, Sun L, Li M Y, Li L J, MacDonald A H and Li X 2016 Nat. Phys. 12 677–682
  * [25] Hao K, Xu L, Wu F, Nagler P, Tran K, Ma X, Schüller C, Korn T, MacDonald A H, Moody G and Li X 2017 2D Mater. 4 025105
  * [26] Moody G, Kavir Dass C, Hao K, Chen C H, Li L J, Singh A, Tran K, Clark G, Xu X, Berghäuser G, Malic E, Knorr A and Li X 2015 Nat. Commun. 6 8315
  * [27] Conway M, Muir J, Earl S, Wurdack M, Mishra R, Tollerud J and Davis J 2022 2D Mater. 9 021001
  * [28] Hao K, Specht J F, Nagler P, Xu L, Tran K, Singh A, Dass C K, Schüller C, Korn T, Richter M, Knorr A, Li X and Moody G 2017 Nat. Commun. 8 15552
  * [29] Ross J S, Wu S, Yu H, Ghimire N J, Jones A M, Aivazian G, Yan J, Mandrus D G, Xiao D, Yao W and Xu X 2013 Nat. Commun. 4 1474
  * [30] Jones A M, Yu H, Ghimire N J, Wu S, Aivazian G, Ross J S, Zhao B, Yan J, Mandrus D G, Xiao D, Yao W and Xu X 2013 Nat. Nanotechnol. 8 634–638
  * [31] Smoleński T, Dolgirev P E, Kuhlenkamp C, Popert A, Shimazaki Y, Back P, Lu X, Kroner M, Watanabe K, Taniguchi T, Esterlis I, Demler E and Imamoğlu A 2021 Nature 595 53–57
  * [32] Popert A, Shimazaki Y, Kroner M, Watanabe K, Taniguchi T, Imamoğlu A and Smoleński T 2022 Nano Lett. 22 7363–7369
  * [33] Efimkin D K and MacDonald A H 2017 Phys. Rev. B 95(3) 035417
  * [34] Chang Y W and Reichman D R 2019 Phys. Rev. B 99(12) 125421
  * [35] Kossacki P, Cibert J, Ferrand D, Merle d’Aubigné Y, Arnoult A, Wasiela A, Tatarenko S and Gaj J A 1999 Phys. Rev. B 60(23) 16018–16026
  * [36] Kossacki P, Płochocka P, Piechal B, Maślana W, Golnik A, Cibert J, Tatarenko S and Gaj J A 2005 Phys. Rev. B 72(3) 035340
  * [37] Hawrylak P 1991 Phys. Rev. B 44(8) 3821–3828
  * [38] Brinkmann D, Kudrna J, Gilliot P, Hönerlage B, Arnoult A, Cibert J and Tatarenko S 1999 Phys. Rev. B 60(7) 4474–4477
  * [39] Koch M, Hellmann R, Bastian G, Feldmann J, Göbel E O and Dawson P 1995 Phys. Rev. B 51(19) 13887–13890
  * [40] Schultheis L, Kuhl J, Honold A and Tu C W 1986 Phys. Rev. Lett. 57(14) 1797–1800
  * [41] Wang H, Ferrio K, Steel D G, Hu Y Z, Binder R and Koch S W 1993 Phys. Rev. Lett. 71(8) 1261–1264
  * [42] Wegener M, Chemla D S, Schmitt-Rink S and Schäfer W 1990 Phys. Rev. A 42(9) 5675–5683
  * [43] Rodek A, Hahn T, Kasprzak J, Kazimierczuk T, Nogajewski K, Połczyńska K E, Watanabe K, Taniguchi T, Kuhn T, Machnikowski P, Potemski M, Wigger D and Kossacki P 2021 Nanophotonics 10 2717–2728
  * [44] Akbari H, Biswas S, Jha P K, Wong J, Vest B and Atwater H A 2022 Nano Lett. 22 7798–7803
  * [45] Savona V and Langbein W 2006 Phys. Rev. B 74 075311
  * [46] Rossi F and Kuhn T 2002 Rev. Mod. Phys. 74(3) 895–950
  * [47] Shah J 1999 Ultrafast spectroscopy of semiconductors and semiconductor nanostructures (Heidelberg: Springer)
  * [48] Langbein W and Patton B 2006 Opt. Lett. 31 1151–1153
  * [49] Niehues I, Schmidt R, Druppel M, Marauhn P, Christiansen D, Selig M, Berghäuser G, Wigger D, Schneider R, Braasch L, Koch R, Castellanos-Gomez A, Kuhn T, Knorr A, Malic E, Rohlfing M, Michaelis de Vasconcellos S and Bratschitsch R 2018 Nano Lett. 18 1751–1757
  * [50] Khatibi Z, Feierabend M, Selig M, Brem S, Linderälv C, Erhart P and Malic E 2018 2D Mater. 6 015015
  * [51] Raja A, Waldecker L, Zipfel J, Cho Y, Brem S, Ziegler J D, Kulig M, Taniguchi T, Watanabe K, Malic E, Heinz T F, Berkelbach T C and Chernikov A 2019 Nat. Nanotechnol. 14 832–837
  * [52] Christiansen D, Selig M, Berghäuser G, Schmidt R, Niehues I, Schneider R, Arora A, Michaelis de Vasconcellos S, Bratschitsch R, Malic E and Knorr A 2017 Phys. Rev. Lett. 119 187402
  * [53] Scarpelli L, Masia F, Alexeev E M, Withers F, Tartakovskii A I, Novoselov K S and Langbein W 2017 Phys. Rev. B 96(4) 045407
  * [54] Zipfel J, Wagner K, Semina M A, Ziegler J D, Taniguchi T, Watanabe K, Glazov M M and Chernikov A 2022 Phys. Rev. B 105(7) 075311
  * [55] Ciulin V, Kossacki P, Haacke S, Ganiére J D, Deveaud B, Esser A, Kutrowski M and Wojtowicz T 2000 Phys. Rev. B 62(24) R16310
  * [56] Kossacki P 2003 Journal of Physics: Condensed Matter 15 R471
  * [57] Mermillod Q, Wigger D, Delmonte V, Reiter D E, Schneider C, Kamp M, Höfling S, Langbein W, Kuhn T, Nogues G and Kasprzak J 2016 Optica 3 377–384
  * [58] Li X, Zhang T, Borca C N and Cundiff S T 2006 Phys. Rev. Lett. 96(5) 057406
  * [59] Singh A, Moody G, Wu S, Wu Y, Ghimire N J, Yan J, Mandrus D G, Xu X and Li X 2014 Phys. Rev. Lett. 112(21) 216804
  * [60] Hamm P and Zanni M 2011 Basics of 2D IR spectroscopy (Cambridge University Press) p 61–87
  * [61] Moody G and Cundiff S T 2017 Adv. Phys.-X 2 641–674
  * [62] Smallwood C L and Cundiff S T 2018 Laser Photonics Rev. 12 1800171
  * [63] Delmonte V, Specht J F, Jakubczyk T, Höfling S, Kamp M, Schneider C, Langbein W, Nogues G, Richter M and Kasprzak J 2017 Phys. Rev. B 96(4) 041124
  * [64] Kasprzak J, Patton B, Savona V and Langbein W 2011 Nature Photon. 5 57–63
  * [65] Singh A, Moody G, Wu S, Wu Y, Ghimire N J, Yan J, Mandrus D G, Xu X and Li X 2014 Phys. Rev. Lett. 112(21) 216804
  * [66] Hao K, Xu L, Nagler P, Singh A, Tran K, Dass C K, Schüller C, Korn T, Li X and Moody G 2016 Nano Lett. 16 5109–5113
  * [67] Tran K, Moody G, Wu F, Lu X, Choi J, Kim K, Rai A, Sanchez D A, Quan J, Singh A, Embley J, Zepeda A, Campbell M, Autry T, Taniguchi T, Watanabe K, Lu N, Banerjee S K, Silverman K L, Kim S, Tutuc E, Yang L, MacDonald A H and Li X 2019 Nature 567 71–75
  * [68] Wilson N P, Lee K, Cenker J, Xie K, Dismukes A H, Telford E J, Fonseca J, Sivakumar S, Dean C, Cao T, Roy X, Xu X and Zhu X 2021 Nat. Mater. 20 1657–1662