# A Semi-Tensor Product based Circuit Simulation for SAT-sweeping

Hongyang Pan1, Ruibing Zhang2, Yinshui Xia2, Lunyao Wang2, Fan Yang1 and Xuan
Zeng1, Zhufei Chu2 1 Fudan University, 2 Ningbo University
{yangfan<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

In recent years, circuit simulators and Boolean satisfiability (SAT) solvers
have been tightly integrated to provide efficient logic synthesis and
verification. Circuit simulation can generate highly expressive simulation
patterns that can either enumerate or filter out most candidates for
synthesis. Subsequently, SAT solvers are employed to check those that remain,
thereby making the logic synthesis process more efficient. This paper
introduces a novel circuit simulator of $k$-input lookup table ($k$-LUT)
networks, based on semi-tensor product (STP). STP-based simulators use
computation of logic matrices, the primitives of logic networks, as opposed to
relying on bitwise logic operations for simulation of $k$-LUT networks.
Experimental results show that our STP-based simulator reduces the runtime by
an average of 7.2$\times$. Furthermore, we integrate this proposed simulator
into a SAT-sweeping engine known as SAT sweeper. Through a combination of
structural hashing, simulation, and SAT queries, SAT sweeper simplifies logic
networks by systematically merging graph vertices from input to output. To
enhance the efficiency, we used STP-based exhaustive simulation, which
significantly reduces the number of false equivalence class candidates,
thereby improving the computational efficiency by reducing the number of SAT
calls required. When compared to the SOTA SAT sweeper, our method demonstrates
an average 35% runtime reduction.

###### Index Terms:

logic synthesis, semi-tensor product of matrices, simulation, SAT-sweeping

##  I Introduction

It is common for logic synthesis and verification applications to integrate
the _Boolean satisfiability_ (SAT) solver with a random or guided simulator.
Currently, logic optimization methods often encode problems as SAT problems to
enhance optimization quality through the use of efficient SAT solvers, where
simulators can assist it in reducing candidates in order to avoid time-
consuming satisfiable runs. Additionally, the SAT solver can find a _counter-
example_ (CE), which the simulator can use to falsify other properties to save
future SAT calls [1]. There are many applications that take advantage of this
integration, such as SAT-sweeping. In SAT-sweeping, the simulator can
pregenerate a set of simulation patterns for a given Boolean network, most
non-equivalence can be effectively eliminated by simply comparing simulation
signatures. After each satisfiable SAT call, it also simulates each CE
immediately in order to disprove as many properties as possible [2, 3, 4].

The efficiency of simulation, encompassing both initial simulation and CE
simulation, plays an important role in SAT-sweeping, with the principal aim of
minimizing the number of expensive NP-hard SAT solver calls. Certain
algorithms incorporate partial simulation as a strategy to expedite simulation
speed. However, in scenarios where the pattern set is insufficiently
expressive, the equivalence classes tend to enlarge, and the refinement via CE
may become computationally more expensive. In practice, the equivalence class
usually contains a few percent of the total gates in a valid merge. Ideally,
we only need to simulate these gates, which must be simulated, with exhaustive
patterns for better runtime.

In this paper, we propose a _$k$ -input lookup table_ ($k$-LUT) circuit
simulator based on _semi-tensor product_ (STP). Modern circuit simulators
support bitwise logic operations, utilizing fast bit-parallel simulation to
enhance efficiency. Nonetheless, simulating $k$-LUT networks poses a unique
challenge, as it does not readily exploit these bit-parallel capabilities.
$k$-LUT simulation must obtain information in turn regarding each simulation
pattern by traversing all nodes in a topological order before computing the
output values of each node. The STP method works by matrices, utilizing logic
matrices for the definition of Boolean variables to prove the basic logic
properties [5]. As primitives in the logic network, logic matrices, which
convert logical reasoning into mathematical computations and preserve
topological information between circuits, are used to simulate the $k$-LUT
network. Motivated by this reasoning, we integrate our STP-based simulator
into the SAT-sweeping engine, where efforts are made in order to run as many
simulations as possible within the given runtime budget, reducing equivalence
classes to avoid unnecessary SAT solver calls.

The main contributions are: (1) a more efficient simulator: compared with
other bitwise logic operation based simulation implementations, the STP-based
simulator does not require specific logic operators, and the output values of
any node can be computed by one matrix pass; (2) the integration of STP-based
simulator into SAT sweeper _& fraig_: we first use SAT-guidance initial
simulations to refine equivalence classes while adding the constant node
substitution, which can significantly reduce the number of false candidates
for merging [6]. Then, map the nodes of non-equivalence classes to $k$-LUTs
and simulate the nodes of equivalence classes with exhaustive simulation
pattern to further refine candidate equivalence classes. In addition, for the
candidate gates to be validated in each valid merge, we consider its
_transitive fanin_ (TFI) to explore the maximum QoR. Experimental results show
that the simulator can reduce the runtime by 7.18$\times$ on average, and the
SAT sweeper can reduce the runtime by 35% on average.

##  II Preliminaries

###  II-A Circuit simulation

Circuit simulation involves visiting nodes in topological order and computing
output values using their input values. The _simulation pattern_ is a
collection of Boolean values assigned to each _primary input_ (PI) of a
network. Practically, multiple simulation patterns can be consolidated by
encoding sequences of Boolean values as machine words, rather than as
individual bit. Modern arithmetic logic units are enabling the computation of
32 or 64 patterns for a node within a single CPU instruction. The _simulation
signature_ of a node is an ordered set of values produced at the node under
each simulation pattern. If a set of simulation patterns covers all possible
combinations of value assignment (a requirement necessitating $2^{k}$ patterns
for $k$ PIs), this set is _exhaustive_ and the simulation signatures are also
known as _truth tables_ (TTs) [7].

Simulation can be executed either globally across the entire network or
locally in a small _window_ that called _partial simulation_. $2^{16}$
patterns are already impractical to handle, however, networks typically
feature a greater number of primary inputs, often exceeding 16. To use the
exhaustive set of patterns, simulation must be restricted to a window
encompassing fewer than 16 (typically within the range of 8 to 10) leaf nodes.

###  II-B Semi-Tensor Product of Matrices and Its Logical Reasoning

This subsection gives a brief review of the STP computation of matrices. We
refer the reader to [8, 9] for more details. The real matrices with $m\times
n$ dimensions are represented by $M^{m\times n}$. Consider two matrices $X\in
M^{m\times n}$ and $Y\in M^{p\times q}$, the STP can produce matrices in any
dimension.

######  Definition 1.

Let $X\in M^{m\times n}$ and $Y\in M^{p\times q}$, the STP of $X$ and $Y$,
denoted by $X\ltimes Y$, is defined as

$X\ltimes Y=(X\otimes I_{t/n})\cdot(Y\otimes I_{t/p})$,

where $\cdot$ represents the common matrix product, $I_{n}$ represents the
identity matrix with dimension $n$, $t$ is the least common multiple of $n$
and $p$, and $\otimes$ is Kronecker product of two arbitrary dimensional
matrices [10].

######  Property 1.

The STP of matrices supports matrix swapping. Let $A$ be a matrix with any
dimensions, if $Z_{r}\in M^{1\times t}$ is a row vector, then $A\ltimes
Z_{r}=Z_{r}\ltimes(I_{t}\otimes A)$. In contrast, if $Z_{c}\in M^{t\times 1}$
is a column vector, then $Z_{c}\ltimes A=(I_{t}\otimes A)\ltimes Z_{c}$.

The matrix form of logic formulas can be used to describe logic
representations in general. We refer to the matrix product as the STP in this
paper and omit the symbol “$\ltimes$” hereinafter. First, we denote the set of
Boolean variables $\mathbb{B}$.

$\mathbb{B}:\left\\{True=\begin{bmatrix}1\\\
0\end{bmatrix},False=\begin{bmatrix}0\\\ 1\end{bmatrix}.\right\\}$ (1)

######  Definition 2.

A $M^{2\times 2^{n}}$ matrix is called a logic matrix if all its columns are
elements in $\mathbb{B}$, where logic matrix $M_{\sigma}$ in which columns are
consistent with the TT (it is read from right to left) of a logic operation
$\sigma$ is called the structural matrix.

######  Property 2.

$a,b\in\mathbb{B}$ and $\sigma$ is an any Boolean operator. The structural
matrix of unary operator “ _not_ ” ($\lnot$) is
$M_{\lnot}=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$. The inversion of variable
$a$ can be converted to matrices multiplication as $\bar{a}=M_{\lnot}a$.
Similarly, for binary operators, the logic representation can be converted as
$a\ \sigma\ b=M_{\sigma}ab$.

Therefore, any Boolean function can be converted into its _STP form_ by
structural matrices, and logic identities can be easily proved using structure
matrices of Boolean operators and STP properties.

######  Example 1.

Prove the logic identity $a\rightarrow b=\bar{a}\lor b$ using the STP based
computation.

###### Proof.

According to the Property 2, the STP form of the left hand side is
$M_{\rightarrow}ab$, while the right hand side is $M_{\lor}(M_{\lnot}a)b$.

$M_{\lor}M_{\lnot}=\begin{bmatrix}1&1&1&0\\\
0&0&0&1\end{bmatrix}\begin{bmatrix}0&1\\\
1&0\end{bmatrix}=\begin{bmatrix}1&0&1&1\\\
0&1&0&0\end{bmatrix}=M_{\rightarrow}$.

Hence, the identity holds. ∎

A key aspect of Boolean function manipulation is the _canonical form_ (CF),
since these functions can be functionally equivalently represented in several
different logic realizations. A canonical form is also available for STP.

######  Property 3.

Any logic expression $\Phi(x_{1},\ldots,x_{n})$ with Boolean variables
$x_{1},\ldots,x_{n}\in\mathbb{B}$ can be computed into a canonical form
$M_{\Phi}$ as

$\Phi(x_{1},\ldots,x_{n})=M_{\Phi}x_{1}\ldots x_{n}$,

where $M_{\Phi}\in M^{2\times 2^{n}}$.

We use an example to explain the STP computation process.

######  Example 2.

There are three persons $a$, $b$, and $c$. They are either honest or liar,
suppose a liar always said a lie and the honest man always told the truth.
Person $a$ said that person $b$ is a liar, person $b$ said person $c$ is a
liar, and person $c$ said that both $a$ and $b$ are liars. Who is/are the
liar(s)?

First, we define logic variable $a$ to indicate person $a$ is honest. Thus
$\bar{a}$ means $a$ is a liar. The definitions also work for Boolean variables
$b$ and $c$. The statements result in the logic expression

$\Phi(a,b,c)=(a\leftrightarrow\bar{b})\land(b\leftrightarrow\bar{c})\land(c\leftrightarrow\bar{a}\land\bar{b}).$
(2)

The STP form of (2) is

$\Phi=M_{\land}^{2}(M_{\leftrightarrow}aM_{\lnot}b)(M_{\leftrightarrow}bM_{\lnot}c)(M_{\leftrightarrow}cM_{\land}M_{\lnot}aM_{\lnot}b)$.

Then, converting the STP form of logic expression into the canonical form
$M_{\Phi}$ as

$\Phi(a,b,c)=M_{\Phi}abc=\begin{bmatrix}0&0&0&0&0&1&0&0\\\
1&1&1&1&1&0&1&1\end{bmatrix}abc.$

If we assume the simulation pattern is 010, that is, $b$ is honest, $a$ and
$c$ are liars,

$a=\begin{bmatrix}0\\\ 1\end{bmatrix},b=\begin{bmatrix}1\\\
0\end{bmatrix},c=\begin{bmatrix}0\\\ 1\end{bmatrix}.$

The STP form $\Phi(a,b,c)$ can be computed, i.e., simulated as

$\displaystyle\Phi(a,b,c)$ $\displaystyle=\begin{bmatrix}0&0&0&0&0&1&0&0\\\
1&1&1&1&1&0&1&1\end{bmatrix}\begin{bmatrix}0\\\
1\end{bmatrix}\begin{bmatrix}1\\\ 0\end{bmatrix}\begin{bmatrix}0\\\
1\end{bmatrix}$ $\displaystyle=\begin{bmatrix}0&1&0&0\\\
1&0&1&1\end{bmatrix}\begin{bmatrix}1\\\ 0\end{bmatrix}\begin{bmatrix}0\\\
1\end{bmatrix}$ $\displaystyle=\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\begin{bmatrix}0\\\ 1\end{bmatrix}=\begin{bmatrix}1\\\
0\end{bmatrix}.$

###  II-C SAT-sweeping

SAT-sweeping is used to detect, prove, and merge (or collect) functionally
equivalent nodes (up to complementation). In SAT-sweeping, two nodes are
checked if they can be merged using SAT [11, 12]. SAT solvers provide a CE in
the event the nodes cannot be merged, which is an input assignment that allows
the two gates to be simulated to have different values. A traditional
implementation of SAT-sweeping would test all possible pairs of nodes. In
order to alleviate this problem, simulation is extensively used in SAT-
sweeping in order to reduce the number of calls to the SAT solver. Using
initial random simulations, nodes can be grouped into equivalence classes,
that is, classes of nodes that always simulate to the same value. In this
case, only calls to SAT are required to prove, or disprove, equivalencies
between gates belonging to the same class. As a result, the number of SAT
queries has already been drastically reduced [13].

##  III STP-based circuit simulator

In this section, we propose a STP-based simulation of $k$-LUT networks. In
modern simulator, the key aspect of simulating $k$-LUTs is applying bitwise
logic operations efficiently. A $k$-LUT takes $k$-input bits, and the
simulator applies bitwise logic operations to its input signals, simulates the
behavior by a predefined TT of LUTs. However, these bitwise operations (AND,
OR, XOR, and NOT) can not provide efficient support for $k$-LUT networks,
making simulation slower. When it comes to STP, any Boolean function can be
easily converted into its $k$-LUT network (bitwise operation is 2-LUT)
represented by matrices, and the simulation is actually matrix multiplication.

###  III-A Circuit Simulation Algorithm

The STP-based simulator can simulate all nodes or some specified nodes. The
former visits all nodes in topological order, while the latter use the cut
algorithm in Section III-B to map these nodes which do not simulated into
$k$-LUT.

The algorithm is shown in Algorithm 1. The input of the algorithm is a $k$-LUT
network _K_ , a simulation pattern set _P_ and the simulation mode _m_. The
simulation mode contains all node simulation ($a$) or specified node
simulation ($s$). The output is the obtained simulation signature _S_
according to the choice of _m_. If the chosen mode is $a$, the simulator will
visit all nodes in the network in a topological order and use their input
value to compute the output value by matrices multiplication(line 2).
Otherwise, when we only need to get the simulation signature of specified
nodes, we perform the simulation by following the steps below. First, we
compute the size _limit_ of a cut based on the number of simulation patterns.
Because it also takes time to compute the TT of cuts, this ensures that the
STP method is more efficient than direct simulation(line 4). Second, we take
nodes in $s$ and _limit_ as the boundary to cut _K_ , so that each cut is a
tree structure with leaf node no larger than _limit_ , and the root node of
each cut is stored in _root_ set(line 5). Then, we use the STP-based matrices
multiplication to compute the TT of all cuts in the _root_ set(line 6).
Finally, each cut node in the root collection is accessed in topological order
and its output is computed based on the input value, as shown in Example 2
(line 7). The algorithm returns the desired simulation signatures.

Input: logic network _K_ , simulation pattern _P_ , _m_(all nodes _a_ or
specified nodes _s_)

Output: simulation signature(_S_)

1if  _m == $a$_ then

2 $S\leftarrow$ sim_all_nodes($K,P$);

3 return $S_{a}$;

4

5else

6 $n\leftarrow P$.size(), _limit_ =_log(n)_ ;

7 $root\leftarrow$ circuit_cut($K,limit,s$);

8 STP_matrices_multiplication($K,root$);

9 $S\leftarrow$sim_nodes($K,P$);

10 return $S_{s}$;

11

12 end if

Algorithm 1 STP-based circuit simulation

###  III-B Cut Algorithm

The conventional method of circuit simulation involves a traversal of all
nodes in the circuit, following a topological order, and computing their
output values based on input values. However, this method becomes redundant
when our objective is to obtain the simulation signature for specific nodes,
rather than all nodes. To address this issue, we introduce a novel cut
algorithm designed to minimize the number of intermediate nodes that require
simulation. This optimization reduces the scale of the simulation pattern
(including only the PI of specific nodes), enables us to swiftly obtain the
simulation signatures of these nodes using exhaustive patterns.

Consider a scenario where a node possesses multiple fan-out connections,
totaling $n$ in number. In the context of simulating all nodes, this specific
node is accessed a total of $n+1$ times: once for computing its output and $n$
additional times for extracting its value. By mapping these multiple fan-out
nodes into a $k$-LUT, we effectively reduce unnecessary accesses.
Consequently, we employ the STP to compute the TTs for each defined cut,
employing nodes requiring simulation as the boundary. This strategic approach
allows the simulator to efficiently acquire the desired node’s value while
incurring minimal computational overhead.

###  III-C Example

Here we get a DAG which havs five PIs ($1,2,3,4,5$) and two POs ($po1,po2$),
as shown in Fig. 1(a). There are six intermediate nodes in the circuit, and
each intermediate node records its TT. For example, the node “6” has two
inputs “1” and “3” and one output “10”, where the TT “0111” indicates that
when the two inputs are assigned values according to the sequence,
respectively “11”, “10”, “01”, “00”, the corresponding values of the output
are “0”, “1”, “1”, “1”, that is, 2-input NAND.

Assume that there are 10 simulation patterns as

$01110010111010011011111001100000000111111010000101,$

where the combination formed by taking out the $i$-th bit of each input
represents the $i$-th simulation pattern. For example, the first simulation
pattern is “01100”.

(a) example logic circuit

(b) cutting algorithm

Figure 1: Illustration of STP-based simulation

If there are two specific nodes, nodes “7” and “8” fall under this category.
We simulate and only get the simulation signature of these two nodes. Firstly,
we need to find the cut that satisfies the leaf node restriction. Since
$3<log_{2}^{10}<4$, we set a limit of 3, meaning that each cut’s input count
should not exceed 3. After cutting, we derive four distinct cuts: (6,10), (7),
(8), and (9,11). For each cut, we proceed to compute its TT, as illustrated in
Fig. 1(b). Finally, all nodes of the root set are sorted according to the
inverse topology order of the graph. We sequentially traverse and simulate
each node within this sorted root set, ultimately leading to the acquisition
of simulation signatures. As an example, nodes “7” and “8” contain three PIs,
namely “2”, “3”, and “4”. The scale of exhaustive patterns of node “7” and “8”
is $2^{2}=4$, $2^{3}=8$, that is less than original 10 patterns. By combining
exhaustive patterns of PIs with the TT of node “7” (0111) and “8” (0111), the
TT of node “7” and “8” can be simulated as follows: 7: 1110, 8: 11110001.

##  IV STP-based SAT-sweeping framework

Input: Network $N$, number $n$

Output: Optimized network $N^{\prime}$

1$n\leftarrow 1000$ ;

2 $S_{e},S_{c}\leftarrow$ SAT_guided_simulation_patterns($N$);

3 _class_ $\leftarrow$ constant_prop($N,S_{c}$), init_equiv_class($N,S_{e}$);

4 _list_ $\leftarrow$ inverse_topo_sort($N$);

5 foreach  _gate $G_{i}\in$ _list_ _ do

6 candidate = $G_{i}$;

7 if  _skip(candidate)_ then

8 continue ;

9

10 end if

11  _class_new_ = _class_($G_{i}$) $\bigcup$ (INV + _class_($\bar{G_{i}}$));

12 sort_topo_order(_class_new_) ;

13 foreach  _gate $G_{j}\in$ _class_new__ do

14 foreach  _$G_{k}\in$ transitive_fanin($G_{j},n$)_ do

15 driver = $G_{k}$;

16 if  _skip(driver, candidate)_ then

17 continue ;

18

19 end if

20 eq = SAT(candidate $\oplus$ driver);

21 if  _eq == unDET_ then

22 mark_dont_touch(candidate);

23 break ;

24

25 end if

26 if  _eq == unSAT_ then

27 $N^{\prime}\leftarrow$ substitute_node($G_{j}$, candidate);

28

29 else

30 CE $\leftarrow$ SAT;

31 STP_simulation($N$);

32 refine_equiv_class(CE, class_new, $N$);

33

34 end if

35

36 end foreach

37

38 end foreach

39

40 end foreach

41return $N^{\prime}$;

Algorithm 2 STP-based SAT-sweeping algorithm

In this section, we integrate our STP-based simulator into the SAT-sweeping
framework, as depicted in Fig. 2, and the algorithm is delineated in Algorithm
2. In essence, our STP-based simulator accomplishes exhaustive simulation by
employing circuit cut algorithm and focusing exclusively on simulating nodes
within equivalence classes. For the SAT solver, we utilize a circuit-based SAT
solver to direct access to the network [14].

As for details, we set the upper limit for the number of nodes that can be
compared within the TFI to 1000 (line 1). Simulation patterns are initially
generated using SAT-guided initial pattern algorithms (line 2). These high-
quality simulation patterns serve as the foundation for computing equivalence
classes, accounting for complementation, and executing the propagation of
constant nodes for substitution (line 3). Next, we arrange the list of gates
that require processing in a reverse topological order, effectively traversing
the circuit from primary outputs to primary inputs (line 4). Subsequently, we
address the sorted gates in sequence. The gate presently under consideration
is denoted as _candidate_ , as it is a potential candidate for removal and
replacement with a preceding gate in the topological order (line 6). To
expedite the process, we promptly examine whether the current candidate should
be skipped by checking for _don’t touch_ conditions (lines 7-9). We treat
equivalence classes for both positive and negative polarity as a single class,
subsequently organizing this general class topologically (lines 10-11). To
maximize the QoR, we consider the TFI cones of each candidate to identify
opportunities for valid merges (lines 12-13). As each member of the
generalized equivalence class is attempted for a merge, it is called _driver_
(line 14). Furthermore, it is essential to verify the conditions of _drivers_
(lines 15-17). The equivalence problem is translated into _Conjunctive Normal
Form_ (CNF) and submitted for resolution to the SAT solver (line 18). In the
event that the solution is _unDET_ , signifying an undetermined outcome, the
_candidate_ is marked as _don’t touch_. This designation greatly enhances
runtime speed and scalability. Subsequently, when we receive the _unSAT_ , we
proceed with node substitution by connecting the fanins of the _candidate_
fanouts to the driver. Upon receiving the _SAT_ , we retrieve the CE from the
solver, employ the STP-based simulator to propagate the CE, and subsequently
refine the equivalence classes based on this information. At the conclusion of
this procedure, any dead nodes are eliminated from consideration.

###  IV-A Refinement using STP-based Simulation

A two-round SAT-guided simulation is first employed to perform tasks such as
engine allocation, computation of candidate equivalence classes, and the
substitution of constant nodes. The first round simulation is to ensure that
gates exhibit simulation signatures characterized by either all zeros or all
ones, which helps reduce network complexity through efficient constant
propagation. Additionally, the second round simulation aims to avoid gates
with only a few ones and the rest zeros, that is, simulation signatures with a
high toggle rate111the toggle rate is the ratio of bit-toggles over the bit-
string length.. The constraints on the characteristics of simulation patterns
can be efficiently formulated as a SAT problem. When the SAT problem is
satisfied, new assignments will be generated at the PIs which satisfy the set
of constraints [6].

Figure 2: The proposed ecosystem.

CE simulation has been employed for the purpose of refining candidate
equivalence classes. However, due to the impractical time demands of
simulations, they often rely on partial simulation patterns. To further
enhance the refinement of candidate equivalence classes, STP-based simulation
exclusively visits nodes belonging to the same equivalence class. For
instance, if an equivalence class comprises only two nodes, “A” and “B”, we
map the remaining nodes into $k$-LUTs, and then simulate the values of “A” and
“B” to disprove their equivalence. Furthermore, during the simulation of each
CE, we first convert nodes not within equivalence classes into $k$-LUTs, and
then simulate candidate nodes to refine the equivalence classes. In contrast
to traditional simulations, STP-based simulations adopt exhaustive simulation
patterns, if the simulation is restricted to a window encompassing fewer than
16 leaf nodes. In terms of practical implementation, the ID of the equivalence
class is stored as an integer array within the _equivalence class manager_.
This information is utilized by the STP-based simulator to exclusively
simulate nodes belonging to the same equivalence class, while the remainder
are mapped to $k$-LUTs.

##  V Experimental results

This section shows the effectiveness of the proposed circuit simulator and
corresponding SAT-sweeping algorithm. The proposed simulator is implemented in
C++ on top of the logic synthesis framework ALSO222Chu Z. ALSO: Advanced logic
synthesis and optimization tool. https://github.com/nbulsi/also, 2022., in
which the source codes are publicly available. All experiments are performed
on a 3.20 GHz Apple M1 CPU with 8GB of main memory.

###  V-A Simulation

First, we evaluate the efficiency of our proposed STP-based simulator on the
EPFL benchmarks suite333EPFL benchmark suite,
https://github.com/lsils/benchmarks by command ‘simulator’ in ALSO. We compare
our proposed simulator (_STP_) with the logic circuit simulation (_Mockturtle_
[15]). Each benchmark is simulated with randomly generated $10^{6}$ simulation
patterns.

TABLE I: Circuit simulation results for EPFL benchmarks.

Benchmark | _Mockturtle_ | _STP_
---|---|---
$T_{A}$(s) |  $T_{L}$(s) |  $T_{A}$(s) | x |  $T_{L}$(s) | x
adder | 1.33 | 12.56 | 1.36 | 0.98 | 0.57 | 22.04
bar | 3.65 | 17.61 | 2.87 | 1.27 | 1.65 | 10.67
div | 65.12 | 414.92 | 69.04 | 0.94 | 56.30 | 7.37
hyp | 236.61 | 1,066.65 | 291.85 | 0.81 | 134.90 | 7.91
log2 | 35.30 | 181.67 | 36.01 | 0.98 | 45.17 | 4.02
max | 3.58 | 33.96 | 2.94 | 1.22 | 5.05 | 6.73
multiplier | 29.89 | 149.81 | 31.20 | 0.96 | 23.72 | 6.32
sin | 5.89 | 31.37 | 6.16 | 0.96 | 7.01 | 4.47
sqrt | 28.21 | 154.44 | 22.03 | 1.28 | 20.68 | 7.47
square | 20.42 | 88.91 | 24.26 | 0.84 | 11.34 | 7.84
artbiter | 12.91 | 67.43 | 14.99 | 0.86 | 21.02 | 3.21
cavlc | 0.77 | 3.09 | 0.81 | 0.95 | 0.53 | 5.83
ctrl | 0.21 | 0.77 | 0.17 | 1.24 | 0.09 | 8.58
dec | 0.28 | 4.98 | 0.62 | 0.45 | 0.90 | 5.54
i2c | 1.61 | 11.58 | 1.45 | 1.11 | 1.55 | 7.47
int2float | 0.28 | 1.47 | 0.28 | 1.00 | 0.18 | 8.19
mem_ctrl | 52.21 | 312.26 | 49.74 | 1.05 | 54.02 | 5.78
priority | 1.19 | 10.03 | 0.95 | 1.26 | 1.17 | 8.58
router | 0.35 | 3.02 | 0.33 | 1.06 | 0.41 | 7.36
voter | 16.28 | 92.69 | 15.99 | 1.02 | 8.55 | 10.84
Geo. | 4.99 | 30.44 | 5.04 |  | 4.24 |
Imp. | 1.0 | 1.0 | 0.99 |  | 7.18 |

  * 1

Imp. : Average Geometric Mean Improvement (old/new).

TABLE II: Comparing the number of SAT calls and the runtime of the SAT
sweepers (runtime in seconds)

Benchmark | Statistics | SAT calls | Total SAT calls | Simulation | Total runtime
---|---|---|---|---|---
PI/PO | Lev | Gate | Result | &fraig | STP | &fraig | STP | &fraig | STP | &fraig | STP | x
6s100 | 127138/97599 | 79 | 636,637 | 627,550 | 2,400 | 160 | 10,085 | 7,845 | 0.30 | 0.63 | 1.73 | 1.35 | 0.78
6s20 | 250/202 | 2,828 | 30,251 | 12,741 | 791 | 45 | 5,613 | 4,867 | 0.08 | 0.20 | 33.35 | 28.92 | 0.87
6s203b41 | 80192/68958 | 65 | 474,322 | 463,531 | 213 | 15 | 5,535 | 5,337 | 0.22 | 0.43 | 1.11 | 1.07 | 0.96
6s281b35 | 268334/177236 | 121 | 2,076,248 | 2,058,408 | 10,908 | 1,148 | 17,308 | 7,548 | 2.69 | 3.56 | 6.84 | 3.98 | 0.58
6s342rb122 | 59253/56839 | 52 | 330,130 | 319,365 | 194 | 9 | 3,224 | 3,039 | 0.06 | 0.19 | 0.24 | 0.23 | 0.94
6s350rb46 | 245680/243400 | 194 | 1,550,412 | 1,545,667 | 163 | 82 | 3,479 | 3,398 | 0.75 | 1.21 | 1.72 | 1.68 | 0.98
6s382r | 106395/104831 | 2,752 | 1,756,654 | 1,704,409 | 1,158 | 74 | 5,911 | 4,827 | 1.87 | 4.08 | 86.36 | 70.52 | 0.82
6s392r | 80920/80151 | 538 | 1,599,275 | 1,583,824 | 583 | 51 | 2,878 | 2,346 | 0.40 | 0.54 | 0.80 | 0.65 | 0.82
beemfwt4b1 | 3851/1595 | 1,212 | 47,368 | 41,580 | 1,009 | 222 | 1,412 | 625 | 0.28 | 0.68 | 16.28 | 7.21 | 0.44
beemfwt5b3 | 7370/3047 | 2,235 | 104,771 | 90,611 | 3,303 | 302 | 3,745 | 744 | 1.63 | 2.48 | 37.25 | 7.40 | 0.20
oski15a07b0s | 8640/4454 | 3,822 | 120,268 | 118,728 | 3,812 | 953 | 4,073 | 1,214 | 2.47 | 3.38 | 17.06 | 5.08 | 0.30
oski2b1i | 26489/13254 | 12,340 | 176,605 | 176,553 | 12,073 | 4,829 | 12,125 | 4,881 | 10.90 | 13.82 | 105.72 | 42.56 | 0.40
b18 | 6630/3322 | 64 | 131,100 | 124,530 | 606 | 17 | 3,926 | 3,337 | 0.05 | 0.25 | 0.40 | 0.34 | 0.85
b19 | 13235/6627 | 77 | 256,503 | 242,499 | 1,283 | 88 | 8,359 | 7,164 | 0.14 | 0.28 | 0.65 | 0.56 | 0.86
leon2 | 298888/291880 | 58 | 789,647 | 787,972 | 1,977 | 92 | 2,898 | 1,013 | 1.55 | 3.63 | 5.11 | 4.79 | 0.94
Geo. |  |  | 335,363 | 305,780 | 1,230 | 119 | 4,943 | 2,965 | 0.53 | 1.04 | 4.83 | 3.15 |
Imp. |  |  | 1.0 | 0.91 | 1.0 | 0.09 | 1.0 | 0.60 | 1.00 | 1.99 | 1.0 | 0.65 |

  * 1

Imp. : Average Geometric Mean Improvement (new/old).

The experimental results, as summarized in Table I, provide valuable insights.
To facilitate comparison, we present the mean simulation time for both the
_And-Inverter Graph_ (AIG), denoted as “$T_{A}$”, and the $6$-LUT networks,
denoted as “$T_{L}$”. Additionally, we quantify the extent of CPU time
acceleration as “x”. To gauge algorithm efficiency, we utilize the geometric
mean, denoted as “Geo.”, and the average geometric mean improvement of
optimized simulation over initialization, expressed as “Imp.”. It is evident
that there is a preference for algorithms capable of achieving superior
benchmark simulation times. In terms of “$T_{A}$”, our method demonstrates
runtime performance comparable to that of _Mockturtle_ , which highlights the
effectiveness of our method in simulating common data structures. For
“$T_{L}$”, the _STP_ exhibits an average runtime improvement of 7.18$\times$
(22.04$\times$ maximum). These findings underscore the enhanced performance of
our proposed simulator in the context of $k$-LUT simulation. This advantage
arises because _Mockturtle_ employs incremental simulation, which accelerates
simulation by selectively re-simulating necessary nodes and re-computing only
the last block of TT. However, it currently lacks support for $k$-LUT
simulation. In $k$-LUT simulation, most simulators are limited to extracting
individual bits of the LUT and simulating them separately.

###  V-B SAT-sweeping

We implemented our method to demonstrate the effectiveness of the proposed SAT
sweeper. Our experiments were conducted using a selected subset of benchmarks
from the HWMCC’15 [16] and IWLS’05[17] benchmark suites. The proposed SAT
sweeper is based on existing engine, _& fraig_ (command ‘&fraig -x’ in
ABC444Mishchenko A. ABC: System for sequential logic synthesis and formal
verification. https://github.com/berkeley-abc/abc, 2022). As of our knowledge,
_& fraig_ represents the most efficient and scalable publicly available SAT
sweeper. All results are veriﬁed by ‘&cec’ command in ABC to ensure functional
correctness.

As shown in Table II, we present the number of PIs and POs, internal AND-nodes
in the original AIG (Gate), logic levels (Lev), and the internal AND-nodes
after SAT-sweeping (Result) in the “Statistics” section. Notably, the number
of Result (9% improvement of Gate) remains consistent across both engines
since they start from the same AIG, with the SAT solver conflict limit
disabled in all runs. In “SAT calls” and “Total SAT calls” sections, we detail
the number of satisfiable SAT calls and Total SAT calls, respectively. We
reduced the number of unsatisfiable SAT calls by using fast $k$-LUT based
exhaustive simulations (restricted to a window encompassing fewer than 16 leaf
nodes), thereby reducing the number of satisfiable solver calls by efficient
equivalence class candidates. Section “Simulation” shows the runtime for
simulation, primarily encompasses initial simulation and CE simulation.
Lastly, the “Total runtime” section provides a direct comparison of the SAT
sweeper’s overall runtime, revealing a substantial 35% reduction when compared
to _& fraig_. While _& fraig_ invests runtime resources in high-quality
initial simulation, our method combines this with an exhaustive CE simulation.
We exhibit almost 2$\times$ increase in simulation time, However, this is
offset by a remarkable 91% reduction in the number of SAT calls. This hybrid
approach not only improves the quality of simulation patterns but also
dramatically reduces the number of SAT calls, ultimately resulting in highly
efficient CE simulations.

##  VI Conclusion

In this paper, we introduce a novel STP-based circuit simulation tailored for
SAT-sweeping. Simulating $k$-LUT networks can be particularly challenging for
conventional fast bitwise simulators available off-the-shelf. STP, however,
provides advantages in the computation of logic matrices and circuit
connectivity. The primary contribution of this paper lies in the introduction
of cut algorithm that substantially enhances simulation speed and promotes
efficient information reuse. When compared to the SOTA simulator, our STP-
based simulation demonstrates an average reduction in CPU time by an average
of 7.18$\times$. Moreover, our simulation approach is deeply integrated into
the SAT sweeper, leading to substantial reductions in runtime for refining
equivalence classes. Consequently, the implemented SAT sweeping approach, on
average, performs 35% reduction in runtimes, without compromising quality.

## References

  * [1] Mishchenko A, Brayton R. “Integrating an AIG package, simulator, and SAT solver,” in _Proc.IWLS_ , 2018: 11-16.
  * [2] Mishchenko A, Chatterjee S, Brayton R, et al. “Improvements to combinational equivalence checking,” in _Proc.ICCAD_ , 2006: 836-843.
  * [3] Zhu Q, Kitchen N, Kuehlmann A, et al. “SAT sweeping with local observability don’t-cares,” in _Proc.DAC_ , 2006: 229-234.
  * [4] Zhang H T, Jiang J H R, Amarú L, et al. “Deep integration of circuit simulator and SAT solver,” in _Proc.DAC_ , 2021: 877-882.
  * [5] Cheng D, Qi H, Xue A. “A Survey on Semi-Tensor Product of Matrices,” in _Journal of Systems Science Complexity_ , 2007, 20: 304-322.
  * [6] Amarú L, Marranghello F, Testa E, et al. “SAT-sweeping enhanced for logic synthesis,” in _Proc.DAC_ , 2020: 1-6.
  * [7] Lee SY, Riener H, Mishchenko A, et al. “A simulation-guided paradigm for logic synthesis and verification,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 2021, 41(8): 2573-2586.
  * [8] Cheng D, Qi H, Zhao Y. “An introduction to semi-tensor product of matrices and its applications,” _World Scientific_ , 2012.
  * [9] Pan HY and Chu ZF, “Exact Synthesis Based on Semi-Tensor Product Circuit Solver,” in _Proc.DATE_ , 2023, pp. 1-6, doi: 10.23919/DATE56975.2023.10137287.
  * [10] Van Loan C F. “The ubiquitous Kronecker product,” _Journal of computational and applied mathematics_ , 2000, 123(1-2): 85-100.
  * [11] Mishchenko A, Chatterjee S, Jiang R, et al. “FRAIGs: A unifying representation for logic synthesis and verification,” _ERL Technical Report_ , 2005.
  * [12] Mishchenko A, Zhang J S, Sinha S, et al. “Using simulation and satisfiability to compute flexibilities in Boolean networks,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 2006, 25(5): 743-755.
  * [13] Lu F, Wang L C, Cheng K T, et al. “A circuit SAT solver with signal correlation guided learning,” in _Proc.DATE_ , 2003: 892-897.
  * [14] Zhang H T, Jiang J H R, Mishchenko A. “A circuit-based SAT solver for logic synthesis,” in _Proc.ICCAD_ , 2021: 1-6.
  * [15] Siang-Yun L. mockturtle: a head-only library for logic synthesis and logic optimization, https://github.com/lsils/mockturtle.
  * [16] HWMCC’15 benchmarks https://fmv.jku.at/hwmcc15/.
  * [17] IWLS 2005 benchmarks https://iwls.org/iwls2005/benchmarks.html.


[1]George Zograf

[1]Timur O. Shegai

1]Department of Physics, Chalmers University of Technology, 412 96, Göteborg,
Sweden

2]Faculty of Physics, University of Warsaw, Pasteura 5, 02-093, Warsaw, Poland

# Combining ultrahigh index with exceptional nonlinearity in resonant
transition metal dichalcogenide nanodisks

<EMAIL_ADDRESS>Alexander Yu. Polyakov  Maria Bancerek  Tomasz J.
Antosiewicz  Betül Küçüköz<EMAIL_ADDRESS>[  [

###### Abstract

Second-order nonlinearity in solids gives rise to a plethora of unique
physical phenomena ranging from piezoelectricity and optical rectification to
optical parametric amplification, spontaneous parametric down-conversion, and
the generation of entangled photon pairs. Monolayer transition metal
dichalcogenides (TMDs), such as MoS2, exhibit one of the highest known second-
order nonlinear coefficients. However, the monolayer nature of these materials
prevents the fabrication of resonant objects exclusively from the material
itself, necessitating the use of external structures to achieve optical
enhancement of nonlinear processes. Here, we exploit the 3R phase of a
molybdenum disulfide multilayer for resonant nonlinear nanophotonics. The lack
of inversion symmetry, even in the bulk of the material, provides a
combination of a massive second-order susceptibility, an extremely high and
anisotropic refractive index in the near-infrared region ($n>$ 4.5), and low
absorption losses, making 3R-MoS2 highly attractive for nonlinear
nanophotonics. We demonstrate this by fabricating 3R-MoS2 nanodisks of various
radii, which support resonant anapole states, and observing substantial ($>$
100-fold) enhancement of second-harmonic generation in a single resonant
nanodisk compared to an unpatterned flake of the same thickness. The
enhancement is maximized at the spectral overlap between the anapole state of
the disk and the material resonance of the second-order susceptibility. Our
approach unveils a powerful tool for enhancing the entire spectrum of optical
second-order nonlinear processes in nanostructured van der Waals materials,
thereby paving the way for nonlinear and quantum high-index TMD-nanophotonics.

###### keywords:

High-index nanophotonics, nonlinear optics, second-harmonic generation,
anapole, TMD-nanophotonics

##  1 Main

Nonlinear optics plays a pivotal role in modern science and technology,
particularly through second-order nonlinearities that enable crucial processes
like second-harmonic generation, optical parametric oscillation, and
parametric down-conversion [1, 2, 3]. These nonlinear phenomena require a
break of inversion symmetry in the bulk of the material. Consequently,
promising material platforms for nonlinear optics applications typically
encompass symmetry-broken dielectrics with low optical absorption.
Specifically, the degree of second-order nonlinearity is quantified by the
second-order nonlinear susceptibility, $\chi^{(2)}$, which gives rise to
nonlinear polarization $P^{(2)}_{i}$ in accordance with:

$P^{(2)}_{i}=\sum_{i,j}\varepsilon_{0}\chi^{(2)}_{ijk}E_{j}E_{k},$ (1)

where $E_{j}$ and $E_{k}$ are the electric field components, and
$\varepsilon_{0}$ is the vacuum permittivity.

Figure 1a illustrates a comparative analysis of commonly used nonlinear
materials alongside crystalline silicon ($c$-Si). The latter, although widely
regarded as a prominent nanophotonics material due to its high refractive
index $n\sim 3.5$ in the lossless near-infrared region [6], lacks $\chi^{(2)}$
due to its centrosymmetry [1]. Lithium niobate (LiNbO3, LN) and barium
titanate (BaTiO3, BTO) are among the most extensively employed nonlinear
materials renowned for their performance in second-order nonlinear processes
and integrated photonics applications [7, 8, 9, 10, 11, 12, 13, 14, 15]. The
advantage of these wide-band semiconductors is their exceptionally low
material losses. In the near-infrared region, LiNbO3 and BaTiO3 support
moderate $\chi^{(2)}$ (largest components) of 50 pm/V [16] and 34 pm/V,
respectively [13]. A significant weakness of these two materials is the
relatively low refractive index in the near-infrared, which is slightly
greater than 2, making the optical contrast between the material and the
substrate dramatically low. Such low contrast, although allowing for them to
perform well in low-loss waveguides and high-quality factor resonators (such
as whispering gallery modes, microrings, etc.) for integrated photonics
applications and discrete optical components [17, 18], makes it difficult to
create compact subwavelength resonant nanostructures supporting optical
nonlinear effects.

Having a refractive index similar to that of $c$-Si, both GaAs and AlGaAs
benefit from a significant value of $\chi^{(2)}$, of the order of 340 pm/V
[16] and 290 pm/V [19], respectively (Figure 1a). These materials have already
demonstrated outstanding performance with record-high second-harmonic
generation (SHG) [19, 20, 21] and spontaneous parametric down-conversion
(SPDC) [22, 23]. However, one of the major disadvantages of using GaAs/AlGaAs
is the complexity and high cost of fabrication, which requires molecular-beam
epitaxy (MBE) techniques. Furthermore, the fabrication process can be
potentially intricate due to the requirement of high-index substrates during
the growth phase, followed by the transfer to low-index substrates [19, 24].

Figure 1: 3R-MoS2 material nonlinear optical properties in the 800 – 1040 nm
range. a) Comparison of linear and nonlinear optical properties for several
selected materials – maximum $|\chi^{(2)}|$ component of each material $vs.$
refractive index $n$ in the near-infrared spectral region. Inset shows the
schematic of the optically non-resonant 3R-MoS2 flake studied in (b-d).
Optical data for various materials are taken from: BaTiO3 [13], $c$-Si [6, 1],
GaAs [1, 16], AlGaAs [19], LiNbO3 from [16], 3R-MoS2 from [25] for refractive
index and the DFT data for nonlinear susceptibility shown in b). b) Colored
diamond scatter plot shows a maximum SHG intensity of the optically non-
resonant 9 nm thick 3R-MoS2 flake on a glass substrate as a function of the
pump wavelength $\lambda_{\mathrm{pump}}$ at 0.3 mW average incident power for
every pump wavelength. Inset shows the optical image of the flake. Cyan solid
line shows $|\chi^{(2)}|$ calculated using DFT. Dashed line corresponds to
normalized $|\chi^{(2)}|^{2}$. c) Power-to-power log-log scale plot of SHG
intensity $vs.$ the pump power. Inset shows polar plots of polarization-
resolved SHG obtained at a 910 nm pump for a 9 nm thick 3R-MoS2 flake (purple)
and MoS2 monolayer (black), respectively. d) SHG spectra of the R-MoS2 flake
for different pump wavelengths $\lambda_{\mathrm{pump}}$ from 800 nm to 1040
nm with 0.3 mW power. $\lambda_{\mathrm{SHG}}$ corresponds to
$\lambda_{\mathrm{pump}}/2$. The color code of the lines corresponds to one of
the diamonds in b).

TMDs have recently emerged as a promising high-index nanophotonics platform in
the near-infrared range, with a significant second-order nonlinear
susceptibility in their atomically thin structures [4, 5, 26, 27, 28, 29].
However, their bulk counterparts are unlikely to exhibit substantial
$\chi^{(2)}$ due to the commonly adopted A1B2 stacking configuration. This is
evident from a vanishing SHG from 2H-MoS2 and similar materials multilayers
containing an even number of layers in the few-layer regime and in the bulk of
the material [4, 5]. Nevertheless, the utilization of nanophotonic resonances
in centrosymmetric multilayer TMD nanostructures has shown a remarkable
enhancement in SHG [30, 31, 32]. Furthermore, the SHG intensity of TMD
monolayers placed on external metasurfaces has been significantly boosted [33,
34]. However, this approach does not allow for implementing an all-TMD
nonlinear nanophotonics concept [35, 36, 37, 38, 39, 40]. As a result, there
is a growing demand for alternative material platforms that can combine high
refractive indices with exceptional nonlinearities.

Unlike the 2H counterpart, the 3R phase of MoS2 does not restore inversion
symmetry, allowing it to possess second-order nonlinearity even in the bulk
crystal. Recently it was shown that 3R-MoS2 has one of the most significant
$\chi^{(2)}$ values, which is useful for SHG [41] and piezoelectricity [42]
applications. Furthermore, the low optical absorption in the near-infrared
range [25] makes 3R-MoS2 advantageous for nonlinear optics applications
compared to lossy Weyl semimetals with giant $\chi^{(2)}$ values [43].
Additionally, 3R-MoS2 supports ultrahigh refractive index ($n>4.5$) [25] and
has prospects for epitaxial wafer-scale fabrication [44, 45]. Therefore,
3R-MoS2 is an excellent candidate for nonlinear all-dielectric nanophotonics
(Figure 1a).

Our density functional theory (DFT) calculations predict that $\chi^{(2)}$ of
3R-MoS2 is particularly high. Specifically, solid lines in Figure 1b
correspond to several selected $\chi^{(2)}_{ijk}$ tensor components (see right
$y$-axis). The cyan-colored line depicts the most dominant ones –
$\chi^{(2)}_{xyx},\chi^{(2)}_{yyy},\chi^{(2)}_{yxx}$, which have pronounced
resonant behavior peaking at $\lambda_{\mathrm{pump}}$ of $\sim$ 950 nm and
reaching values as high as $\sim$ 800 pm/V. Therefore, 3R-MoS2 exhibits
superior performance compared to conventional nonlinear materials, as
summarized in Figure 1a, in terms of both second-order nonlinearity and
refractive index. The measured SHG as a function of the pump wavelength for a
9 nm thin, optically non-resonant flake qualitatively agrees with the DFT
calculations and peaks around 910 – 920 nm. Moreover, an even better agreement
is met when the experimental SHG intensity is plotted against
$|\chi^{(2)}_{xyx}|^{2}$, since $I_{\mathrm{SHG}}\propto|\chi^{(2)}|^{2}$. We
attribute the observed bell-shaped SHG spectrum to a material resonance of
nonlinearity in 3R-MoS2. This resonance is likely of excitonic origin, similar
to previously observed SHG in monolayer MoS2 [4] and WSe2 [26].

To verify the second-order nature of the observed nonlinear process, we plot
the power-to-power dependence of SHG at the optimal 910 nm pump. Figure 1c
shows the incident pump power versus SHG intensity in a log-log scale. The
slope of the fitting function is $2.02\pm 0.02$ with a 95% confidence,
manifesting the second-order nonlinear process. The damage threshold for this
flake occurs at laser powers between 50 and 72 mW, which is in agreement with
previous reports [41]. Note that the laser used in this experiment has a
repetition rate of 80 MHz and a pulse duration of $\sim$ 100 fs, implying that
an averaged power of 50 mW corresponds to a peak power of $\sim$
6.25$\times$103 W. The latter is nearly 500 times higher than the damage
threshold of AlGaAs nanodisks on SiO2/ITO [19], making 3R-MoS2 attractive for
high-power nonlinear optics applications.

Furthermore, we performed polarization-resolved SHG measurements on the same
flake (optical image is shown as an inset in Figure 1b). The result displays a
6-fold symmetry profile, as shown in Figure 1c by a purple polar plot. For
comparison, we also provide the SHG data of a MoS2 monolayer, shown in black.
The monolayer signal is more than 10 times weaker than that of the 9 nm thick
3R-MoS2 flake measured under the same conditions. Additionally, we measured
the SHG spectra as a function of excitation wavelength in the 800 – 1040 nm
range at a fixed pump power to demonstrate the precise correspondence between
the pump and the emission wavelength, namely
$\lambda_{\mathrm{SHG}}=\lambda_{\mathrm{pump}}/2$, Figure 1d. The color codes
used in Figure 1d are the same as in Figure 1b. Remarkably, the SHG spectra
are intense and lack any residual or background signal, allowing us to use an
alternative SHG detection scheme – an avalanche photodiode (APD) – which is
advantageous for extended SHG mapping experiments discussed below.

Optically-resonant nanostructures. We successfully optimized the optical pump
conditions to achieve the maximum SHG intensity in optically non-resonant
3R-MoS2 flakes, relying solely on the resonant properties of the material. In
addition, the DFT calculations were performed using the bulk model of 3R-MoS2,
which excluded any consideration of optical resonances. Our current objective
is to go beyond this and engineer the resonant optical response, in
combination with the material properties, to further enhance the efficiency of
SHG. The core concept of this work is to ensure that the engineered resonant
optical response of the nanostructures aligns spectrally with the material
resonances of the 3R-MoS2, creating an effective overlap between them. Indeed,
the SHG can be dramatically enhanced by using high-$Q$ optical resonators. If
the resonance is tuned to the pump wavelength, the SHG output power scales as
$P^{2\omega}\propto\left[Q_{1}P^{\omega}\right]^{2}$ [19, 21, 32], where
$P^{\omega}$ is the pump power, $\omega$ is the pump photon frequency, and
$Q_{1}$ – the quality factor of the resonance at the pump wavelength. Here, we
do not account for the optical mode at the SHG wavelength due to the high
optical losses of 3R-MoS2 in that range.

We proceed with the design and fabrication of 3R-MoS2 subwavelength optical
resonators. Due to challenges in fabricating spherical nanoresonators capable
of hosting optical Mie modes [46, 47], TMD nanodisks are preferred for their
highest possible symmetry [35, 48]. Specifically, TMD nanodisks are
insensitive to incident linear polarization and can support optical modes
across visible and near-infrared ranges. Therefore, in this work, we focus on
fabricating and studying optically resonant subwavelength nanodisks. These
nanodisks can support anapole states in the 880 – 960 nm region, which enables
efficient spectral overlap with the material-specific $\chi^{(2)}$ resonance.
The anapole states can substantially enhance the SHG signal by capturing
optical excitations within the disk instead of radiating them out to the far
field [30, 49, 20, 50].

Figure 2: Linear optical properties of 3R-MoS2 nanodisks. a) Numerical
calculation of $|E|^{2}$, characterizing the stored electromagnetic energy
inside the 3R-MoS2 nanodisk, in normalized units. A1, A2, A3 – anapole-like
states. Inset shows the schematics of the calculation geometry. D1 and D2
label nanodisks with resonantly enhanced stored electromagnetic energy. b)
Numerical calculation of the elastic light scattering spectra by individual
3R-MoS2 nanodisks as a function of the nanodisk’s radius in dark-field
geometry using conditions similar to experimental ones. c) Lilac
($R_{\mathrm{disk}}$ = 350 nm) and mint (190 nm) lines correspond to D2 and D1
numerical (dashed) and experimental (solid) dark-field scatterings. Inset
shows SEM images of the two selected nanodisks. The scale bar is 500 nm. Side
panels show the corresponding near-field distributions calculated at a 920 nm
pump. d) The color map shows experimental dark-field scattering spectra for
3R-MoS2 nanodisks in a logarithmic scale. Each scattering spectrum obtained
for a particular radius $R_{\mathrm{disk}}$ is normalized to its maximum
value. The white dashed lines indicate the dispersion of anapole states with
the nanodisk’s radius.

We start with a numerical design of an appropriate nanodisk geometry. Our
analysis involves studying resonantly enhanced electromagnetic energy stored
within the volume of the nanodisk, namely $u_{\mathrm{NP}}\propto|E|^{2}$,
where $|E|^{2}=EE^{*}$, and $E^{*}$ stands for complex conjugation of the
electric field. Figure 2a shows the colormap of normalized $|E|^{2}$ in the
volume of an individual 3R-MoS2 nanodisk of 65 nm height and a range of radii
$R_{\mathrm{disk}}$ from 130 nm to 500 nm under normal plane-wave excitation
from the bottom (glass substrate) side. Such excitation geometry is chosen
because it resembles the experimental SHG setup (discussed below), which
allows for improved scattering signal collection due to a higher $k$-vector in
the glass substrate. The choice of a 3R-MoS2 nanodisk height of $h=65$ nm is
based on our numerical analysis, which indicates that the optimal SHG response
occurs within the range of heights $h\approx$ 60 – 90 nm for a free-standing
nanodisk. Furthermore, our numerical analysis allows for using directional
excitation and specific apertures for collecting the scattered light. One can
see that the nanodisks host three resonant features (denoted as A1, A2, A3)
that show enhanced $|E|^{2}$ in the spectral region of the highest nonlinear
response. The region below 700 nm demonstrates no resonant features due to
high optical losses. We mark the optimal disks with points D1 and D2,
corresponding to the highest stored electromagnetic energy at the wavelength
of the material $\chi^{(2)}$ resonance ($\sim$ 920 nm). These disks are
promising for subsequent SHG experiments and their theoretically predicted
parameters are $R_{\mathrm{disk}}$ = 190 nm and 350 nm, respectively.

In order to assign the resonant features A1, A2, A3 to optical anapoles, we
first calculate the scattering response of the nanodisks. Figure 2b shows a
colormap of the calculated scattering cross-section for various 3R-MoS2
nanodisk sizes. Note that in this case, the light is incident on the sample
from the top, to mimic the experimental dark-field (DF) scattering setup,
discussed below. The scattering cross-section at each specific nanodisk radius
is normalized to unity to enhance the visual representation of the data.
Indeed, the A1, A2, and A3 features demonstrate anapole-like behavior, which
is manifested in a significant drop in their scattering signal. The scattering
spectra of D1 and D2 nanodisks (mint and lilac color dashed lines) alongside
their calculated electromagnetic near-field maps are shown in Figure 2c. The
near-field maps were evaluated at the minima of the scattering spectra
(labeled as D1 and D2) and show that the electric fields are confined mostly
within the nanodisk volume and exhibit typical anapole profiles. We thus
conclude that within the studied parameter range, the appearance of optical
anapoles in 3R-MoS2 nanodisks is feasible [35].

Based on the theoretical predictions, we fabricated 65 nm high 3R-MoS2
nanodisks with radii in the range of 130 – 500 nm. Such a design results in a
number of resonant features in the 880 – 960 nm spectral range, making it
highly promising for simultaneous enhancement of SHG through the combined
effect of optical and material resonances.

We performed morphological characterization of the fabricated nanodisks by
scanning electron microscopy (SEM) and atomic-force microscopy. Figure 2c
depicts experimental DF scattering of two individual nanodisks with
experimental radii of $\sim$ 225 nm (mint solid) and $\sim$ 380 nm (lilac
solid). These spectra are plotted against the theoretical 190 nm and 350 nm
radii disks (corresponding colors – dashed line). The agreement between theory
and experiment is good, with an important remark that there is a slight offset
between the sizes used in the calculation model and the actual experimental
sizes. This difference may arise due to the disk’s geometry deviating from a
perfect cylinder and having a slight frustum shape (bottom and top radii are
slightly different). The insets in Figure 2c depict the corresponding SEM
images of experimentally fabricated nanodisks, demonstrating this point.
Taking this into account, the experimental DF scattering response is slightly
shifted in comparison to the numerical design, as one can see from Figure 2d.
In Figure 2d, we highlight the anapole-like dips A1-3 observed in experimental
DF scattering spectra and their dispersion with the size of the nanodisks.

Figure 3: Second-harmonic generation in 3R-MoS2 nanodisks. a) Schematic of the
experiment: an individual 3R-MoS2 nanodisk is excited with a near-infrared
femtosecond laser and efficiently converts the pump light into SHG due to a
combination of the anapole state hosted by the nanodisk and a substantial
$\chi^{(2)}$ nonlinearity. A quarter of the nanodisk has been removed to
visualize the near-field distribution at the anapole state. b) Numerical
calculation of the SHG intensity of 3R-MoS2 nanodisks with varying disk radii
on a glass substrate at different pump wavelengths. The SHG intensity is
obtained as total radiation generated at
$\lambda_{\mathrm{SHG}}$=$\lambda_{\mathrm{pump}}/2$. c) Experimentally
measured SHG color map of individual 3R-MoS2 nanodisks on a glass substrate at
different pump wavelengths. The intensity of the SHG signal is obtained at
$\lambda_{\mathrm{pump}}/2$ for every nanodisk. d) SHG spectra measured at 910
nm pump for an unpatterned 3R-MoS2 flake (red line), resonant nanodisks: D1
(mint) and D2 (lilac), and an arbitrary off-resonant nanodisk D3 (black). Note
that the $y$-axis is shown in a logarithmic scale.

Resonantly-enhanced SHG in individual 3R-MoS2 nanodisks. We now turn our
attention to SHG of individual 3R-MoS2 nanodisks. Specifically, we study their
nonlinear optical properties in a configuration schematically shown in Figure
3a. The SHG intensity calculation for nanodisks adopting this excitation-
collection geometry is shown in Figure 3b. A Comsol Multiphysics numerical
model was used in two steps: i) the near-infrared wavelength excitation
$\lambda_{\mathrm{pump}}$ of the nanodisk followed by an extraction of the
near-field data; ii) evaluation of the nonlinear polarization of the disk at
$\lambda_{\mathrm{pump}}/2$ using near-fields from i) weighted by the
$\chi^{(2)}$ tensor, and, finally, SHG at the $\lambda_{\mathrm{pump}}/2$.
Remarkably, the areas of the highest SHG intensity in the spectra appear at
the overlap between the $\chi^{(2)}$ tensor resonance (Figure 1b) and in good
agreement with the A1, A2, and A3 curves predicted in the linear response
(Figure 2a).

We further proceed to SHG measurements. Figure 3c shows SHG maps obtained for
nanodisk radii in the 130 – 500 nm range and pump wavelengths 800 – 1040 nm
(fixed power of $\sim$ 300 $\mu$W). The data was recorded by scanning the
sample with respect to the focused laser spot over the entire nanodisk array
(5 $\mu$m disk-to-disk separation) at a given pump wavelength using a piezo
stage. The nonlinear signal was collected using an APD. As evaluated by the
stored electromagnetic energy and linear DF scattering analysis (Figure 2),
the nanodisks with experimental radii of $\sim$ 225 nm (D1) and $\sim$ 380 nm
(D2) support anapole-like behavior at the wavelength at which we observe
maximum experimental $\chi^{(2)}$ (see Figure 1). These disks, therefore, are
expected to display the highest SHG intensity peaking at $\sim$ 910 – 920 nm
$\lambda_{\mathrm{pump}}$, which matches the condition of the highest
experimental value of $\chi^{(2)}$ for the material. Remarkably, the
experimentally measured anapole-enhanced SHG for a 380 nm nanodisk (D2) at a
910 nm pump provides close to 2 orders ($\sim$ 80-fold) of magnitude increase
in the SHG signal compared to the arbitrary off-resonant nanodisk at this
wavelength (labeled as D3 with $R_{\mathrm{disk}}~{}\sim$ 310 nm). Moreover,
the difference between the simultaneously optically and materially resonant D2
nanodisk pumped at $\lambda_{\mathrm{pump}}=910~{}\mathrm{nm~{}}$ and an off-
resonant (both optically and materially) one can be significantly higher and
approaches $\sim 4$ orders of magnitude (limited by detection resolution,
$\sim$ 5000-fold). An SHG enhancement exceeding 2-orders of magnitude is
reached between the resonant anapole nanodisk and the parent flake from which
the disk was fabricated. However, in the case of the unpatterned flake, the
volume of the emitting SHG material is likely much larger due to the
diffraction limit of the pump beam. Remarkably, the highest SHG intensity
areas spectrally match with the DF scattering anapole dips confirming the
hypothesis that the SHG signal is dramatically enhanced by optical anapole
states. Furthermore, the agreement with the calculations of the SHG response
of the nanodisks on a glass substrate shown in Figure 3b is also remarkable. A
slight mismatch between calculations and experiments is mainly due to two
reasons: i) the DFT-calculated $\chi^{(2)}$ tensor, used for numerical
calculations, peaks around 950 – 960 nm pump wavelength, whereas the
experimental one is around 910 – 920 nm (Figure 1b), and, ii) the SEM images
demonstrate that the fabricated nanodisks slightly deviate from the perfect
cylindrical shape, which affects their optical response.

To further improve our understanding of the SHG enhancement mechanism, we have
performed an explicit head-to-head comparison of the SHG signals from the two
most resonant disks D1, D2, the non-resonant disk D3 (at the same
$\lambda_{\mathrm{pump}}$ = 910 nm wavelength, but without anapole
enhancement) and the unpatterned parent flake. Additionally, we collected the
SHG signal of these disks under off-material resonance conditions under 820
and 1020 nm pumping. These measurements show SHG signals slightly above the
background noise level. Such a comparison allows us to roughly estimate the
relative contributions of material and optical enhancements to the SHG signal.
Specifically, our analysis reveals that the optical anapole enhancement is
responsible for about 80-fold enhancement (from the comparison between D1, D2,
and D3). The material resonance provides an additional 100-fold increase at
the 910 nm pump wavelength when compared to 820 to 1020 nm pumps (from Figure
1b). Jointly, the optical and material resonance grant SHG signal tuning in an
enhancement range approaches 4 orders of magnitude within the near-IR spectral
range ($\sim$ 5000-fold). This is illustrated in Figure 3d, which depicts the
SHG spectra (log scale) obtained using a 910 nm pump wavelength at low pump
power 0.3 mW and just 1 s collection time. The two most resonant disks – D2
(lilac line) and D1 (mint line) possess an SHG signal that is more than 2
orders of magnitude stronger than the unpatterned 3R-MoS2 parent flake (red),
whereas the off-resonant nanodisk D3 is only slightly stronger than the
unpatterned flake. One of the most off-resonant disks in which an SHG signal
can be detected with the spectrometer at the same experimental conditions is
the D1 disk pumped at 820 nm. As a consequence of being driven outside the
anapole state and at the edge of the $\chi^{(2)}$ resonance, its SHG spectrum
at 410 nm peaks slightly above the noise floor and is approximately 1 and 2
orders of magnitude weaker than, respectively, both the unpatterned flake and
optically off-resonant disk D3 when driven at the material resonance (see
Figure 3d).

In conclusion, we demonstrated that the SHG emission can be dramatically
enhanced in bulk $\chi^{(2)}$ 3R-MoS2 material via combining material and
nanophotonic resonances. The material response was optimized by exciting the
resonance of the $\chi^{(2)}$ at $\sim$ 910 nm pump wavelength, whereas
nanophotonic resonances were engineered through appropriate nanodisk
dimensions supporting the anapole state in the wavelength range overlapping
with the material resonance. Such an approach provides more than 3 orders of
magnitude enhancement of SHG signal for the most resonant nanodisks _vs._ the
most off-resonant case ($\sim$ 5000-fold) and more than 2 orders of magnitude
enhancement compared to the unpatterned flake of the same thickness ($\sim$
400-fold). We anticipate that by manipulating the thickness of the 3R-MoS2
flake, one can potentially achieve even higher enhancement of SHG by
engineering nanophotonic resonances. Moreover, we envision that nanostructured
3R-MoS2 holds great potential as a versatile platform for a wide range of all-
TMD nonlinear nanophotonics applications beyond SHG. This includes various
$\chi^{(2)}$ nonlinear processes such as the generation of entangled photon
pairs [23], optical parametric amplification, electro-optical effects, and
more, wherein the resonator and nonlinear medium are integrated within the
same nanostructured object, resonant metasurface, or nanophotonic circuit.

## Acknowledgments

G.Z., A.P., B.K., and T.O.S acknowledge funding from the Swedish Research
Council (VR Miljö project, grant No: 2016-06059 and VR project, grant No:
2017-04545), the Knut and Alice Wallenberg Foundation (grant No: 2019.0140),
Chalmers Area of Advance Nano, 2D-TECH VINNOVA competence center (Ref.
2019-00068) and Olle Engkvist foundation (grant No: 211-0063). M.B. and T.J.A.
acknowledge support from the Polish National Science Center via the project
2019/34/E/ST3/00359. This work was performed in part at Myfab Chalmers and at
the Chalmers Material Analysis Laboratory, CMAL. Calculations were partially
done at the ICM, UW (#GC84-51).

## Conflicts of interests

The authors declare no conflicts of interest.

## References

  *   * [1] Boyd, R. W.  _Nonlinear optics_ (Academic press, Orlando, 2020).
  * [2] Dousse, A. _et al._ Ultrabright source of entangled photon pairs.  _Nature_ 466, 217–220 (2010).
  * [3] Lu, X. _et al._ Chip-integrated visible–telecom entangled photon pair source for quantum communication.  _Nat. Phys._ 15, 373–381 (2019).
  * [4] Kumar, N. _et al._ Second harmonic microscopy of monolayer MoS2.  _Phys. Rev. B_ 87, 161403 (2013).
  * [5] Li, Y. _et al._ Probing symmetry properties of few-layer MoS2 and h-BN by optical second-harmonic generation.  _Nano Lett._ 13, 3329–3333 (2013).
  * [6] Aspnes, D. E. & Studna, A.  Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV.  _Phys. Rev. B_ 27, 985 (1983).
  * [7] Carletti, L. _et al._ Second harmonic generation in monolithic lithium niobate metasurfaces.  _Opt. Express_ 27, 33391–33398 (2019).
  * [8] Solntsev, A. S., Kumar, P., Pertsch, T., Sukhorukov, A. A. & Setzpfandt, F.  LiNbO3 waveguides for integrated spdc spectroscopy.  _APL Photonics_ 3, 021301 (2018).
  * [9] Fedotova, A. _et al._ Second-harmonic generation in resonant nonlinear metasurfaces based on lithium niobate.  _Nano Lett._ 20, 8608–8614 (2020).
  * [10] Okoth, C., Cavanna, A., Santiago-Cruz, T. & Chekhova, M.  Microscale generation of entangled photons without momentum conservation.  _Phys. Rev. Lett._ 123, 263602 (2019).
  * [11] Santiago-Cruz, T. _et al._ Photon pairs from resonant metasurfaces.  _Nano Lett._ 21, 4423–4429 (2021).
  * [12] Pohl, D. _et al._ An integrated broadband spectrometer on thin-film lithium niobate.  _Nat. Photonics_ 14, 24–29 (2020).
  * [13] Karvounis, A., Timpu, F., Vogler-Neuling, V. V., Savo, R. & Grange, R.  Barium titanate nanostructures and thin films for photonics.  _Adv. Opt. Mater._ 8, 2001249 (2020).
  * [14] Timpu, F., Sergeyev, A., Hendricks, N. R. & Grange, R.  Second-harmonic enhancement with Mie resonances in perovskite nanoparticles.  _ACS Photonics_ 4, 76–84 (2017).
  * [15] Savo, R. _et al._ Broadband Mie driven random quasi-phase-matching.  _Nat. Photonics_ 14, 740–747 (2020).
  * [16] Shoji, I., Kondo, T., Kitamoto, A., Shirane, M. & Ito, R.  Absolute scale of second-order nonlinear-optical coefficients.  _J. Opt. Soc. Am. B_ 14, 2268–2294 (1997).
  * [17] Fedotova, A. _et al._ Lithium niobate meta-optics.  _ACS Photonics_ 9, 3745–3763 (2022).
  * [18] Zhang, M. _et al._ Broadband electro-optic frequency comb generation in a lithium niobate microring resonator.  _Nature_ 568, 373–377 (2019).
  * [19] Koshelev, K. _et al._ Subwavelength dielectric resonators for nonlinear nanophotonics.  _Science_ 367, 288–292 (2020).
  * [20] Timofeeva, M. _et al._ Anapoles in free-standing III–V nanodisks enhancing second-harmonic generation.  _Nano Lett._ 18, 3695–3702 (2018).
  * [21] Frizyuk, K.  Second-harmonic generation in dielectric nanoparticles with different symmetries.  _J. Opt. Soc. Am. B_ 36, F32–F37 (2019).
  * [22] Marino, G. _et al._ Spontaneous photon-pair generation from a dielectric nanoantenna.  _Optica_ 6, 1416–1422 (2019).
  * [23] Santiago-Cruz, T. _et al._ Resonant metasurfaces for generating complex quantum states.  _Science_ 377, 991–995 (2022).
  * [24] Camacho-Morales, R. _et al._ Nonlinear generation of vector beams from AlGaAs nanoantennas.  _Nano Lett._ 16, 7191–7197 (2016).
  * [25] Munkhbat, B., Wróbel, P., Antosiewicz, T. J. & Shegai, T. O.  Optical constants of several multilayer transition metal dichalcogenides measured by spectroscopic ellipsometry in the 300–1700 nm range: High index, anisotropy, and hyperbolicity.  _ACS Photonics_ 9, 2398–2407 (2022).
  * [26] Wang, G. _et al._ Giant enhancement of the optical second-harmonic emission of WSe2 monolayers by laser excitation at exciton resonances.  _Phys. Rev. Lett._ 114, 097403 (2015).
  * [27] Trovatello, C. _et al._ Optical parametric amplification by monolayer transition metal dichalcogenides.  _Nat. Photonics_ 15, 6–10 (2021).
  * [28] Khan, A. R. _et al._ Optical harmonic generation in 2D materials.  _Adv. Funct. Mater._ 32, 2105259 (2022).
  * [29] Kumar, P. _et al._ Light–matter coupling in large-area van der Waals superlattices.  _Nat. Nanotechnol._ 17, 182–189 (2022).
  * [30] Busschaert, S. _et al._ Transition metal dichalcogenide resonators for second harmonic signal enhancement.  _ACS Photonics_ 7, 2482–2488 (2020).
  * [31] Nauman, M. _et al._ Tunable unidirectional nonlinear emission from transition-metal-dichalcogenide metasurfaces.  _Nat. Commun._ 12, 5597 (2021).
  * [32] Popkova, A. A. _et al._ Nonlinear exciton-mie coupling in transition metal dichalcogenide nanoresonators.  _Laser Photonics Rev._ 16, 2100604 (2022).
  * [33] Hu, G. _et al._ Coherent steering of nonlinear chiral valley photons with a synthetic Au–WS2 metasurface.  _Nat. Photonics_ 13, 467–472 (2019).
  * [34] Bernhardt, N. _et al._ Quasi-bic resonant enhancement of second-harmonic generation in WS2 monolayers.  _Nano Lett._ 20, 5309–5314 (2020).
  * [35] Verre, R. _et al._ Transition metal dichalcogenide nanodisks as high-index dielectric mie nanoresonators.  _Nat. Nanotechnol._ 14, 679–683 (2019).
  * [36] Munkhbat, B. _et al._ Transition metal dichalcogenide metamaterials with atomic precision.  _Nat. Commun._ 11, 4604 (2020).
  * [37] Ling, H., Li, R. & Davoyan, A. R.  All van der Waals integrated nanophotonics with bulk transition metal dichalcogenides.  _ACS Photonics_ 8, 721–730 (2021).
  * [38] Sung, J. _et al._ Room-temperature continuous-wave indirect-bandgap transition lasing in an ultra-thin WS2 disk.  _Nat. Photonics_ 16, 792–797 (2022).
  * [39] Munkhbat, B., Küçüköz, B., Baranov, D. G., Antosiewicz, T. J. & Shegai, T. O.  Nanostructured transition metal dichalcogenide multilayers for advanced nanophotonics.  _Laser Photonics Rev._ 17, 2200057 (2023).
  * [40] Zotev, P. G. _et al._ Van der Waals materials for applications in nanophotonics.  _Laser Photonics Rev._ n/a, 2200957 (2023).
  * [41] Xu, X. _et al._ Towards compact phase-matched and waveguided nonlinear optics in atomically layered semiconductors.  _Nat. Photonics_ 16, 698–706 (2022).
  * [42] Dong, Y. _et al._ Giant bulk piezophotovoltaic effect in 3R-MoS2.  _Nat. Nanotechnol._ 18, 36–41 (2023).
  * [43] Wu, L. _et al._ Giant anisotropic nonlinear optical response in transition metal monopnictide Weyl semimetals.  _Nat. Phys._ 13, 350–355 (2017).
  * [44] Li, T. _et al._ Epitaxial growth of wafer-scale molybdenum disulfide semiconductor single crystals on sapphire.  _Nat. Nanotechnol._ 16, 1201–1207 (2021).
  * [45] Liu, L. _et al._ Uniform nucleation and epitaxy of bilayer molybdenum disulfide on sapphire.  _Nature_ 605, 69–75 (2022).
  * [46] Tselikov, G. I. _et al._ Transition metal dichalcogenide nanospheres for high-refractive-index nanophotonics and biomedical theranostics.  _Proc. Natl Acad. Sci. U.S.A._ 119, e2208830119 (2022).
  * [47] Chernikov, A. S. _et al._ Tunable optical properties of transition metal dichalcogenide nanoparticles synthesized by femtosecond laser ablation and fragmentation.  _J. Mater. Chem. C_ 11, 3493–3503 (2023).
  * [48] Green, T. D. _et al._ Optical material anisotropy in high-index transition metal dichalcogenide mie nanoresonators.  _Optica_ 7, 680–686 (2020).
  * [49] Baranov, D. G., Verre, R., Karpinski, P. & Kall, M.  Anapole-enhanced intrinsic raman scattering from silicon nanodisks.  _ACS Photonics_ 5, 2730–2736 (2018).
  * [50] Ushkov, A. A. _et al._ Anapole states and scattering deflection effects in anisotropic van der waals nanoparticles.  _Phys. Rev. B_ 106, 195302 (2022).
  * [51] Taghizadeh, A., Thygesen, K. S. & Pedersen, T. G.  Two-Dimensional Materials with Giant Optical Nonlinearities near the Theoretical Upper Limit.  _ACS Nano_ 14, 7155–7167 (2021).