#  Variational Quantum Support Vector Machine based on $\Gamma$ matrix
expansion
and Variational Universal-Quantum-State Generator

Motohiko Ezawa  Department of Applied Physics, University of Tokyo, Hongo
7-3-1, 113-8656, Japan

###### Abstract

We analyze a binary classification problem by using a support vector machine
based on variational quantum-circuit model. We propose to solve a linear
equation of the support vector machine by using a $\Gamma$ matrix expansion.
In addition, it is shown that an arbitrary quantum state is prepared by
optimizing a universal quantum circuit representing an arbitrary $U(2^{N})$
based on the steepest descent method. It may be a quantum generalization of
Field-Programmable-Gate Array (FPGA).

Introduction

Quantum computation is a hottest topic in contemporary physicsFeynman ; DiVi ;
Nielsen . An efficient application of quantum computations is machine
learning, which is called quantum machine learningLloyd ; Schuld ; Biamonte ;
Wittek ; Harrow ; Wiebe ; Reben ; ZLi ; SchuldB ; Hav ; Lamata ; Cong . A
support vector machine is one of the most fundamental algorithms for machine
learningVap ; Noble ; Suy , which classifies data into two classes by a
hyperplane. The optimal hyperplane is determined by an associated linear
equation $F|\psi_{\text{in}}\rangle=|\psi_{\text{out}}\rangle$, where $F$ and
$|\psi_{\text{out}}\rangle$ are given. A quantum support vector machine solves
this linear equation by a quantum computerZhaokai ; Hav ; Reben . Usually, the
linear equation is solved by the Harrow-Hassidim-Lloyd (HHL) algorithmHHL .
However, this algorithm requires many quantum gates. Thus, the HHL algorithm
is hard to be executed by using a near-term quantum computer. Actually, this
algorithm has experimentally been verified only for two and three qubitsXCai ;
Barz ; JPan . In addition, it requires a unitary operator to execute
$e^{iFt}$, which is quite hard to be implemented.

The number of qubits in current quantum computers is restricted. Variational
quantum algorithms are appropriate for these small-qubit quantum computers,
which use both quantum computers and classical computers. Various methods have
been proposed such as Quantum Approximate Optimization Algorithm (QAOA)QAOA ,
variational eigenvalue solverPeru , quantum circuit learningMitarai  and
quantum linear solverPrie ; Xxu . We use wave functions with variational
parameters in QAOA, which are optimized by minimizing the expectation value of
the Hamiltonian. A quantum circuit has variational parameters in quantum
circuit learningMitarai , which are optimized by minimizing a certain cost
function. A quantum linear solver solves a linear equation by variational
ansatzPrie ; Xxu . The simplest method of the optimization is a steepest-
descent method.

In this paper, we present a variational method for a quantum support vector
machine by solving an associated linear equation based on variational quantum
circuit learning. We propose a method to expand the matrix $F$ by the $\Gamma$
matrices, which gives simple quantum circuits. We also propose a variational
method to construct an arbitrary state by using a universal quantum circuit to
represent an arbitrary unitary matrix $U(2^{N})$. We prepare various internal
parameters for a universal quantum circuit, which we optimize by minimizing a
certain cost function. Our circuit is capable to determine the unitary
transformation $U$ satisfying
$U|\psi_{\text{initial}}\rangle=|\psi_{\text{final}}\rangle$ with arbitrary
given states $|\psi_{\text{initial}}\rangle$ and
$|\psi_{\text{final}}\rangle$. It will be a quantum generalization of field-
programmable-gate array (FPGA), which may execute arbitrary outputs with
arbitrary inputs.

Risults

Support vector machine.

A support vector machine (SVM) is a computer algorithm that learns by examples
to assign labels to objects. It is a typical method to solve a binary-
classification problemVap . A simplest example reads as follows. Suppose that
there are red and blue points whose distributions are almost separated into
two dimensions. We classify these data points into two classes by a line, as
illustrated in Fig.1.

In general, $M$ data points are spattered in $D$ dimensions, which we denote
$\boldsymbol{x}_{j}$, where $1\leq j\leq M$. The problem is to determine a
hyperplane,

$\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}=0,$ (1)

separating data into two classes with the use of a support vector machine. We
set

$\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}>0$ (2)

for red points and

$\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}<0$ (3)

for blue points. These conditions are implemented by introducing a function

$f\left(\boldsymbol{x}\right)=\text{sgn}\left(\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}\right),$
(4)

which assigns $f\left(\boldsymbol{x}\right)=1$ to red points and
$f\left(\boldsymbol{x}\right)=-1$ to blue points. In order to determine
$\omega_{0}$ and $\boldsymbol{\omega}$ for a given set of data
$\boldsymbol{x}_{j}$, we introduce real numbers $\alpha_{j}$ by

$\boldsymbol{\omega}=\sum_{j=1}^{M}\alpha_{j}\boldsymbol{x}_{j}.$ (5)

A support vector machine enables us to determine $\omega_{0}$ and $\alpha_{j}$
by solving the linear equation

$F\left(\begin{array}[]{c}\omega_{0}\\\ \alpha_{1}\\\ \vdots\\\
\alpha_{M}\end{array}\right)=\left(\begin{array}[]{c}0\\\ y_{1}\\\ \vdots\\\
y_{M}\end{array}\right),$ (6)

where $y_{i}=f(x_{i})=\pm 1$, and $F$ is a $(M+1)\times(M+1)$ matrix given by

$F=\left(\begin{array}[]{cccc}0&1&\cdots&1\\\ 1&&&\\\
\vdots&&K+I_{M}/\gamma&\\\ 1&&&\end{array}\right).$ (7)

Here,

$K_{ij}=\boldsymbol{x}_{i}\cdot\boldsymbol{x}_{j},$ (8)

is a Kernel matrix, and $\gamma$ is a certain fixed constant which assures the
existence of the solution of the linear equation (6) even when the red and
blue points are slightly inseparable. Note that $\gamma\rightarrow\infty$
corresponds to the hard margin condition. Details of the derivation of Eq.(6)
are given in Method A.

Figure 1: (a) Binary classification of red and blue points based on a quantum
support vector machine with soft margin. A magenta (cyan) line obtained by an
exact solution (variational method). (b) Evolution of the cost function. The
vertical axis is the log${}_{10}E_{\text{cost}}$. The horizontal axis is the
variational step number. We have used $r=2$, $\xi_{1}=0.001$ and
$\xi_{2}=0.0005$ and $\gamma=1$. We have runed simulations ten times.

Quantum linear solver based on $\Gamma$ matrix expansion.

We solve the linear equation (6) by a quantum computer. In general, we solve a
linear equation

$F\left|\psi_{\text{in}}\right\rangle=c\left|\psi_{\text{out}}\right\rangle,$
(9)

for an arbitrary given non-unitary matrix $F$ and an arbitrary given state
$\left|\psi_{\text{out}}\right\rangle$. Here, the coefficient $c$ is
introduced to preserve the norm of the state, and it is given by

$c=\sqrt{\left\langle\psi_{\text{in}}\right|F^{\dagger}F\left|\psi_{\text{in}}\right\rangle}.$
(10)

The HHL algorithmHHL  is a most famous algorithm to solve this linear equation
by a quantum computer. We first construct a Hermitian matrix by

$H=\left(\begin{array}[]{cc}0&F\\\ F^{\dagger}&0\end{array}\right).$ (11)

Then, a unitary matrix associated with $F$ is uniquely obtained by $e^{iHt}$.
Nevertheless, it requires many quantum gates. In addition, it is a nontrivial
problem to implement $e^{iHt}$.

Recently, variational methods have been proposedPrie  to solve the linear
equation (9). In one of the methods, the matrix $F$ is expanded in terms of
some unitary matrices $U_{j}$ as

$F=\sum_{j=0}^{2^{N}-1}c_{j}U_{j}.$ (12)

In general, a complicated quantum circuit is necessary to determine the
coefficient $c_{j}$.

We start with a trial state $|\tilde{\psi}_{\text{in}}\rangle$ to determine
the state $|\psi_{\text{in}}\rangle$. Application of each unitary matrix to
this state is efficiently done by a quantum computer,
$U_{j}|\tilde{\psi}_{\text{in}}\rangle=|\tilde{\psi}_{\text{out}}^{\left(j\right)}\rangle$,
and we obtain

$F|\tilde{\psi}_{\text{in}}\rangle=\sum_{j=0}^{2^{N}-1}c_{j}U_{j}|\tilde{\psi}_{\text{in}}\rangle=\sum_{j=0}^{2^{N}-1}c_{j}|\tilde{\psi}_{\text{out}}^{\left(j\right)}\rangle\equiv
c|\tilde{\psi}_{\text{out}}\rangle,$ (13)

where $|\tilde{\psi}_{\text{out}}\rangle$ is an approximation of the given
state $\left|\psi_{\text{out}}\right\rangle$. We tune a trial state
$|\tilde{\psi}_{\text{in}}\rangle$ by a variational method so as to minimize
the cost functionPrie

$E_{\text{cost}}\equiv
1-\left|\langle\tilde{\psi}_{\text{out}}|\psi_{\text{out}}\rangle\right|^{2},$
(14)

which measures the similarity between the approximate state
$|\tilde{\psi}_{\text{out}}\rangle$ and the state
$\left|\psi_{\text{out}}\right\rangle$ in (9). We have $0\leq
E_{\text{cost}}\leq 1$, where $E_{\text{cost}}=0$ for the exact solution. The
merit of this cost function is that the inner product is naturally calculated
by a quantum computer.

Let the dimension of the matrix $F$ be $2^{N}$. It is enough to use $N$
satisfying $2^{N-1}<D\leq 2^{N}$ without loss of generality by adding trivial
$2^{N}-D$ components to the linear equation. We propose to expand the matrix
$F$ by the gamma matrices $\Gamma_{j}$ as

$F=\sum_{j=0}^{2^{N}-1}c_{j}\Gamma_{j},$ (15)

with

$\Gamma_{j}=\bigotimes_{\beta=1}^{N}\sigma_{\alpha}^{\left(\beta\right)},$
(16)

where $\alpha$ $=0,x,y$ and $z$.

The merit of our method is that it is straightforward to determine $c_{j}$ by
the well-known formula

$c_{j}=\text{Tr}\left[\Gamma_{j}F\right].$ (17)

In order to construct a quantum circuit to calculate $c_{j}$, we express the
matrix $F$ by column vectors as

$F=\left\\{\left|f_{0}\right\rangle,\cdots,\left|f_{2^{N}-1}\right\rangle\right\\}.$
(18)

We have $\left(\left|f_{q-1}\right\rangle\right)_{p}=F_{pq}$, where subscript
$p$ denotes the $p$-th component of $\left|f_{q-1}\right\rangle$. Then $c_{j}$
is given by

$c_{j}=\sum_{q=0}^{2^{N}-1}\left(\Gamma_{j}\left|f_{q}\right\rangle\right)_{q}=\sum_{q=0}^{2^{N}-1}\left\langle\\!\left\langle
q\right|\right.\Gamma_{j}\left|f_{q}\right\rangle,$ (19)

where the subscript $q$ denotes the ($q+1$)-th component of
$\Gamma_{j}\left|f_{q}\right\rangle$. We have introduced a notation
$\left|q\right\rangle\\!\rangle\equiv|n_{1}n_{2}\cdots n_{N}\rangle$ with
$n_{i}=0,1$, where $q$ is the decimal representation of the binary number
$n_{1}n_{2}\cdots n_{N}$. See explicit examples for one and two qubits in
Method B.

The state $\left.\left|q\right\rangle\\!\right\rangle\equiv|n_{1}n_{2}\cdots
n_{N}\rangle$ is generated as follows. We prepare the NOT gates
$\sigma_{x}^{\left(i\right)}$ for the $i$-th qubit if $n_{i}=1$. Using all
these NOT gates we define

$U_{X}^{\left(q\right)}=\bigotimes\limits_{n_{i}=1}\sigma_{x}^{\left(i\right)}.$
(20)

We act it on the initial state $\left|0\right\rangle\\!\rangle$ and obtain

$U_{X}^{\left(q\right)}\left|0\right\rangle\\!\rangle=\left|q\right\rangle\\!\rangle.$
(21)

Next, we construct a unitary gate $U_{f_{q}}$ generating
$\left|f_{q}\right\rangle$,

$U_{f_{q}}\left|0\right\rangle\\!\rangle=\left|f_{q}\right\rangle.$ (22)

We will discuss how to prepare $U_{f_{q}}$ by a quantum circuit soon later;
See Eq.(28). By using these operators, $c_{j}$ is expressed as

$c_{j}=\sum_{q=0}^{2^{N}-1}\left\langle\\!\left\langle
0\right|\right.U_{X}^{\left(q\right)}\Gamma_{j}U_{f_{q}}\left|0\right\rangle\\!\rangle,$
(23)

which can be executed by a quantum computer. We show explicit examples in
Fig.2.

Once we have $c_{j}$, the final state is obtained by applying $\Gamma_{j}$ to
$|\tilde{\psi}_{\text{in}}\rangle$ and taking sum over $j$, which leads to

$|\tilde{\psi}_{\text{out}}\rangle=F|\tilde{\psi}_{\text{in}}\rangle=\sum_{j=0}^{2^{N}-1}c_{j}\Gamma_{j}|\tilde{\psi}_{\text{in}}\rangle.$
(24)

The implementation of the $\Gamma$ matrix is straightforward in quantum
circuit, because the $\Gamma$ matrix is composed of the Pauli sigma matrices,
as shown in Fig.2.

Figure 2: Quantum circuits determining $c_{j}$. We show an example with (a)
$\Gamma_{yx0}=\sigma_{y}\otimes\sigma_{x}\otimes\sigma_{0}$.
$U_{X}^{\left(6\right)}\left.\left|0\right\rangle\\!\right\rangle=\sigma_{x}^{\left(1\right)}\sigma_{x}^{\left(2\right)}\left|000\right\rangle=\left|110\right\rangle=\left.\left|6\right\rangle\\!\right\rangle$
and (b) $\Gamma_{xyz}=\sigma_{x}\otimes\sigma_{y}\otimes\sigma_{z}$.
$U_{X}^{\left(5\right)}\left|0\right\rangle\\!\rangle=\sigma_{x}^{\left(1\right)}\sigma_{x}^{\left(3\right)}\left|000\right\rangle=\left|101\right\rangle=\left.\left|5\right\rangle\\!\right\rangle$.

We may use the steepest descent method to find an optimal trial state
$|\tilde{\psi}_{\text{in}}\rangle$ closest to the state
$|\psi_{\text{in}}\rangle$. We calculate the difference of the cost function
$\Delta E_{\text{cost}}$ when we slightly change the trial state
$|\tilde{\psi}_{\text{in}}(t)\rangle$ at step $t$ by the amount of
$\Delta|\tilde{\psi}_{\text{in}}(t)\rangle$ as

$\Delta E_{\text{cost}}\equiv
E_{\text{cost}}\left(|\tilde{\psi}_{\text{in}}(t)\rangle+\Delta|\tilde{\psi}_{\text{in}}(t)\rangle\right)-E_{\text{cost}}\left(|\tilde{\psi}_{\text{in}}(t)\rangle\right)\simeq\frac{\Delta
E_{\text{cost}}}{\Delta|\tilde{\psi}_{\text{in}}(t)\rangle}\Delta|\tilde{\psi}_{\text{in}}(t)\rangle.$
(25)

We explain how to construct $|\tilde{\psi}_{\text{in}}(t)\rangle$ by a quantum
circuit soon later; See Eq.(28). Then, we renew the state as

$|\tilde{\psi}_{\text{in}}(t)\rangle\rightarrow|\tilde{\psi}_{\text{in}}(t)\rangle-\eta_{t}\frac{\Delta
E_{\text{cost}}}{\Delta|\tilde{\psi}_{\text{in}}(t)\rangle}\Delta|\tilde{\psi}_{\text{in}}(t)\rangle,$
(26)

where we use an exponential function for $\eta_{t}$,

$\eta_{t}=\xi_{1}e^{-\xi_{2}t}.$ (27)

We choose appropriate constants $\xi_{1}$ and $\xi_{2}$ for an efficient
search of the optimal solution, whose explicit examples are given in the
caption of Fig.2. We stop the renewal of the variational step when the
difference $\Delta|\tilde{\psi}_{\text{in}}(t)\rangle$ becomes sufficiently
small, which gives the optimal state of the linear equation (9).

Variational universal-quantum-state generator.

In order to construct the trial state $|\tilde{\psi}_{\text{in}}(t)\rangle$,
it is necessary to prepare an arbitrary state $\left|\psi\right\rangle$ by a
quantum circuit. Alternatively, we need such a unitary transformation $U$ that

$U\left|0\right\rangle\\!\rangle=\left|\psi\right\rangle.$ (28)

It is known that any unitary transformation is done by a sequential
application of the Hadamard, the $\pi/4$ phase-shift and the CNOT gatesDeutsch
; Dawson . Indeed, an arbitrary unitary matrix is decomposable into a
sequential application of quantum gatesDeutsch ; Dawson , each of which is
constructed as a universal quantum circuit systematicallyKraus ; Vidal ; Motto
; Shende ; Vatan ; Sousa . Universal quantum circuits have so far been
demonstrated experimentally for two and three qubitsHanne ; DiCarlo ; Qiang ;
Roy .

We may use a variational method to construct $U$ satisfying Eq.(28). Quantum
circuit learning is a variational methodMitarai , where angle variables
$\theta_{i}$ are used as variational parameters in a quantum circuit $U$, and
the cost function is optimized by tuning $\theta_{i}$. We propose to use a
quantum circuit learning for a universal quantum circuit. We show that an
arbitrary state $\left|\psi\left(\theta_{i}\right)\right\rangle$ can be
generated by tuning $U\left(\theta_{i}\right)$ starting from the initial state
$\left|0\right\rangle\\!\rangle$ as

$U\left(\theta_{i}\right)\left|0\right\rangle\\!\rangle=\left|\psi\left(\theta_{i}\right)\right\rangle.$
(29)

We adjust $\theta_{i}$ by minimizing the cost function

$E_{\text{cost}}\left(\theta_{i}\right)\equiv
1-\left|\left\langle\psi\left(\theta_{i}\right)\left|\psi\right\rangle\right.\right|^{2},$
(30)

which is the same as that of the variational quantum support vector machine.
We present explicit examples of universal quantum circuits for one, two and
three qubits in Method C.

Figure 3: Evolution of the cost function for (a) two qubits and (b) three
qubits. The vertical axis is the log${}_{10}E_{\text{cost}}$. The horizontal
axis is the number of variational steps. We use $c_{1}=0.005$ and
$c_{2}=0.005$ for both the two- and three-qubit universal quantum circuits. We
prepare random initial and final states, where we have runed simulations ten
times.

Quantum Field-Programmable-Gate Array.

We next consider a problem to find a unitary transformation $U_{\text{ini-
fin}}$ which maps an arbitrary initial state
$\left|\psi_{\text{initial}}\right\rangle$ to an arbitrary final state
$\left|\psi_{\text{final}}\right\rangle$,

$U_{\text{ini-
fin}}\left|\psi_{\text{initial}}\right\rangle=\left|\psi_{\text{final}}\right\rangle.$
(31)

Since we can generate an arbitrary unitary matrix as in Eq.(28), it is
possible to generate such matrices $U_{\text{ini}}$ and $U_{\text{fin}}$ that

$U_{\text{ini}}\left|0\right\rangle\\!\rangle=\left|\psi_{\text{initial}}\right\rangle,\qquad
U_{\text{fin}}\left|0\right\rangle\\!\rangle=\left|\psi_{\text{final}}\right\rangle.$
(32)

Then, Eq.(31) is solved as

$U_{\text{fin}}=U_{\text{ini-fin}}U_{\text{ini}},$ (33)

since $U_{\text{ini-
fin}}\left|\psi_{\text{initial}}\right\rangle=U_{\text{ini-
fin}}U_{\text{ini}}\left|0\right\rangle\\!\rangle=\left|\psi_{\text{final}}\right\rangle=U_{\text{fin}}\left|0\right\rangle\\!\rangle$.

An FPGA is a classical integrated circuit, which can be programmable by a
customer or a designer after manufacturing in a factory. An FPGA executes any
classical algorithms. On the other hand, our variational universal quantum-
state generator creates an arbitrary quantum state. We program by using the
variational parameters $\theta_{i}$. In this sense, the above quantum circuit
may be considered as a quantum generalization of FPGA, which is a quantum FPGA
(q-FPGA).

We show explicitly how the cost function is renewed for each variational step
in the case of two- and three-qubit universal quantum circuits in Fig.3, where
we have generated the initial and the final states randomly. We optimize 15
parameters $\theta_{i}$ for two-qubit universal quantum circuits and 82
parameters $\theta_{i}$ for three-qubit universal quantum circuits. We find
that $U_{\text{ini-fin}}$ is well determined by variational method as in
Fig.3.

Variational quantum support vector machine.

We demonstrate a binary classification problem in two dimensions based on the
support vector machine. We prepare a data set, where red points have a
distribution around $\left(r\cos\Theta,r\sin\Theta\right)$ with variance $r$,
while blue points have a distribution around
$\left(-r\cos\Theta,-r\sin\Theta\right)$ with variance $r$. We assume the
Gaussian normal distribution. We choose $\Theta$ randomly. We note that there
are some overlaps between the red and blue points, which is the soft margin
model.

As an example, we show the distribution of red and blue points and the lines
obtained by the variational method marked in cyan and by the direct solution
of (6) marked in magenta in Fig.1. They agrees well with one another, where
both of the lines well separate red and blue points. We have prepared 31 red
points and 32 blue points, and used six qubits.

Discussion

We have proposed that the matrix $F$ is efficiently inputted into a quantum
computer by using the $\Gamma$-matrix expansion method. There are many ways to
use a matrix in a quantum computer such as linear regression and principal
component analysis. Our method will be applicable to these cases.

Although it is possible to obtain the exact solution for the linear equation
by the HHL algorithm, it requires many gates. On the other hand, it is often
hard to obtain the exact solution by variational methods since trial functions
may be trapped to a local minimum. However, this problem is not serious for
the machine learning problem because it is more important to obtain an
approximate solution efficiently rather than an exact solution by using many
gates. Indeed, our optimized hyperplane also well separates red and blue
points as shown in Fig.1(a).

In order to classify $M$ data, we need to prepare $\log_{2}M$ qubits. It is
hard to execute a large number of data points by current quantum computers.
Recently, it is shown that electric circuits may simulate universal quantum
gatesEzawaUniv ; EzawaDirac ; LCBit  based on the fact that the Kirchhoff law
is rewritten in the form of the Schrödinger equationEzawaSch . Our variational
algorithm will be simulated by using them.

Methods

A: Support vector machine.

A support vector machine is an algorithm for supervised learningVap ; Noble ;
Suy . We first prepare a set of training data, where each point is marked
either in red or blue. Then, we determine a hyperplane separating red and blue
points. After learning, input data are classified into red or blue by
comparing the input data with the hyperplane. The support vector machine
maximizes a margin, which is a distance between the hyperplane and data
points. If red and blue points are perfectly separated by the hyperplane, it
is called a hard margin problem [Fig.4(a)]. Otherwise, it is called a soft
margin problem [Fig.4(b)].

We minimize the distance $d_{j}$ between a data point $\boldsymbol{x}_{j}$ and
the hyperplane given by

$d_{j}=\frac{\left|\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right|}{\left|\boldsymbol{\omega}\right|}.$
(34)

We define support vectors $\boldsymbol{x}$ as the closest points to the
hyperplane. There is such a vector in each side of the hyperplane, as shown in
Fig.4(a). This is the origin of the name of the support vector machine.
Without loss of generality, we set

$\left|\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}\right|=1$ (35)

for the support vectors, because the hyperplane is present at the equidistance
of two closest data points and because it is possible to set the magnitude of
$\left|\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}\right|$ to be $1$ by
scaling $\boldsymbol{\omega}$ and $\omega_{0}$. Then, we maximize the distance

$d=\frac{\left|\boldsymbol{\omega}\cdot\boldsymbol{x}+\omega_{0}\right|}{\left|\boldsymbol{\omega}\right|}=\frac{1}{\left|\boldsymbol{\omega}\right|},$
(36)

which is identical to minimize $\left|\boldsymbol{\omega}\right|$.

Figure 4: Illustration of the hyperplane and the support vector. Two support
vectors are marked by red and blue squares. (a) Hard margin where red and blue
points are separated perfectly, and (b) soft margin where they are separated
imperfectly.

First, we consider the hard margin problem, where red and blue points are
perfectly separable. All red points satisfy
$\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}>1$ and all blue points
satisfy $\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}<-1$. We
introduce variables $y_{j}$, where $y_{j}=1$ for red points and $y_{j}=-1$ for
blue points. Using them, the condition is rewritten as

$\left(\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)y_{j}\geq
1$ (37)

for each $j$. The problem is reduced to find the minimum of
$\left|\boldsymbol{\omega}\right|^{2}$ under the above inequalities. The
optimization under inequality conditions is done by the Lagrange multiplier
method with the Karush-Kuhn-Tucker conditionKKT . It is expressed in terms of
the Lagrangian as

$L\left(\boldsymbol{\omega},\omega_{0},\boldsymbol{\alpha}\right)=\frac{1}{2}\left|\boldsymbol{\omega}\right|^{2}-\sum_{j}\beta_{j}[\left(\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)y_{j}-1],$
(38)

where $\beta_{j}$ are Lagrange multipliers to ensure the constraints.

For the soft margin case, we cannot separate two classes exactly. In order to
treat this case, we introduce slack variables $\xi_{j}$ satisfying

$\left(\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)y_{j}\geq
1-\xi_{j},\qquad\xi_{j}\geq 0$ (39)

and redefine the cost function as

$E_{\text{cost}}=\frac{1}{2}\left|\boldsymbol{\omega}\right|^{2}+\gamma\sum_{j=1}^{M}\xi_{j}^{2}.$
(40)

Here, $\gamma=\infty$ corresponds to the hard margin. The second term
represents the penalty for some of data points to have crossed over the
hyperplane. The Lagrangian is modified as

$L\left(\boldsymbol{\omega},\omega_{0},\xi_{i},\boldsymbol{\beta}\right)=\frac{1}{2}\left|\boldsymbol{\omega}\right|^{2}+\gamma\sum_{j=1}^{M}\xi_{j}^{2}-\sum_{j=1}^{M}\left[\left(\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)\beta_{j}y_{j}-\left(1-\xi_{i}\right)\right].$
(41)

The stationary points are determined by

$\displaystyle\frac{\partial L}{\partial\boldsymbol{\omega}}=$
$\displaystyle\boldsymbol{\omega}-\sum_{j=1}^{M}\beta_{j}y_{j}\boldsymbol{x}_{j}=0,$
(42) $\displaystyle\frac{\partial L}{\partial\omega_{0}}=$
$\displaystyle-\sum_{j=1}^{M}\beta_{j}y_{j}=0,$ (43)
$\displaystyle\frac{\partial L}{\partial\xi_{j}}=$
$\displaystyle\gamma\xi_{j}-\beta_{j}=0,$ (44) $\displaystyle\frac{\partial
L}{\partial\beta_{j}}=$
$\displaystyle\left(\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)y_{j}-\left(1-\xi_{i}\right)=0.$
(45)

We may solve these equations to determine $\boldsymbol{\omega}$ and $\nu_{j}$
as

$\boldsymbol{\omega}=\sum_{j=1}^{M}\beta_{j}y_{j}\boldsymbol{x}_{j},$ (46)

from (42), and

$\xi_{j}=\beta_{j}/\gamma$ (47)

from (44). Inserting them into (45), we find

$y_{j}\sum_{i=1}^{M}\left(\beta_{i}y_{i}\boldsymbol{x}_{i}\cdot\boldsymbol{x}_{j}+\omega_{0}\right)-\left(1-\beta_{j}/\gamma\right)=0.$
(48)

Since $y_{j}^{2}=1$, it is rewritten as

$\omega_{0}+\sum_{i=1}^{M}\left(\boldsymbol{x}_{i}\cdot\boldsymbol{x}_{j}+\delta_{ij}/\gamma\right)\beta_{i}y_{i}=y_{j}.$
(49)

Since $\beta_{j}$ appears always in a pair with $y_{j}$, we introduce a new
variable defined by

$\alpha_{j}=\beta_{j}y_{j},$ (50)

and we define the Kernel matrix $K_{ij}$ as

$K_{ij}=\boldsymbol{x}_{i}\cdot\boldsymbol{x}_{j}.$ (51)

Then, $\omega_{0}$ and $\alpha_{j}$ are obtained by solving linear equations

$\displaystyle\sum_{i=1}^{M}\alpha_{j}=$ $\displaystyle 0,$ (52)
$\displaystyle\omega_{0}+\sum_{i=1}^{M}\left(\boldsymbol{x}_{i}\cdot\boldsymbol{x}_{j}+\delta_{ij}/\gamma\right)\alpha_{i}=$
$\displaystyle y_{j},$ (53)

which are summarized as

$\left(\begin{array}[]{cccc}0&1&\cdots&1\\\ 1&&&\\\ \vdots&&K+I_{M}/\gamma&\\\
1&&&\end{array}\right)\left(\begin{array}[]{c}\omega_{0}\\\ \alpha_{1}\\\
\vdots\\\ \alpha_{M}\end{array}\right)=\left(\begin{array}[]{c}0\\\ y_{1}\\\
\vdots\\\ y_{M}\end{array}\right),$ (54)

which is Eq.(6) in the main text. Finally, $\boldsymbol{\omega}$ is determined
by

$\boldsymbol{\omega}=\sum_{j=1}^{M}\alpha_{j}\boldsymbol{x}_{j}.$ (55)

Once the hyperplane is determined, we can classify new input data into red if

$\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}>0$ (56)

and blue if

$\boldsymbol{\omega}\cdot\boldsymbol{x}_{j}+\omega_{0}<0.$ (57)

Thus, we obtain the hyperplane for binary classification.

B: $\Gamma$ matrix expansion.

We explicitly show how to calculate $c_{j}$ in (17) based on the $\Gamma$
matrix expansion for the one and two qubits.

One qubit:

We show an explicit example of the $\Gamma$-matrix expansion for one qubit.
Ome qubit is represented by a $2\times 2$ matrix,

$F=\left(\begin{array}[]{cc}F_{11}&F_{12}\\\ F_{21}&F_{22}\end{array}\right).$
(58)

The column vectors are explicitly given by

$\displaystyle\left|f_{1}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{11}\\\
F_{21}\end{array}\right)=F_{11}\left|0\right\rangle+F_{21}\left|1\right\rangle,\qquad$
(61) $\displaystyle\left|f_{2}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{12}\\\
F_{22}\end{array}\right)=F_{12}\left|0\right\rangle+F_{22}\left|1\right\rangle.$
(64)

The coefficient $c_{j}$ in (17) is calculated as

$c_{j}=\text{Tr}\left[\sigma_{j}F\right]=\langle
0|\sigma_{j}\left|f_{1}\right\rangle+\left\langle
1\right|\sigma_{j}\left|f_{2}\right\rangle=\sum_{p=0,1}\langle
p|\sigma_{j}\left|f_{p}\right\rangle=\sum_{p=0,1}\langle
0|U_{X}^{\left(p\right)}\sigma_{j}U_{f_{p}}\left|0\right\rangle.$ (65)

Two qubits:

Next, we show an explicit example of the $\Gamma$-matrix expansion for two
qubits. Two qubits are represented by a $4\times 4$ matrix,

$F=\left(\begin{array}[]{cccc}F_{11}&F_{12}&F_{13}&F_{14}\\\
F_{21}&F_{22}&F_{23}&F_{24}\\\ F_{31}&F_{32}&F_{33}&F_{34}\\\
F_{41}&F_{42}&F_{43}&F_{44}\end{array}\right).$ (66)

The column vectors are explicitly given by

$\displaystyle\left|f_{1}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{11}\\\ F_{21}\\\ F_{31}\\\
F_{41}\end{array}\right)=F_{11}\left|00\right\rangle+F_{21}\left|01\right\rangle+F_{31}\left|10\right\rangle+F_{41}\left|11\right\rangle,$
(71) $\displaystyle\left|f_{2}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{12}\\\ F_{22}\\\ F_{32}\\\
F_{42}\end{array}\right)=F_{12}\left|00\right\rangle+F_{22}\left|01\right\rangle+F_{32}\left|10\right\rangle+F_{42}\left|11\right\rangle,$
(76) $\displaystyle\left|f_{3}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{13}\\\ F_{23}\\\ F_{33}\\\
F_{43}\end{array}\right)=F_{13}\left|00\right\rangle+F_{23}\left|01\right\rangle+F_{33}\left|10\right\rangle+F_{43}\left|11\right\rangle,$
(81) $\displaystyle\left|f_{4}\right\rangle=$
$\displaystyle\left(\begin{array}[]{c}F_{14}\\\ F_{24}\\\ F_{34}\\\
F_{44}\end{array}\right)=F_{14}\left|00\right\rangle+F_{24}\left|01\right\rangle+F_{34}\left|10\right\rangle+F_{44}\left|11\right\rangle.$
(86)

The coefficient $c_{j}$ in (17) is calculated as

$c_{j}=\text{Tr}\left[\Gamma_{j}F\right]=\left\langle
00\right|\Gamma_{j}\left|f_{1}\right\rangle+\left\langle
01\right|\Gamma_{j}\left|f_{2}\right\rangle+\left\langle
10\right|\Gamma_{j}\left|f_{3}\right\rangle+\left\langle
11\right|\Gamma_{j}\left|f_{4}\right\rangle=\sum_{p=0}^{3}\langle\\!\left\langle
p\right|\Gamma_{j}\left|f_{p}\right\rangle=\sum_{p=0}^{3}\langle\\!\left\langle
0\right|U_{X}^{\left(p\right)}\Gamma_{j}U_{f_{p}}\left|0\right\rangle\\!\rangle.$
(87)

C: Universal quantum circuits.

Angle variables are used as variational parameters in a universal quantum
circuit learning. We present examples for one, two and three qubits.

One-qubit universal quantum circuit:

The single-qubit rotation gates are defined by

$\displaystyle R\left(\theta,\phi\right)=$
$\displaystyle\exp\left[-i\theta\left(\sigma_{x}\cos\phi+\sigma_{y}\sin\phi\right)/2\right],$
(88) $\displaystyle R_{z}\left(\phi_{z}\right)=$
$\displaystyle\exp\left[-i\sigma_{z}\phi_{z}/2\right].$ (89)

The one-qubit universal quantum circuit is constructed as

$U^{\left(1\right)}\left(\theta,\phi,\phi_{z}\right)=R\left(\theta,\phi\right)R_{z}\left(\phi_{z}\right)=\left(\begin{array}[]{cc}e^{-i\phi_{z}/2}\cos\frac{\theta}{2}&-ie^{i\left(\phi_{z}/2-\phi\right)}\sin\frac{\theta}{2}\\\
-ie^{-i\left(\phi_{z}/2-\phi\right)}\sin\frac{\theta}{2}&e^{i\phi_{z}/2}\cos\frac{\theta}{2}\end{array}\right).$
(90)

We show a quantum circuit in Fig.5(a). There are three variational parameters.

It is obvious that an arbitrary state is realized starting from the state
$\left|0\right\rangle$ as

$U\left(1\right)\left(\begin{array}[]{c}1\\\
0\end{array}\right)=\left(\begin{array}[]{c}e^{-i\phi_{z}/2}\cos\frac{\theta}{2}\\\
-ie^{-i\left(\phi_{z}/2-\phi\right)}\sin\frac{\theta}{2}\end{array}\right).$
(91)

Figure 5: Universal quantum circuits for (a) one, (b) two and (c) three
qubits.

Two-qubit universal quantum circuit:

The two-qubit universal quantum circuit is constructed asHanne

$\displaystyle U\left(2\right)\equiv$
$\displaystyle\left[U^{\left(1\right)}\left(\theta_{A},\phi_{A},\phi_{z,A}\right)\otimes
U^{\left(1\right)}\left(\theta_{B},\phi_{B},\phi_{z,B}\right)\right]U_{G}\left[R\left(\theta_{E},0\right)\otimes
R\left(\frac{3\pi}{2},0\right)\right]U_{G}\left[R\left(\theta_{F},\frac{\pi}{2}\right)\otimes
R\left(\frac{3\pi}{2},\theta_{G}\right)\right]$ $\displaystyle
U_{G}\left[U^{\left(1\right)}\left(\theta_{C},\phi_{C},\phi_{z,C}\right)\otimes
U^{\left(1\right)}\left(\theta_{D},\phi_{D},\phi_{z,D}\right)\right],$ (92)

where the entangling two-qubit gate is defined byHanne

$U_{G}=e^{-i\pi/4}\exp\left[\frac{i\pi}{4}\sigma_{z}\otimes\sigma_{z}\right].$
(93)

The two-qubits universal quantum circuit contains 15 variational parameters.
We show a quantum circuit in Fig.5(b).

Three-qubit universal quantum circuit:

The three-qubit universal quantum circuit is constructed as

$U\left(3\right)\equiv\left[U_{A}^{\left(2\right)}\otimes
U_{A}^{\left(1\right)}\right]U_{A}\left(3\right)\left[U_{B}^{\left(2\right)}\otimes
U_{B}^{\left(1\right)}\right]U_{C}\left(3\right)\left[U_{C}^{\left(2\right)}\otimes
U_{C}^{\left(1\right)}\right]U_{B}\left(3\right)\left[U_{D}^{\left(2\right)}\otimes
U_{D}^{\left(1\right)}\right],$ (94)

where $U_{A}^{\left(1\right)}$, $U_{B}^{\left(1\right)}$,
$U_{C}^{\left(1\right)}$, and $U_{D}^{\left(1\right)}$ are one-qubit universal
quantum circuits, while $U_{A}^{\left(2\right)}$, $U_{B}^{\left(2\right)}$,
$U_{C}^{\left(2\right)}$, and $U_{D}^{\left(2\right)}$ are two-qubit universal
quantum circuit and

$\displaystyle U_{A}\left(3\right)=$
$\displaystyle\exp\left[i\left(\theta_{xxz}^{A}\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{z}+\theta_{yyz}^{A}\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{z}+\theta_{zzz}^{A}\sigma_{z}\otimes\sigma_{z}\otimes\sigma_{z}\right)\right],$
(95) $\displaystyle U_{B}\left(3\right)=$
$\displaystyle\exp\left[i\left(\theta_{xxz}^{B}\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{z}+\theta_{yyz}^{B}\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{z}+\theta_{zzz}^{B}\sigma_{z}\otimes\sigma_{z}\otimes\sigma_{z}\right)\right],$
(96) $\displaystyle U_{C}\left(3\right)=$
$\displaystyle\exp\left[i\left(\theta_{xxx}^{C}\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{x}+\theta_{yyx}^{C}\sigma_{y}\otimes\sigma_{y}\otimes\sigma_{x}+\theta_{zzx}^{C}\sigma_{z}\otimes\sigma_{z}\otimes\sigma_{x}+\theta_{00x}^{C}\sigma_{0}\otimes\sigma_{0}\otimes\sigma_{x}\right)\right].$
(97)

Eplicit quantum circuits for $U_{A}\left(3\right)$, $U_{B}\left(3\right)$ and
$U_{C}\left(3\right)$ are shown in Ref.Sousa . The three-qubits universal
quantum circuit contains 82 variational parameters. We show a quantum circuit
in Fig.5(c).

Multi-qubit universal quantum circuit:

General multi-qubit universal quantum circuit is constructed in Ref.Motto .
The minimum numbers of variational parameters are $4^{N}-1$ for $N$-qubit
unicersal quantum circuits. However, we need more variational parameters in
the currently known algorithm for $N\geq 3$.

## References

  * (1) R. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
  * (2) D. P. DiVincenzo, Science 270, 255 (1995).
  * (3) M. Nielsen and I. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2016, p. 189; ISBN 978-1-107-00217-3.
  * (4) S. Lloyd, M. Mohseni and P. Rebentrost, arXiv:1307.0411
  * (5) M. Schuld, I. Sinayskiy, F. Petruccione, Contemporary Physics. 56, 172 (2014).
  * (6) P. Wittek, Quantum Machine Learning: What Quantum Computing Means to Data Mining Academic Press (2014)
  * (7) J. Biamonte, Nature. 549, 195 (2017).
  * (8) A. W. Harrow, A. Hassidim and S. Lloyd, Phys. Rev. Lett. 103 150502 (2009)
  * (9) N. Wiebe, D. Braun and S. Lloyd, Phys. Rev. Lett. 109, 050505 (2012)
  * (10) P. Rebentrost, M. Mohseni, S. Lloyd, Phys. Rev. Lett. 113, 130503 (2014)
  * (11) Z. Li, X. Liu, N. Xu and J. Du, Phys. Rev. Lett. 114, 140504 (2015).
  * (12) M. Schuld and N. Killoran, Phys. Rev. Lett. 122, 040504 (2019).
  * (13) V. Havlicek, A. D. Corcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, J. M. Gambetta, Nature. vol. 567, pp. 209 (2019)
  * (14) I. Cong, S. Choi and M. D. Lukin, Nature Physics 15, 1273 (2019).
  * (15) L. Lamata, Mach. Learn.: Sci. Technol. 1, 033002 (2020).
  * (16) V. Vapnik and A. Lerner, Automation and Remote Control, 24, (1963)
  * (17) W. S. Noble, Nature Biotechnology 24, 1565 (2006)
  * (18) J. A. K. Suykens and J. Vandewalle, Neural Processing Letters 9, 293 (1999)
  * (19) L. Zhaokai, L. Xiaomei, X. Nanyang, D. jiangfeng, Phys. Rev. Lett. 114, 140504 (2015)
  * (20) A. W. Harrow, A. Hassidim, S. Lloyd, Phys. Rev. Lett. 15, 150502 (2009)
  * (21) X.-D. Cai, C. Weedbrook, Z.-E. Su, M.-C. Chen, Mile Gu, M.-J. Zhu, L. Li, N.-L. Liu, C.-Y. Lu, J.-W. Pan, Phys. Rev. Lett. 110, 230501 (2013)
  * (22) S. Barz, I. Kassal, M. Ringbauer, Y. O. Lipp, B. Dakic, A. Aspuru-Guzik, P. Walther, Scientific Reports 4,115 (2014)
  * (23) J. Pan, Y. Cao, X. Yao, Z. Li, C. Ju, X. Peng, S. Kais and J. Du, Phys. Rev. A 89, 022313 (2014)
  * (24) E. Farhi, J. Goldstone, S. Gutmann, MIT-CTP/4610
  * (25) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, J. L. O’Brien, Nature Communications, 5:4213, (2014)
  * (26) K. Mitarai, M. Negoro, M. Kitagawa and K. Fujii, Phys. Rev. A 98, 032309 (2018)
  * (27) C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, P. J. Coles, LA-UR-19-29101
  * (28) X. Xu, J. Sun, S. Endo, Y. Li, S. C. Benjamin, X. Yuan, arXiv:1909.03898
  * (29) D. Deutsch, Proceedings of the Royal Society A. 400, 97 (1985)
  * (30) C. M. Dawson and M. A. Nielsen arXiv:quant-ph/0505030.
  * (31) B. Kraus and J.I. Cirac, Phys. Rev. A 63, 062309 (2001)
  * (32) G. Vidal, C. M. Dawson, Phys. Rev. A 69, 010301 (2004)
  * (33) M. Mottonen, J. J. Vartiainen, V. Bergholm, M. M. Salomaa, Phys. Rev. Lett. 93, 130502 (2004)
  * (34) V. V. Shende, I. L. Markov, S. S. Bullock, Physical Review A 69, 062321 (2004).
  * (35) F. Vatan and C. P. Williams, arXiv:quant-ph/0401178
  * (36) P. B. M. Sousa, R. V. Ramos, arXiv:quant-ph/0602174
  * (37) D. Hanneke, J. P. Home, J. D. Jost, J. M. Amini, D. Leibfried, D. J. Wineland, Nature Physics 6, 13 (2009)
  * (38) L. DiCarlo, J. M. Chow, J. M. Gambetta, Lev S. Bishop, B. R. Johnson, D. I. Schuster, J. Majer, A. Blais, L. Frunzio, S. M. Girvin, R. J. Schoelkopf, Nature 460, 240 (2009)
  * (39) X. Qiang, X. Zhou, J. Wang, C. M. Wilkes, T. Loke, S. O’Gara, L. Kling, G. D. Marshall, R. Santagati, T. C. Ralph, J. B. Wang, J. L. O’Brien, M. G. Thompson, J. C. F. Matthews, Nature Photonics 12, 534 (2018)
  * (40) T. Roy, S. Hazra, S. Kundu, M. Chand, M. P. Patankar and R. Vijay, Phys. Rev. Applied 14, 014072 (2020)
  * (41) M. Ezawa, Phys. Rev. Research 2, 023278 (2020).
  * (42) M. Ezawa, J. Phys. Soc. Jpn. 89, 124712 (2020)
  * (43) M. Ezawa, cond-mat/arXiv:2012.06124
  * (44) M. Ezawa, Phys. Rev. B 100, 165419 (2019).
  * (45) H. W. Kuhn and A. W. Tucker, Proceedings of 2nd Berkeley Symposium. Berkeley: University of California Press. pp. 481–492.

Acknowledgements

The author is very much grateful to E. Saito and N. Nagaosa for helpful
discussions on the subject. This work is supported by the Grants-in-Aid for
Scientific Research from MEXT KAKENHI (Grants No. JP17K05490 and No.
JP18H03676). This work is also supported by CREST, JST (JPMJCR16F1 and
JPMJCR20T2).

Author contributions

M.E. conceived the idea, performed the analysis, and wrote the manuscript.

Additional information

Competing financial and non-financial interests: The author declares no
competing financial and non-financial interests.


# Direct view of phonon dynamics in atomically thin MoS2

Tristan L Britt<EMAIL_ADDRESS>[  Qiuyang Li  [  Laurent P.
René de Cotret  [  Nicholas Olsen  [  Martin Otto  Syed Ali Hassan  [  Marios
Zacharias  Department of Mechanical and Materials Science Engineering, Cyprus
University of Technology, P.O. Box 50329, Limassol 3603, Cyprus  Fabio Caruso
Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-
Universität zu Kiel, Leibnizstraße 15, Kiel 24118, Germany  Xiaoyang Zhu  [
Bradley J Siwick  [<EMAIL_ADDRESS>
###### Abstract

Transition metal dichalcogenide monolayers and heterostructures are highly
tunable material systems that provide excellent models for physical phenomena
at the two-dimensional (2D) limit. While most studies to date have focused on
electrons and electron-hole pairs, phonons also play essential roles. Here, we
apply ultrafast electron diffraction and diffuse scattering to directly
quantify, with time and momentum resolution, electron-phonon coupling (EPC) in
monolayer molybdenum disulfide (MoS2) and phonon transport from the monolayer
to a silicon nitride (Si3N4) substrate. Optically generated hot carriers
result in a profoundly anisotropic distribution of phonons in the monolayer on
the $\sim$ 5 ps time scale. A quantitative comparison with ab-initio ultrafast
dynamics simulations reveals the essential role of dielectric screening in
weakening EPC. Thermal transport from the monolayer to the substrate occurs
with the phonon system far from equilibrium. While screening in 2D is known to
strongly affect equilibrium properties, our findings extend this understanding
to the dynamic regime.

###### keywords:

ultrafast, monolayer, transition metal dichalcogenide, phonon, coupling

McGillPhysics]Department of Physics, Center for the Physics of Materials,
McGill University, 3600 rue Université, Montréal, Québec H3A 2T8, Canada
NYUChem]Department of Chemistry, Columbia University, New York City, New York
10027, USA McGillPhysics]Department of Physics, Center for the Physics of
Materials, McGill University, 3600 rue Université, Montréal, Québec H3A 2T8,
Canada NYUChem]Department of Chemistry, Columbia University, New York City,
New York 10027, USA McGillPhysics]Department of Physics, Center for the
Physics of Materials, McGill University, 3600 rue Université, Montréal, Québec
H3A 2T8, Canada Institut FOTON - UMR 6082, Univ Rennes, INSA Rennes, 20 Avenue
des Buttes de Coesmes, Rennes 35700, France NYUChem]Department of Chemistry,
Columbia University, New York City, New York 10027, USA
McGillPhysics]Department of Physics, Center for the Physics of Materials,
McGill University, 3600 rue Université, Montréal, Québec H3A 2T8, Canada
Department of Chemistry, McGill University, 801 rue Sherbrooke Ouest,
Montréal, Québec H3A 0B8, Canada TMDC, DFT, BTE, EPC, UV, NIR

##  1 Introduction

Figure 1:  (a) Schematic of the pump-probe UED and UEDS experiment. The MoS2
monolayer-Si3N4 (abbreviated Si:N) heterostructure is photoexcited (pumped) by
an ultrashort laser pulse (400nm, 3.1 eV, fluence of 1.8 mJ-cm-2) at
$t=t_{0}$. After a controllable delay time $\tau$, an ultrashort electron
pulse probes the sample at normal incidence (in $\hat{z}$), and transmission
elastic/inelastic (phonon-diffuse) scattering is imaged at the detector
(right). A representative Brillouin zone of 1L-MoS2 with the Bragg peak
$\mathbf{H}$, the scattering vector $\mathbf{Q}$, and reduced scattering
vector $\mathbf{q}$ are indicated. (b) Schematic of the photodoped carrier
populations formed via photoexcitation (left) and after rapid carrier
thermalization (right). The color bar represents electron occupation.
Conduction band electrons appear yellow and valence band holes appear blue.
Hot carriers in the right panel are described by a Fermi-Dirac distribution at
an elevated charge carrier temperature $T_{0}^{\text{cc}}\simeq 1700$K.

2D Transition-metal dichalcogenides (TMDs) are currently subject of intense
research due to their exciting electronic and opto-electronic properties.
These arise from the highly correlated nature of their spin1, 2, valley3, 4,
electronic5, and vibrational6 degrees of freedom which can vary dramatically
with the number of layers, most notably with the breaking of spatial inversion
symmetry in monolayers. In molybdenum disulfide (MoS2) monolayers (1L), strong
spin-orbit coupling additionally 7 allows for direct control of the spin and
valley degrees of freedom 3, 8 and the possibility of exciting chiral phonons
9, 10. The coupling of monolayers in heterostructures 11, moiré superlattices
12, 13 or to underlying substrates 14 and modified dielectric environments 15
have opened further avenues for wide tunability in properties and control of
correlated phases.

Many effects in these (and other) 2D materials arise specifically from the
interplay between electronic and vibrational degrees of freedom, known as
electron-phonon coupling (EPC)16. EPC is a central concern in materials
physics generally, providing the fundamental origin of phenomena as diverse as
conventional superconductivity, charge-density waves and soft-mode phase
transitions, playing an important role in determining both charge and thermal
transport properties. However, unlike their 3D counterparts where EPC is a
fundamental property of the materials, EPC should be uniquely sensitive to the
environment of a 2D material where the Coulomb potential is known to be poorly
screened within the 2D layer 17.

Here we show that ultrafast electron diffraction (UED) combined with the
recently developed technique of ultrafast electron diffuse spectroscopy (UEDS)
can provide a momentum-resolved picture of the nonequilibrium phonon dynamics
that follow photocarrier generation in monolayer MoS2. This electron analog of
ultrafast diffuse x-ray scattering 18, 19 has previously provided momentum-
resolved information on inelastic electron-phonon scattering 20, 21, 22, 23,
soft phonon physics, 24 and polaron formation 25 in bulk materials. The
extension of this technique to monolayers is possible due to advancements in
instrumentation and the $>10^{5}$ enhancement in atomic scattering cross
section for electrons compared to x-rays. This provides access to novel
phenomena, including effects of the local dielectric environment on EPC within
the monolayer.

The technique also provides access to the microscopic, nonequilibrium phonon
population distribution at all momenta during phonon (heat) transport across
the monolayer-substrate interface. Comparing the experimental results with
first-principles calculations of the coupled electron-phonon dynamics26, 27
provides validation of this approach for the same in-monolayer processes
probed experimentally. Overall, the current work demonstrates a new approach
to probe EPC and nanoscale phonon transport within a monolayer/substrate
heterostructure that is broadly applicable to monolayer heterostructures and
moiré superlattices.

The opto-electronic properties, exciton physics 28 and ultrafast carrier
dynamics in monolayer TMDs have been extensively studied by both conventional
and ultrafast spectroscopies 29, 30, 31, 32, 33. To date, experimental
investigations of EPC effects and nonequilibrium phonon dynamics in monolayer
and few-layer TMDs have been limited to what can be learned about zone-
centered optical phonons via conventional and time-resolved Raman spectroscopy
34, 35, 36 or averaged atomic mean-squared displacements via UED37, 38, 39.
The momentum selectivity of phonon emission and the time-dependent anisotropic
phonon populations, expected throughout the Brillouin zone (BZ) of TMDs26, 21,
is hidden from view with these techniques. By contrast, UEDS provides
momentum-resolved information on EPC and its modification due to (i) the
dielectric environment, (ii) nonequilibrium phonon relaxation in 1L-MoS2
(including anharmonic decay within the monolayer), and (iii) phonon transport
to the underlying substrate as we show here.

In these pump-probe experiments (Fig 1a), the excitation photon energy of
3.1eV (400 nm) is resonant with vertical excitation near zone center,
photodoping electrons (holes) into the conduction (valence) bands as shown in
Fig 1b (left panel). These nonequilibrium charge carriers rapidly relax
towards the band edges ($K$ and $Q$ points) in momentum and energy via
electron-electron scattering 26 as indicated in Fig 1b (right panel). Carrier
cooling proceeds primarily through inelastic electron-phonon scattering, which
is expected to be strongly anisotropic in momentum due to the structure of the
valence and conduction band valleys (Fig. 2). It is this electron-phonon
coupling (from the perspective of the phonon system) that is followed with
time and mode resolution with UEDS through the elastic and inelastic phonon
scattered electron intensity measured as a function of scattering vector,
$I(\mathbf{Q},\tau)$ [Å2] (see Fig 1). See Methods for details.

Figure 2: The structure of the valence and conduction bands of 1L-MoS2
indicating available intra-valley and inter-valley inelastic phonon scattering
pathways for phonon-mediated carrier cooling. $\Gamma$ (black), $K$ (dark
blue), $M$ (light blue), and $Q$ (purple) phonon scattering processes are
shown for electron (hole) relaxation. Umklapp processes are not shown. Inset
in red is the $\Gamma$-centered rhomboidal primitive unit cell in reciprocal
space. For the conduction band, values are given relative to $\Delta$, the
Kohn-Sham direct band-gap. Figure 3: Photocarrier-phonon equilibration in
1L-MoS2 as measured via Bragg peak Debye-Waller dynamics. (a) Relative change
in (300) Bragg peak intensity following photoexcition (blue circles). First-
principles calculation of free-standing 1L-MoS2 (red curve), and including the
dielectric environment provided by the Si:N substrate (orange curve). The
best-fit (dashed lines) are given for long times, described by the 1D heat
kernel (see text). The inset highlights the long-time behavior of each signal
and the blue band represents the 1$\sigma$ uncertainty bound on the data
points. (b) Increase in MSD extracted from Bragg intensity by Eqn 1. Red,
orange, and dashed curves as in panel (a). (c) Average phonon mode temperature
(momentum integrated) from ab-initio simulations including the Si:N dielectric
environment.

##  2 Dielectric Screening of EPC in the 1L-MoS2 /Si:N Heterostructure

Following photoexcitation, all measured 1L-MoS2 Bragg peaks are suppressed 40,
41 due to the Debye-Waller effect42, 43 (Fig 3a). The relative peak intensity
is directly related to the increase in in-plane atomic mean-squared
displacements (MSD) $\langle u^{2}\rangle(\tau)$ given by:

$\langle u^{2}\rangle(\tau)-\langle u^{2}_{0}\rangle=$
$-\frac{3}{4\pi^{2}}\frac{\ln\\{\widetilde{I}_{0}(\mathbf{H}_{ij},\tau)/\widetilde{I}_{0}(\mathbf{H}_{ij},\tau<\tau_{0})\\}}{|\mathbf{H}_{ij}|^{2}}$
(1)

where $\langle u^{2}_{0}\rangle$ is the average of the equilibrium in-plane
MSD tensor, and $\mathbf{H}_{ij}$ is the scattering vector of the 1L-MoS2
Bragg reflection with Miller index ${(i,j)}$. The measured changes in MSD
following photoexcitation is shown in Fig 3b.

The transient rise in MSD provides an average measure of the rate at which
photocarrier excitation energy is transferred to phonons in the monolayer,
informing on carrier-phonon equilibration through EPC. We compare these
measurements directly with the results obtained by combining together ab-
initio calculations of nonequilibrium dynamics26 and scattering intensities27,
44 (see SI) in Figs. 3a-b. The ultrafast dynamics simulations for a free-
standing monolayer film are shown in red, predicting a much higher rate of MSD
increase than that observed in the experimental data. When the dielectric
environment provided by Si:N is included using a semi-infinite slab model with
no free parameters, there is quantitative agreement between the rise is MSD
measured and that predicted within experimental uncertainties.

The dielectric screening provided by Si:N in the heterostructure renormalizes
the EPC matrix elements according to $\tilde{g}_{mn\nu}({\bf
k,q})=g_{mn\nu}({\bf k,q})/\varepsilon^{\rm sub}_{\infty}$, where
$g_{mn\nu}({\bf k,q})$ is the coupling matrix element of a free-standing
monolayer (see SI 3.2.3). Here, the high-frequency dielectric constant of the
semi-infinite slab is the average permittivity between ambient vacuum
conditions and the Si:N substrate, namely
$\varepsilon^{\mathrm{sub}}_{\infty}=(1+\varepsilon_{\infty})/2$, where
$\varepsilon_{\infty}=7.8$ is dielectric

constant of bulk Si:N45, 46. The average value of these matrix elements for
the LA and LO ($E^{{}^{\prime}}$) phonon modes in free-standing 1L-MoS2 at $K$
are 19 meV and 23 meV respectively. These are each reduced by a factor of
$\varepsilon^{\mathrm{sub}}_{\infty}=4.4$ due to the presence of the Si:N
dielectric environment, providing a quantitative explanation for the
approximately order-of-magnitude reduction in the rate of photocarrier-phonon
energy transfer compared to predictions for a free-standing film (SI, Table
S1). The inclusion of the renormalized EPC matrix elements produces
quantitative agreement between first-principles and experimental results in
both amplitude and rate of change.

Figure 4: Momentum-resolved phonon re-equilibration dynamics (a) All-phonon
differential diffuse scattering pattern of 1L-MoS2 calculated from first
principles as $\Delta I=I(\mathbf{Q},\tau=5\text{
ps})-I(\mathbf{Q},T=300\text{ K})$. Inset (upper left) is the thermal
differential diffuse scattering pattern calculated as $\Delta
I=I(\mathbf{Q},T=380\text{ K})-I(\mathbf{Q},T=300\text{ K})$. The temperature
of 380 K corresponds to an effective lattice temperature as extracted from the
observed MSD at $\tau=5$ ps (see text), shown on the same color scale. Black
hexagons indicate BZ boundaries. Regions for which data is shown in (b-d) are
indicated with the matching color. (b-d) The relative change in diffuse
intensity at reduced scattering vectors (b) $K$, (c) $M$, (d) LA phonons at
${Q}$, and for (e) TA phonons at ${Q}$. Signals are obtained by integration
over the colored regions in (a), as well as over every visible BZ (see SI).
Acoustic signals are extracted by integrating over the segmented annuli given
in (a), with LA and TA distinction possible due to phonon scattering selection
rules (SI, Fig. S5). Red and orange curves as in Fig 3.

The excitation energy has equilibrated between carriers and phonons in 1L-MoS2
by $\tau\sim 10$ ps, as shown by the peak value of MSD in Fig 3b. The roll-
over and decay of MSD for $\tau\ >10$ ps indicates the reduction in
vibrational energy in the monolayer due to phonon transport into the Si:N
substrate (Fig. 3c). The reported rise-time for the MSD following
photoexcitation (i.e. the observed Debye-Waller decay of Bragg peak
intensities) in the case of 1L-MoS2 on sapphire is slower than reported
here37, 38, while the opposite is true for 1L-MoS2 on SiO239. This is as
expected based on the dielectric constants of these substrate materials,
providing additional evidence in support of our conclusions regarding
dielectric screening of the electron-phonon interaction in 1L-MoS2 .

Since photoexcitation was effectively uniform over the probed 250 $\mu$m
region ($<$ 10% variation), thermal transport from the photoexcited monolayer
into the Si:N substrate is well described as one-dimensional (1D) on
picosecond timescales; a delta-function heat impulse (following
photoexcitation of the monolayer) diffuses into the substrate. In this
geometry, with these initial conditions, an analytical solution for the
temperature of the monolayer is given by the Green’s function of the 1D heat
equation. This model provides an excellent fit of the data for long times in
Fig 3, far superior than a single exponential decays, yielding a thermal
conductivity $\kappa=313\pm 4$ Wm-1K-1 which is in reasonable agreement with
simulated and experimental values in the $\bf{c}$-axis of Si:N47, 48. We
conclude the thermal boundary resistance between 1L-MoS2 and Si:N in the
heterostructure is small.

##  3 Mode-Resolved Nonequilibrium Phonon Dynamics in 1L-MoS2

The increase in MSD determined from the Debye-Waller suppression of Bragg peak
intensities does not uniquely define the microscopic state of the phonon
system. By contrast, UEDS measurements together with ab-initio simulations
yield unprecedented details of mode- and momentum-resolved nonequilibriumn
phonon population distributions in the monolayer that underlie the changes in
MSD observed via the Bragg peak dynamics. The transient UEDS signals from
1L-MoS2 following photoexcitation show the strongest increases at the $K$, $M$
and $Q$ points of the BZ (Fig 4). The time- and momentum-resolved phonon
excitation dynamics at each of these points in the BZ are in good agreement
with our ab-initio predictions provided the effects of Si:N substrate
dielectric screening are included (Fig 4b-e).

These combined UEDS and first-principles analyses show that the nonequilibrium
state of the phonon system several picoseconds after photoexcitation is
profoundly anisotropic in momentum. This anisotropy is primarily determined by
the momentum-dependent electron-phonon interaction strength and the available
inelastic electron-phonon scattering pathways that are open to the hot
carriers. These pathways are constrained by the electronic band structure and
carrier distribution (Fig 2) and explain the observed and computed momentum-
dependent phonon heating dynamics. Within 5 - 10 ps, the carrier and phonon
systems in 1L-MoS2 have equilibrated with respect to the partition of
excitation energy, but the phonon system remains profoundly out of equilibrium
internally.

Previous work has demonstrated the possibility of defining a time-dependent
effective phonon temperature, $T_{\text{eff}}(\tau)$, that corresponds to the
observed MSD using the model49

$\langle
u^{2}\rangle(\tau)=\frac{3\hbar}{2M}\int\limits_{0}^{\infty}\coth\bigg{(}\frac{\hbar\omega}{2k_{B}T_{\text{eff}}(\tau)}\bigg{)}\frac{F(\omega)}{\omega}d\omega$
(2)

where $F(\omega)$ is the phonon density of states (DOS). However, such a
$T_{\text{eff}}(\tau)$ provides a misleading view of the nonequilibrium state
of the phonon system during carrier-phonon equilibration. This is illustrated
in Fig 4a where the nonequilibrium phonon-diffuse differential scattering
intensity at 5 ps is compared with a thermalized phonon-diffuse differential
intensity distribution at $T_{\text{eff}}(5\text{ps})$ = 380K (inset), the
effective temperature determined by Eqn 2 and the measured MSD at 5 ps (Fig
3). The phonon population distribution in 1L-MoS2 is still profoundly
nonthermal and not well described by an effective temperature.

Further relaxation of these anisotropic nonequilibrium phonons in 1L-MoS2
involves coupling processes internal to the monolayer and heat transfer
between the monolayer and Si:N substrate in the heterostructure. These
distinct processes are both resolved by these measurements. In Fig 5, the
diffuse intensity dynamics at $K$ out to 150 ps are compared against the MSD
dynamics extracted from the Bragg peaks, whose $\sim$ 50 ps decay time (single
exponential fit) indicates the cooling rate of the monolayer to the underlying
substrate. The observed decay of diffuse intensity at $K$ is in poor agreement
with these MSD dynamics, indicating a different process is involved. The
single exponential decay time constant determined for dynamics at $K$ is 25
ps, twice as rapid as the MSD dynamics, but in good agreement with the ab-
initio anharmonic decay rate of $E^{{}^{\prime}}$ optical phonons at $K$ (22
ps) to which UEDS is most sensitive (SI, Fig S8).

Figure 5: Optical phonon anharmonic decay vs monolayer cooling. $K$ valley
optical phonon diffuse scattering (blue) compared to MSD decay (green). The
single exponential decay of the diffuse intensity at $K$ (25 ps) is in good
agreement with the anharmonic decay rate of $E^{{}^{\prime}}$ optical phonons
computed via first-principles calculations (see SI Sec 3.3).

The observed decay of the MSD is, however, in reasonable agreement with the
decay in phonon-diffuse scattering measured for both the mid-BZ LA and TA
modes whose heating dynamics are shown in Fig 4d-e (SI, Table S1). For times
$<$ 30 ps, this sub-nanoscale phonon transport across the 1L-MoS2 /Si:N
heterostructure interface involves a profoundly nonequilibrium population of
phonons in the monolayer, including a much higher occupancy of high-wavevector
acoustic phonons than would be expected of a thermalized distribution. These
conditions are the result of phonon (heat) transport across the monolayer-
substrate interface occurring on similar timescales as the nonequilibrium
phonon relaxation within the monolayer itself. UEDS provides a direct,
momentum-resolved window on sub-nanoscale phonon transport in this far-from-
equilibrium regime.

##  4 Conclusion

These results provide a time- and momentum-resolved view of electron-phonon
coupling, anharmonic phonon decay and thermal transport in a monolayer 2D
semiconductor supported on a dielectric substrate. We show that the dielectric
environment provided by Si:N leads to strong renormalization of the electron-
phonon coupling in the monolayer. This finding fills a critical gap in our
understanding of how poor Coulomb screening at the 2D limit results in a range
of physical properties sensitive to the local dielectric environment. First-
principles ultrafast dynamics simulations using a recently developed framework
are in excellent agreement with these measurements. Combined, these approaches
provide a momentum-resolved protocol which can yield details of coupling
dynamics in 2D material systems and their heterostructures.

##  5 Methods

The total scattered intensity can be decomposed into

$I_{\text{all}}(\mathbf{Q},\tau)=I_{0}(\mathbf{Q},\tau)+I_{1}(\mathbf{Q},\tau)+\cdots$
(3)

The zeroth-order term $I_{0}(\mathbf{Q},\tau)$ is the elastic Bragg scattering
and the first-order contribution $I_{1}(\mathbf{Q},\tau)$ is the inelastic
single-phonon diffuse scattering that is the primary focus of this work.
Adopting phonon normal mode coordinates gives:

$I_{1}(\mathbf{Q},t)\propto\sum_{\nu}\underbrace{\frac{n_{\nu}(\mathbf{q},t)+1/2}{\omega_{\nu}(\mathbf{q},t)}}_{\textstyle\begin{array}[]{c}|a_{\mathbf{q}\nu}|^{2}\\\
\end{array}}\big{|}\mathcal{F}_{1\nu}(\mathbf{Q},t)\big{|}^{2}$ (4)

where the label $\nu$ indicates the specific phonon mode, $\bf{Q}$ is the
electron scattering vector, $\bf{q}$ is the reduced phonon wavevector (i.e.
$\bf{q}$ = $\bf{Q}$ \- $\bf{H}$, where $\bf{H}$ is the closest Bragg peak),
$a_{\bf{q}\nu}$ is the vibrational amplitude of mode $\nu$, $n_{\nu}$ is the
mode-resolved occupancy with energy $\hbar\omega_{\nu}$, and
$\mathcal{F}_{1\nu}$ are known as the one-phonon structure factors. $I_{1}$
provides momentum-resolved information on the nonequilibrium distribution of
phonons across the entire BZ, since it depends only on phonon modes with
wavevector $\bf{q}$ (Fig 1a). The $\mathcal{F}_{1\nu}$ are geometrical weights
that describe the relative strength of scattering from different phonon modes
and depend sensitively on the atomic polarization vectors
$\left\\{\bf{e}_{\bf{q}\nu\kappa}\right\\}$50. Most importantly,
$\mathcal{F}_{1\nu}\left(\bf{Q}\right)$ are relatively large when the phonon
mode $\nu$ is polarized parallel to the reduced scattering vector $\bf{q}$.
These phonon scattering selection rules mean that $\mathcal{F}_{1\nu}$ for the
out-of-plane (Z-polarized) modes and the optical modes of
$E^{{}^{\prime\prime}}$ symmetry are very weak in the geometry of these
experiments (SI, Fig S8). These experiments primarily probe the
$\mathbf{q}$-dependent population dynamics of the $E^{{}^{\prime}}$ optical
and LA/TA modes. Terms of higher-order than $I_{1}$ represent multi-phonon
scattering which have lower cross-sections and do not contribute significantly
to the interpretation of the 1L-MoS2 signals reported here27, 44. The 1L-MoS2
/Si:N specimens used in these experiments provide two distinct contributions
to $I_{\text{all}}$ that are both evident in Fig 1: (i) elastic scattering
from the amorphous Si:N substrate layer which is distributed as diffuse rings,
and (ii) the Bragg and phonon-diffuse scattering from the 1L-MoS2 . The
qualitatively different character of these signals makes the amorphous Si:N
contribution to the scattering signals,
$I^{\text{sub}}(\mathbf{Q},\tau)=I^{\text{sub}}(|\mathbf{Q}|)$, easily
subtracted from the data as a background. See SI for details.

The techniques of Liu et al were used to generate the 1L-MoS2 sample onto the
supporting Si:N substrate 51. A 150 nm-thick Au film was deposited onto a Si
wafer (from Nova Electronic Materials) with e-beam evaporation (0.05 nm/s).
Polyvinylpyrrolidone (PVP) solution (from Sigma Aldrich, mw 40000, 10% wt in
ethanol/acetonitrile wt 1/1) was spin-coated on the top of the Au film (1500
rpm, acceleration 500 rpm/s, 2 min) and then heated at 150 °C for 5 min. A
piece of thermal release tape was used to pick up PVP/Au film and was then
gently pressed onto a freshly cleaved bulk MoS2 single crystal (from HQ
graphene). The thermal release tape was gently lifted up and a MoS2 monolayer
is left, attached on the Au. The MoS2 monolayer was then transferred by
pressing the PVP/Au film with MoS2 monolayer onto a 30 nm thick amorphous
silicon nitrite (Si:N) TEM window (Norcada). The thermal release tape was then
removed by heating at 130 °C. The PVP layer was removed by dissolving in
deionized (DI) water for 2 h. The monolayer covered by Au was rinsed with
acetone and cleaned by O2 plasma for 3 min. After removing the Au cover with a
gold etchant solution (2.5g I2 and 10g KI in 100 mL DI water), the monolayer
on substrate was rinsed with DI water and isopropanol, and dried under N2
flow.

Electronic and vibrational properties of monolayer MoS2 have been obtained
from density functional and density functional perturbation theory,
respectively, and were used to calculate electron-phonon coupling matrix
elements via the Quantum Espresso52 and EPW53 codes. From here, the time-
propagation of the time-dependent Boltzmann equation (see SI Sec 3.2.1) yields
time-, momentum-, and phonon mode-resolved occupations, which completely
describe the nonthermal vibrational state. These occupations are used to
compute the dynamic all-phonon diffuse scattering intensity (including Bragg
scattering) within the Laval-Born-James theory27 as implemented in EPW (see SI
Sec 3.2.2).

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC), the Fonds de Recherche du Québec–Nature et
Technologies (FRQNT), the Canada Foundation for Innovation (CFI) and McGill
Fessenden Professorship.

B.J.S. and X.Y.Z conceived the research. XYZ acknowledges support for sample
preparation by the Materials Science and Engineering Research Center (MRSEC)
through US-NSF grant DMR-2011738. N.O and Taketo Handa, supervised by X.Y.Z,
performed PL measurements verifying the number of layers of the sample. T.L.B.
performed the experiments with assistance from L.P.R.de C., M.O., S.A.H., and
computed the phonon dispersion, equilibrium one-phonon structure factors, and
electronic band structure. The samples were prepared by Q.L. First-principles
calculations based on the TDBE were performed by F.C., who acknowledges
funding from the Deutsche Forschungsgemeinschaft (DFG) - Projektnummer
443988403. The dynamic diffuse scattering intensities and MSD were computed by
M.Z., who acknowledges financial support from the program
META$\Delta$I$\Delta$AKT$\Omega$P of the Cyprus University of Technology. We
acknowledge that the results of this research have been achieved using the
DECI resource Prometheus at CYFRONET in Poland with support from the PRACE
aisbl.

Experimental procedures, sample characterization, data preparation, and ab-
initio methodology are given in the supplementary information. Raw data is
available upon reasonable request with permission from the author(s).

## References

  * Zhou et al. 2019 Zhou, B. T.; Taguchi, K.; Kawaguchi, Y.; Tanaka, Y.; Law, K. T. Spin-orbit coupling induced valley Hall effects in transition-metal dichalcogenides. _Communications Physics_ 2019, _2_ , 26
  * Molina-Sánchez et al. 2013 Molina-Sánchez, A.; Sangalli, D.; Hummer, K.; Marini, A.; Wirtz, L. Effect of spin-orbit interaction on the optical spectra of single-layer, double-layer, and bulk MoS2. _Phys. Rev. B_ 2013, _88_ , 045412
  * Xiao et al. 2012 Xiao, D.; Liu, G.-B.; Feng, W.; Xu, X.; Yao, W. Coupled Spin and Valley Physics in Monolayers of ${\mathrm{MoS}}_{2}$ and Other Group-VI Dichalcogenides. _Phys. Rev. Lett._ 2012, _108_ , 196802
  * Beyer et al. 2019 Beyer, H.; Rohde, G.; Grubišić Čabo, A.; Stange, A.; Jacobsen, T.; Bignardi, L.; Lizzit, D.; Lacovig, P.; Sanders, C. E.; Lizzit, S.; Rossnagel, K.; Hofmann, P.; Bauer, M. 80% Valley Polarization of Free Carriers in Singly Oriented Single-Layer ${\mathrm{WS}}_{2}$ on Au(111). _Phys. Rev. Lett._ 2019, _123_ , 236802
  * Qiu et al. 2013 Qiu, D. Y.; da Jornada, F. H.; Louie, S. G. Optical Spectrum of ${\mathrm{MoS}}_{2}$: Many-Body Effects and Diversity of Exciton States. _Phys. Rev. Lett._ 2013, _111_ , 216805
  * Yao et al. 2008 Yao, W.; Xiao, D.; Niu, Q. Valley-dependent optoelectronics from inversion symmetry breaking. _Phys. Rev. B_ 2008, _77_ , 235406
  * Kadantsev and Hawrylak 2012 Kadantsev, E. S.; Hawrylak, P. Electronic structure of a single MoS2 monolayer. _Solid State Communications_ 2012, _152_ , 909–913
  * Song and Dery 2013 Song, Y.; Dery, H. Transport Theory of Monolayer Transition-Metal Dichalcogenides through Symmetry. _Phys. Rev. Lett._ 2013, _111_ , 026601
  * Zhang and Niu 2015 Zhang, L.; Niu, Q. Chiral phonons at high-symmetry points in monolayer hexagonal lattices. _Physical review letters_ 2015, _115_ , 115502
  * Zhu et al. 2018 Zhu, H.; Yi, J.; Li, M.-Y.; Xiao, J.; Zhang, L.; Yang, C.-W.; Kaindl, R. A.; Li, L.-J.; Wang, Y.; Zhang, X. Observation of chiral phonons. _Science_ 2018, _359_ , 579–582
  * Balch et al. 2020 Balch, H. B.; Evans, A. M.; Dasari, R. R.; Li, H.; Li, R.; Thomas, S.; Wang, D.; Bisbey, R. P.; Slicker, K.; Castano, I., et al. Electronically Coupled 2D Polymer/MoS2 Heterostructures. _Journal of the American Chemical Society_ 2020, _142_ , 21131–21139
  * Zhang et al. 2018 Zhang, N.; Surrente, A.; Baranowski, M.; Maude, D. K.; Gant, P.; Castellanos-Gomez, A.; Plochocka, P. Moiré intralayer excitons in a MoSe2/MoS2 heterostructure. _Nano letters_ 2018, _18_ , 7651–7657
  * Maity et al. 2022 Maity, I.; Mostofi, A. A.; Lischner, J. Chiral valley phonons and flat phonon bands in moire materials. _Physical Review B_ 2022, _105_ , L041408
  * Trainer et al. 2020 Trainer, D. J.; Wang, B.; Bobba, F.; Samuelson, N.; Xi, X.; Zasadzinski, J.; Nieminen, J.; Bansil, A.; Iavarone, M. Proximity-induced superconductivity in monolayer MoS2. _ACS nano_ 2020, _14_ , 2718–2728
  * Borghardt et al. 2017 Borghardt, S.; Tu, J.-S.; Winkler, F.; Schubert, J.; Zander, W.; Leosson, K.; Kardynał, B. E. Engineering of optical and electronic band gaps in transition metal dichalcogenide monolayers through external dielectric screening. _Physical review materials_ 2017, _1_ , 054001
  * Giustino 2017 Giustino, F. Electron-phonon interactions from first principles. _Rev. Mod. Phys._ 2017, _89_ , 015003
  * Zhu et al. 2015 Zhu, X.; Monahan, N. R.; Gong, Z.; Zhu, H.; Williams, K. W.; Nelson, C. A. Charge Transfer Excitons at van der Waals Interfaces. _Journal of the American Chemical Society_ 2015, _137_ , 8313–8320, PMID: 26001297
  * Trigo et al. 2013 Trigo, M.; Fuchs, M.; Chen, J.; Jiang, M.; Cammarata, M.; Fahy, S.; Fritz, D. M.; Gaffney, K.; Ghimire, S.; Higginbotham, A., et al. Fourier-transform inelastic X-ray scattering from time-and momentum-dependent phonon–phonon correlations. _Nature Physics_ 2013, _9_ , 790–794
  * Teitelbaum et al. 2021 Teitelbaum, S. W.; Henighan, T. C.; Liu, H.; Jiang, M. P.; Zhu, D.; Chollet, M.; Sato, T.; Murray, É. D.; Fahy, S.; O’Mahony, S., et al. Measurements of nonequilibrium interatomic forces using time-domain x-ray scattering. _Physical Review B_ 2021, _103_ , L180101
  * Stern et al. 2018 Stern, M. J.; René de Cotret, L. P.; Otto, M. R.; Chatelain, R. P.; Boisvert, J.-P.; Sutton, M.; Siwick, B. J. Mapping momentum-dependent electron-phonon coupling and nonequilibrium phonon dynamics with ultrafast electron diffuse scattering. _Phys. Rev. B_ 2018, _97_ , 165416
  * Waldecker et al. 2017 Waldecker, L.; Bertoni, R.; Hübener, H.; Brumme, T.; Vasileiadis, T.; Zahn, D.; Rubio, A.; Ernstorfer, R. Momentum-Resolved View of Electron-Phonon Coupling in Multilayer ${\mathrm{WSe}}_{2}$. _Phys. Rev. Lett._ 2017, _119_ , 036803
  * Chase et al. 2016 Chase, T.; Trigo, M.; Reid, A. H.; Li, R.; Vecchione, T.; Shen, X.; Weathersby, S.; Coffee, R.; Hartmann, N.; Reis, D. A.; Wang, X. J.; Dürr, H. A. Ultrafast electron diffraction from non-equilibrium phonons in femtosecond laser heated Au films. _Applied Physics Letters_ 2016, _108_ , 041909
  * Seiler et al. 2021 Seiler, H.; Zahn, D.; Zacharias, M.; Hildebrandt, P.-N.; Vasileiadis, T.; Windsor, Y. W.; Qi, Y.; Carbogno, C.; Draxl, C.; Ernstorfer, R.; Caruso, F. Accessing the Anisotropic Nonthermal Phonon Populations in Black Phosphorus. _Nano Letters_ 2021, _21_ , 6171–6178, PMID: 34279103
  * Otto et al. 2021 Otto, M. R.; Pöhls, J.-H.; René de Cotret, L. P.; Stern, M. J.; Sutton, M.; Siwick, B. J. Mechanisms of electron-phonon coupling unraveled in momentum and time: The case of soft phonons in TiSe2. _Science Advances_ 2021, _7_
  * René de Cotret et al. 2022 René de Cotret, L. P.; Otto, M. R.; Pöhls, J.-H.; Luo, Z.; Kanatzidis, M. G.; Siwick, B. J. Direct visualization of polaron formation in the thermoelectric SnSe. _Proceedings of the National Academy of Sciences_ 2022, _119_
  * Caruso 2021 Caruso, F. Nonequilibrium Lattice Dynamics in Monolayer MoS2. _The Journal of Physical Chemistry Letters_ 2021, _12_ , 1734–1740, PMID: 33569950
  * Zacharias et al. 2021 Zacharias, M.; Seiler, H.; Caruso, F.; Zahn, D.; Giustino, F.; Kelires, P. C.; Ernstorfer, R. Efficient first-principles methodology for the calculation of the all-phonon inelastic scattering in solids. _Physical review letters_ 2021, _127_ , 207401
  * Wang et al. 2018 Wang, G.; Chernikov, A.; Glazov, M. M.; Heinz, T. F.; Marie, X.; Amand, T.; Urbaszek, B. Colloquium: Excitons in atomically thin transition metal dichalcogenides. _Reviews of Modern Physics_ 2018, _90_ , 021001
  * Liu et al. 2019 Liu, F.; Ziffer, M. E.; Hansen, K. R.; Wang, J.; Zhu, X. Direct Determination of Band-Gap Renormalization in the Photoexcited Monolayer ${\mathrm{MoS}}_{2}$. _Phys. Rev. Lett._ 2019, _122_ , 246803
  * Lee et al. 2021 Lee, W.; Lin, Y.; Lu, L.-S.; Chueh, W.-C.; Liu, M.; Li, X.; Chang, W.-H.; Kaindl, R. A.; Shih, C.-K. Time-resolved ARPES Determination of a Quasi-Particle Band Gap and Hot Electron Dynamics in Monolayer MoS2. _Nano Letters_ 2021, _21_ , 7363–7370
  * He et al. 2020 He, M.; Rivera, P.; Van Tuan, D., et al. Valley phonons and exciton complexes in a monolayer semiconductor. _Nature Communications_ 2020, _11_ , 618
  * Li et al. 2015 Li, Y.; Dong, N.; Zhang, S.; Zhang, X.; Feng, Y.; Wang, K.; Zhang, L.; Wang, J. Giant two-photon absorption in monolayer MoS2. _Laser & Photonics Reviews_ 2015, _9_ , 427–434
  * Kime et al. 2017 Kime, G.; Leontiadou, M. A.; Brent, J. R.; Savjani, N.; O’Brien, P.; Binks, D. Ultrafast Charge Dynamics in Dispersions of Monolayer MoS2 Nanosheets. _The Journal of Physical Chemistry C_ 2017, _121_ , 22415–22421
  * Sohier et al. 2019 Sohier, T.; Ponomarev, E.; Gibertini, M.; Berger, H.; Marzari, N.; Ubrig, N.; Morpurgo, A. F. Enhanced electron-phonon interaction in multivalley materials. _Physical Review X_ 2019, _9_ , 031019
  * Quan et al. 2021 Quan, J.; Linhart, L.; Lin, M.-L.; Lee, D.; Zhu, J.; Wang, C.-Y.; Hsu, W.-T.; Choi, J.; Embley, J.; Young, C., et al. Phonon renormalization in reconstructed MoS 2 moiré superlattices. _Nature materials_ 2021, 1–6
  * Ferrante et al. 2022 Ferrante, C.; Battista, G. D.; López, L. E. P.; Batignani, G.; Lorchat, E.; Virga, A.; Berciaud, S.; Scopigno, T. Picosecond energy transfer in a transition metal dichalcogenide-graphene heterostructure revealed by transient Raman spectroscopy. _Proceedings of the National Academy of Sciences_ 2022, _119_ , e2119726119
  * Mannebach et al. 2015 Mannebach, E. M.; Li, R., et al. Dynamic Structural Response and Deformations of Monolayer MoS2 Visualized by Femtosecond Electron Diffraction. _Nano Letters_ 2015, _15_ , 6889–6895, PMID: 26322659
  * He et al. 2020 He, X.; Chebl, M.; Yang, D.-S. Cross-Examination of Ultrafast Structural, Interfacial, and Carrier Dynamics of Supported Monolayer MoS2. _Nano Letters_ 2020, _20_ , 2026–2033, PMID: 32031381
  * Luo et al. 2021 Luo, D.; Tang, J.; Shen, X.; Ji, F.; Yang, J.; Weathersby, S.; Kozina, M. E.; Chen, Z.; Xiao, J.; Ye, Y., et al. Twist-Angle-Dependent Ultrafast Charge Transfer in MoS2-Graphene van der Waals Heterostructures. _Nano Letters_ 2021, _21_ , 8051–8057
  * Grosso and Parravicini 2014 Grosso, G.; Parravicini, G. P. In _Solid State Physics (Second Edition)_ , second edition ed.; Grosso, G., Parravicini, G. P., Eds.; Academic Press: Amsterdam, 2014; pp 437–482
  * Fultz and Howe 2008 Fultz, B.; Howe, J. M. _Transmission Electron Microscopy and Diffractometry of Materials_ ; Springer Berlin Heidelberg: Berlin, Heidelberg, 2008; pp 163–221
  * Debye 1913 Debye, P. Interferenz von Röntgenstrahlen und Wärmebewegung. _Annalen der Physik_ 1913, _348_ , 49–92
  * Waller 1923 Waller, I. Zur Frage der Einwirkung der Wärmebewegung auf die Interferenz von Röntgenstrahlen. _Zeitschrift für Physik_ 1923, _17_ , 398–408
  * Zacharias et al. 2021 Zacharias, M.; Seiler, H.; Caruso, F.; Zahn, D.; Giustino, F.; Kelires, P. C.; Ernstorfer, R. Multiphonon diffuse scattering in solids from first principles: Application to layered crystals and two-dimensional materials. _Phys. Rev. B_ 2021, _104_ , 205109
  * Warlimont 2005 Warlimont, H. In _Springer Handbook of Condensed Matter and Materials Data_ ; Martienssen, W., Warlimont, H., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2005; pp 431–476
  * Balland and Glachant 1999 Balland, B.; Glachant, A. In _New Insulators, Devices and Radiation Effects_ ; Barbottin, G., Vapaille, A., Eds.; Instabilities in Silicon Devices; North-Holland, 1999; Vol. 3; pp 3–144
  * Hu et al. 2020 Hu, F.; Xie, Z.-P.; Zhang, J.; Hu, Z.-L.; An, D. Promising high-thermal-conductivity substrate material for high-power electronic device: silicon nitride ceramics. _Rare Metals_ 2020, _39_ , 463–478
  * Hirosaki et al. 2002 Hirosaki, N.; Ogata, S.; Kocer, C.; Kitagawa, H.; Nakamura, Y. Molecular dynamics calculation of the ideal thermal conductivity of single-crystal $\alpha$\- and $\beta-{\mathrm{Si}}_{3}{\mathrm{N}}_{4}$. _Phys. Rev. B_ 2002, _65_ , 134110
  * Waldecker et al. 2016 Waldecker, L.; Bertoni, R.; Ernstorfer, R.; Vorberger, J. Electron-Phonon Coupling and Energy Flow in a Simple Metal beyond the Two-Temperature Approximation. _Phys. Rev. X_ 2016, _6_ , 021003
  * René de Cotret et al. 2019 René de Cotret, L. P.; Pöhls, J.-H.; Stern, M. J.; Otto, M. R.; Sutton, M.; Siwick, B. J. Time- and momentum-resolved phonon population dynamics with ultrafast electron diffuse scattering. _Phys. Rev. B_ 2019, _100_ , 214115
  * Liu et al. 2020 Liu, F.; Wu, W.; Bai, Y.; Chae, S. H.; Li, Q.; Wang, J.; Hone, J.; Zhu, X.-Y. Disassembling 2D van der Waals crystals into macroscopic monolayers and reassembling into artificial lattices. _Science_ 2020, _367_ , 903–906
  * Giannozzi et al. 2017 Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. _Journal of Physics: Condensed Matter_ 2017, _29_ , 465901
  * Poncé et al. 2016 Poncé, S.; Margine, E.; Verdi, C.; Giustino, F. EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions. _Computer Physics Communications_ 2016, _209_ , 116–133


All of these claims in Section 2 are proved based  on   (<ref>),   (<ref>), and some properties of weighted Sobolev spaces $H^{\gamma}_{p,\Theta}(\cD)$ taken e.g from [25].  Since these properties in [25] hold true on arbitrary domains, the exactly same proofs of Section 2 work with $\cD$ replaced by $\cO$.

In the corresponding theorem to Theorem <ref> we do not need the additional condition for the initial condition since our domain $\domain$ is bounded and hence the distance to the boundary is bounded. For the comparison with conic domain case, see Remark <ref>.

For $m=1,\ldots,M$, let $\kappa_m$ denote the interior angle at the vertex $p_m$, and denote

\begin{equation}\label{kappas}
\kappa_0:=\max_{1\leq m\leq M}\kappa_m.
\end{equation}

Also, for each $m$,  let $\cD_m$ denote the conic domain in $\bR^2$ such that

\cO \cap B_{\varepsilon}(p_m) \cap \{p_m+x: x\in \cD_m\} = \cO \cap B_{\varepsilon}(p_m)

for all sufficiently small $\varepsilon>0$. Denote

$$\lambda^{\pm}_{c,\cL,\cO}:=\min_{m}\lambda^{\pm}_{c,\cL, \cD_m} \quad \text{if $\cL$ is non-random}

$$\lambda_{c,\cO}(\nu_1,\nu_2):=\min_{m}\lambda_c(\nu_1,\nu_2,\cD_m)\quad \text{ if $\cL$ is random}.

In Theorem <ref> below, we pose  the condition

\begin{equation}
\label{theta poly}
p(1-\lambda^+_{c,\cL,\cO})<\theta< p(1+\lambda^-_{c,\cL,\cO})
\end{equation}

if $\cL$ is non-random, and

\begin{equation}
\label{theta poly2}
p(1-\lambda_{c,\cO}(\nu_1,\nu_2))<\theta <p(1+\lambda_{c,\cO}(\nu_1,\nu_2))
\end{equation}

if $\cL$ is random.  By Proposition <ref>,

\lambda^{\pm}_{c,\cL,\cO} \geq \lambda_{c,\cO}(\nu_1,\nu_2) \geq  \sqrt{\frac{\nu_1}{\nu_2}}\cdot\frac{\pi}{\kappa_0}

with $\kappa_0$ from (<ref>).
Hence,  the following is a   sufficient condition for (<ref>) and (<ref>) :

\begin{equation}
\label{concrete}
\end{equation}

Here are our main results on polygonal domains.

Let $p\in[2,\infty)$, $\gamma \geq -1$,  and Assumption <ref> hold.  Also assume that

\begin{equation}
\label{theta application}
\end{equation}

and condition (<ref>) holds if $\cL$ is non-random, condition (<ref>) holds if $\cL$ is random.
Then for given $f^0\in\bK^{\gamma \vee 0}_{p,\theta+p,\Theta+p}(\domain,T)$, $\tbf=(f^1,\cdots,f^d) \in\bK^{\gamma +1}_{p,\theta,\Theta}(\domain,T,d)$,  $g\in\bK^{\gamma+1}_{p,\theta,\Theta}(\domain,T,\ell_2)$, and $u_0\in\bU^{\gamma+2}_{p,\theta,\Theta}(\cO)$, the equation

\begin{equation}\label{stochastic parabolic equation polygon}
d u =\left(\cL u+f^0+\sum_{i=1}^d f^i_{x^i}\right)dt +\sum^{\infty}_{k=1} g^kdw_t^k,\quad   t\in(0,T]\quad\,; \quad u(0,\cdot)=u_0
\end{equation}

admits  a unique solution $u$  in the class $\cK^{\gamma+2}_{p,\theta,\Theta}(\domain,T)$.
Moreover, the estimate

\begin{eqnarray}
\|u\|_{\cK^{\gamma+2}_{p,\theta,\Theta}(\domain,T)}
&\leq& N\big(\|f^0\|_{\bK^{\gamma\vee 0}_{p,\theta+p,\Theta+p}(\domain,T)}
&&\quad\quad+\|u_0\|_{\bU^{\gamma+2}_{p,\theta,\Theta}(\cO)}\big)\label{main estimate-2}
\end{eqnarray}

holds with a constant $N=N(\domain,p,\gamma,\nu_1,\nu_2,\theta,\Theta,T)$.

Since $d=2$ in this section, the range of $\Theta$ in (<ref>) coincides with $(d-1,d-1+p)$ which we have kept throughout this article.

Let $p\geq 2$, $\theta, \Theta\in \bR$ and $u\in \cK^{\gamma+2}_{p,\theta,\Theta}(\cO,T)$.

(i) If $\gamma+2-\frac{d}{p}\geq n+\delta$, where $n\in \bN_0$ and $\delta\in (0,1]$, then for any $k\leq n$,

|\rho^{k-1+\frac{\Theta}{p}} \rho^{(\theta-\Theta)/p}_{\circ} D^{k}u(\omega,t,\cdot)|_{C(\cO)}+
[\rho^{n-1+\delta+\frac{\Theta}{p}} \rho^{(\theta-\Theta)/p}_{\circ} D^{k} u(\omega,t,\cdot)]_{C^{\delta}(\cO)}<\infty

holds for almost all $(\omega,t)$.  In particular,

\begin{equation}
\label{eqn 8.10.21p}
|u(\omega,t,x)|\leq N(\omega,t) \rho^{1-\frac{\Theta}{p}}(x) \rho^{(-\theta+\Theta)/p}_{\circ}(x).
\end{equation}

(ii) Let

2/p<\alpha<\beta\leq 1, \quad \gamma+2-\beta-d/p \geq m+\varepsilon,

where $m\in \bN_0$ and $\varepsilon\in (0,1]$.  Put $\eta=\beta-1+\Theta/p$. Then for any $k\leq m$,

\begin{eqnarray}
\label{eqn 8.11.11p}
&&\bE \sup_{t,s\leq T}  \frac {\left|\rho^{\eta+k}  \rho^{(\theta-\Theta)/p}_{\circ} \left(D^ku(t)-D^ku(s)\right)\right|^p_{C(\cO)}}
{|t-s|^{p\alpha/2-1}}<\infty, \\
&& \bE \sup_{t,s\leq T}  \frac {\left[\rho^{\eta+m+\varepsilon}  \rho^{(\theta-\Theta)/p}_{\circ} \left(D^mu(t)-D^mu(s)\right)\right]^p_{C^{\varepsilon}(\cO)}}
{|t-s|^{p\alpha/2-1}} <\infty.  \label{eqn 8.11.21p}
\end{eqnarray}

The claims follow from the corresponding results of (<ref>) and (<ref>) mentioned in Theorem <ref>.

For the proof of  Theorem <ref>, we first prove the following estimate.

Let Assumptions in Theorem <ref> hold. Then there exists a constant $N=N(d,p,\theta,\Theta,\nu_1,\nu_2,\cO,T)$ such that  the a priori estimate

\begin{eqnarray}
\|u\|_{\cK^{\gamma+2}_{p,\theta,\Theta}(\domain,T)}&\leq& N\big(\|f^0\|_{\bK^{\gamma\vee 0}_{p,\theta+p,\Theta+p}(\domain,T)}
&&\quad\quad+ \|u_0\|_{\bU^{\gamma+2}_{p,\theta,\Theta}(\cO)}\big)\label{polygon a priori}
\end{eqnarray}

holds provided  that a solution $u\in \cK^{\gamma+2}_{p,\theta,\Theta}(\cD,T)$ to equation (<ref>)  exists.

First, choose a sufficiently small constant $r>0$ such that each $B_{3r}(p_m)$ contains only one vertex $p_m$ and intersects with only two edges for each $m=1,\,\ldots,\,M$. Then we choose a function $\xi\in\cC_c^{\infty}(\bR^2)$ satisfying

\begin{align*}
1_{B_r(0)}(x)\leq \xi(x)\leq 1_{B_{2r}(0)}(x)\quad\text{for all }x\in\bR^2.
\end{align*}

Let $\xi_m(x):=\xi(x-p_m)$ and $\xi_0:=1-\sum_{m=1}^M\xi_m$.
By the choice of $r$ and $\xi$, $supp(\xi_m)$s are disjoint and hence $0\leq \xi_0\leq 1$.
Moreover, $\xi_0(x)=1$ if $\rho_{V}(x)>2r$.

For $m=1,\ldots,M$, let $\cD_m$ be the angular (conic) domain centered at $p_m$ with interior angle $\kappa_m$ such that $\cD_m\cap B_{3r}(p_m)=\cO\cap B_{3r}(p_m)$.

Now let $G$ be a $C^1$-domain in $\cO$ such that

\xi_0(x)=0\quad\text{for }x\in\cO\setminus G\quad\text{and}\quad \inf_{x\in G} \rho_{\circ}(x)\geq c>0\text{ with a constant}\; c.

Then, due to the choices of  $\xi_m$ and  $\cD_m$ $(m=1,\ldots,M)$,    (<ref>) and (<ref>)  together easily yield

\begin{align*}
\|\xi_m v\|^p_{K^n_{p,\theta,\Theta}(\cO)} \sim\|\xi_m v\|_{K^n_{p,\theta,\Theta}(\cD_m)},\quad m=1,\ldots,M.
\end{align*}

for any $\theta,\Theta\in\bR$, $n\in\{0,1,2,\ldots\}$, and $v\in K^n_{p,\theta,\Theta}(\cO)$.  Similarly,

\|\xi_0 v\|^p_{K^n_{p,\theta,\Theta}(\cO)}\sim \int_{G} |\xi_0v|^p \rho^{\Theta-d}dx\sim \|\xi_0 v\|^p_{H^n_{p,\Theta}(G)},

and the same relations hold for $\ell_2$-valued functions. Denote

\bH^{\gamma}_{p,\Theta}(G,T):=L_p(\Omega\times (0,T], \cP; H^{\gamma}_{p,\Theta}(G)),

\bH^{\gamma}_{p,\Theta}(G,T,\ell_2):=L_p(\Omega\times (0,T], \cP; H^{\gamma}_{p,\Theta}(G;\ell_2)).

Then, the above observations  in particular imply

\begin{align}\label{partition-of-unity.eq}
\|v\|_{\bK^n_{p,\theta,\Theta}(\cO,T)} \sim \left(\|\xi_0 v\|_{\bH^n_{p,\Theta}(G,T)}+\sum_{m=1}^M\|\xi_m v\|_{\bK^n_{p,\theta,\Theta}(\cD_m,T)}\right)
\end{align}

for any $v\in \bK^n_{p,\theta,\Theta}(\cO,T)$, where $n\in \{0,1,2,\cdots\}$.

Now, for each $m=1,\ldots,M$ we define $u_m:=\xi_mu$. Then, since $\gamma+2\geq 1$,  $u_m$ belongs to $\bK^{1}_{p,\theta-p,\Theta-p}(\cD_m,T)$. Also, $\xi_0 u$ belongs to $\bH^{1}_{p,\Theta-p}(G,T)$.  Note that each $u_m$ satisfies

\begin{equation}\label{equation for m}
d(u_m)=\left(\cL u_m+f^0_m+\sum_{i=1}^d (f^i_m)_{x^i}\right)dt+\sum_k g^{k}_m dw_t^k,\quad t\in (0,T]
\end{equation}

in the sense of distributions on  $\cD_m$ with the initial condition $
u_m(0,\cdot)=\xi_m u_0$ and $\xi_0 u$ satisfies

\begin{equation}\label{equation for m=0}
d(\xi_0 u)=\left(\cL (\xi_0 u)+f^0_0+\sum_{i=1}^d (f^i_0)_{x^i}\right)dt+\sum_k g^{k}_0 dw_t^k,\quad t\in (0,T]
\end{equation}

in the sense of distributions on $G$ with the initial condition $w(0,\cdot)=\xi_0 u_0$,  where

\begin{equation}\label{f g for m}
f^0_m=f^0\xi_m-\sum_{i=1}^d f^i (\xi_m)_{x^i}-u \cL(\xi_{m}), \quad f^i_m=2\sum_{j=1}^d a^{ij}u\,(\xi_{m})_{x^j},\quad
g_m=\, g\xi_m
\end{equation}

for $m=0,1,2,\ldots,M$.

Since $supp(\xi_m) \subset \overline{B_{2r}(p_m)}$ and $(\xi_m)_x=0$ on a neighborhood of $p_m$ for $m=1,\ldots,M$, we have

\begin{align*}
\|u (\xi_m)_x\|_{\bL_{p,\theta,\Theta}(\cO,t)}+\|u(\xi_m)_{xx}\|_{\bL_{p,\theta+p,\Theta+p}(\cO,t)}\leq N\|u\|_{\bL_{p,\theta,\Theta}(\cO,t)}
\end{align*}

for $t\leq T$, where $N$ depends only on $\cO,\,p,\,\theta$ and $\Theta$.

Hence, for $m=1,\,\ldots,\,M$, by Theorems <ref> and <ref>, which our range of $\theta$ allows us to use,   we have for any $t\leq T$,

\begin{align*}
&\quad \quad \|\xi_m u\|_{\bK^1_{p,\theta-p,\Theta-p}(\cD_m,t)}
\\
& \leq N\big(\|f^0_m\|_{\bL_{p,\theta+p,\Theta+p}(\cD_m,t)}+ \sum_{i=1}^d \|f^i_m\|_{\bL_{p,\theta,\Theta}(\cD_m,t)}+\|g_m\|_{\bL_{p,\theta,\Theta}(\cD_m,t,\ell_2)}\nonumber\\
&\hspace{1cm}+\|\xi_m u_0\|_{\bU^1_{p,\theta,\Theta}(\cD_m)}\big)
\\
& \leq N\big(\|u\|_{\bL_{p,\theta,\Theta}(\cO,T)}+\|f^0\|_{\bL_{p,\theta+p,\Theta+p}(\cO,T)}
+\|\tbf\|_{\bL_{p,\theta,\Theta}(\cO,T,d)} +\|g\|_{\bL_{p,\theta,\Theta}(\cO,T,\ell_2)}\\
\end{align*}

For $m=0$, by  <cit.> (or <cit.>),
we  have

\begin{align*}
&\quad \quad\quad \|\xi_0 u\|_{\bH^1_{p,\Theta-p}(G,t)}\\
&\leq N\left(\|f^0_0\|_{\bL_{p,\Theta+p}(G,t)}+ \sum_{i=0}^d \|f^i_0\|_{\bL_{p,\Theta}(G,t)}+\|g_0\|_{\bL_{p,\Theta}(G,t,\ell_2)}+\|\xi_0u_0\|_{L_p(\Omega;H^{1-2/p}_{p,\Theta+2-p}(G))}\right)
\\
&\leq N\big(\|u\|_{\bL_{p,\theta,\Theta}(\cO,T)}+\|f^0\|_{\bL_{p,\theta+p,\Theta+p}(\cO,T)}
+\|\tbf\|_{\bL_{p,\theta,\Theta}(\cO,T,d)} +\|g\|_{\bL_{p,\theta,\Theta}(\cO,T,\ell_2)}\\
&\hspace{1cm}+ \|u_0\|_{\bU^1_{p,\theta,\Theta}(\cO)}\big).
\end{align*}

Summing up over all $m=0,\,\ldots,\,M$ and using (<ref>), for each $t\leq T$, we have

\begin{align*}
\\
&\leq N \big(\|u\|_{\bL_{p,\theta,\Theta}(\cO,t)}+\|f^0\|_{\bL_{p,\theta+p,\Theta+p}(\cO,T)}
+\|\tbf\|_{\bL_{p,\theta,\Theta}(\cO,T)} +\|g\|_{\bL_{p,\theta,\Theta}(\cO,T,\ell_2)}\\
&\hspace{1cm}  \|u_0\|_{\bU^1_{p,\theta,\Theta}(\cO)}\big).
\end{align*}

Using this and  the polygonal versions of (<ref>) and (<ref>),  which mentioned in Theorem <ref>,  we get, for each $t\leq T$,

\begin{align*}
&\quad\quad \|u\|^p_{\cK^1_{p,\theta,\Theta}(\cO,t)}\\
\leq &\,N\int^t_0\|u\|^p_{\cK^1_{p,\theta,\Theta}(\cO,s)}ds
\\
\end{align*}

Applying Gronwall's inequality, we further obtain

\begin{align*}
&\quad\quad \|u\|_{\cK^1_{p,\theta,\Theta}(\cO,T)}\\
&\leq N\left(\|f^0\|_{\bL_{p,\theta+p,\Theta+p}(\cO,T)}+\|\tbf\|_{\bL_{p,\theta,\Theta}(\cO,T)}+\|g\|_{\bL_{p,\theta,\Theta}(\cO,T,\ell_2)}+\|u_0\|_{\bU^1_{p,\theta,\Theta}(\cO)}\right).
\end{align*}

This and the polygonal version of Lemma <ref>, which is mentioned in Theorem <ref>,  yield a priori estimate (<ref>).
The lemma is proved.

The following is a $C^1$-domain version of Lemma <ref>. We use it in the  proof of Theorem <ref> below.

Let $G$ be a bounded $\cC^1$ domain in $\bR^d$ and let $p_j\in[2,\infty)$, $\Theta_j\in (d-1,d-1+p_j)$ for $j=1,2$. Assume that $u\in \bH^1_{p_1,\Theta_1-p_1}(G,T)$

\begin{align*}
du=(\cL u+f^0+\sum_{i=1}^d f^i_{x^i})\,dt+\sum_kg^kdw^k_t,\quad t\in(0,T]\quad
\end{align*}

in the sense of distributions on $G$ with the initial condition $u(0,\cdot)=u_0(\cdot)$ and $f^0$, $f^i$ ($i=1,2,\ldots,$), $g$, $u_0$ satisfying

f^0\in \bL_{p_j,\Theta_j+p_j}(G,T)\cap \bL_{p_j,d+p_j}(G,T),\quad f^i \in \bL_{p_j,\Theta_j}(G,T)\cap \bL_{p_j,d}(G,T), \, i=1,\cdots,d,

g\in \bL_{p_j,\Theta_j}(G,T,\ell_2)\cap \bL_{p_j,d}(G,T,\ell_2),

\quad u_0\in L_p(\Omega,\rF_0;H^{1-2/p_j}_{p_j,\Theta_j+2-p_j}(G))\cap  L_p(\Omega,\rF_0;H^{1-2/p_j}_{p_j,d+2-p_j}(G))

for both $j=1$ and $j=2$.
Then $u$ belongs to $\bH^1_{p_2,\Theta_2-p_2}(G,T)$.

See <cit.>. We remark that  only $\Delta$ is considered in [2], however the proof of   <cit.> works for general case without any changes since the proof depends only on   <cit.> (or <cit.>),
which involves operators having coefficients measurable in $(\omega,t)$ and continuous in $x$.

We recall $d=2$ in this section.

Due to Lemma <ref>, we only need to prove the existence result. Furthermore, relying on standard approximation argument, we may assume

f^0\in\bK^{\infty}_c(\cO,T),\quad \tbf\in\bK^{\infty}_c(\cO,T,2), \quad  g\in\bK^{\infty}_c(\cO,T,\ell_2)\quad u_0\in \bK^{\infty}_c(\cO).

Considering $u-u_0$ as usual, we may assume $u_0\equiv 0$.
Also, note that $g^k=0$ for all large $k$ (say, for all $k> N$), and each $g^k$ is of the type $\sum_{j=1}^{n(k)} 1_{(\tau^k_{j-1},\tau^k_j]}(t) h^{kj}(x)$, where $\tau^k_j$ are bounded stopping times and $h^{jk}\in C^{\infty}_c(\cO)$.
Thus the function $v$ defined by

v(t,x):=\sum_{k=1}^{\infty}\int^t_0 g^k dw^k_s=\sum_{k\leq N} \sum_{j\leq n(k)} \left(w^k_{\tau^k_j \wedge t}-w^k_{\tau^k_{j-1} \wedge t}\right) h^{kj}(x)

is  infinitely differentiable in $x$ and vanishes near the boundary of $\cO$. Consequently $v$ belongs to $\cK^{\nu+2}_{p,\theta,\Theta}(\cO,T)$ for any $\nu, \theta, \Theta \in \bR$ as we consult with Definition <ref>.  Now, $u$ satisfies equation (<ref>) if and only if $\bar{u}:=u-v$ satisfies

d\bar{u}=(\cL\bar{u}+\bar{f}^0+\sum_{i=1}^2 f^i_{x^i})dt, \quad t\in(0,T]\quad;\quad \bar{u}(0,\cdot)= 0,

where $\bar{f}^0=f^0+\cL v$.  Hence, considering $\bar{f}^0$ in place of $f^0$, to prove the existence we may further assume $g=0$.

Then, by the  classical results without weights  for $p=2$,  (see, e.g.  [27] or <cit.>),
there exists a  solution $u$ in $\cK^1_{2,2,2}(\cO,T)$  to equation (<ref>), which now is simplified as

u_t=\cL u+f^0+\sum_{i=1}^2 f^i_{x^i}, \quad t\in(0,T]\quad; \quad  u(0,\cdot)=0.

By Theorem A in [1] (or see  estimate (2.11) and proof of Theorem 2.4 in  [12] for more detail), for any $r>4$, we have

\begin{equation}
\label{bound}
\bE \sup_{t,x} |u(t,x)|^p \leq N \bE \|\,|f^0|+|\tbf|\,\|^p_{L_r((0,T]\times \cO))} <\infty.
\end{equation}

Now we prove $u\in \cK^1_{p,\theta,\Theta}(\cO,T)$ using Lemma <ref> and  Lemma <ref> along with $u\in \cK^1_{2,2,2}(\cO,T)$.   Define $u_m:=\xi_m u$ in the same way we did in  the proof of Lemma <ref>. Then $\xi_m u$ satisfies (<ref>) in the sense of distributions on  $\cD_m$ for $m=1,\ldots,M$ and $\xi_0u$ satisfies (<ref>)
on $G$ for $m=0$ with the same $f^0_m,\,f^i_m,\, \xi_mu_0$ as in (<ref>).
Note that since $f^0,\tbf$ are bounded and $f^0,\tbf, (\xi_m)_x, (\xi_m)_{xx}$ vanish near vertices,  we have for any $\theta\in \bR$, $q\geq 2$ and $1<\Theta<1+q$,

\begin{eqnarray*}
&&\|f^0_m\|^q_{\bL_{q,\theta+q,\Theta+q}(\cO,T)} +\sum_{i=1}^2\|f^i_m\|^q_{\bL_{q,\theta,\Theta}(\cO,T)}\\
&\leq& N \bE \int^T_0 \int_{\cO}(1+|u|^q) \rho^{\Theta-2}dx\leq N \left(\int_{\cO}\rho^{\Theta-2}dx\right) \bE \sup_{t,x}(1+|u|^q)<\infty.
\end{eqnarray*}

For the last inequality we used (<ref>) and the fact $\Theta-2>-1$.
Hence, $f^0_m, \,f^i_m$ along with $\xi_mu_0$ satisfy assumptions in Lemma <ref> and  Lemma <ref>. Consequently  $\xi_m u \in \bK^1_{p,\theta-p,\Theta-p}(\cD_m,T)$ as $\xi_m u \in \cK^1_{p,\theta,\Theta}(\cD_m,T)$ for $m=1,2,\cdots, M$ and $\xi_0 u \in \bH^1_{p,\Theta-p}(G,T)$.  These and  (<ref>)  with $n=1$ yield $u\in \bK^1_{p,\theta-p,\Theta-p}(\cO,T)$ and in turn $u\in \cK^1_{p,\theta,\Theta}(\cO,T)$.

Finally, the polygonal version of Lemma <ref>, prescribed in Theorem <ref>,  proves that  the solution $u$ found above actually belongs to the space $u\in \cH^{\gamma+2}_{p,\theta,\Theta}(\cO,T)$.  The theorem is proved.

[1]
D. G. Aronson, Non-negative solutions of linear parabolic equations,
Ann. Scuola Norm. Sup. Pisa 22 (1968), no. 3, 607-694.

[2]
P. A. Cioica, K. Kim, and K. Lee,  On the regularity of the stochastic heat equation on polygonal domains in $\mathbb{R}^2$, J. Differential Equations 267 (2019), 6447–6479.

[3]
P. A. Cioica, K. Kim, K. Lee, and F. Lindner,  An $L_p$-estimate for the stochastic heat equation on an angular domain in $\mathbb{R}^2$, Stoch. Partial Differ. Equ. Anal. Comput. 6 (2018), no.  1, 45–72.

[4]
F. Flandoli, Dirichlet boundary value problem for stochastic parabolic equations: compatibility relation and regularity of solutions,
Stochastics Stochastis Rep., 29 (1990), no. 3, 331–357.

[5]
J. V. Beervan, M. Veraar and L. Weis, Stochastic Maximal $L^p$-regularity,
Annals of probab., 40 (2012), no. 2, 788–812.

[6]
D. Kim, Elliptic equations with nonzero boundary conditions in weighted Sobolev spaces,
J. Math. Anal. Appl. 337 (2008), 1465–1479.

[7]
K. Kim, K. Lee and J. Seo, A refined Green's function estimate of the time measurable parabolic operators with conic domains, to appear in Potential Anal.  https://doi.org/10.1007/s11118-020-09886-w.

[8]
K. Kim, K. Lee and J. Seo, A weighted sobolev regularity theory of the parabolic equations with measurable coefficients on conic domains in $\bR^d$, to appear in J. Differential Equations, arXiv:2103.10049.

[9]
K. Kim,  A weighted Sobolev space theory of parabolic stochastic PDEs on non-smooth domains,
Journal Theoret. Probab.,
27 (2014), no. 1, 107–136.

[10]
K. Kim,  On $L_p$-theory of stochastic partial differential equations of divergence form in $C^1$ domains,
Probab. Theory Relat. Fields, 130 (2004), 473-492.

[11]
K. Kim,  On stochastic partial differential equations with variable coefficients in $C^1$ domains,
Stochastic Processes and their Applications.,
112 (2004), no. 2, 261–283.

[12]
K. Kim, $L_q(L_p)$ theory and Hölder estimates for parabolic SPDEs,
Stoch. Processes and their applications, 114 (2004), 313-330.

[13]
K.Kim and N. V. Krylov, On the Sobolev space theory of parabolic and elliptic equations in $C^1$ domains, SIAM J. Math. Anal.,
36 (2004), no. 2, 618–642.

[14]
K. Kim and N. V. Krylov, On SPDEs with variable coefficients in one space dimension,  Potential Anal., 21 (2004), 209-239.

[15]
V. A. Kozlov, A. Nazarov,  The Dirichlet problem for non-divergence parabolic equations
with discontinuous in time coefficients in a wedge, Math. Nachr.  287
(2014), no. 10, 1142–1165.

[16]
N. V. Krylov, Some properties of traces for stochastic and deterministic parabolic weighted Sobolev spaces, J. Functional Anal. 183 (2001), 1–41.

[17]
N. V. Krylov, An analytic approach to SPDEs, in Stochastic partial differential equations: six perspectives, Math. Surveys Monogr., vol. 64, R. Carmona and B. Rozovskii, eds., AMS, Providence, RI, (1999), 185–-242.

[18]
N. V. Krylov, Some properties of weighted spaces in $\bR^d_+$,
Ann. Scuola Norm. Sup. Pisa Cl. Sci., 28 (1999), no. 4, 675–693.

[19]
N. V. Krylov and S. V. Lototsky, A Sobolev space theory of SPDEs with constant coefficients in a half space, SIAM J. Math. Anal., 31 (1999), no. 1, 19–33.

[20]
N. V. Krylov and S. V. Lototsky, A Sobolev space theory of SPDEs with constant coefficients in a half line, SIAM J. Math. Anal., 30 (1999), no. 2, 298–325.

[21]
N. V. Krylov, Weighted Sobolev spaces and Laplace's equation and the heat equations in a half space,
Commu. in Partial Differential Equations, 24 (1999), no 9&10, 1611–1653.

[22]
N. V. Krylov, On $L_p$-theory of stochastic partial differential equations,
SIAM J. Math. Anal., 27 (1996), 313-340.

[23]
A. Kufner,  Weighted Sobolev spaces, John
Wiley and Sons Inc, (1984).

[24]
S. V. Lototsky,
Linear stochastic parabolic equations, degenerating on the
boundary of a domain, Electronic Journal of Probability,
6 (2001), no. 24, 1–14.

[25]
S. V. Lototsky,
Sobolev spaces with weights in domains and boundary value problems for degenerate elliptic equations,
Methods Appl. 1 (2000), no. 1, 195–204.

[26]
A. Nazarov,
$L_p$-estimates for a solution to the Dirichlet problem and to the Neumann problem for the heat equation in
a wedge with edge of arbitrary codimension,  Probl. Mat. Anal. 22 (2001), 126–159 (Russian); English transl.: J. Math.
Sci. 106 (2001), no.3, 2989–3014.

[27]
B. L. Rozovsky,  Stochastic Evolution Systems. Linear Theory and Applications to Nonlinear Filtering, translated from the Russian by A. Yarkho,  Mathematics and its Applications (Soviet Series) 35, Kluwer Academic Publishers Group, Dordrecht, (1990).

[28]
V. A. Solonnikov, ${L}_p$-estimates for solutions of the heat
equation in a dihedral angle, Rend. Mat. Appl. (7) 21 (2001),
no. 1, 1–15.