# Excitonic insulator states in molecular functionalized atomically-thin
semiconductors

Dominik Christiansen1 Malte Selig1 Mariana Rossi2,3 Andreas Knorr1 1Institut
für Theoretische Physik, Nichtlineare Optik und Quantenelektronik, Technische
Universität Berlin, 10623 Berlin, Germany  2Max Planck Institute for the
Structure and Dynamics of Matter, 22761 Hamburg, Germany  3Fritz Haber
Institute of the Max Planck Society, 14195 Berlin, Germany

###### Abstract

The excitonic insulator is an elusive electronic phase exhibiting a correlated
excitonic ground state. Materials with such a phase are expected to have
intriguing properties such as excitonic high-temperature superconductivity.
However, compelling evidence on the experimental realization is still missing.
Here, we theoretically propose hybrids of two-dimensional semiconductors
functionalized by organic molecules as prototypes of excitonic insulators,
with the exemplary candidate WS2-F6TCNNQ. This material system exhibits an
excitonic insulating phase at room temperature with a ground state formed by a
condensate of interlayer excitons. To address an experimentally relevant
situation, we calculate the corresponding phase diagram for the important
parameters: temperature, gap energy, and dielectric environment. Further, to
guide future experimental detection, we show how to optically characterize the
different electronic phases via far-infrared to terahertz (THz) spectroscopy.

Introduction: The excitonic insulator (EI) is a charge neutral, strongly
interacting insulating phase that arises from spontaneous formation of
excitons. The EI was first predicted theoretically in 1965 by L.V. Keldysh
Keldysh and Kopaev (1965) and systematically investigated by W. Kohn Jérome et
al. (1967). The insulating phase is anticipated to appear in semiconductors at
thermodynamic equilibrium, as long as the predicted exciton binding energy
exceeds the band gap. It presents an interesting platform for realizing many-
body ground states of condensed bosons in solids. Although the concept has
been known for almost 60 years, to date compelling experimental evidence of
the excitonic insulator is still missing. Since the EI phase is predicted to
host many novel properties, such as superfluidity Gupta et al. (2020);
Perfetto and Stefanucci (2021), excitonic high-temperature superconductivity
Parmenter and Henson (1970); Wang et al. (2019a), and exciton condensation,
breakthroughs in finding this new class of insulators has attracted great
attention over the last decades Varsano et al. (2017); Wang et al. (2019a);
Perfetto et al. (2019, 2020); Perfetto and Stefanucci (2020); Jia et al.
(2020).

There are a few materials, which are suspected to possess EI ground states in
a solid state, namely 1T-TiSe2, Ta2NiSe5 or TmSe0.45Te0.55 Cercellier et al.
(2007); Monney et al. (2010); Wakisaka et al. (2009); Lu et al. (2017);
Werdehausen et al. (2018); Bucher et al. (1991); Bronold and Fehske (2006);
Fukutani et al. (2021). However, it has been difficult to establish whether
the EI state has been realized, because the expected electronic phase
transition is accompanied by a structural phase transition, which makes it
difficult to distinguish between EI and normal insulator or Peierls transition
Rossnagel et al. (2002); Zenker et al. (2014); Jiang et al. (2018); Volkov et
al. (2021). Just recently, the discussion started about EIs in TMDC
heterostructures Du et al. (2017); Ataei et al. (2021); Ma et al. (2021);
Zhang et al. (2021).

Figure 1: (a) The HIOS consisting of a monolayer WS2 and a layer of F6-TCNNQ
molecules with level alignment in the free-particle picture. (b)-(d) The
interlayer gap dispersion of the dashed boxed for the three possible
electronic phases and the character of their optical excitation.

In this Letter, we present a blueprint for the first realization of an EI
based on the interlayer excitons of hybrid inorganic-organic systems (HIOS).
HIOS is a growing field with increasing technological importance because it
combines the best of two worlds: The strong light-matter
interaction/tunability of the transition energies of organic molecules with
the high carrier mobility of inorganic semiconductors Sun et al. (2019); Song
et al. (2017); Koch (2007); Fahlman et al. (2019). For the construction of an
EI, the low dielectric constant of the molecular lattices Torabi et al. (2015)
and the strong localization of their electrons, due to an infinite effective
mass, are excellent conditions for large exciton binding energies of HIOS
interface excitons. In particular, the functionalization of atomically-thin
semiconductors with organic molecules allows to choose a material combination
with a band gap in an appropriate range Groom et al. (2016); Stuke et al.
(2020) and the spatial indirect character of the exciton allows for static
dipoles and thus unambiguous fine tuning via static Stark shifts Lorchat et
al. (2021). A the same time, the energy level tunability of the molecular
layer has advantages over TMDC bilayer or heterobilayers, where the EI phase
is still under discussion Du et al. (2017); Ataei et al. (2021); Ma et al.
(2021); Zhang et al. (2021).

Figure 1(a) depicts the investigated HIOS with interlayer dispersion close to
the K-point. The relevant electronic interlayer transition for the EI built up
occurs between TMDC valence and molecular conduction band highlighted by the
dashed line box. Due to the direct-gap character of the heterostructure
dispersion, we also circumvent the formation of Peierls charge density waves.
Depending on the interlayer band gap (adjustable also by tuning of an applied
voltage) three different electronic phases can be expected: semi-metal,
semiconductor, and EI, cp. Fig. 1(b)-(d). We will show that all three phases
can be distinguished by far-infrared/THz absorption: The semi-metal phase
exhibits a vanishing or negative band gap with a free electron gas as ground
state, cf. Fig. 1(b). The optical response will then be determined by the
Drude model known from metals Drude (1900). The semiconducting phase, cf. Fig.
1(c), is determined by an exciton binding energy smaller than the free-
particle band gap resulting in a filled valence band as ground state. The
optical response is given as Lorentz response, however due to the interlayer
character of the transition with considerably small oscillator strength Gillen
(2021); Deilmann and Thygesen (2018). In contrast, the excitonic insulating
phase, cf. Fig. 1(d), occurs if the expected exciton binding energy is larger
than the band gap. We will exploit that the far-infrared response of such EIs
are characterized by transitions in the exciton ladder Kira and Koch (2006);
Bhattacharyya et al. (2014) of condensed ground state excitons to higher bound
states.

Figure 2: (a) Numerical solution of the ordering parameter
$\Delta_{\mathbf{k}}$ and gap dispersion $\Sigma_{\mathbf{k}}$ as function of
wave number at room temperature. Both quantities constitute the Bogoliubov
dispersion $E_{\mathbf{k}}$. (b) The EI exhibits a non-vanishing coherence and
an intrinsic inversion. The coherence and occupations peak in the pocket of
the Bogoliubov dispersion.

In the following, even without applied electric field for energy band fine
tuning, we predict an EI from a hybrid structure of a WS2 monolayer and a
self-assembled layer of F6-TCNNQ molecules, cf. Fig. 1(a). The organic
molecules form a periodic two-dimensional lattice, where we investigate the
realistic ratio of one molecule per 16 WS2 unit cells due to the weak
molecule-WS2 interactionPark et al. (2021). From first-principles calculations
we find a type-II band alignment heterostructure with an interlayer gap
between WS2 valence band and F6-TCNNQ lowest unoccupied molecular orbital of
$E_{G}=0.12\,\mathrm{eV}$. The Kohn-Sham energy levels of F6-TCNNQ and WS2
were obtained with the range-separated HSE06 hybrid functional Krukau et al.
(2006), as implemented in the FHI-aims program Blum et al. (2009); Ren et al.
(2012); Levchenko et al. (2015) and using standard „intermediate“ settings
Park et al. (2021); Wang et al. (2019b). The naively calculated interlayer
exciton binding energy, using well documented methods Malic et al. (2018),
with completely filled valence band as ground state is
$E_{B}=0.14\,\mathrm{eV}$. This is on the same order of magnitude as the band
gap of 0.12 eV, indicating the possibility of forming a strongly correlated
insulating phase to minimize energy. Based on this observation, we develop in
the following the description of the EI for the introduced hybrid and
calculate its optical response. In particular, we show that the excitonic
condensate can be stable up to room temperature under the proper choice of
substrate.

Ground state of EI: To calculate the EI ground state, we diagonalize the
field-independent part of the mean-field Hamiltonian (cf. appendix) by a
Bogoliubov transformation Bogoljubov et al. (1958); Valatin (1958) with
$\alpha^{\dagger}_{\mathbf{k}}=u^{*}_{\mathbf{k}}v^{\dagger}_{\mathbf{k}}-w^{*}_{\mathbf{k}}c^{\dagger}_{\mathbf{k}}$
and
$\beta^{\dagger}_{\mathbf{k}}=w_{\mathbf{k}}v^{\dagger}_{\mathbf{k}}+u_{\mathbf{k}}c^{\dagger}_{\mathbf{k}}$
Jérome et al. (1967); Keldysh and Kozlov (1968); Comte and Nozieres (1982).
The fermionic operators $\alpha^{\dagger}_{\mathbf{k}}$ and
$\beta^{\dagger}_{\mathbf{k}}$ create an electron in a linear combination of
valence $v^{(\dagger)}_{\mathbf{k}}$ and conduction
$c^{(\dagger)}_{\mathbf{k}}$ bands in analogy to the Bogoliubov particle
operators from BCS superconductivity theory. The diagonalized Hamiltonian
reads
$H=\sum_{n=\\{\alpha,\beta\\},\mathbf{k}}E_{n,\mathbf{k}}n^{\dagger}_{\mathbf{k}}n^{\mathstrut}_{\mathbf{k}}$
with the hybridized bands
$E_{\alpha/\beta,\mathbf{k}}=(\tilde{\varepsilon}_{c,\mathbf{k}}+\tilde{\varepsilon}_{v,\mathbf{k}})/2\mp\sqrt{\Sigma^{2}_{\mathbf{k}}+\Delta^{2}_{\mathbf{k}}}$.
The excitation spectrum of the new quasi-particles corresponds to the
Bogoliubov dispersion
$E_{\mathbf{k}}=\sqrt{\Sigma^{2}_{\mathbf{k}}+\Delta^{2}_{\mathbf{k}}}$, where
the gap dispersion is defined as
$\Sigma_{\mathbf{k}}=(\tilde{\varepsilon}_{c,\mathbf{k}}-\tilde{\varepsilon}_{v,\mathbf{k}})/2$.
The quantity $\Delta_{\mathbf{k}}$ is determined via the transcendental gap
equation Kozlov and Maksimov (1965); Sabio et al. (2010); Stroucken et al.
(2011)

$\displaystyle\Delta_{\mathbf{k}}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{k^{\prime}}}V_{\mathbf{k-k^{\prime}}}\frac{\Delta_{\mathbf{k}^{\prime}}}{\sqrt{\Sigma^{2}_{\mathbf{k}^{\prime}}+\Delta^{2}_{\mathbf{k}^{\prime}}}}\left(f_{\alpha,\mathbf{k}^{\prime}}-f_{\beta,\mathbf{k}^{\prime}}\right)$
(1) $\displaystyle 2\Sigma_{\mathbf{k}}$
$\displaystyle=\varepsilon_{c,\mathbf{k}}-\varepsilon_{v,\mathbf{k}}+\sum_{\mathbf{k}^{\prime}}\mathcal{V}^{\text{mol}}_{\mathbf{k-k^{\prime}}}\left(u^{2}_{\mathbf{k}^{\prime}}f_{\beta,\mathbf{k}^{\prime}}+w^{2}_{\mathbf{k}^{\prime}}f_{\alpha,\mathbf{k}^{\prime}}\right)$
$\displaystyle-\sum_{\mathbf{k}^{\prime}}\mathcal{V}^{\text{WS}_{2}}_{\mathbf{k-k^{\prime}}}\left(u^{2}_{\mathbf{k}^{\prime}}f_{\alpha,\mathbf{k}^{\prime}}+w^{2}_{\mathbf{k}^{\prime}}f_{\beta,\mathbf{k}^{\prime}}\right)\;.$
(2)

The coupled system of equations determines $\Delta_{\mathbf{k}}$,
$u_{\mathbf{k}}$, and $w_{\mathbf{k}}$ via temperature and band gap. The
quantity $\Delta_{\mathbf{k}}$ can be identified as ordering parameter
determining the phase of the heterostructure. A finite value accounts for a
finite probability to create electron-hole pairs, which designates the
excitonic instability. In the limit $\Delta_{\mathbf{k}}\rightarrow 0$ the
dipolar excitonic insulator Zimmermann and Schindler (2007); Schindler and
Zimmermann (2008); Gu et al. (2021) converts to the conventional phases of
semiconductor or semi-metal depending on the band gap. In Eq. (1) and (2) we
defined the matrix element
$\mathcal{V}^{l}_{\mathbf{k-k^{\prime}}}=V_{\mathbf{0}}^{l}-V_{\mathbf{k-k^{\prime}}}^{l}-V_{\mathbf{0}}+V_{\mathbf{k-k^{\prime}}}$
with intralayer Coulomb potential $V^{l}_{\mathbf{k}}$
($l=\\{\text{molecule,WS}_{2}\\}$), and interlayer potential
$V_{\mathbf{k}}$Ovesen et al. (2019) (cf. appendix). Equation (1) and (2)
depend on the occupations of the hybridized bands $f_{n,\mathbf{k}}=\langle
n^{\dagger}_{\mathbf{k}}n^{\mathstrut}_{\mathbf{k}}\rangle$ with
$n=\\{\alpha,\beta\\}$. The temperature $T=0\,\mathrm{K}$ limit of the gap
equation, i.e. $f_{\alpha,\mathbf{k}}-f_{\beta,\mathbf{k}}=1$, can also be
obtained from a minimization of the energy Sabio et al. (2010); Stroucken et
al. (2011). Clearly, the magnitude of the ordering parameter
$\Delta_{\mathbf{k}}$ depends on band gap and temperature, which enters via
the Fermi functions $f_{\alpha/\beta,\mathbf{k}}$ of the hybridized bands
$E_{\alpha/\beta,\mathbf{k}}$. The chemical potential is chosen such that the
charge density is a conserved quantity as function of temperature and band
gap. In our case, we consider a charge neutral structure, i.e. the density of
holes in WS2 equals the density of electron in the molecule. Because the
ordering parameter enters also in the gap dispersion $\Sigma_{\mathbf{k}}$ via
the Fermi functions $f_{n,\mathbf{k}}(E_{n,\mathbf{k}})$, both quantities have
to be solved simultaneously. If not stated otherwise we use an hBN substrate
($\epsilon=4.5$) entering the Coulomb potential and a surrounding of air
($\epsilon=1$) for the numerical evaluation.

Figure 3: (a) Phase diagram of a WS2-F6TCNNQ stack on hBN and air surrounding.
The EI phase is seperated from the semi-metal by the semiconducting phase. The
numbers denote the position for the absorption spectra in Fig. 4. Number 3
corresponds to zero applied field. (b) Phase diagram of the WS2-F6TCNNQ stack
on hBN and altering supstrate, characterized by the dielectric constant.
Environments with low dielectric constant are favourable to stabilize the
excitonic insulating phase at higher temperature.

Figure 2(a) displays the numerical solution of $\Delta_{\mathbf{k}}$ and
$\Sigma_{\mathbf{k}}$ at room temperature, and the resulting Bogoliubov
dispersion $E_{\mathbf{k}}$ with the original band gap of
$E_{G}=0.12\,\mathrm{eV}$. The ordering parameter $\Delta_{\mathbf{k}}$ is
symmetric to $k=0$ and displays a monotonous decrease with the wave number,
comparable to the gap function of s-wave superconductors Bardeen et al.
(1957); Tinkham (2004); Reyes et al. (2016). Together with
$\Sigma_{\mathbf{k}}$ it yields a sombrero-like Bogoliubov dispersion. We
stress that Fig. 2(a) displays the results for 300 K already suggesting that
the excitonic condensate is stable up to room temperature. Finally, Fig. 2(a)
displays also the gap dispersion $\Sigma_{\mathbf{k}}$. For the discussion of
the gap dispersion, it is convenient to introduce the spontaneously forming
inversion $f^{(0)}_{\mathbf{k}}\equiv
f^{(0)}_{v,\mathbf{k}}-f^{(0)}_{c,\mathbf{k}}=\Sigma_{\mathbf{k}}/\sqrt{\Sigma^{2}_{\mathbf{k}}+\Delta^{2}_{\mathbf{k}}}(f_{\alpha,\mathbf{k}}-f_{\beta,\mathbf{k}})$
and the macroscopic coherence $p^{(0)}_{\mathbf{k}}=\langle
v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}\rangle^{(0)}=\frac{1}{2}\Delta_{\mathbf{k}}/\sqrt{\Sigma^{2}_{\mathbf{k}}+\Delta^{2}_{\mathbf{k}}}(f_{\alpha,\mathbf{k}}-f_{\beta,\mathbf{k}})$
forming without external source and plotted in Fig. 2(b). Here, in contrast to
a conventional semiconductor ($\Delta_{\mathbf{k}}=0$) for
$\Delta_{\mathbf{k}}\neq 0$ the ground state coherence $p^{(0)}_{\mathbf{k}}$
has a non-vanishing value and the occupation inversion deviates from unity
clarifying why $\Delta_{\mathbf{k}}$ is referred to as ordering parameter. Due
to the Hartree-Fock renormalizations the gap dispersion $\Sigma_{\mathbf{k}}$
turns negative, cf. Fig. 2(a). This leads to an intrinsic inversion
$f^{(0)}_{\mathbf{k}}<0$ close to the band extremum. From Fig. 2(b) we see
that the ground state coherence and occupations peak within the pockets of the
Bogoliubov dispersion, which we can identify with the Fermi wave number
$k_{F}$.

Phase diagram: Depending on the two external parameters temperature and band
gap, tunable by static electric fields Deilmann and Thygesen (2018); Lorchat
et al. (2021), affecting the ordering parameter $\Delta_{\mathbf{k}}\neq 0$ we
can expect three different electronic phases: EI, semiconductor, and semi-
metal, cf. Fig. 1(b)-(d). While the EI phase is present for a finite ordering
parameter $\Delta_{\mathbf{k}}$, the other two phases reveal a vanishing
$\Delta_{\mathbf{k}}$. However, semi-metal and semiconductor can be
distinguished via the inversion $f^{(0)}_{\mathbf{k}}$: While for a
conventional semiconductor the inversion is one in the ground state, the value
is smaller than one for a semi-metal reflecting the presence of a free
electron gas. Figure 3(a) shows the calculated phase diagram of a WS2-F6TCNNQ
stack on hBN substrate as function of temperature and band gap. As guidance we
include the coexistence lines between the different phases. We find that the
excitonic insulating phase is stable for temperatures up to 350 K underlining
the ability of our proposed structure as high-temperature EI. The EI phase
appears in the band gap range of 0.04 eV to 0.38 eV. We see that by increasing
the band gap the heterostructure enters its semiconducting phase. When
decreasing the band gap we approach the semi-metal limit. Interestingly, we
find no coexistence line between EI and semi-metal, but the heterostructure
traverses the semiconducting phase again. This results from a fast decrease of
the ordering parameter, vanishing prior to a negative band gap. The
corresponding separating area between EI and semi-metal could be understood as
an excited semiconducting phase ($\Delta_{\mathbf{k}}=0$,
$E_{\mathbf{k}}=\tilde{\varepsilon}_{c,\mathbf{k}}-\tilde{\varepsilon}_{v,\mathbf{k}}>0$
as in the semiconducting phase but with $f_{c,\mathbf{k}}\neq 0$,
$f^{(0)}_{\mathbf{k}}<1$). To tune the band gap over the full energy range,
which is discussed in Fig. 3(a), a maximum applied voltage of 0.5 V between
molecule and inorganic semiconductor is necessary. Finally, we see a
sublimation line at high enough temperature – a direct transition from
semiconductor to semi-metal. The billiard balls in Fig. 3(a) are discussed
later together with Fig. 4.

We stress that the phase diagram strongly depends on the dielectric
environment influencing the interlayer Coulomb potential Ovesen et al. (2019):
In Fig. 3(b) we show the phases for the fixed original band gap of
$E_{G}=0.12\,\mathrm{eV}$ with altering temperature and dielectric
surrounding. The WS2-F6TCNNQ stack is still placed on top of an hBN substrate,
as in all calculations before, but the dielectric constant of a top layer is
changed from vacuum to hBN encapsulation: The stability of the condensate for
higher temperature rapidly decreases with increasing dielectric constant.
Obviously, a dielectric environment with low mean dielectric constant is
favourable to stabilize the EI phase at high temperatures.

Figure 4: Absorption spectrum at 100 K along the dotted line in Fig. 3(a). The
Drude response in the semi-metal phase is modulated by a $p$-exciton Rydberg
series for the EI phase in the terahertz regime. In the semiconducting phase
the terahertz response vanishes.

Optical response: Based on the ground state calculations, we can now calculate
the frequency-dependent absorption coefficient via a self-consistent coupling
of material and wave equation Knorr et al. (1996); Malic and Knorr (2013);
Katsch and Knorr (2020). The linear absorption with respect to the dynamical
field $\mathbf{E}(t)$ is determined by the susceptibility described by the
microscopic polarization $p_{\mathbf{k}}$ and occupations $f_{\mathbf{k}}$.
Both quantities are expanded up to first order in the exciting electric field
Stroucken et al. (2011); Stroucken and Koch (2015). The initial conditions
arise from the ground state (initially calculated $p^{(0)}_{\mathbf{k}}$ and
$f^{(0)}_{\mathbf{k}}$). The dynamical correction to first order in
$\mathbf{E}(t)$ is denoted by $p^{(1)}_{\mathbf{k}}$ and
$f^{(1)}_{\mathbf{k}}$. The Heisenberg equation of motion for the microscopic
polarization $p^{(1)}_{\mathbf{k}}$ can be diagonalized exploiting the
Bogoliubov-Wannier equation Stroucken and Koch (2015). After diagonalization,
based on a 1$s$ ground state, it describes far-infrared transitions to higher
lying excitonic states $\mu$ with their excitonic energy $E_{\mu}$ and wave
function $\varphi_{\mu,\mathbf{k}}$. Consequently, the microscopic
polarization can be expressed as excitonic polarization
$p_{\mu}=\sum_{\mathbf{k}}\varphi^{*}_{\mu,\mathbf{k}}p^{(1)}_{\mathbf{k}}$
Kira and Koch (2006); Katsch et al. (2020); Dong et al. (2021a) (details in
appendix). A direct solution of the Heisenberg equations of motion for
excitonic polarization and excited occupation in frequency domain yields the
susceptibility, detectable in far-infrared/THz experiments:

$\displaystyle\chi(\omega)=-\frac{1}{\epsilon_{0}}\sum_{\mu}\frac{\mathbf{d}_{\mu}\otimes\mathbf{d}_{\mu}+\mathbf{j}_{\mu}\otimes\mathbf{j}_{\mu}}{\hbar\omega-
E_{\mu}+i\gamma}+\frac{e_{0}^{2}}{\epsilon_{0}\hbar}\sum_{\mathbf{k}}\frac{\mathbf{v}_{\mathbf{k}}\otimes\nabla_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}}{\omega^{2}+i\gamma\omega/\hbar}$
(3)

with the particle velocity $\mathbf{v}_{\mathbf{k}}=\hbar\mathbf{k}/m$,
elementary charge $e_{0}$, and valence electron mass $m$. The excitonic
interband matrix element reads
$\mathbf{d}_{\mu}=\mathbf{d}\sum_{\mathbf{k}}f^{(0)}_{\mathbf{k}}\varphi^{*}_{\mu,\mathbf{k}}$
with electronic dipole moment $\mathbf{d}$. The second summand is driven by
the excitonic intraband matrix element
$\mathbf{j}_{\mu}=e_{0}\cdot\sum_{\mathbf{k}}\varphi^{*}_{\mu,\mathbf{k}}\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}$.
Both contributions exhibit resonances at the exciton energy $E_{\mu}$. For an
EI with $s$-symmetric ground state the first excited exciton is of
$p$-symmetry, that we can expect the 1$s$-2$p$ transition as the energetically
lowest resonance Kira and Koch (2006); Bhattacharyya et al. (2014). For a
$p$-excited state ($\mu=p$ in Eq. (3)), the interband source vanishes because
of its uneven parity. In contrast, the intraband source $\mathbf{j}_{\mu}$ is
finite. When entering the semiconducting phase, the Bogoliubov-Wannier
equation yields as lowest excited state 1$s$ excitons. In this phase, the
optical source is of interband nature from a fully occupied valence band as
ground state, cf. Fig. 1(c). In contrast the intraband source vanishes due to
symmetry reasons. Additionally, now holds $f^{(0)}_{v,\mathbf{k}}=1$ (full
valence band) that the third term in Eq. (3), which resembles to the
conductivity tensor of a plasma, vanishes. Finally, for zero or negative band
gap the heterostructure is in its semi-metallic phase and the optical matrix
elements $\mathbf{d}_{\mu}$ and $\mathbf{j}_{\mu}$ are zero due to the
corresponding ground state coherence and inversion. The optical response is
now solely described by the third term in Eq. (3). We can perform a partial
integration to bring the third term in the susceptibility into the form
$\chi_{\text{Drude}}(\omega)=\omega_{pl}^{2}/(\omega^{2}+i\gamma\omega/\hbar)$,
which corresponds to the plasma response in a Drude model for a free electron
gas, cf. Fig. 1(b). It stems from intraband transitions of the microscopic
occupation with the plasma frequency
$\omega^{2}_{pl}=e^{2}_{0}\sum_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}/\epsilon_{0}m$.
Also for the EI phase the Drude response is present, since the ground state
valence band occupation $f^{(0)}_{v,\mathbf{k}}$ exhibits a Fermi edge. In
order to account for the broadening of the response, we include a
phenomenological dephasing $\gamma$ Selig et al. (2016); Christiansen et al.
(2017); Helmrich et al. (2021). A detailed investigation of the exciton
lifetime in the insulating phase would require an evaluation of exciton-phonon
interaction into the HIOS Bloch equations.

Figure 4 displays the calculated absorption along the enumerated line in Fig.
3(a). In the semi-metal phase (dot 1) with overlapping valence and conduction
band we observe the well-known Drude response for a free electron gas. Due to
the infinite mass of the conduction band electrons in the molecules the
response stems solely from the half filled WS2 valence band. Opening the gap
(dot 2) an additional feature arises, which stems from the intraband matrix
element in Eq. (3). The rising macroscopic coherence $p^{(0)}_{\mathbf{k}}$ is
convoluted with the excitonic wave function, which interpolates between $s$
and $p$ state. When entering the EI phase (dot 3) the Drude response is
modulated by a $p$-excitonic Rydberg series in the far-infrared/THz regime
stemming from the transitions from 1$s$ ground state to $p$-excited states.
Since the optical source of the observed 1$s$-2$p$ transition is of intraband
nature, it has sufficient oscillator strength (dot 3 and 4) to be observed in
optical experiments. The oscillator strength decreases as a function of the
exciton number due to a decrease of the integral in Eq. (3). For a further gap
opening we observe a blue shift of the exciton and a decreasing oscillator
strength since the ordering parameter and the connected ground state functions
$f^{(0)}_{\mathbf{k}}$ and $p_{\mathbf{k}}^{(0)}$ decrease. Finally, for a
band gap of 0.4 eV (dot 5) the wave function of the lowest excited state
changes to $s$-symmetry, that the intraband source vanishes, marking the
transition to the semiconducting phase. The interlayer exciton is now driven
by the interband dipole matrix element. Because of the large detuning of
intra- and interlayer exciton, the hybridization is small yielding an
extremely small electronic dipole element Lorchat et al. (2021); Deilmann and
Thygesen (2018); Gillen (2021). Also the Drude response vanishes due to
$f^{(0)}_{\mathbf{k}}=1$. Therefore, the interlayer exciton is not observable
in absorption, the far-infrared/THz response vanishes, and the heterostructure
becomes transparent for these wavelength. For the situation at room
temperature a similar picture emerges characterized by a generally weaker
oscillator strength since the ground states acting as sources are less
populated.

Conclusion: In summary, we propose HIOS as a candidate for the realization of
an EI exploiting interlayer excitons. For this, the molecule decoration of an
atomically-thin semiconductor is chosen in a way that the interlayer band gap
lies in the range of the binding energy of the interlayer excitons. Additional
static fields can be used for fine tuning or to induce different electronic
phases. The occurring optical far-infrared/THz response can be used to
characterize the electronic phases: While the conventional interlayer
semiconductor exhibits a $s$-like Rydberg series with minimal oscillator
strength, the EI’s Rydberg series has $p$-character with strong oscillator
strength including a Drude-like response from the ground state occupation. We
expect that our results trigger new studies on high-temperature exciton
condensation and an experimental realization of the predicted room temperature
EI could lead to new optoelectronic applications.

We thank Florian Katsch, Manuel Katzer, Lara Greten (TU Berlin) and Kiril
Bolotin (FU Berlin) for fruitful discussions. We acknowledge financial support
from the Deutsche Forschungsgemeinschaft (DFG) through SFB 951 (D.C., M.S.,
A.K.) Projektnummer 182087777. D.C. thanks the graduate school Advanced
Materials (SFB 951) for support.

##  Appendix A Hamiltonian

Due to the periodic arrangement of the molecules, the molecular operators of
the lowest unoccupied orbital $c$ of molecule $\nu$ can be transformed to a
Bloch basis
$c_{\mathbf{k}}=\sum_{\nu}\exp(-i\mathbf{k}\cdot\mathbf{R}_{\nu})c_{\nu}/\sqrt{N_{\text{m}}}$
with $N_{\text{m}}$ the number of molecular unit cells and wave vector
$\mathbf{k}$ Slater (1934); Kaplan and Ruvinskiǐ (1976); Specht et al. (2018).
At the same time, valence band electrons in the TMDC layer are described by
operators $v_{\mathbf{k}}$. In this notation, we can construct a many-particle
Hamiltonian, which is parameterized from electronic structure ab initio
calculations:

$\displaystyle H$
$\displaystyle=\sum_{\lambda,\mathbf{k}}\varepsilon_{\lambda,\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}}+ie_{0}\sum_{\lambda,\mathbf{k}}\mathbf{E}(t)\cdot(\nabla_{\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}})\lambda^{\mathstrut}_{\mathbf{k}}+\hbar\sum_{\mathbf{k}}\Omega_{\mathbf{k}}(t)\left(v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}+c^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}\right)+\frac{1}{2}\sum_{\begin{subarray}{c}\lambda,\lambda^{\prime},\nu,\nu^{\prime}\\\
\mathbf{k,k^{\prime},q}\end{subarray}}V^{\lambda\nu\nu^{\prime}\lambda^{\prime}}_{\mathbf{k,k^{\prime},q}}\leavevmode\nobreak\
\lambda^{\dagger}_{\mathbf{k+q}}\nu^{\dagger}_{\mathbf{k^{\prime}-q}}\nu^{\prime\mathstrut}_{\mathbf{k^{\prime}}}\lambda^{\prime\mathstrut}_{\mathbf{k}}.$
(4)

The first term describes the free kinetic energy. The single-particle energies
are $\varepsilon_{c,\mathbf{k}}=E_{G}$ and
$\varepsilon_{v,\mathbf{k}}=\hbar^{2}\mathbf{k}^{2}/2m$ with effective hole
mass $m$. The second and third term describe intra- and interband transitions,
respectively. The last term includes Coulomb interaction. To account for
inter- and intraband renormalization effects and exciton formation we can
restrict the band indices to the combinations of: all band indices correspond
to valence or conduction band and to even numbers of valence and conduction
band indices. Uneven numbers of band indices would describe Meitner-Auger
processes. However, since valence and conduction band are in different layers,
the Meitner-Auger process would require a large wave function overlap to
significantly contribute and can therefore be neglected Dong et al. (2021b).
Together with a random phase approximation for electronic occupations
$\langle\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}^{\prime}}\rangle\rightarrow\langle\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}^{\prime}}\rangle\delta_{\mathbf{k,k^{\prime}}}$
we obtain the Hartree-Fock Hamiltonian

$\displaystyle H$
$\displaystyle=\sum_{\mathbf{k}}\left(\varepsilon_{c,\mathbf{k}}+\sum_{\mathbf{k^{\prime}}}V^{\text{mol}}_{\mathbf{0}}f_{c,\mathbf{k^{\prime}}}+\sum_{\mathbf{k^{\prime}}}V_{\mathbf{0}}f_{v,\mathbf{k^{\prime}}}-\sum_{\mathbf{k^{\prime}}}V^{\text{mol}}_{\mathbf{k-k^{\prime}}}f_{c,\mathbf{k^{\prime}}}-\sum_{\mathbf{k^{\prime}}}V_{\mathbf{k-k^{\prime}}}f_{v,\mathbf{k^{\prime}}}\right)c^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}$
$\displaystyle+\sum_{\mathbf{k}}\left(\varepsilon_{v,\mathbf{k}}+\sum_{\mathbf{k^{\prime}}}V^{\text{WS}_{2}}_{\mathbf{0}}f_{v,\mathbf{k^{\prime}}}+\sum_{\mathbf{k^{\prime}}}V_{\mathbf{0}}f_{c,\mathbf{k^{\prime}}}-\sum_{\mathbf{k^{\prime}}}V^{\text{WS}_{2}}_{\mathbf{k-k^{\prime}}}f_{v,\mathbf{k^{\prime}}}-\sum_{\mathbf{k^{\prime}}}V_{\mathbf{k-k^{\prime}}}f_{c,\mathbf{k^{\prime}}}\right)v^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}$
$\displaystyle+ie_{0}\mathbf{E}(t)\cdot\sum_{\lambda,\mathbf{k}}\nabla_{\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}}$
$\displaystyle-\sum_{\mathbf{k}}\left(\sum_{\mathbf{k}^{\prime}}V_{\mathbf{k-k^{\prime}}}p^{*}_{\mathbf{k^{\prime}}}-\mathbf{d}_{\mathbf{k}}\cdot\mathbf{E}(t)\right)v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}-\sum_{\mathbf{k}}\left(\sum_{\mathbf{k}^{\prime}}V_{\mathbf{k-k^{\prime}}}p_{\mathbf{k^{\prime}}}-\mathbf{d}^{*}_{\mathbf{k}}\cdot\mathbf{E}(t)\right)c^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}$
(5)
$\displaystyle=\sum_{\lambda,\mathbf{k}}\tilde{\varepsilon}_{\lambda,\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}}-\sum_{\mathbf{k}}\Delta_{\mathbf{k}}\left(v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}+c^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}\right)+ie_{0}\sum_{\lambda,\mathbf{k}}\mathbf{E}(t)\cdot(\nabla_{\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}})\lambda^{\mathstrut}_{\mathbf{k}}+\hbar\sum_{\mathbf{k}}\Omega_{\mathbf{k}}(t)\left(v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}+c^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}\right).$
(6)

which reproduces the semiconductor Bloch equations in the Hartree-Fock limit,
including effects as exciton formation or band gap renormalization. The
single-particle energies $\tilde{\varepsilon}_{\lambda,\mathbf{k}}$ are
renormalized by Hartree-Fock contributions of intra- and interlayer Coulomb
interaction, which contain the carrier occupation
$f_{\lambda,\mathbf{k}}=\langle\lambda^{\dagger}_{\mathbf{k}}\lambda^{\mathstrut}_{\mathbf{k}}\rangle$
in TMDC and molecule layer and
$\Delta_{\mathbf{k}}=\sum_{\mathbf{k}^{\prime}}V_{\mathbf{k-k^{\prime}}}p_{\mathbf{k}^{\prime}}$
accounts for the formation of bound excitons and corresponds to the built up
of a macroscopic coherence in the EI state. The microscopic transitions are
defined as $p_{\mathbf{k}}=\langle
v^{\dagger}_{\mathbf{k}}c^{\mathstrut}_{\mathbf{k}}\rangle$ and
$V_{\mathbf{k}}$ denotes the interlayer Coulomb potential Ovesen et al.
(2019). The last two terms include light-matter interaction consisting of
intra- and interband transitions. The latter are determined by the Rabi
frequency
$\Omega_{\mathbf{k}}=\mathbf{d}_{\mathbf{k}}\cdot\mathbf{E}(t)/\hbar$ with
electronic dipole moment $\mathbf{d}_{\mathbf{k}}$ and electric field
$\mathbf{E}(t)$. In this manuscript, we use the band gap energy as parameter,
tunable by a static external additional out-of-plane electric field via the
Stark effect Lorchat et al. (2021). This tuning enables different electronic
phases of the heterostructure. Since the excitonic insulator phase results
from a spontaneous formation of excitons, we first focus on the ground state
calculation in absence of an external exciting optical field.

The Hamiltonian is then diagonalized with a Bogoliubov transformation. The new
ground state is given by Jérome et al. (1967)

$\displaystyle|\tilde{\Psi}_{0}\rangle=\Pi_{\mathbf{k}}\left(u^{*}_{\mathbf{k}}-w^{*}_{\mathbf{k}}c^{\dagger}_{\mathbf{k}}v^{\mathstrut}_{\mathbf{k}}\right)|\Psi_{0}\rangle=\Pi_{\mathbf{k}}\alpha^{\dagger}_{\mathbf{k}}|0\rangle$
(7)

with $|\Psi_{0}\rangle=\Pi_{\mathbf{k}}v^{\dagger}_{\mathbf{k}}|0\rangle$ as
the conventional semiconducting ground state constructed from the vacuum state
$|0\rangle$. The coherence factors $|w_{\mathbf{k}}|^{2}$ and
$|u_{\mathbf{k}}|^{2}$ describe the probabilities that the pair state is
occupied or unoccupied, respectively. The diagonalization yields the coherence
factors

$\displaystyle|u_{\mathbf{k}}|^{2}$
$\displaystyle=\frac{1}{2}\left(1+\frac{\Sigma_{\mathbf{k}}}{\sqrt{\Sigma_{\mathbf{k}}^{2}+\Delta_{\mathbf{k}}^{2}}}\right),\quad\text{and}\quad|w_{\mathbf{k}}|^{2}=\frac{1}{2}\left(1-\frac{\Sigma_{\mathbf{k}}}{\sqrt{\Sigma_{\mathbf{k}}^{2}+\Delta_{\mathbf{k}}^{2}}}\right)\;.$
(8)

Finally, we define the Coulomb potential

$\displaystyle V_{\mathbf{q}}^{ll^{\prime}}$
$\displaystyle=\frac{e_{0}^{2}}{2\epsilon_{0}A|\mathbf{q}|\epsilon^{ll^{\prime}}_{\mathbf{q}}};\quad\epsilon_{\mathbf{q}}^{ll^{\prime}}=\left\\{\begin{array}[]{ll}\epsilon_{\mathbf{q}},&l\neq
l^{\prime}\\\ \epsilon^{i}_{\mathbf{q}},&l=l^{\prime}\equiv
i\end{array}\right.$ (11)

with $i=\\{0,1\\}$ for molecular and TMDC layer. The dielectric functions read

$\displaystyle\epsilon_{\mathbf{q}}=\kappa
g^{0}_{|\mathbf{q}|}g^{1}_{|\mathbf{q}|}f_{|\mathbf{q}|}\quad\text{and}\quad\epsilon^{i}_{\mathbf{q}}$
$\displaystyle=\frac{\kappa
g^{1-i}_{|\mathbf{q}|}f_{|\mathbf{q}|}}{\cosh(\delta_{1-i}{|\mathbf{q}|}/2)h^{i}_{|\mathbf{q}|}}$
(12)

with the abbreviations

$\displaystyle f_{\mathbf{q}}$
$\displaystyle=1+\frac{1}{2}\left((\frac{\kappa_{0}}{\kappa}+\frac{\kappa}{\kappa_{0}})\tanh(\delta_{0}|\mathbf{q}|)+(\frac{\kappa_{1}}{\kappa}+\frac{\kappa}{\kappa_{1}})\tanh(\delta_{1}|\mathbf{q}|)+(\frac{\kappa_{0}}{\kappa_{1}}+\frac{\kappa_{1}}{\kappa_{0}})\tanh(d_{0}|\mathbf{q}|)\tanh(\delta_{1}|\mathbf{q}|)\right)$
(13) $\displaystyle h_{\mathbf{q}}^{i}$
$\displaystyle=1+\frac{\kappa}{\kappa_{i}}\tanh(\delta_{i}|\mathbf{q}|)+\frac{\kappa}{\kappa_{1-i}}\tanh(\delta_{1-i}|\mathbf{q}|/2)+\frac{\kappa_{i}}{\kappa_{1-i}}\tanh(\delta_{i}|\mathbf{q}|)\tanh(\delta_{1-i}|\mathbf{q}|/2)$
(14) $\displaystyle g_{\mathbf{q}}^{i}$
$\displaystyle=\frac{\cosh(\delta_{i}|\mathbf{q}|)}{\cosh(\delta_{1-i}|\mathbf{q}|/2)}\left(1+\frac{\kappa}{\kappa_{i}}\tanh(\delta_{i}|\mathbf{q}|/2)\right).$
(15)

The parameters are
$\kappa_{i}=\sqrt{\epsilon_{\parallel}^{i}\epsilon_{\perp}^{i}}$ and $\kappa$
for the dielectric background,
$\alpha_{i}=\sqrt{\epsilon^{i}_{\parallel}/\epsilon^{i}_{\perp}}$,
$\delta_{i}=\alpha_{i}d_{i}$ with the layer thickness $d_{i}$.

##  Appendix B Equations of motion

The macroscopic polarization is defined as $\mathbf{P}(t)=-\delta
H_{lm}/\delta\mathbf{E}(t)$ with the light-matter Hamiltonian $H_{lm}$.
Already in excitonic basis, the macroscopic polarization reads

$\displaystyle\mathbf{P}(t)=-\sum_{\mu,\mathbf{k}}\left(f^{(0)}_{\mathbf{k}}\varphi_{\mu,\mathbf{k}}\mathbf{d}_{\mathbf{k}}+ie_{0}\varphi^{*}_{\mu,\mathbf{k}}\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}\right)P_{\mu}+ie_{0}\sum_{\lambda,\mathbf{k}}[\nabla_{\mathbf{k}}\lambda^{\dagger}_{\mathbf{k}}]\lambda^{\mathstrut}_{\mathbf{k}}$
(16)

from the macroscopic polarization we can calculate the optical current

$\displaystyle\mathbf{j}(t)$
$\displaystyle=-\sum_{\mu,\mathbf{k}}\left(\mathbf{d}_{\mathbf{k}}\varphi_{\mu,\mathbf{k}}f^{(0)}_{\mathbf{k}}+ie_{0}\varphi_{\mu,\mathbf{k}}\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}\right)\frac{d}{dt}P_{\mu}(t)+e_{0}\sum_{\mathbf{k}}\mathbf{v}_{\mathbf{k}}F^{(1)}_{v,\mathbf{k}}(t)$
(17) $\displaystyle\mathbf{j}(\omega)$
$\displaystyle=-\sum_{\mu,\mathbf{k}}\left(\mathbf{d}_{\mathbf{k}}\varphi_{\mu,\mathbf{k}}f^{(0)}_{\mathbf{k}}+ie_{0}\varphi_{\mu,\mathbf{k}}\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}\right)i\omega
P_{\mu}(\omega)+e_{0}\sum_{\mathbf{k}}\mathbf{v}_{\mathbf{k}}F^{(1)}_{v,\mathbf{k}}(\omega)$
(18)

where we introduced the particle velocity
$\mathbf{v}_{\mathbf{k}}=\hbar\mathbf{k}/m$ and Fourier transformed to get
from Eq. (17) to Eq. (18). Since the molecule electrons are infinitely heavy
they do not contribute to the current that only the valence band electrons
$F_{v,\mathbf{k}}^{(1)}$ are relevant. We derive the equations of motion for
the optical excitations, which read

$\displaystyle i\hbar\frac{d}{dt}p_{\mathbf{k}}$
$\displaystyle=\left(2\Sigma_{\mathbf{k}}+ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}\right)p_{\mathbf{k}}-\Delta_{\mathbf{k}}f_{\mathbf{k}}+\hbar\Omega_{\mathbf{k}}f_{\mathbf{k}}$
(19) $\displaystyle i\hbar\frac{d}{dt}f_{\mathbf{k}}$
$\displaystyle=2i\Im\textnormal{m}\left(\Delta_{\mathbf{k}}p_{\mathbf{k}}\right)+e_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}f_{\mathbf{k}}$
(20)

with
$2\Sigma_{\mathbf{k}}=\tilde{\varepsilon}_{c,\mathbf{k}}-\tilde{\varepsilon}_{v,\mathbf{k}}$
and inversion $f_{\mathbf{k}}=f_{v,\mathbf{k}}-f_{c,\mathbf{k}}$. Then, we
expand the polarization and density into orders of the exciting electric field

$\displaystyle p_{\mathbf{k}}$
$\displaystyle=p^{(0)}_{\mathbf{k}}+p^{(1)}_{\mathbf{k}}+\mathcal{O}(2)$ (21)
$\displaystyle f_{\mathbf{k}}$
$\displaystyle=f^{(0)}_{\mathbf{k}}+f^{(1)}_{\mathbf{k}}+\mathcal{O}(2).$ (22)

Staying in the first order of the electric field the HIOS Bloch equation reads

$\displaystyle i\hbar\frac{d}{dt}p^{(1)}_{\mathbf{k}}$
$\displaystyle=2\Sigma^{(0)}_{\mathbf{k}}p^{(1)}_{\mathbf{k}}-\Delta^{(1)}_{\mathbf{k}}f^{(0)}_{\mathbf{k}}+2\Sigma^{(1)}_{\mathbf{k}}p^{(0)}_{\mathbf{k}}-\Delta^{(0)}_{\mathbf{k}}f^{(1)}_{\mathbf{k}}+\hbar\Omega_{\mathbf{k}}f^{(0)}_{\mathbf{k}}+ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}$
(23) $\displaystyle i\hbar\frac{d}{dt}f^{(1)}_{\mathbf{k}}$
$\displaystyle=2i\Im\left(\Delta^{(0)}_{\mathbf{k}}p^{(1)}_{\mathbf{k}}+\Delta^{(1)}_{\mathbf{k}}p^{(0)}_{\mathbf{k}}\right)+ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}f^{(0)}_{\mathbf{k}}.$
(24)

Since $p^{(0)}_{\mathbf{k}}$ and $f^{(0)}_{\mathbf{k}}$ describe the ground
state their dynamics vanish. The first two terms of Eq. (23) correspond to the
semiconductor Bloch equation. The second two terms are new and couple ground
state distributions and excited quantities. The last two terms describe the
optical excitation with interband source with included Pauli blocking Selig et
al. (2019) and intraband source, respectively. The coupling induced by the
third and fourth term can be resolved via the transformation Stroucken and
Koch (2015)

$\displaystyle P^{(1)}_{\mathbf{k}}$
$\displaystyle=\frac{E_{\mathbf{k}}+\Sigma^{(0)}_{\mathbf{k}}}{E_{\mathbf{k}}}p^{(1)}_{\mathbf{k}}-\frac{E_{\mathbf{k}}-\Sigma^{(0)}_{\mathbf{k}}}{E_{\mathbf{k}}}p^{*(1)}_{\mathbf{k}}-\frac{\Delta^{(0)}_{\mathbf{k}}}{E_{\mathbf{k}}}f^{(1)}_{\mathbf{k}}$
(25) $\displaystyle F^{(1)}_{\mathbf{k}}$
$\displaystyle=\frac{\Delta^{(0)}_{\mathbf{k}}}{E_{\mathbf{k}}}\left(p^{(1)}_{\mathbf{k}}+p^{*(1)}_{\mathbf{k}}\right)+\frac{\Sigma^{(0)}_{\mathbf{k}}}{E_{\mathbf{k}}}f^{(1)}_{\mathbf{k}}.$
(26)

We obtain

$\displaystyle i\hbar\frac{d}{dt}P^{(1)}_{\mathbf{k}}$
$\displaystyle=2E_{\mathbf{k}}P^{(1)}_{\mathbf{k}}-\sum_{\mathbf{k^{\prime}}}V_{\mathbf{k-k^{\prime}}}P^{(1)}_{\mathbf{k}^{\prime}}+\hbar\Omega_{\mathbf{k}}f^{(0)}_{\mathbf{k}}+ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}$
(27) $\displaystyle i\hbar\frac{d}{dt}F^{(1)}_{\mathbf{k}}$
$\displaystyle=ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}f^{(0)}_{\mathbf{k}}$
(28)

We can identify the Bogoliubov-Wannier equation

$\displaystyle
2E_{\mathbf{k}}\varphi_{\mu,\mathbf{k}}-\sum_{\mathbf{k^{\prime}}}V_{\mathbf{k-k^{\prime}}}\varphi_{\mu,\mathbf{k}^{\prime}}=E_{\mu}\varphi_{\mu,\mathbf{k}}$
(29)

which yields the excitonic equations

$\displaystyle i\hbar\frac{d}{dt}P^{(1)}_{\mu}$
$\displaystyle=E_{\mu}P^{(1)}_{\mu}+\sum_{\mathbf{k}}f^{(0)}_{\mathbf{k}}\varphi^{*}_{\mu,\mathbf{k}}\hbar\Omega_{\mathbf{k}}+ie_{0}\mathbf{E}(t)\cdot\sum_{\mathbf{k}}\varphi^{*}_{\mu,\mathbf{k}}\nabla_{\mathbf{k}}p^{(0)}_{\mathbf{k}}$
(30) $\displaystyle i\hbar\frac{d}{dt}F^{(1)}_{\mathbf{k}}$
$\displaystyle=ie_{0}\mathbf{E}(t)\cdot\nabla_{\mathbf{k}}f^{(0)}_{\mathbf{k}}$
(31)

Equation (30) can be Fourier transformed an inserted into Eq. (18). For the
occupations we Fourier transform Eq. (31) with phenomenological dephasing and
insert into Eq. (18), which yields for the last term

$\displaystyle\mathbf{j}(\omega)=-e_{0}\sum_{\mathbf{k}}\mathbf{v}_{\mathbf{k}}\frac{ie_{0}\mathbf{E}(t)\nabla_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}}{\hbar\omega+i\gamma}=\left(-\frac{ie_{0}^{2}}{\hbar}\sum_{\mathbf{k}}\frac{\mathbf{v}_{\mathbf{k}}\otimes\nabla_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}}{\omega+i\gamma/\hbar}\right)\cdot\mathbf{E}(t)$
(32)

Comparing with the definition of the current
$\mathbf{j}(\omega)=-i\omega\epsilon_{0}\chi(\omega)\mathbf{E}(\omega)$ we can
identify the scalar susceptibility

$\displaystyle\chi(\omega)=-\frac{1}{\epsilon_{0}}\sum_{\mu}\frac{\mathbf{d}_{\mu}\otimes\mathbf{d}_{\mu}+\mathbf{j}_{\mu}\otimes\mathbf{j}_{\mu}}{\hbar\omega-
E_{\mu}+i\gamma}+\frac{e_{0}^{2}}{\epsilon_{0}\hbar}\sum_{\mathbf{k}}\frac{\mathbf{v}_{\mathbf{k}}\otimes\nabla_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}}{\omega^{2}+i\omega\gamma/\hbar}.$
(33)

In case that the electronic phase has a valence occupation
$f^{(0)}_{v,\mathbf{k}}$ with Fermi edge, we can partially integrate the last
term in Eq. (32) and find

$\displaystyle\left(-\frac{ie_{0}^{2}}{\hbar}\sum_{\mathbf{k}}\frac{\mathbf{v}_{\mathbf{k}}\otimes\nabla_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}}{\omega+i\gamma/\hbar}\right)\cdot\mathbf{E}(t)=\frac{ie_{0}^{2}}{m}n_{el}\frac{1}{\omega+i\gamma/\hbar}\mathbf{E}(\omega)$
(34)

where defined the carrier number
$n_{el}=\sum_{\mathbf{k}}f^{(0)}_{v,\mathbf{k}}$. Again comparing with the
definition of the optic current and assuming a perpendicular excitation we
find the scalar susceptibility

$\displaystyle\chi(\omega)=-\frac{1}{\epsilon_{0}}\sum_{\mu}\frac{|d_{\mu}|^{2}+|j_{\mu}|^{2}}{\hbar\omega-
E_{\mu}+i\gamma}+\frac{\omega^{2}_{pl}}{\omega^{2}+i\omega\gamma/\hbar}$ (35)

with the plasma frequency $\omega^{2}_{pl}=e_{0}^{2}n_{el}/\epsilon_{0}m$.

##  Appendix C Optical selection rules of excitonic insulator

Investigating the gap equation Eq. (1) of the main text we see that it
corresponds to the Bogoliubov-Wannier equation with vanishing exciton binding
energy. This suggests that also the ground state polarization could be
projected onto the wave function acting as solution for $E_{\mu}=0$:
$p^{(0)}_{\mathbf{k}}=\sum_{\nu,E_{\nu}=0}\varphi_{\nu,\mathbf{k}}p^{(0)}_{\nu}$.
Then the momentum-gradient in the intraband matrix element acts onto the wave
function. For a $s$-type ground state the angular derivative vanishes.
Together with an analytical treatment of the angle-sum we obtain for the
intraband source

$\displaystyle
j_{\mu}=e_{0}\pi\mathbf{e}\cdot\sum_{\nu,k}\varphi_{\mu,k}^{*}\partial_{k}\varphi_{\nu,k}\left(\begin{array}[]{c}1\\\
\pm i\end{array}\right)$ (38)

where the sign stands for $\mu=p_{+}$ and $\mu=p_{-}$ final states. These two
final states exhibits a circular dichroismic selection rule, comparable to KK-
and K-K--excitons in monolayer TMDCs.

## References

  * Keldysh and Kopaev (1965) L. Keldysh and Y. V. Kopaev, Soviet Physics Solid State, USSR 6, 2219 (1965).
  * Jérome et al. (1967) D. Jérome, T. Rice, and W. Kohn, Physical Review 158, 462 (1967).
  * Gupta et al. (2020) S. Gupta, A. Kutana, and B. I. Yakobson, Nature communications 11, 1 (2020).
  * Perfetto and Stefanucci (2021) E. Perfetto and G. Stefanucci, Physical Review B 103, L241404 (2021).
  * Parmenter and Henson (1970) R. Parmenter and W. Henson, Physical Review B 2, 140 (1970).
  * Wang et al. (2019a) R. Wang, O. Erten, B. Wang, and D. Xing, Nature communications 10, 1 (2019a).
  * Varsano et al. (2017) D. Varsano, S. Sorella, D. Sangalli, M. Barborini, S. Corni, E. Molinari, and M. Rontani, Nature communications 8, 1 (2017).
  * Perfetto et al. (2019) E. Perfetto, D. Sangalli, A. Marini, and G. Stefanucci, Phys. Rev. Materials 3, 124601 (2019).
  * Perfetto et al. (2020) E. Perfetto, A. Marini, and G. Stefanucci, Phys. Rev. B 102, 085203 (2020).
  * Perfetto and Stefanucci (2020) E. Perfetto and G. Stefanucci, Phys. Rev. Lett. 125, 106401 (2020).
  * Jia et al. (2020) Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. Jäck, S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak, et al., arXiv preprint arXiv:2010.05390 (2020).
  * Cercellier et al. (2007) H. Cercellier, C. Monney, F. Clerc, C. Battaglia, L. Despont, M. Garnier, H. Beck, P. Aebi, L. Patthey, H. Berger, et al., Physical review letters 99, 146403 (2007).
  * Monney et al. (2010) C. Monney, E. Schwier, M. G. Garnier, N. Mariotti, C. Didiot, H. Cercellier, J. Marcus, H. Berger, A. Titov, H. Beck, et al., New Journal of Physics 12, 125019 (2010).
  * Wakisaka et al. (2009) Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, Physical review letters 103, 026402 (2009).
  * Lu et al. (2017) Y. Lu, H. Kono, T. Larkin, A. Rost, T. Takayama, A. Boris, B. Keimer, and H. Takagi, Nature communications 8, 1 (2017).
  * Werdehausen et al. (2018) D. Werdehausen, T. Takayama, M. Höppner, G. Albrecht, A. W. Rost, Y. Lu, D. Manske, H. Takagi, and S. Kaiser, Science advances 4, eaap8652 (2018).
  * Bucher et al. (1991) B. Bucher, P. Steiner, and P. Wachter, Physical review letters 67, 2717 (1991).
  * Bronold and Fehske (2006) F. X. Bronold and H. Fehske, Physical Review B 74, 165107 (2006).
  * Fukutani et al. (2021) K. Fukutani, R. Stania, C. I. Kwon, J. S. Kim, K. J. Kong, J. Kim, and H. W. Yeom, arXiv preprint arXiv:2105.11101 (2021).
  * Rossnagel et al. (2002) K. Rossnagel, L. Kipp, and M. Skibowski, Physical Review B 65, 235101 (2002).
  * Zenker et al. (2014) B. Zenker, H. Fehske, and H. Beck, Physical Review B 90, 195118 (2014).
  * Jiang et al. (2018) Z. Jiang, Y. Li, S. Zhang, and W. Duan, Physical Review B 98, 081408 (2018).
  * Volkov et al. (2021) P. A. Volkov, M. Ye, H. Lohani, I. Feldman, A. Kanigel, and G. Blumberg, arXiv preprint arXiv:2104.07032 (2021).
  * Du et al. (2017) L. Du, X. Li, W. Lou, G. Sullivan, K. Chang, J. Kono, and R.-R. Du, Nature communications 8, 1 (2017).
  * Ataei et al. (2021) S. S. Ataei, D. Varsano, E. Molinari, and M. Rontani, Proceedings of the National Academy of Sciences 118 (2021).
  * Ma et al. (2021) L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak, and J. Shan, Nature 598, 585 (2021).
  * Zhang et al. (2021) Z. Zhang, E. C. Regan, D. Wang, W. Zhao, S. Wang, M. Sayyad, K. Yumigeta, K. Watanabe, T. Taniguchi, S. Tongay, et al., arXiv preprint arXiv:2108.07131 (2021).
  * Sun et al. (2019) J. Sun, Y. Choi, Y. J. Choi, S. Kim, J.-H. Park, S. Lee, and J. H. Cho, Advanced Materials 31, 1803831 (2019).
  * Song et al. (2017) Z. Song, T. Schultz, Z. Ding, B. Lei, C. Han, P. Amsalem, T. Lin, D. Chi, S. L. Wong, Y. J. Zheng, et al., ACS nano 11, 9128 (2017).
  * Koch (2007) N. Koch, ChemPhysChem 8, 1438 (2007).
  * Fahlman et al. (2019) M. Fahlman, S. Fabiano, V. Gueskine, D. Simon, M. Berggren, and X. Crispin, Nature Reviews Materials 4, 627 (2019).
  * Torabi et al. (2015) S. Torabi, F. Jahani, I. Van Severen, C. Kanimozhi, S. Patil, R. W. Havenith, R. C. Chiechi, L. Lutsen, D. J. Vanderzande, T. J. Cleij, et al., Advanced Functional Materials 25, 150 (2015).
  * Groom et al. (2016) C. R. Groom, I. J. Bruno, M. P. Lightfoot, and S. C. Ward, Acta Crystallographica Section B: Structural Science, Crystal Engineering and Materials 72, 171 (2016).
  * Stuke et al. (2020) A. Stuke, C. Kunkel, D. Golze, M. Todorović, J. T. Margraf, K. Reuter, P. Rinke, and H. Oberhofer, Scientific data 7, 1 (2020).
  * Lorchat et al. (2021) E. Lorchat, M. Selig, F. Katsch, K. Yumigeta, S. Tongay, A. Knorr, C. Schneider, and S. Höfling, Physical Review Letters 126, 037401 (2021).
  * Drude (1900) P. Drude, Annalen der physik 306, 566 (1900).
  * Gillen (2021) R. Gillen, physica status solidi (b) , 2000614 (2021).
  * Deilmann and Thygesen (2018) T. Deilmann and K. S. Thygesen, Nano letters 18, 2984 (2018).
  * Kira and Koch (2006) M. Kira and S. Koch, Progress in quantum electronics 30, 155 (2006).
  * Bhattacharyya et al. (2014) J. Bhattacharyya, S. Zybell, F. Eßer, M. Helm, H. Schneider, L. Schneebeli, C. Böttge, B. Breddermann, M. Kira, S. W. Koch, et al., Physical Review B 89, 125313 (2014).
  * Park et al. (2021) S. Park, H. Wang, T. Schultz, D. Shin, R. Ovsyannikov, M. Zacharias, D. Maksimov, M. Meissner, Y. Hasegawa, T. Yamaguchi, et al., Advanced Materials , 2008677 (2021).
  * Krukau et al. (2006) A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, The Journal of chemical physics 125, 224106 (2006).
  * Blum et al. (2009) V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009).
  * Ren et al. (2012) X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, and M. Scheffler, New Journal of Physics 14, 053020 (2012).
  * Levchenko et al. (2015) S. V. Levchenko, X. Ren, J. Wieferink, R. Johanni, P. Rinke, V. Blum, and M. Scheffler, Computer Physics Communications 192, 60 (2015).
  * Wang et al. (2019b) H. Wang, S. V. Levchenko, T. Schultz, N. Koch, M. Scheffler, and M. Rossi, Advanced Electronic Materials 5, 1800891 (2019b).
  * Malic et al. (2018) E. Malic, M. Selig, M. Feierabend, S. Brem, D. Christiansen, F. Wendler, A. Knorr, and G. Berghäuser, Physical Review Materials 2, 014002 (2018).
  * Bogoljubov et al. (1958) N. Bogoljubov, V. V. Tolmachov, and D. Širkov, Fortschritte der physik 6, 605 (1958).
  * Valatin (1958) J. Valatin, Il Nuovo Cimento (1955-1965) 7, 843 (1958).
  * Keldysh and Kozlov (1968) L. Keldysh and A. Kozlov, Sov. Phys. JETP 27, 521 (1968).
  * Comte and Nozieres (1982) C. Comte and P. Nozieres, Journal de Physique 43, 1069 (1982).
  * Kozlov and Maksimov (1965) A. Kozlov and L. Maksimov, Sov. Phys. JETP 21, 790 (1965).
  * Sabio et al. (2010) J. Sabio, F. y. Sols, and F. Guinea, Physical Review B 82, 121413 (2010).
  * Stroucken et al. (2011) T. Stroucken, J. Grönqvist, and S. Koch, Physical Review B 84, 205445 (2011).
  * Zimmermann and Schindler (2007) R. Zimmermann and C. Schindler, Solid state communications 144, 395 (2007).
  * Schindler and Zimmermann (2008) C. Schindler and R. Zimmermann, Physical Review B 78, 045313 (2008).
  * Gu et al. (2021) J. Gu, L. Ma, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, arXiv preprint arXiv:2108.06588 (2021).
  * Ovesen et al. (2019) S. Ovesen, S. Brem, C. Linderälv, M. Kuisma, T. Korn, P. Erhart, M. Selig, and E. Malic, Communications Physics 2, 1 (2019).
  * Bardeen et al. (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Physical review 108, 1175 (1957).
  * Tinkham (2004) M. Tinkham, Introduction to superconductivity (Courier Corporation, 2004).
  * Reyes et al. (2016) D. Reyes, M. A. Continentino, C. Thomas, and C. Lacroix, Annals of Physics 373, 257 (2016).
  * Knorr et al. (1996) A. Knorr, S. Hughes, T. Stroucken, and S. Koch, Chemical physics 210, 27 (1996).
  * Malic and Knorr (2013) E. Malic and A. Knorr, Graphene and carbon nanotubes: ultrafast optics and relaxation dynamics (John Wiley & Sons, 2013).
  * Katsch and Knorr (2020) F. Katsch and A. Knorr, Physical Review X 10, 041039 (2020).
  * Stroucken and Koch (2015) T. Stroucken and S. W. Koch, Journal of Physics: Condensed Matter 27, 345003 (2015).
  * Katsch et al. (2020) F. Katsch, D. Christiansen, R. Schmidt, S. M. de Vasconcellos, R. Bratschitsch, A. Knorr, and M. Selig, Physical Review B 102, 115420 (2020).
  * Dong et al. (2021a) S. Dong, M. Puppin, T. Pincelli, S. Beaulieu, D. Christiansen, H. Hübener, C. W. Nicholson, R. P. Xian, M. Dendzik, Y. Deng, et al., Natural Sciences , e10010 (2021a).
  * Selig et al. (2016) M. Selig, G. Berghäuser, A. Raja, P. Nagler, C. Schüller, T. F. Heinz, T. Korn, A. Chernikov, E. Malic, and A. Knorr, Nature communications 7, 1 (2016).
  * Christiansen et al. (2017) D. Christiansen, M. Selig, G. Berghäuser, R. Schmidt, I. Niehues, R. Schneider, A. Arora, S. M. de Vasconcellos, R. Bratschitsch, E. Malic, et al., Physical review letters 119, 187402 (2017).
  * Helmrich et al. (2021) S. Helmrich, K. Sampson, D. Huang, M. Selig, K. Hao, K. Tran, A. Achstein, C. Young, A. Knorr, E. Malic, et al., Physical Review Letters 127, 157403 (2021).
  * Slater (1934) J. C. Slater, Reviews of Modern Physics 6, 209 (1934).
  * Kaplan and Ruvinskiǐ (1976) I. Kaplan and M. Ruvinskiǐ, Soviet Journal of Experimental and Theoretical Physics 44, 1127 (1976).
  * Specht et al. (2018) J. F. Specht, E. Verdenhalven, B. Bieniek, P. Rinke, A. Knorr, and M. Richter, Physical Review Applied 9, 044025 (2018).
  * Dong et al. (2021b) S. Dong, S. Beaulieu, M. Selig, P. Rosenzweig, D. Christiansen, T. Pincelli, M. Dendzik, J. D. Ziegler, J. Maklar, R. P. Xian, et al., arXiv preprint arXiv:2108.06803 (2021b).
  * Selig et al. (2019) M. Selig, F. Katsch, R. Schmidt, S. M. de Vasconcellos, R. Bratschitsch, E. Malic, and A. Knorr, Physical Review Research 1, 022007 (2019).


Classical models of computation have been successful in capturing the very essence of individual computing devices. Although they are useful to understand computability power and limitations in the small, such models are not suitable to study large-scale complex computations. Accordingly, plenty of formalisms have been proposed in the last half century as an attempt to raise the level of abstraction, with the aim of describing not only a single computing device but interactions among a collection of them. In this paper, we encompass such formalisms into a common framework which we refer to as Models of High-Level Computation. We particularly discuss the semantics, some of the key properties, paradigms and future directions of such models.

§ INTRODUCTION

Classical models of computation formally emerged during the first half of the century 20th, as an attempt to capture the very essence of computation in the form of information processing. Under the Church-Turing thesis, these models have been successful to describe any possible effective procedure construed as an algorithmic process. Although they can be used to describe complex computations, classical models fail to provide expressive means to do so, apart from not being suitable to accept information streams at computation-time. Describing complex computations has become essential nowadays since the scale and complexity of computing systems is exponentially increasing (especially in the realm of distributed computing). Accordingly, plenty of Models of High-Level Computation (MHCs) have been proposed over the last half century so as to capture the essence of complex computations. In this paper, we describe the meaning, key properties, current paradigms and future directions of MHCs.

The rest of the paper is structured as follows. Section <ref> describes the general notion of MHCs. Section <ref> sketches the design of a simple MHC with the aim of discussing the key elementary properties of such models. Section <ref> presents a preliminary classification of MHCs by considering three major paradigms. Section <ref> outlines the conclusions and future directions of the research area this paper introduces.

§ WHAT IS A MODEL OF HIGH-LEVEL COMPUTATION?

A MHC raises the level of abstraction of its classical, low-level counterpart (e.g., Turing Machines) by providing a birds-eye-view of multiple computing devices rather than focusing on a single one.[We should not confuse high-level with higher-order computations [1]. Higher-order computability is a well-established field in theoretical computer science, whose aim is to study computations that receive and produce other computations.] As individual devices are treated as modular black boxes, their internal details are irrelevant. What matters is how to glue together multiple (low- or even high-level) computations in order to form more complex ones. Thus, a MHC induces modularity and describes how computing devices interact with the aim of processing information that may be initially encoded, that can continuously be produced by an exogenous entity or any combination thereof.[By exogenous we mean an entity that is out of the scope of the interacting computing devices, e.g., a human operator.] That is, a MHC can be open or closed.[Classical, low-level models of computation are inherently closed because they have a fixed input before computing (e.g., on a tape) so outputs can only be read after termination. Open models of computation contrast with this behaviour by allowing the processing of data streams from the external world while computing. An example is a process that waits for a human operator's input before proceeding, just as in the original Turing c-machines [2].]

By the above, it is evident that a MHC needs to intrisically separate computation from interaction. Therefore, a so-called composition mechanism needs to be used. A composition mechanism particularly specifies how computing devices interact from a high-level perspective. An interaction always defines control flow either implicitly or explicitly. It is explicit when there is a clear construct specifying the order of invocation of computing devices, and it is implicit otherwise. Apart from control, an interaction can optionally define implicit or explicit data flow to establish a data exchange scheme among computing devices.

Rather than describing how a single computing device produces an output from a given input, a MHC gives a birds-eye-view on multiple interacting computing devices. For this, a MHC provides formal rules to determine the moment in which a each participant device should compute and how devices must interact.

Typically, computing devices facilitate composition by providing interfaces that serve as endpoints to entry or exit some internal computational structure. Enabling these endpoints is what allows treating computing devices as high-level, modular black boxes which can be composed algebraically or non-algebraically. Only algebraic composition enables compositionality, a necessary property to ensure closure by the provision of higher-order operators that receive a number of computing devices and produce a composite one which can, in turn, be used as an operand to construct even more complex computing devices (see Figure <ref>).

[draw=black] (3.8,0) rectangle ++(2.2,2);
[align=center] at (4.9,1)Composition

at (0,1.7)(Composite) Computing Device $1$;
[label=$\vdots$] at (0,0.55) ;
at (0,0.3)(Composite) Computing Device $M$;
[->] (3,1.7) – (3.8,1.7);
[->] (3,0.3) – (3.8,0.3);

at (9.4,1)Composite Computing Device;
[->] (6,1) – (6.8,1);

Conceptual representation of an algebraic composition operator that takes $M$ computing devices and produces a composite one. Each device operand can be either composite or non-composite.

A MHC (i) can be open or closed, (ii) induces modularity by black-boxing computing devices (with clear interfaces), (iii) separates computation from interaction, (iv) always defines explicit or implicit control flow, (v) optionally defines implicit or explicit data flow, and (vi) enables interaction among computing devices via a  composition mechanism (which can be algebraic or non-algebraic).

§ DESIGNING A SIMPLE MHC

For the sake of argument (even if not computationally powerful), in this section we design a simple MHC to demonstrate the key properties discussed in Section <ref>. For this, we consider computing devices in the form of Nondeterministic Finite Automata (NFAs) as well as composition operators for concatenating ($\oplus$) and parallelising ($\otimes$) — see Definition <ref>. The behaviour of $\oplus$ and $\otimes$ is directly derivable from the respective proofs of closure under concatenation and under union of regular languages [3], as shown by Definitions <ref> and <ref>.

Let $\mathbb{A}$ be the universe of NFAs. A NFA ${N \in \mathbb{A}}$ is a tuple ${(\Sigma,S,s,\delta,F)}$ where $\Sigma$ is a finite set of input symbols which always contains the empty string $\epsilon$, $S$ is a finite set of states, ${s \in S}$ is called an initial state, ${\delta:S \times \Sigma \rightarrow \mathcal{P}(S)}$ is a transition function and ${F \subseteq S}$ is an empty or nonempty set of final states. We say that $N$ accepts a string ${w=x_1 x_2 \ldots x_n}$ over $\Sigma$ if there is a sequence ${(p_0,p_1,\ldots,p_n) \in S^n}$ satisfying the following conditions:

* ${p_0=s}$,

* ${p_{i+1} \in \delta(p_i,x_{i+1})}$ for ${i=0,\ldots,n-1}$ and

* ${p_n \in F}$.

The concatenation operator ${\oplus: \mathbb{A} \times \mathbb{A} \rightarrow \mathbb{A}}$ is a binary function that receives two NFAS, ${N_1=(\Sigma_1,S_1,s_1,\delta_1,F_1)}$ and ${N_2=(\Sigma_2,S_2,s_2,\delta_2,F_2)}$, and produces a concatenative NFA ${N_1 \oplus N_2=(\Sigma_0,S_0,s_0,\delta_0,F_0)}$ where ${\Sigma_0=\Sigma_1\cup\Sigma_2}$, ${S_0=S_1\cup S_2}$, $s_0$ is the initial state $s_1$ of $N_1$, ${F_0=F_2}$ and:

\[
\delta_0(s,x) =
\begin{cases}
\delta_1(s,x) & \text{ if } s \in S_1 \text{ and } s \notin F_1 \\
\delta_1(s,x) & \text{ if } s \in F_1 \text{ and } x\neq\epsilon \\
\delta_1(s,x)\cup\{s_2\} & \text{ if } s \in F_1 \text{ and } x=\epsilon \\
\delta_2(s,x) & \text{ if } s \in S_2
\end{cases}
\]

for any state ${s \in S_0}$ and any symbol ${x \in \Sigma_0}$.

The parallelising operator ${\otimes: \mathbb{A} \times \mathbb{A} \rightarrow \mathbb{A}}$ is a binary function that receives two NFAS, ${N_1=(\Sigma_1,S_1,s_1,\delta_1,F_1)}$ and ${N_2=(\Sigma_2,S_2,s_2,\delta_2,F_2)}$, and produces a parallel NFA ${N_1 \otimes N_2=(\Sigma_0,S_0,s_0,\delta_0,F_0)}$ where ${\Sigma_0=\Sigma_1\cup\Sigma_2}$, ${S_0=S_1\cup S_2 \cup\{s_0\}}$, ${F_0=F_1\cup F_2}$ and:

\[
\delta_0(s,x) =
\begin{cases}
\delta_i(s,x) & \text{ if } s \in S_i \text{ for } i=1,2 \\
\{s_1,s_2\} & \text{ if } s=s_0 \text{ and } x=\epsilon \\
\emptyset & \text{ if } s=s_0 \text{ and } x \neq \epsilon
\end{cases}
\]

for any state ${s \in S_0}$ and any symbol ${x \in \Sigma_0}$.

After defining the algebraic composition operators for our MHC, let us now consider the NFA $N_1$ depicted in Figure <ref> which recognises the language ${L_1=\{w\}}$ where $w$ is a string over $\{a,b\}$ containing the symbol “b" at the third position (from right to left) such as $abaabaa$ [3]. Also, consider the NFA $N_2$ depicted in Figure <ref> which recognises the language ${L_2=\{a^m b^n\}}$ for ${m>0}$ and ${n \geq 0}$.

[draw=black!38!white] (-1.1,-0.6) rectangle ++(7.2,2.2);
[state, initial] (p0) $p_0$;
[state, right of=p0] (p1) $p_1$;
[state, right of=p1] (p2) $p_2$;
[state, accepting, right of=p2] (p3) $p_3$;
(p0) edge[loop above] node$a,b$ (p0)
(p0) edge[above] node$b$ (p1)
(p1) edge[above] node$a,b$ (p2)
(p2) edge[above] node$a,b$ (p3);

NFA $N_1$.

[draw=black!38!white] (-1.1,-0.6) rectangle ++(3.6,2.2);

[state, initial] (q0) $q_0$;
[state, accepting, right of=q0] (q1) $q_1$;
(q0) edge[loop above] node$a$ (q0)
(q0) edge[above] node$a$ (q1)
(q1) edge[loop above] node$b$ (q1);

NFA $N_2$.

NFAs serving as building blocks to construct more complex NFAs.

Treating $N_1$ and $N_2$ as high-level computations (by just considering their respective initial and final states), enable us to compose them into more complex automata via the operators described in Definitions <ref> and <ref>. Particularly, $N_1$ and $N_2$ can be composed into the concatenative NFA ${N_1 \oplus N_2}$ to accept strings in the language ${\{wa^mb^n\}}$. The structure of such a composite is depicted in Figure <ref>. Figure <ref> shows that $N_1$ and $N_2$ can also be composed into the parallel NFA $N_1 \otimes N_2$ to recognising the language ${L_1 \cup L_2}$. As it is well-known that regular languages are closed under concatenation and union, both ${N_1 \oplus N_2}$ and ${N_1 \otimes N_2}$ can in turn be inductively composed into more complex NFAs using the same composition operators. That is, our MHC satisfies closure. It also satisfies the closed-world property since it is well-known that every (classical) NFA does not admit external data streams while computing.

[draw=black!38!white] (-1.6,-1.1) rectangle ++(12.2,3.2);

[draw=black!38!white] (-1.1,-0.6) rectangle ++(7.2,2.2);
[state, initial] (p0) $p_0$;
[state, right of=p0] (p1) $p_1$;
[state, right of=p1] (p2) $p_2$;
[state, right of=p2] (p3) $p_3$;
(p0) edge[loop above] node$a,b$ (p0)
(p0) edge[above] node$b$ (p1)
(p1) edge[above] node$a,b$ (p2)
(p2) edge[above] node$a,b$ (p3);

[draw=black!38!white] (6.8,-0.6) rectangle ++(3.3,2.2);
[state, right of=p3, xshift=0.4cm] (q0) $q_0$;
[state, accepting, right of=q0] (q1) $q_1$;
(q0) edge[loop above] node$a$ (q0)
(q0) edge[above] node$a$ (q1)
(q1) edge[loop above] node$b$ (q1);

(p3) edge[above] node$\epsilon$ (q0);

Concatenative NFA ${N_1 \oplus N_2}$ (as per Definition <ref>).

[state, initial] (r0) $r_0$;
[draw=black!38!white] (-1.1,-2.55) rectangle ++(9.55,5.9);

[draw=black!38!white] (1.2,0.6) rectangle ++(6.7,2.2);
[state, right of=r0, yshift=1.2cm] (p0) $p_0$;
[state, right of=p0] (p1) $p_1$;
[state, right of=p1] (p2) $p_2$;
[state, accepting, right of=p2] (p3) $p_3$;
(p0) edge[loop above] node$a,b$ (p0)
(p0) edge[above] node$b$ (p1)
(p1) edge[above] node$a,b$ (p2)
(p2) edge[above] node$a,b$ (p3);

[draw=black!38!white] (1.2,-2) rectangle ++(3.1,2.2);
[state, below of=p0, yshift=-0.8cm] (q0) $q_0$;
[state, accepting, right of=q0] (q1) $q_1$;
(q0) edge[loop above] node$a$ (q0)
(q0) edge[above] node$a$ (q1)
(q1) edge[loop above] node$b$ (q1);

(r0) edge[left] node$\epsilon$ (p0)
(r0) edge[left] node$\epsilon$ (q0);

Parallel NFA ${N_1 \otimes N_2}$ (as per Definition <ref>).

Composing the NFAs of Figure <ref> in two different ways. An enclosing box means that the internal details of an automaton are abstracted away. Only initial and final states are visible from a high-level perspective. That is, both $N_1$ and $N_2$ are treated as black boxes.

Control flow between composed NFAs is implicitly defined in the added state transitions of a composite NFA. For instance, the transition ${\delta(p_3,\epsilon)=\{q_0\}}$ of ${N_1\oplus N_2}$ serves to implicitly pass control from $N_1$ to $N_2$, whereas the transition ${\delta(r_0,\epsilon)=\{p_0,q_0\}}$ of ${N_1 \otimes N_2}$ defines implicit control for the simultaneous activation of $N_1$ and $N_2$.

Although data flow is not directly supported by our MHC, it is important to analyse how data is processed within a composite automaton. In the case of a concatenative NFA, an input string is esentially computed in two different chunks. In a parallel NFA, there is no data partitioning but each composed NFA computes the whole input independently. For instance, if $aabaaaab$ needs to be processed by ${N_1\otimes N_2}$ then both $N_1$ and $N_2$ would compute and reject on $aabaaaab$ simultaneously. Injecting the same input into ${N_1\oplus N_2}$ would lead to $N_1$ passing control to $N_2$ just after implicitly determining that $aabaa$ is in $L_1$. After receiving control, $N_2$ would just compute and accept on $aab$. If we were to extend our MHC with the notion of explicit data flows (e.g., by introducing direct message passing), ${N_1\oplus N_2}$ could sequentially pass $aab$ from $N_1$ to $N_2$. By the same token, we could also enable data replication within ${N_1\otimes N_2}$ to make sure that both $N_1$ and $N_2$ receive a copy of the whole input before computing.

§ WHAT CLASSES OF MHCS DO CURRENTLY EXIST?

In this section, we classify MHCs into three different classes: state-oriented, data-oriented and control-oriented. For each of them, we provide a brief description and present a few examples.

We devise three major classes of MHCs, namely state-oriented, data-oriented and control-oriented.

State-oriented MHCs are built upon traditional Finite State Machines (FSMs) to enable the description of complex computing systems. The most representative models in this class are the so-called Communicating Finite State Machines (CFSMs) [4] and Hierarchical Finite State Machines (HFSMs) [5].[Note that there also is a hybrid variant called Communicating Hierarchical State Machines [6].] CFSMs introduce “send" and “receive" operations in their semantics to enable concurrent interactions among distinct FSMs via (unbounded) First-In-First-Out (FIFO) channels. The overall behaviour of a system is described through traces of configurations each of which is a tuple of both the state of the FSMs involved and the content of FIFO channels. Extensions of CFSMs include open CFSMs [7] and parameterised CFSMs [8]. The other subclass of state-oriented MHCs is that of HFSMs which allow the specification of nested FSMs in the form of so-called superstates. Due to their simplicity to modelling complex high-level computations, HFSMs have influenced the design of several industry-oriented formalisms such as the UML Superstructure specification which, in turn, has become the de facto standard for modelling complex software systems. In Section <ref>, we designed a simple MHC which can be classified as a HFSM. Other examples pertaining to this class are Hierarchical Featured State Machines [9] and scope-dependent HFSMs [10].

Data-oriented MHCs originated in late sixtees as an attempt to model concurrent computations in the form of direct message passing among a collection of interacting computing devices (usually referred to as actors). In these models, control flow is not explicit but implicit in the collaborative exchange of data. Here, data-driven computations are typically expressed as directed graphs in which nodes and edges denote data processing actors and explicit communication channels for directionally exchanging data streams, respectively [11].[Some data-oriented MHCs (e.g., Token Flows [12]) introduce special actors to enable control constructs such as branching.] Perhaps the most representative MHCs pertaining to this class are Kahn Process Networks (KPNs) [13] and the actor model [14]. On the one hand, KPNs define denotational semantics for modelling (high-level) concurrent actors that communicate through FIFO queues. Although operational semantics have been proposed to explicitly describe the flow of data values [11], it has been showed that such a semantics coincide with the original denotational notion of KPNs [15]. The operational idea is that there is a static network of actors in which every actor consumes input tokens, perfoms some computation on those tokens and produces output tokens that can potentially be consumed by other actors. The actor model does not operate on a static network of data processing actors, but it allows actors to dynamically create other actors, send messages to other actors and decide how to handle data. Other data-oriented MHCs that have been proposed over time (many of them variants of KPNs) include Synchronous Data Flow [16], dataflow process networks [17], reactive process networks [18] and synchronous blocks [19]. Interestingly, there are some MHCs (e.g., Reo [20]) which do not allow direct data passing among computing devices, but data exchanges are “coordinated exogenously" via well-defined algebraic dataflow structures.

Control-oriented MHCs enable the specification of explicit order in which individual computing devices need to be invoked so that interacting devices behave as passive units of computation. Examples of MHCs belonging to this class are exogenous connectors for encapsulated components [21], Behaviour Trees [22] and workflow control flow models [23]. Particularly, exogenous connectors allow the formation of hierarchical control flow structures that define (exogenous) coordination for the invocation of different computing devices. A few extensions to this model have been proposed over time (e.g., [24, 25]). Behaviour trees are similar to exogenous connectors in the sense that hierarchical control flow structures are formed. The difference lies in the operational semantics of control flow coordination. Over time, several extensions and enhacements of Behaviour Trees have been proposed (see [26]), mainly to modelling modular Robot behaviour. Apart from exogenous connectors and Behaviour Trees, we also have workflow control flow models which define formal rules to govern the computation of workflows. Here, control flow specifies the order in which workflow activities are activated. A workflow activity is a fundamental unit of computation which can either be indivisible or contain other activities. Specific examples of workflow control flow models include  Workflow Nets [27], formal BPEL process models [28], YAWL [29], among others.[A workflow engine is an abstract machine able to receive a workflow control flow model $M$ and a workflow specification $S$, before computing the activities in $S$ according to the rules imposed by $M$.]

It is important to note that the three classes of MHCs we just considered are not exhaustive but indicative of the vast range of MHCs that have been proposed over time in the need of describing complex computations. For instance, we did not consider process calculi [30] which are KPN-influenced models that borrow properties from both data- and control-oriented MHCs, in order to enable the description of inter-process communication through well-defined algebraic laws for equational reasoning. In the last years, there has been an increasing tendency to move towards describing high-level processes via algebras over operads of wiring diagrams [31] which provide formal constructs to reasoning about functional and concurrent computations in an intuitive yet rigorous manner, typically in the language of symmetric monoidal categories wherein morphisms express high-level computations that can graphically be depicted as string diagrams [32].[String diagrams have their roots in the Penrose graphical notation for tensor networks [33].] We believe that such category-theoretic operadic models are collectively forging the foundations of promising compositional MHCs. Examples within this paradigm include the calculus of signal flows [34], the resource calculus [35] and the zx-calculus [36]. Although not operadic in nature, other formalisms deserving attention are the so-called component models [37] since many of them form MHCs themselves. For instance, Reo [20] and exogenous connectors [21] are component models that yield data- and control-oriented MHCs, respectively, as mentioned previously.

§ CONCLUSIONS AND FUTURE DIRECTIONS

In this paper, we introduced a new class of models of computation referred to as MHCs. Contrary to their low-level counterpart, the purpose of a MHC is to define interactions among diverse computing devices through a certain composition mechanism, so that interaction results from composition (not the other way round). As interactions occur outside the internal structure of each composed computing device, a MHC separates computation from interaction and induces modularity by treating computing devices as black bloxes. Composition can be done either algebraically or non-algebraically in order to specify (explicit/implicit) control flow and, optionally, (explicit/implicit) data flow.

In Section <ref>, we presented a closed MHC for facilitating the interaction of NFAs through an algebraic composition mechanism that defines implicit control flow only. Through this example, we demonstrated the modularity property that permeates the notion of MHCs. Although our examplary MHC satisfies closure over NFAs, there are other MHCs that do not necessarily compose state machines so as to deal with the well-known state explosion problem, e.g., [20, 21]. There even are MHCs that provide mechanisms to sequentially compose Turing Machines [38]. In Section <ref>, we classified MHCs into three major classes: state-oriented, data-oriented and control-oriented. Although our review is introductory only, it can serve as a starting point to devise a more complete classification scheme (or even a whole taxonomy). For a more detailed analysis and comparison of MHC classes, one can consider the properties presented in Section <ref> and beyond such as determinism vs non-determinism or algebraic versus non-algebraic.

We believe that algebraically composing computing devices is of paramount importance to tame the complexity of interactions, especially as the size of computing devices becomes larger and larger. Algebraic composition can also be beneficial to compositionally verify certain properties such as termination and reachability. Nevertheless, there are number of directions that need to be addressed before unleashing the full potential of MHCs, including (i) compositional semantics for concurrent high-level computations, (ii) design of programming languages built upon the notion of MHCs and (iii) classification of expressive power (especially of open MHCs), just to name a few directions.

Moving up the ladder of abstraction from low- to high-level computations resembles the paradigm shift from low- to high-level programming languages. It is also analogous to the change in perspective from concrete mathematical structures to high-level ones (as in Category Theory). The paradigm shift we present in this paper clearly indicates that raising the level of abstraction is inevitably fundamental to deal with the intrinsic complexity that surround us. Accordingly, we envision that MHCs will play a crucial role in modelling complex, large-scale computing systems (or systems of systems) across various domains in the coming years. We invite the theoretical computing community to continue exploring and expanding the promising and intriguing frontiers of what we call Models of High-Level Computation.

[1]
John Longley and Dag Normann.
Higher-Order Computability.
Springer, Berlin; Heidelberg, 1st edition, 2015.

[2]
A. M. Turing.
On Computable Numbers, with an Application to the
Proceedings of the London Mathematical Society,
s2-42(1):230–265, 1937.

[3]
Michael Sipser.
Introduction to the Theory of Computation.
Cengage Learning, Boston, MA, 3rd edition, 2013.

[4]
Daniel Brand and Pitro Zafiropulo.
On Communicating Finite-State Machines.
Journal of the ACM, 30(2):323–342, 1983.

[5]
David Harel.
Statecharts: a visual formalism for complex systems.
Science of Computer Programming, 8(3):231–274, 1987.

[6]
Rajeev Alur, Sampath Kannan, and Mihalis Yannakakis.
Communicating Hierarchical State Machines.
In Jirí Wiedermann, Peter van Emde Boas, and Mogens Nielsen,
editors, Automata, Languages and Programming, Lecture Notes in
Computer Science, pages 169–178, Berlin, Heidelberg, 1999. Springer.

[7]
Franco Barbanera, Ugo de'Liguoro, and Rolf Hennicker.
Connecting open systems of communicating finite state machines.
Journal of Logical and Algebraic Methods in Programming,
109:1–34, 2019.

[8]
Benedikt Bollig.
Logic for communicating automata with parameterized topology.
In Proceedings of the Joint Meeting of the 23rd EACSL
Annual Conference on Computer Science Logic (CSL) and the 29th
Annual ACM/IEEE Symposium on Logic in Computer Science
(LICS), pages 1–10, New York, NY, USA, 2014. ACM.

[9]
Vanderson Hafemann Fragal, Adenilso Simao, and Mohammad Reza Mousavi.
Hierarchical featured state machines.
Science of Computer Programming, 171:67–88, 2019.

[10]
Salvatore La Torre, Margherita Napoli, Mimmo Parente, and Gennaro Parlato.
Verification of scope-dependent hierarchical state machines.
Information and Computation, 206(9):1161–1177, 2008.

[11]
Jack B. Dennis.
First version of a data flow procedure language.
In B. Robinet, editor, Proceedings of Symposium on
Programming, pages 362–376, Paris, France, 1974.

[12]
J.T. Buck and E.A. Lee.
Scheduling dynamic dataflow graphs with bounded memory using the
token flow model.
In IEEE International Conference on Acoustics, Speech,
and Signal Processing, volume 1, pages 429–432, 1993.

[13]
G. Kahn.
The Semantics of a Simple Language for Parallel
North-Holland Publishing Co., 1974.

[14]
Carl Hewitt.
Viewing control structures as patterns of passing messages.
Artificial Intelligence, 8(3):323–364, 1977.

[15]
Edward A. Lee and Eleftherios Matsikoudis.
The semantics of dataflow with firing.
In Gordon Plotkin, Gérard Huet, Jean-Jacques Lévy, and Yves Bertot,
editors, From Semantics to Computer Science: Essays in Honour
of Gilles Kahn, pages 71–94. Cambridge University Press, Cambridge,

[16]
E.A. Lee and D.G. Messerschmitt.
Synchronous data flow.
Proceedings of the IEEE, 75(9):1235–1245, 1987.

[17]
E.A. Lee and T.M. Parks.
Dataflow process networks.
Proceedings of the IEEE, 83(5):773–801, 1995.

[18]
Marc Geilen and Twan Basten.
Reactive process networks.
In Proceedings of the 4th ACM international conference on
Embedded software, pages 137–146, New York, NY, USA, 2004. ACM.

[19]
Stephen A. Edwards and Edward A. Lee.
The semantics and execution of a synchronous block-diagram language.
Science of Computer Programming, 48(1):21–42, 2003.

[20]
Farhad Arbab.
Reo: a channel-based coordination model for component composition.
Mathematical Structures in Computer Science, 14(3):329–366,

[21]
Kung-Kiu Lau, Mario Ornaghi, and Zheng Wang.
A Software Component Model and Its Preliminary
In Frank S. de Boer, Marcello M. Bonsangue, Susanne Graf, and
Willem-Paul de Roever, editors, Formal Methods for Components and
Objects, Lecture Notes in Computer Science, pages 1–21, Berlin,
Heidelberg, 2006. Springer.

[22]
Michele Colledanchise and Petter Ögren.
Behavior Trees in Robotics and AI: An Introduction.
CRC Press, Boca Raton London New York, 1st edition edition, 2018.

[23]
Nick Russell, W. M. P. van der Aalst, and Arthur H. M. ter Hofstede.
Workflow Patterns: The Definitive Guide.
The MIT Press, 2016.

[24]
Kung-Kiu Lau and Mario Ornaghi.
Control Encapsulation: A Calculus for Exogenous Composition
of Software Components.
In Grace A. Lewis, Iman Poernomo, and Christine Hofmeister, editors,
Component-Based Software Engineering, Lecture Notes in
Computer Science, pages 121–139, Berlin, Heidelberg, 2009. Springer.

[25]
Damian Arellanes, Kung-Kiu Lau, and Rizos Sakellariou.
Decentralized Data Flows for the Functional Scalability of
Service-Oriented IoT Systems.
The Computer Journal, 66(6):1477–1506, 2023.

[26]
Matteo Iovino, Edvards Scukins, Jonathan Styrud, Petter Ögren, and Christian
A survey of Behavior Trees in robotics and AI.
Robotics and Autonomous Systems, 154:104096, 2022.

[27]
W. M. P. Van der Aalst.
The application of Petri-nets to workflow management.
Journal of Circuits, Systems and Computers, 8(1):21–66, 1998.

[28]
Chun Ouyang, Eric Verbeek, Wil M. P. van der Aalst, Stephan Breutel, Marlon
Dumas, and Arthur H. M. ter Hofstede.
Formal semantics and analysis of control flow in WS-BPEL.
Science of Computer Programming, 67(2):162–198, 2007.

[29]
W. M. P. Van der Aalst and A. H. M. ter Hofstede.
YAWL: yet another workflow language.
Information Systems, 30(4):245–275, 2005.

[30]
J. C. M. Baeten, T. Basten, and M. A. Reniers.
Process Algebra: Equational Theories of Communicating
Cambridge Tracts in Theoretical Computer Science. Cambridge
University Press, Cambridge, 2009.

[31]
Donald Yau.
Operads of Wiring Diagrams.
Springer, New York, NY, 1st edition, 2018.

[32]
Robin Piedeleu and Fabio Zanasi.
An Introduction to String Diagrams for Computer Scientists,

[33]
Roger Penrose.
Applications of Negative Dimensional Tensors.
Combinatorial mathematics and its Applications, pages 221–244,

[34]
Filippo Bonchi, Paweł Sobociński, and Fabio Zanasi.
The Calculus of Signal Flow Diagrams I: Linear relations
on streams.
Information and Computation, 252:2–29, 2017.

[35]
Filippo Bonchi, Joshua Holland, Robin Piedeleu, Paweł Sobociński, and Fabio
Diagrammatic algebra: from linear to concurrent systems.
Proceedings of the ACM on Programming Languages, 3(POPL):1–28,

[36]
Bob Coecke and Ross Duncan.
Interacting quantum observables: categorical algebra and
New Journal of Physics, 13(4):1–85, 2011.

[37]
Kung-Kiu Lau and Simone Di Cola.
An Introduction to Component-based Software
World Scientific, Singapore, 1st edition, 2017.

[38]
Dina Q. Goldin, Scott A. Smolka, Paul C. Attie, and Elaine L. Sonderegger.
Turing machines, transition systems, and interaction.
Information and Computation, 194(2):101–128, 2004.