# Deep Learning of ab initio Hessians for Transition State Optimization

Eric C.-Y. Yuan*1,3, Anup Kumar*4, Xingyi Guan1,3, Eric D. Hermes5, Andrew S.
Rosen6,7,
Judit Zádor8, Teresa Head-Gordon1-3, Samuel M. Blau4

###### Abstract

Identifying transition states – saddle points on the potential energy surface
connecting reactant and product minima – is central to predicting kinetic
barriers and understanding chemical reaction mechanisms. In this work, we
train an equivariant neural network potential, NewtonNet, on an ab initio
dataset of thousands of organic reactions from which we derive the analytical
Hessians from the fully differentiable machine learning (ML) model. By
reducing the computational cost by several orders of magnitude relative to the
Density Functional Theory (DFT) ab initio source, we can afford to use the
learned Hessians at every step for the saddle point optimizations. We have
implemented our ML Hessian algorithm in Sella, an open source software package
designed to optimize atomic systems to find saddle point structures, in order
to compare transition state optimization against quasi-Newton Hessian updates
using DFT or the ML model. We show that the full ML Hessian robustly finds the
transition states of 240 unseen organic reactions, even when the quality of
the initial guess structures are degraded, while reducing the number of
optimization steps to convergence by 2–3$\times$ compared to the quasi-Newton
DFT and ML methods. All data generation, NewtonNet model, and ML transition
state finding methods are available in an automated workflow.

1Kenneth S. Pitzer Theory Center and Department of Chemistry,

2Departments of Bioengineering and Chemical and Biomolecular Engineering,

University of California, Berkeley, CA, 94720 USA

3Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,
CA, 94720 USA

4Energy Technologies Area, Lawrence Berkeley National Laboratory, Berkeley,
CA, 94720 USA

5Quantum-Si, Branford, CT, 06405 USA

6Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,
CA, 94720 USA

7Department of Materials Science and Engineering, University of California,
Berkeley, CA 94720, USA 8Combustion Research Facility, Sandia National
Laboratories, Livermore, CA, 94550, USA

*Authors contributed equally

Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Computational identification of transition states (TSs) on the quantum
mechanical potential energy surface (PES) is central to predicting reaction
barriers and understanding chemical reactivity.1, 2 The height of the barrier
exponentially impacts the reaction rate coefficient via the Eyring equation,
and the geometric character of the metastable state is informative about the
kinetic mechanism, making TSs key to describing a broad range of chemical
kinetic outcomes for enzymes, next-generation synthetic catalysts, batteries,
and conformational changes of molecules and materials.

Transition states are first-order saddle points, and locating them via
mathematical optimization is particularly challenging on high dimensional PESs
relevant in complex molecular systems. Locating an equilibrium geometry, i.e.,
a local minimum on the PES, can be found using bracketing methods based on
function evaluations (0th order)3, 4 and with methods that use gradient
information, such as steepest descent5 or conjugate gradient6 (1st order).
Under the quadratic approximation, identification of the transition state can
be more robustly achieved in fewer steps by 2nd order methods7 using the
Hessian matrix, whose elements $H_{ij}$ are defined as the second derivative
of the energy $E$ with respect to atomic positions $R_{i}$ and $R_{j}$.
However, a metastable first-order saddle point is characterized by a single
negative Hessian eigenvalue, and hence 2nd order methods are indispensable to
optimize for molecular TS energies and geometries.1 The Newton-Raphson (NR)
method and its variants, including restricted and augmented methods like the
trust radius method (TRM) and the rational function optimization (RFO) method,
select the displacement vector $\Delta R^{(k)}$ (or in internal coordinates)
at step $k$ using the inverse Hessian $[H^{(k)}]^{-1}$ and the gradient
$g^{(k)}$.8, 9, 10 The geometry direct inversion of the iterative subspaces
(GDIIS) method and its variants similarly utilize Hessians and gradients to
define a search space for optimization.11

However, the evaluation of analytical Hessians for ab initio methods such as
Density Functional Theory (DFT) requires solving coupled-perturbed equations,
which scales one power of system size $N$ higher than the energy or the
gradient and thus can be prohibitively expensive. Consequently, several TS
optimization approaches have been proposed that construct cheaper approximate
Hessians using only gradient information to avoid expensive Hessian
calculations, in general referred to as quasi-Newton (QN) methods.12, 13, 14,
15, 16, 17, 18 Arguably the most widely used Hessian approximation for
minimization is arguably the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method,
where the Hessian is iteratively updated using a rank-2 matrix generated from
the displacements and gradients. Such Hessian updates, however, are positive
definite by design and thus cannot be applied to TS searches. Instead, methods
such as symmetric rank-one (SR1) or Murtagh–Sargent,19 Powell–symmetric-
Broyden (PSB), Murtagh–Sargent–Powell (MSP), Bofill, and TS-BFGS methods are
developed for an indefinite approximate Hessian.20, 21, 22, 23, 24 On a
complex PES, an optimization step can displace the molecule from the preceding
quadratic region such that the resulting updated approximate Hessian quickly
diverges from the true one and thus requires expensive reconstruction.10
Alternatively, double- and single-ended interpolation methods such as nudged
elastic band (NEB),25, 26 quadratic synchronous transit (QST),27 and growing
string method (GSM)28 have been well-established in the field in the past few
decades, and are often combined with QN updates. Despite all these efforts, TS
optimization still requires significant user involvement and relies on trial
and error when robust Hessian information is absent. However, if the full
Hessian is available at every optimization step, concerns regarding the
quality of the initial Hessian and subsequent updates become much less of a
problem when determining a TS.

The recent development of deep learning models for the PES provides a whole
new possibility for acquiring and applying the Hessian in chemically relevant
tasks.29, 30, 31, 32, 33 Intuitively, the power of a fully differentiable
machine learning (ML) force field does not stop at forces or gradients but
also broadly applies to second (and higher order) derivative properties such
as the Hessian matrix $H_{ij}$. In this case, it is possible to calculate
Hessians analytically by automatic differentiation, by finite differences
using gradients from the machine learning model, or by estimating Hessians
using first order information as per the Davidson procedure.1 For example,
such an idea has recently been explored using Gaussian process regression,
where an ML PES was locally trained on semiempirical energies, forces, and
optionally Hessians and used to estimate the updated Hessians.34, 35 Yet, the
high memory demand using kernel-based methods can significantly reduce the
applicability on all but small systems, and the semi-empirical level of theory
can be deficient for reliable chemistry.

In this work, we train an equivariant message-passing neural network (eMPNN),
NewtonNet,31 on an augmented version of the Transition-1X (T1x) dataset,36 a
new benchmark dataset containing $\sim$10 million configurations generated by
the NEB method on $\sim$10 thousand gas-phase organic reactions evaluated with
DFT. Although the full training data is comprised of only energies and
gradients of the molecular configurations, with no Hessians provided, the
whole neural network is fully differentiable such that we can generate the
Hessian $H_{ij}$ through back propagation. We then apply the ML Hessians to TS
optimization on an independent data set of 240 organic reactions previously
proposed by Hermes and co-workers,37, 38 and which are outside of the training
set. We have adapted the Sella code38 to read in full ML Hessians to perform
TS optimizations for these reactions and utilize the same code and
optimization settings in order to compare against QN Hessian optimization with
either ML or DFT.

We find that incorporating explicit Hessians from the NewtonNet ML model into
TS optimization yields a 2–3$\times$ reduction in search steps compared to
approximate Hessian methods, demonstrating a remarkable efficiency improvement
by ensuring higher-confidence search directions that are closer to the optimal
path. The more accurate description of the Hessian also leads to improved
robustness against structural perturbation such that the TS optimization is
less reliant on a good initial guess. With our deep learning model, the
Hessian calculation is over 1000$\times$ faster than the corresponding ab
initio calculation and is consistently more robust in finding TSs than quasi-
Newton methods using the ML or DFT potential energy surface. The combination
of greater efficiency, reduced reliance on good initial guesses, and robust TS
convergence for unseen reactions unveils a new realm of opportunities to
utilize full Hessians for TS optimizations with appropriately constructed data
sets of complex reactive chemistry.

##  2 Results

###  2.1 Machine Learned Prediction of DFT Hessians

Figure 1(a) shows the NewtonNet eMPNN model in which the DFT-computed
molecular energy $E$ is predicted by transforming and aggregating atomic
features $a_{i}$ that accumulate local chemical environmental information from
spatial neighbors $a_{j}$ and interatomic distances $R_{ij}$ through message
passing layers.31 The molecular energy $E$ is then differentiated with respect
to the atomic positions $R_{i}$ to predict atomic forces $F_{i}$ or gradients
$g_{i}$, but of relevance here is that it can be auto-differentiated twice to
obtain $H_{ij}$. We have demonstrated that the energies and forces can be
predicted with excellent accuracy across a whole range of chemistry including
small organic molecules31 as well as for methane and hydrogen combustion, even
with a limited amount of training examples.31, 39

Figure 1: The NewtonNet model and its performance on the ANI-1 and T1-x data
sets. (a) the eMPNN designed for 3D molecular graphs with nodes $\\{Z_{i}\\}$
and edges $\\{R_{ij}\\}$ to predict molecular energies $E$ and atomic forces
$F_{i}$.31 In this work, we further differentiate the network to derive
$H_{i}j$ for TS optimization tasks. (b) The distribution of atomic pairwise
distances, $R$, relative to equilibrium bond distances, $R_{0}$, among
datasets we used for training, where the T1x data set provides more data in
the TS region, and further augmented with the ANI-1x data to add data
corresponding to bond compression. The predicted (c) potential energies and
(d) forces along the reaction coordinate for an unseen reaction for the pre-
trained and fine-tuned model. A comprehensive statistical analysis of energy
and force prediction errors along the reaction coordinates for 1248 unseen
test reactions is summarized in Supplementary Fig. S1 for the pre-trained and
fine-tuned models. Details of the training protocols are described in Methods
and Supplementary Fig. S2 and S3.

Like all ML potentials, the quality of the learned PES and its derivative
properties depends on the availability of relevant training data. Our ML model
is pre-trained on the ANI-1 dataset, which contains more than 20 million off-
equilibrium conformations of small organic molecules up to 8 heavy atoms and
is evaluated with the $\omega$B97X density functional40 and 6-31G* basis
set.41 Figure 1(b) demonstrates that the original ANI-1 dataset is mostly
composed of near-equilibrium geometries41, 42 and that the reaction pathways
are notably undersampled around the metastable states of the reactions.36 As a
result, the pre-trained ML model predicts the energies and forces accurately
(with respect to the underlying DFT data) at the reactant and product states
but fails significantly around the TS (Figure1(c,d)).

Hence, it is fortunate to have the T1x dataset,36 which is a new benchmark for
TS-related ML tasks, containing 9,644,740 molecular configurations generated
by NEB from 10,073 organic reactions, at a level of DFT commensurate with the
ANI-1 data. This new data better represents the entire reaction pathway as
seen in Figure 1(b) and allows us to fine-tune the pre-trained model. As seen
in Figure 1(c,d) when using the fine-tuned model, the reaction barrier is more
smoothly interpolated between NEB images, the false identification of an
energy maximum has been eliminated, and the atomic forces are correctly
predicted to be negligible as expected of a first-order saddle point.

Due to the high cost of the Hessian calculation and storage, the training
datasets we use do not contain ab initio Hessian reference samples. Despite
the lack of such training examples, a prediction of the atomic forces from an
ML model that is continuous and smooth strongly suggests the possibility of
achieving Hessian predictions without explicit training on such tasks. Based
on this assumption, one approach is a finite-difference Hessian estimation
that can be easily realized using the gradients predicted by our model by
stepping along each Cartesian axis. However, an analytical gradient of the
first derivative is more cost-effective than a finite-difference method, and
such a gradient can be performed as long as the neural network is at least
twice differentiable. In this regard, the NewtonNet model is designed using
sigmoid linear unit (SiLU)43 and a polynomial cutoff function44 and is
therefore $C^{2}$ continuous. Utilizing the automatic differentiation in our
neural network, the forces and Hessians can be analytically acquired by back
propagation.

Leveraging the smoothness of the deep learning PES after fine-tuning, Figure 2
shows that reasonably accurate Hessian predictions of the reference DFT model
can be acquired on molecular TSs when compared to the pre-trained model. The
root mean square error (RMSE) of the predictions of eigenvalues is
quantitatively quite good given the range of negative to positive values, and
the eigenvector similarity between the ML and DFT models also improves
significantly after fine-tuning the ML model with additional T1x data. It is
worth noting that the majority of the error arises from the positive
eigenspace with a constant 20% underestimation of Hessian eigenvalues from
DFT. This can be understood in part from Figure 1(b) that the T1x dataset we
use for training is biased toward weaker bond and greater anharmonicity, such
that lower apparent bond strength and force constants will be observed and
modeled. We therefore augmented the T1x dataset by selecting 1,232,469
molecules from the ANI-1x dataset42 that share the same chemical formula with
the T1x dataset but that exhibit compressed chemical bonds, which further
improves accuracy of the Hessian predictions (Figure 1(c)). We note that some
of this improvement may simply be due to the addition of more data, but it
does alleviate some of the underprediction of the positive eigenvalues. Most
importantly, the predicted Hessians by our fine-tuned model using the aug-T1x
data have very accurate leftmost eigenvalue and eigenvector, which are the
most critical ingredients in TS optimization and iterative Hessian
diagonalization.45, 10

Figure 2: Pre-trained and fine-tuned NewtonNet eMPNN performance on Hessian
prediction of the test set. (a) The pre-trained model accurately predicts
Hessians at reactant R and product P minima geometries but fails dramatically
at TSs. (b) The fine-tuned model using the T1x data significantly improves the
accuracy at TSs but with notable underestimation of Hessian eigenvalues. (c)
Augmenting the T1x dataset with compressed bond configurations creates a more
balanced training data and improves the overall performance. More
comprehensive comparisons of the pre-trained and fine-tuned ML prediction
accuracy for Hessians is provided in Supplementary Figs. S4-S8.

###  2.2 Transition State Optimization using Machine Learned Hessians

The fine-tuned NewtonNet model for predicting TS properties is subsequently
employed in practical TS optimization scenarios involving new reactions
independent of the augmented ANI-1X/T1x training and test data. These include
hydrogen migration reactions, endo- and exo-cyclization, generalized Korcek
step 2 reactions, retro-ene reactions, and reverse 1,2 and 1,3 insertions (see
Supplementary Figs. S9–S12); given the training data these involve only
closed-shell molecules. We focus on TS optimization for these unseen reactions
in order to compare a traditional QN method that approximates Hessians using
gradient information from DFT calculations or ML predictions versus a full
explicit ML Hessian used at every step.

We interfaced our fine-tuned NewtonNet model with Sella, a state-of-the-art
open source TS geometry optimizer.38 In Sella, the interconversion between the
Cartesian coordinates and the redundant internal coordinates is automatically
handled, and the Hessian is iteratively diagonalized for the leftmost
eigenvector10 used in the geodesic saddle point optimization.46 In order to
start the TS optimization for the 240 Sella benchmark reactions, we generated
initial guesses with KinBot using reaction templates,37 where each template
also defines the intended reactant and product end states for a given
reaction. We employed restricted step partitioned rational function
optimization (RS-PRFO)47, 48, 49 for the TS optimizations, with dynamically
adjusted step sizes determined by evaluating the confidence of each step (see
Methods for details). After TS optimization, we follow the intrinsic reaction
coordinate (IRC) from the optimized TS structure to find the minimum energy
path that connects the reactant and product; the robustness of the TS
optimization methods is quantified by comparing the intended reactions and the
predicted reactions. The complete list of found transition states of the 240
predicted benchmark reactions is summarized in Supplementary Table S1.

Figure 3: Efficiency improvement using full-Hessian TS optimization compared
to the quasi-Newton approach. (a) The full-Hessian TS optimization requires
50% fewer steps to reach convergence than the QN approximate-Hessian approach,
using the identical NewtonNet force field on the same reactions. (b) The
improved efficiency of the full-Hessian TS optimization comes from both more
confident steps (top) and more direct paths (bottom) to converge. In this
efficiency comparison, exploration steps of initial Hessian construction or
Hessian re-construction for QN restarts have been excluded, whether using DFT
or NewtonNet for gradient calculations.

Figure 3 shows how optimization efficiency is dramatically improved by
providing full explicit Hessians at every optimization step, which is now
affordable relative to DFT as illustrated in Supplementary Figure S13. In
Figure 3(a) we find that the number of steps required to converge to a TS can
be reduced by 2$\times$ of that required by the QN approach using the ML (or
DFT; see Supplementary Figure S14). The trend is notably non-linear, and full-
Hessian optimization is even more advantageous for challenging tasks that
require larger numbers of optimization steps. If the iterative diagonalization
steps for initial Hessian construction and Hessian reconstruction when the QN
approximation breaks down are included, a reduction of close to 3$\times$ of
the required steps is observed when considering these gradient calls. Greater
optimization efficiency can stem from two factors: increased confidence at
each optimization step to increase the step size and a more optimal overall
path of optimization. As shown in Figure 3(b), both factors improve the
efficiency of the TS optimization with full analytical Hessians. In particular
use of the full analytical Hessians exhibits an increase in confidence as
measured by the RS-PRFO method, which allows for increased step size on
average. A shorter, closer to optimal optimization path also plays a smaller
but significant role with the analytical Hessians compared to the QN
approaches whether on the ML or DFT PES.

Figure 4: The quality of optimized TSs using NewtonNet. (a) NewtonNet predicts
reactions that match the intended reactions with higher success rates, both
with and without full Hessians, compared to DFT. The full Hessian also finds
$\sim$10% more TSs which involve chemical reactions as opposed to
conformational changes. (b) The value of the full Hessians over the
approximate QN approaches is apparent when the quality of the initial TS
guesses deteriorates. The QN convergence decays with additional noise to the
guess structures, while the full Hessian convergence is more robust to
perturbations. (c) Comparing whether the left most frequency found on the ML
PES is also a negative frequency on the DFT PES using the ML geometry. (d)
Reoptimizing the ML transition state structure on the DFT surface demonstrates
superior performance for 2-end matches and identifying chemical reactions
compared with starting from the original KinBot initial guess.

Next we consider the robustness of TS optimization using the fully analytical
Hessians versus approximate Hessians (using ML or DFT) by comparing the
intended reactions from KinBot to those predicted from the IRC after TS
optimization. As shown in Figure 4(a), the NewtonNet full or QN Hessian
yielded the intended IRC reactant and product endpoints (2-end match) more
often then the QN Hessian from DFT. In some cases only a 1-end match is found
because the predicted product is more stable, and more chemically plausible,
than the intended product from KinBot. We also characterize the types of TSs
found in Figure 4(a), in which the graphs of the reactants and products are
inequivalent for chemical reactions, whereas isomorphic endpoints indicate a
conformational TS; in either case, both converge to a first-order saddle
point. All KinBot initial guess geometries are intended to yield chemical
reaction TSs, and we see that employing full Hessians in TS optimization
yields $\sim$10% more TSs that involve chemical reactions as opposed to
conformational changes compared to the approximate QN method.

The unreliability of a generated initial guess structure can lead to poor
convergence or inaccurate predictions. Therefore, we consider a measure of
robustness in which the TS optimization must recover from a poor starting
structure, which we analyze by systematically introducing random structural
perturbations to the guess structures generated from KinBot. Figure 4(b)
demonstrates the robustness of NewtonNet’s full-Hessian optimization,
maintaining consistent performance and even exhibiting slight improvements as
noise levels increase up until 2–3 pm, unlike the QN methods. In contrast, the
performance using approximated QN Hessians rapidly decays, even using DFT,
highlighting the importance of accurate Hessians for robust TS optimization.

A final important metric of robustness is whether NewtonNet predicted TS
structures have negative eigenvalues on the DFT PES and/or improve the outcome
of TS optimization on the DFT PES. In Figure 4(c) we compare the vibrational
frequencies of the NewtonNet optimized TS saddle point structures, using both
the full Hessian and the QN Hessian, and use those structures as input to
calculate the frequencies from the DFT Hessian. We then identify the DFT
frequency mode that corresponds to the negative frequency mode from NewtonNet.
For $\sim$96$\%$ of the reactions, NewtonNet predicts highly accurate
imaginary frequencies, regardless of full or QN Hessians. However, there are
10 cases for the QN-ML saddle points, and 11 for the full-ML Hessian saddle
points, which have positive DFT Hessian eigenvalues. We further reoptimized
all 240 reactions using the full Hessian ML TS structures as initial guesses
on the DFT surface (Figure 4(d)). In comparison with the optimization outcomes
when starting with the original KinBot guess structures, we see an increase in
both 2-end matches and chemical transformations under the DFT reoptimization.
Thus, starting from ML optimized TS structures we find overall improved
solutions on the DFT surface.

##  3 Conclusion

We have presented a highly generalizable approach for predicting ab initio
Hessians using machine learning based solely on energy and gradient data, and
only requiring the property of second-order differentiability. Although Denzel
and co-workers found that the feasibility of predicting Hessians using ML
without access to explicit Hessian training data was poor,35 our study shows
that solely relying on energy and force data using a well-trained ML model can
efficiently and accurately predict the Hessian for reactive systems. We
attribute our contrasting conclusions to the sufficiently high-quality PES
obtained through the deep equivariant message-passing neural network and the
mathematical relationship between the potential energy and its derivatives
using a ML model that is $C^{2}$ continuous. To our knowledge, none of the
widely used chemical datasets currently include Hessian information, and it is
unlikely that such datasets will become available in the foreseeable future
due to the high cost of generating ab initio Hessians. Thus, it is good news
that models trained on energies and forces are sufficient to derive meaningful
ML Hessians.

The ability to generate high-quality explicit Hessians by deep learning models
obviates the complexities and assumptions associated with standard TS
optimization approaches, in which 240 new reactions never seen in the training
data are predicted with greater efficiency, accuracy, and robustness compared
to QN ML or DFT. The implementation in Sella to utilize the ML Hessians incurs
minimal additional computational overhead and requires no model retraining.
Using NewtonNet, the Hessians can be calculated at least three orders of
magnitude faster than the DFT calculation while use of the full ML Hessian
takes 2–3$\times$ fewer steps in the TS optimization compared to the QN
approach.

This work emphasizes new TS state search methodology and hence used available
training data from ANI-1, ANI-1X, and T1x that is specialized for applications
involving reactive molecular organic systems but at a low-level of DFT and
basis set quality. Hence, for accurate predictions it would be desirable to
recalculate these data sets at a higher level of theory, likely better density
functionals and certainly larger basis sets, in order to predict quantitative
barriers. Of course, with appropriate new data sets, we can generalize the ML-
Hessian approach for many practical applications for TS optimization in many
areas of chemical and material sciences. This is made possible by the tight
software integration we have developed for NewtonNet with Sella, and workflows
which could be generalized to other relevant systems.

We envision several new areas of TS optimization using the ML approach we have
described here. For example, the TS discovery for a reaction using methods
such as NEB,25, 26 QST,27 and GSM28 often turn to local TS optimization
methods when a reasonable approximate structure has been obtained.10 We showed
that our ML transition state structures do improve the DFT optimizations.
Therefore, a stepwise procedure that integrates mathematical optimization and
machine learning to collaboratively achieve TS discovery for chemical
reactions should be viable. In addition, since we can accurately calculate
vibrational frequencies at the optimized TS structures, molecular free energy
can be efficiently derived using the harmonic approximation. Further, since
the quadratic correction can be applied at each step throughout the
optimization process, it is possible to optimize TS structures on the free
energy surface rather than PES, as is extending the feasibility of variational
transition state theory (VTST) to larger systems where ab initio vibrational
analysis becomes impractical.50

##  4 Methods

###  4.1 Data preparation

The T1x dataset for training is split in two different ways to assess the
accuracy of the model predictions, illustrated in Supplementary Figure S2. The
original splitting in the literature is based on molecular compositions. All
geometries with the same formula (equilibrium and non-equilibrium) are part of
the same training, validation, and test set, leading to minimum data leakage.
The error from this splitting can be regarded as the worst-case uncertainty
estimation of an unseen configuration for the application on the real system
in Figure 1 and 2. On the other hand, we wish to maximize the chemical
knowledge from the dataset learned by our model, so a more conventional
splitting among molecular conformations is devised to learn the PES of all
reactions. Hence in the design of the test set, we ensured that no reaction
had both reactant and product pairs found in the training set.

The reference DFT Hessians are performed using Q-Chem 6.0.0,51 using the
$\omega$B97X functional40 and 6-31G* basis set52 in order to maintain
compatibility with the T1x dataset.36 The eigenvalues are assigned based on
the cosine similarity between the predicted and reference eigenvectors as a
linear sum assignment problem using the Jonker–Volgenant algorithm.53

The performance of the ML-Hessian and ML- and DFT-QN TS optimizations is
evaluated by the Sella benchmark dataset.38 The dataset contains 500 small
organic molecules between 7 and 25 atoms in configurations that approximate
the TS geometries the across reaction families, among which the 240 reactions
are closed shell. We regenerate 240 of them with further refinement, in
contrast to the original data, in that the guess structures in our work are
constrained minima on ab initio PES instead of saddles on semi-empirical
PES.36 We also inject Gaussian noise up to 50 pm directly onto the atomic
positions in the Cartesian coordinate of the initial guesses in order to
degrade them for the purpose of understanding a given methods ability to still
find the TS.

###  4.2 NewtonNet Model and training details

The NewtonNet model with three message passing layers is trained using the
same architecture as described in Reference 31. Each node encodes an atomic
environment into 128 features initialized by atom types $Z_{i}$, and each edge
encodes an interatomic distance $R_{ij}$ in 20 radial basis functions with
polynomial cutoff of 5 Å.54, 44 The node features are equivariantly updated
with messages from neighboring nodes and edges. The molecular energy
$\tilde{E}$ is the sum of all atomic energies $\tilde{E}_{i}$,55

$\tilde{E}=\sum_{i^{\prime}}^{A}\tilde{E}_{i^{\prime}}(\\{Z_{i}\\},\\{R_{ij}\\})$
(1)

where atomic energies $\tilde{E}_{i}$ are predicted from the node features at
the final layer and $A$ is the total number of atoms. The predicted atomic
forces $\tilde{F}_{i}$ are calculated as the first derivative of the molecular
energy $\tilde{E}$ with respect to atomic positions $R_{i}$,56

$\tilde{F}_{i}=-\nabla_{i}\tilde{E}=-\frac{\partial\tilde{E}}{\partial R_{i}}$
(2)

and the predicted atomic Hessians $\tilde{H}_{ij}$ are further calculated as
the second analytical derivatives of the energy,

$\tilde{H}_{ij}=\frac{\partial^{2}\tilde{E}}{\partial R_{j}\partial
R_{i}}=\frac{\partial^{2}\tilde{E}}{\partial R_{i}\partial R_{j}}$ (3)

However, only the energy $\tilde{E}$ and forces $\tilde{F}_{i}$ are trained in
the loss function $\mathscr{L}$,

$\mathscr{L}=\frac{\lambda_{E}}{M}\sum_{m}^{M}(\tilde{E}_{m}-E_{m})^{2}+\frac{\lambda_{F}}{M}\sum_{m}^{M}\frac{1}{3A_{m}}\sum_{i}^{A_{m}}||\tilde{F}_{mi}-F_{mi}||^{2}$
(4)

where $M$ is the total number of molecular graphs, which is 8 million for
training, 1 million for validation, and 1 million for testing.

We use mini-batch gradient descent algorithm with a batch size of 100 to
minimize the loss function using the Adam optimizer57 with initial learning
rate of $10^{-4}$ and decay rate of 0.7 on plateau. Fully connected neural
networks with sigmoid linear unit (SiLU) nonlinearity43 for all functions were
used throughout the message passing layer. The application of smooth
activation functions like SiLU is critical because the network has to be at
least twice differentiable for Hessian calculations. We take $\lambda_{E}=1$
and $\lambda_{F}=20$ in the loss function in Equation 4 to put more emphasis
on forces for derivative properties, and additional L2 regularization of
$10^{-5}$ is applied on all trainable parameters to further smooth out the
potential energy surface. Layer normalization58 on the atomic features at
every message passing layer is applied for the stability of training. An
ensemble of four models is trained on each splitting manner to ensure the
reproducibility and reliability of the prediction. Outlier among the 4
predictions is removed if its absolute difference from the closest number
compared to the difference from farthest number is larger than the 95%
confidence limit of the Dixon Q’s test.39

###  4.3 Transition state optimization

For the transition state calculations, we use the Quantum Accelerator
(QuAcc),59 a Python package for high-throughput computations with an easy-to-
use interface for Atomic Simulation Environment (ASE)60 optimizers. We utilize
Sella38 as the ASE optimizer for TS and intrinsic reaction coordinate (IRC)
calculations. A new feature to Sella is introduced in this work to provide an
external Hessian matrix at each optimization step. Details of the Sella
settings are described in the Supplementary Information and we refer readers
to the Sella paper for more information.10 The comparison of reactants and
products are based on graph isomorphism. Molecular connectivity graphs are
created using Open Babel61 with atom indexing and compared using the VF2
algorithm.62 The optimization path length is calculated in the Cartesian
coordinate with the Kabsch algorithm.63 The path length in Figure 3(b) only
accounts ones with 2-end matches.

##  5 DATA AVAILABILITY

All data 64 including initial transition state guess structures, optimized
transition states, and corresponding reactants and products with their
coordinates of geometry, energy, forces and hessians are available at
https://doi.org/10.6084/m9.figshare.25356616. Source data for Figures 1-4 is
available with this manuscript.

##  6 CODE AVAILABILITY

The codebase is comprised of several publicly available packages and tools
that contribute to the project. Sella38 is publicly accessible at
https://github.com/zadorlab/sella and comes with comprehensive documentation.
NewtonNet65, another integral part of the project, is also publicly available
at https://github.com/THGLab/NewtonNet/tree/v1.0.0. The recipes implemented in
QuAcc for NewtonNet and Q-Chem, utilizing Sella as the ASE optimizer for
transition state and IRC calculations, are publicly accessible and accompanied
by thorough documentation.59 Wrappers to QuAcc and Q-Chem, employed for
running calculations and storing data in the MongoDB database, are available
for NewtonNet at https://github.com/Quantum-Accelerators/quacc. The full
workflow and the analysis scripts66, responsible for generating molecular
graphs, retrieving data from the MongoDB database, and performing graph
isomorphisms to analyze reactions, are available at
https://github.com/THGLab/MLHessian-TSopt/tree/main. This comprehensive
summary provides insights into the availability of the codebase for potential
readers and collaborators.

##  7 ACKNOWLEDGMENTS

S.M.B., A. K., E.C.-Y.Y., and T.H-G. thank the LBL Laboratory Director
Research Development program for the work supporting transition state methods.
E.C.-Y.Y., X.G. and T.H-G. thank the CPIMS program, Office of Science, Office
of Basic Energy Sciences, Chemical Sciences Division of the U.S. Department of
Energy under Contract DE-AC02-05CH11231 for support of the machine learning.
A.S.R. acknowledges support via a Miller Research Fellowship from the Miller
Institute for Basic Research in Science, University of California, Berkeley.
This work used computational resources provided by the National Energy
Research Scientific Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility operated under Contract DE-AC02-05CH11231, and
the Lawrencium computational cluster resource provided by the IT Division at
the Lawrence Berkeley National Laboratory (Supported by the Director, Office
of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy
under Contract No. DE-AC02-05CH11231).

##  8 AUTHOR CONTRIBUTIONS STATEMENT

S.M.B. and T.H.G. designed the project. E.C.-Y.Y. and X.G. carried out the
NewtonNet training and Hessian calculations, E.H., A.K. and S.M.B. interfaced
NewtonNet with the Sella software package, A.K., S.M.B., and A.S.R.
implemented workflows, A.K. and S.M.B. executed TS workflows, J.Z. ran KinBot
to generate initial guess structures and reference endpoints, and E.C.-Y.Y.
and T.H.G. wrote the paper. All authors discussed the results and made
comments and edits to the manuscript.

##  9 COMPETING INTERESTS STATEMENT

The authors declare no competing interests.

## References

  * Davidson 1975 Davidson, E. R. The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. _J. Comp. Phys._ 1975, _17_ , 87–94.
  * Amos and Rice 1989 Amos, R. D.; Rice, J. E. Implementation of analytic derivative methods in quantum chemistry. _Comp. Phys. Rep._ 1989, _10_ , 147–187.
  * Nelder and Mead 1965 Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. _Comput. J._ 1965, _7_ , 308–313.
  * Brent 1973 Brent, R. _Algorithms for Minimization without Derivatives_ ; Prentice-Hall: Englewood Cliffs, NJ, 1973.
  * Debye 1909 Debye, P. Näherungsformeln für die Zylinderfunktionen für große Werte des Arguments und unbeschränkt veränderliche Werte des Index. _Math. Ann._ 1909, _67_ , 535–558.
  * Hestenes and Stiefel 1952 Hestenes, M. R.; Stiefel, E. Methods of conjugate gradients for solving linear systems. _J. Res. Natl. Bureau Stand._ 1952, _49_ , 409–435.
  * Nocedal and Wright 2006 Nocedal, J.; Wright, S. J. _Numerical Optimization_ , 2nd ed.; Springer: New York, NY, USA, 2006.
  * Baker 1986 Baker, J. An algorithm for the location of transition states. _J. Comp. Chem._ 1986, _7_ , 385–395.
  * Schlegel 2011 Schlegel, H. B. Geometry optimization. _WIREs: Comput. Mol. Sci._ 2011, _1_ , 790–809.
  * Hermes et al. 2019 Hermes, E. D.; Sargsyan, K.; Najm, H. N.; Zádor, J. Accelerated Saddle Point Refinement through Full Exploitation of Partial Hessian Diagonalization. _J. Chem. Theo. Comput._ 2019, _15_ , 6536–6549.
  * Császár and Pulay 1984 Császár, P.; Pulay, P. Geometry optimization by direct inversion in the iterative subspace. _J. Mol. Struct._ 1984, _114_ , 31–34.
  * Schlegel 1982 Schlegel, H. B. Optimization of equilibrium geometries and transition structures. _J. Comp. Chem._ 1982, _3_ , 214–218.
  * Schlegel 1984 Schlegel, H. B. Estimating the hessian for gradient-type geometry optimizations. _Theo. Chim. Acta_ 1984, _66_ , 333–340.
  * Schlegel 1987 Schlegel, H. B. _Advances in Chemical Physics_ ; John Wiley & Sons, Ltd, 1987; p 249–286.
  * Fischer and Almlof 1992 Fischer, T. H.; Almlof, J. General methods for geometry and wave function optimization. _J. Phys. Chem._ 1992, _96_ , 9768–9774.
  * Lindh et al. 1995 Lindh, R.; Bernhardsson, A.; Karlstrom, G.; Malmqvist, P.-A. On the use of a Hessian model function in molecular geometry optimizations. _Chem. Phys. Lett._ 1995, _241_ , 423–428.
  * Jensen 2003 Jensen, F. Using force fields methods for locating transition structures. _J. Chem. Phys._ 2003, _119_ , 8804–8808.
  * Chantreau Majerus et al. 2021 Chantreau Majerus, R.; Robertson, C.; Habershon, S. Assessing and rationalizing the performance of Hessian update schemes for reaction path Hamiltonian rate calculations. _J. Chem. Phys._ 2021, _155_ , 204112.
  * Byrd et al. 1996 Byrd, R. H.; Khalfan, H. F.; Schnabel, R. B. Analysis of a Symmetric Rank-One Trust Region Method. _SIAM J. Opt._ 1996, _6_ , 1025–1039.
  * Dennis Jr and Schnabel 1996 Dennis Jr, J. E.; Schnabel, R. B. _Numerical methods for unconstrained optimization and nonlinear equations_ ; SIAM, 1996.
  * Anglada and Bofill 1998 Anglada, J. M.; Bofill, J. M. How good is a Broyden–Fletcher–Goldfarb–Shanno-like update hessian formula to locate transition structures? Specific reformulation of Broyden–Fletcher–Goldfarb–Shanno for optimizing saddle points. _J. Comp. Chem._ 1998, _19_ , 349–362.
  * Fletcher 2000 Fletcher, R. _Practical methods of optimization_ ; John Wiley and Sons, 2000\.
  * Bofill 2003 Bofill, J. M. Remarks on the updated Hessian matrix methods. _Int. J. Quant. Chem._ 2003, _94_ , 324–332.
  * Hratchian and Schlegel 2005 Hratchian, H. P.; Schlegel, H. B. Using Hessian Updating To Increase the Efficiency of a Hessian Based Predictor-Corrector Reaction Path Following Method. _J. Chem. Theo. Comput._ 2005, _1_ , 61–9.
  * Jonsson et al. 1998 Jonsson, H.; Mills, G.; Jacobsen, K. W. _Classical and Quantum Dynamics in Condensed Phase Simulations_ ; World Scientific, 1998; pp 385–404.
  * Henkelman et al. 2000 Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. _J. Chem. Phys._ 2000, _113_ , 9901–9904.
  * Govind et al. 2003 Govind, N.; Petersen, M.; Fitzgerald, G.; King-Smith, D.; Andzelm, J. A generalized synchronous transit method for transition state location. _Computational Materials Science_ 2003, _28_ , 250–258.
  * Peters et al. 2004 Peters, B.; Heyden, A.; Bell, A. T.; Chakraborty, A. A growing string method for determining transition states: Comparison to the nudged elastic band and string methods. _J. Chem. Phys._ 2004, _120_ , 7877–7886.
  * Zhang et al. 2020 Zhang, Y.; Wang, H.; Chen, W.; Zeng, J.; Zhang, L.; Wang, H.; E, W. DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models. _Comp. Phys. Comm._ 2020, _253_ , 107206.
  * Bac et al. 2022 Bac, S.; Patra, A.; Kron, K. J.; Mallikarjun Sharada, S. Recent Advances toward Efficient Calculation of Higher Nuclear Derivatives in Quantum Chemistry. _J. Phys. Chem. A_ 2022, _126_ , 7795–7805.
  * Haghighatlari et al. 2022 Haghighatlari, M.; Li, J.; Guan, X.; Zhang, O.; Das, A.; Stein, C. J.; Heidar-Zadeh, F.; Liu, M.; Head-Gordon, M.; Bertels, L.; Hao, H.; Leven, I.; Head-Gordon, T. NewtonNet: a Newtonian message passing network for deep learning of interatomic potentials and forces. _Dig. Disc._ 2022, _1_ , 333–343.
  * Duan et al. 2023 Duan, C.; Du, Y.; Jia, H.; Kulik, H. J. Accurate transition state generation with an object-aware equivariant elementary reaction diffusion model. _Nature Comp. Sci._ 2023, _3_ , 1045–1055.
  * Kim et al. 2024 Kim, S.; Woo, J.; Kim, W. Y. Diffusion-based generative AI for exploring transition states from 2D molecular graphs. _Nature Comm._ 2024, _15_ , 341.
  * Denzel and Kästner 2018 Denzel, A.; Kästner, J. Gaussian Process Regression for Transition State Search. _J. Chem. Theo. Comput._ 2018, _14_ , 5777–5786.
  * Denzel and Kästner 2020 Denzel, A.; Kästner, J. Hessian Matrix Update Scheme for Transition State Search Based on Gaussian Process Regression. _J. Chem. Theo. Comput._ 2020, _16_ , 5083–5089.
  * Schreiner et al. 2022 Schreiner, M.; Bhowmik, A.; Vegge, T.; Busk, J.; Winther, O. Transition1x - a dataset for building generalizable reactive machine learning potentials. _Scientific Data_ 2022, _9_ , 779.
  * Van de Vijver and Zádor 2020 Van de Vijver, R.; Zádor, J. KinBot: Automated stationary point search on potential energy surfaces. _Comp. Phys. Comm._ 2020, _248_ , 106947.
  * Hermes et al. 2022 Hermes, E. D.; Sargsyan, K.; Najm, H. N.; Zádor, J. Sella, an Open-Source Automation-Friendly Molecular Saddle Point Optimizer. _J. Chem. Theo. Comput._ 2022, _18_ , 6974–6988.
  * Guan et al. 2023 Guan, X.; Heindel, J. P.; Ko, T.; Yang, C.; Head-Gordon, T. Using machine learning to go beyond potential energy surface benchmarking for chemical reactivity. _Nature Comp. Sci._ 2023, _3_ , 965–974.
  * Chai and Head-Gordon 2008 Chai, J.-D.; Head-Gordon, M. Systematic optimization of long-range corrected hybrid density functionals. _J. Chem. Phys._ 2008, _128_ , 084106\.
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. _Chem. Sci._ 2017, _8_ , 3192–3203.
  * Smith et al. 2018 Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. _J. Chem. Phys._ 2018, _148_ , 241733.
  * Elfwing et al. 2018 Elfwing, S.; Uchibe, E.; Doya, K. Sigmoid-weighted linear units for neural network function approximation in reinforcement learning. _Neural Networks_ 2018, _107_ , 3–11.
  * Gasteiger et al. 2022 Gasteiger, J.; Groß, J.; Günnemann, S. Directional Message Passing for Molecular Graphs. 2022, arXiv:2003.03123 [physics, stat].
  * Sharada et al. 2014 Sharada, S. M.; Bell, A. T.; Head-Gordon, M. A finite difference Davidson procedure to sidestep full ab initio hessian calculation: Application to characterization of stationary points and transition state searches. _J. Chem. Phys._ 2014, _140_ , 164115.
  * Hermes et al. 2021 Hermes, E. D.; Sargsyan, K.; Najm, H. N.; Zádor, J. Geometry optimization speedup through a geodesic approach to internal coordinates. _J. Chem. Phys._ 2021, _155_ , 094105.
  * Banerjee et al. 1985 Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. Search for stationary points on surfaces. _J. Phys. Chem._ 1985, _89_ , 52–57.
  * Anglada and Bofill 1997 Anglada, J. M.; Bofill, J. M. A reduced-restricted-quasi-Newton–Raphson method for locating and optimizing energy crossing points between two potential energy surfaces. _J. Comp. Chem._ 1997, _18_ , 992–1003.
  * Besalú and Bofill 1998 Besalú, E.; Bofill, J. M. On the automatic restricted-step rational-function-optimization method. _Theo. Chem. Acc._ 1998, _100_ , 265–274.
  * Bao and Truhlar 2017 Bao, J. L.; Truhlar, D. G. Variational transition state theory: theoretical framework and recent developments. _Chem. Soc. Rev._ 2017, _46_ , 7548–7596.
  * Epifanovsky et al. 2021 Epifanovsky, E. et al. Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package. _J. Chem. Phys._ 2021, _155_ , 084801.
  * Ditchfield et al. 2003 Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules. _J. Chem. Phys._ 2003, _54_ , 724–728.
  * Crouse 2016 Crouse, D. F. On implementing 2D rectangular assignment algorithms. _IEEE Transactions on Aerospace and Electronic Systems_ 2016, _52_ , 1679–1696.
  * Schütt et al. 2021 Schütt, K. T.; Unke, O. T.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. _Proceedings of the 38th International Conference on Machine Learning_ 2021, _139_ , 9377–9388.
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. _Phys. Rev. Lett._ 2007, _98_ , 146401.
  * Chmiela et al. 2018 Chmiela, S.; Sauceda, H. E.; Müller, K.-R.; Tkatchenko, A. Towards exact molecular dynamics simulations with machine-learned force fields. _Nature Comm._ 2018, _9_ , 3887.
  * Kingma and Ba 2017 Kingma, D. P.; Ba, J. Adam: A Method for Stochastic Optimization. 2017, arXiv:1412.6980 [cs].
  * Ba et al. 2016 Ba, J. L.; Kiros, J. R.; Hinton, G. E. Layer Normalization. 2016, arXiv:1607.06450 [cs, stat].
  * Rosen 2023 Rosen, A. quacc – The Quantum Accelerator. 2023; https://github.com/Quantum-Accelerators/quacc.
  * Larsen et al. 2017 Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. _J. Phys.: Cond. Matt._ 2017, _29_ , 273002.
  * O’Boyle et al. 2011 O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. _J. Cheminform._ 2011, _3_ , 33.
  * Foggia et al. 2001 Foggia, P.; Sansone, C.; Vento, M. An Improved Algorithm for Matching Large Graphs. 2001.
  * Kabsch 1976 Kabsch, W. A solution for the best rotation to relate two sets of vectors. _Acta Cryst. Sec. A_ 1976, _32_ , 922–923.
  * Yuan et al. 2024 Yuan, E.; Kumar, A.; Guan, X.; Hermes, E. D.; Rosen, A.; Zador, J.; Blau, S.; Head-Gordon, T. Data for deep learning of Ab Initio Hessians for transition state optimization. 2024; https://doi.org/10.6084/m9.figshare.25356616.v1.
  * Yuan et al. 2024 Yuan, E.; Haghighatlari, M.; Rosen, A. S.; Guan, N. X.; Li, J.; Head-Gordon, T. THGLab/newtonnet: V0.1.0. 2024.
  * Yuan et al. 2024 Yuan, E.; Kumar, A.; Guan, X.; Hermes, E. D.; Rosen, A.; Zador, J.; Blau, S.; Head-Gordon, T. THGLab/mlhessian-tsopt: V1.0.1. 2024; https://zenodo.org/doi/10.5281/zenodo.10799018.


These authors contributed equally to this work. These authors contributed
equally to this work. [1,2]James W McIver []Corresponding author

1]Max Planck Institute for the Structure and Dynamics of Matter, Hamburg,
Germany 2]Department of Physics, Columbia University, New York, NY, USA
3]Institut für Theorie der Statistischen Physik, RWTH Aachen University and
JARA-Fundamentals of Future Information Technology, Aachen, Germany

4]Research Center for Electronic and Optical Materials, National Institute for
Materials Science, Tsukuba, Japan 5]Research Center for Materials
Nanoarchitectonics, National Institute for Materials Science, Tsukuba, Japan
6]Center for Computational Quantum Physics, Simons Foundation Flatiron
Institute, New York, NY, USA 7]CNano-BioSpectroscopy Group, Departamento de
Fisica de Materiales, Universidad del País Vasco, San Sebastián, Spain

8]Institute for Theoretical Physics and Bremen Center for Computational
Materials Science, University of Bremen, Bremen, Germany

9]Laboratoire de Physique de l’Ecole Normale Supérieure, Université PSL, CNRS,
Sorbonne Université, Université Paris-Citéo, Paris, France 10]Institut
Universitaire de France, France

# Cavity electrodynamics of van der Waals heterostructures

Gunda Kipp  Hope M Bretscher  Benedikt Schulte  Dorothee Herrmann  Kateryna
Kusyak  Matthew W Day  Sivasruthi Kesavan  Toru Matsuyama  Xinyu Li  Sara
Maria Langner  Jesse Hagelstein  Felix Sturm  Alexander M Potts  Christian J
Eckhardt  Yunfei Huang  Kenji Watanabe  Takashi Taniguchi  Angel Rubio  Dante
M Kennes  Michael A Sentef  Emmanuel Baudin  Guido Meier  Marios H Michael  *
[  [  [  [  [  [  [  [  [  [

###### Abstract

Van der Waals (vdW) heterostructures host many-body quantum phenomena that can
be tuned in situ using electrostatic gates. These gates are often
microstructured graphite flakes that naturally form plasmonic cavities,
confining light in discrete standing waves of current density due to their
finite size. Their resonances typically lie in the GHz - THz range,
corresponding to the same $\text{\,}\mathrm{\SIUnitSymbolMicro eV}$ \- meV
energy scale characteristic of many quantum effects in the materials they
electrically control. This raises the possibility that built-in cavity modes
could be relevant for shaping the low-energy physics of vdW heterostructures.
However, capturing this light-matter interaction remains elusive as devices
are significantly smaller than the diffraction limit at these wavelengths,
hindering far-field spectroscopic tools. Here, we report on the sub-wavelength
cavity electrodynamics of graphene embedded in a vdW heterostructure plasmonic
microcavity. Using on-chip THz spectroscopy, we observed spectral weight
transfer and an avoided crossing between the graphite cavity and graphene
plasmon modes as the graphene carrier density was tuned, revealing their
ultrastrong coupling. Our findings show that intrinsic cavity modes of
metallic gates can sense and manipulate the low-energy electrodynamics of vdW
heterostructures. This opens a pathway for deeper understanding of emergent
phases in these materials and new functionality through cavity control.

##  1 Main

A vast array of low-energy quantum phenomena have been observed in vdW
heterostructures, constructed by mechanically stacking and patterning
exfoliated two dimensional (2D) materials. Examples include superconductivity
[1, 2], correlated insulating states [3], polaritonic phenomena [4, 5, 6, 7,
8], 2D magnetism [9] and both integer and fractional quantum anomalous Hall
effects [10, 11, 12, 13]. These effects are highly sensitive to their
electromagnetic environment [14] and can be tuned in situ using electrostatic
gates, which control the material’s carrier density and electronic structure.
We show here that these gates, typically patterned from graphite flakes for
ease and improved device quality [15], play an additional role in vdW
heterostructures. Due to their sub-wavelength size, they naturally form
plasmonic cavities. The edges of the sample set up standing waves of current
density, confining light deep in the near-field [16, 17, 18, 19, 20]. The
plasmonic cavity resonances of prototypical graphite gates ($\sim
25\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ \-
$2.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$) coincidentally fall on the same
energy scale of many phenomena in vdW heterostructures ($\sim$0.1 -
$10\text{\,}\mathrm{meV}$). This raises the question: can discrete modes of
light confined in a graphite gate modify or even control the electrodynamics
of a vdW heterostructure?

Engineering the properties of quantum materials using the enhanced light-
matter interaction provided by a cavity has gained traction as a route to new
functionality [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36]. When the cavity coupling strength $g$ approaches the cavity resonance
frequency $\nu_{0}$, collective matter modes, such as phonons, plasmons or
magnons, are transformed into light-matter hybrids, whose formation can
significantly impact the macroscopic properties of the material itself [21,
23]. As the interaction becomes stronger and begins to reach the so-called
ultrastrong coupling regime, fewer photons, or even only photon vacuum
fluctuations, are required to create new thermodynamic ground states [37, 38].
Pioneering experiments have demonstrated an enhancement of ferromagnetism
[39], a change in the critical temperature of a metal to insulator transition
[40] and the modification of quantum Hall states in an empty cavity [41].

Fig. 1: Plasmonic microcavity setup for ultrastrong coupling of collective
modes in 2D materials

Figure 1: a Schematic of a 2D material embedded in a sub-wavelength
($W\ll\lambda_{\text{THz}}$) plasmonic cavity formed by a few-nm graphite
flake and surrounding dielectric environment. A collective mode of a 2D
material (red), such as a graphene plasmon studied in this work, can hybridize
with the plasmonic bare cavity mode (blue). b The frequency of the screened
plasmonic mode in the 2D material (light magenta line) is tuned with carrier
density into resonance with the screened bare cavity mode, $\nu_{0}$ (light
blue line). The cavity and plasmon modes hybridize due to coupling mediated by
unscreened currents between the metal traces and an avoided crossing appears
[42], with an energy splitting of twice the coupling strength, $g$. When
$\eta$, the ratio of $g/\nu_{0}$, exceeds the value of 0.1, the system is in
the ultrastrong interaction regime.

The diverse phenomena found in vdW heterostructures should be particularly
sensitive to the plasmonic cavity modes [43, 44, 45] naturally provided by the
metallic gates. Fluctuations are increased due to their reduced
dimensionality, preempting the formation of long-range order [46, 47].
Additionally, the low energy scale of emergent properties makes it easier to
attain ultrastrong light matter interaction, where qualitatively new
properties could be induced purely by coupling to vacuum fluctuations.
However, developing routes for control requires demonstrating and
understanding the interaction between cavity modes and a micron-sized vdW
heterostructure. Therefore, spectral probes on the meV energy scale are
needed, which are experimentally challenging to generate [48]. Furthermore,
typical device sizes are orders of magnitude smaller than the wavelength of
THz light ($\sim 300\text{\,}\mathrm{\SIUnitSymbolMicro m}$), precluding far-
field measurements.

Here, we report on the cavity electrodynamics of monolayer graphene embedded
in a vdW heterostructure plasmonic microcavity, probed using on-chip THz
spectroscopy. We found that the plasmonic collective oscillations in graphene
hybridized with the plasmonic cavity mode, native to a microstructured
graphite electrostatic gate (see Fig. 1). We observed spectral weight transfer
from the graphite cavity mode to the graphene modes and quantified the
normalized coupling strength, $\eta=g/\nu_{0}$, accessing the evasive
ultrastrong light-matter interaction regime ($\eta>0.1$), where cavity
engineering becomes possible. We identify the coupling mechanism and provide
generalizable design principles for enhancing or minimizing this
hybridization, delivering a platform for either engineering low-energy
dynamics or spectroscopically probing the undisturbed many-body ground states
of quantum materials.

Our findings reveal the general relevance of cavities in vdW heterostructures
for tailoring macroscopic properties. This suggests the necessity for
understanding the role that the built-in cavity modes play in the known ground
state properties of these materials, which could then be applied to develop
deterministic cavity control of electronic phases. Such a route heralds the
possibility of collective quantum phenomena and new functionality, such as
Bose-Einstein condensation of plasmons [49, 50], single photon detection in
the THz [51] or plasmonic quantum networks [52]. By delivering a cavity
platform on a chip that enables deterministic tuning of light-matter
interactions, spectral read-out of coupling strength, and simultaneous DC
transport to capture the electronic properties of an undriven cavity, this
work provides a generalizable, experimental route for investigation and cavity
control of emergent phases in vdW heterostructures.

##  2 Cavity design and spectral readout

The cavities reported here combine gate-tunable vdW heterostructure devices of
typical dimensions ($\sim 10\times$$10\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}$) with ultrafast optoelectronic circuitry [53]. This circuitry
confines light to the near-field, where $(\nu\lambda_{\text{THz}})/c\ll 1$
($c$ is the speed of light), and enables spectral-readout from 0.1 -
$1\text{\,}\mathrm{THz}$ [54, 55, 18, 56, 17, 57, 53, 58, 59, 60]. Each cavity
was constructed using a precision-cut graphite flake, encapsulated in hBN and
placed on a sapphire substrate. The hBN preserves the intrinsic properties,
like high mobility, of the 2D material from the deleterious effect of the
environment [61], and electrically isolates it from the circuitry. These
spectrometers feature in situ referencing circuitry, allowing the measurement
of the time-domain cavity and reference THz pulses on separate transmission
line arms (Fig. 2a-b). The Fourier transforms of the recorded time-domain
signals were then used to calculate the transmission coefficient, which in
turn was used for computing the real and imaginary parts of the graphite
cavity conductivity, plotted in Fig. 2c [42].

Fig. 2: On-chip Thz spectroscopy of plasmonic microcavity modes

Figure 2: a Circuit design for on-chip THz spectroscopy of vdW heterostructure
microcavities: Two identical THz pulses (dark blue) are emitted from a
voltage-biased THz antenna and coupled into two symmetric coplanar strip
transmission line arms. The THz pulse interacts with the graphite cavity (dark
$\rightarrow$ light blue) and is measured using a photoconductive switch
triggered by a delayed probe beam on the left. The signal can be compared to
the reference signal measured on the detector on the right. b Experimentally-
measured THz transients of a typical bare cavity (without an additional 2D
material) and reference as a function of time delay. For geometric parameters,
see [42]. The detected cavity signal is significantly modified as compared to
the reference trace. c Real and imaginary parts of the near-field cavity
conductivity obtained from the data in b.

The experimental bare graphite cavity response in Fig. 2c exhibits a resonance
at $\nu_{0}$ that is well described by a Lorentz oscillator model and results
from the excitation of a 2D plasmonic standing wave, which is predominantly
confined beneath the metal traces. The electric field of the THz pulse
generates a current in the graphite that reflects at discontinuities in the
dielectric environment, particularly at the sample boundaries, forming a
standing wave of current density ($j_{m=0}$ in Fig. 2b, inset). The wavelength
and thus momentum of the confined plasmonic resonance is determined by the
geometry and dielectric environment of the heterostructure [16, 17, 18, 19,
20]. We developed an analytical model to describe the plasmonic modes of the
bare cavity and their conductivity, whose results match with both experiment
and finite element simulations that capture the full electromagnetic field
propagation through the cavity [42, 62]. Thus, the resonance frequency of the
cavity mode can be systematically controlled by the microstructured cavity
design [42].

These cavities and associated circuitry provide multiple key functions. First,
the circuitry allows for the contact-less, direct measurement of the complex
cavity conductivity on the length-scales of typical vdW heterostructure
dimensions and spanning the “THz gap”, impossible to access with far-field
techniques [48]. This allows the extraction of parameters such as linewidths,
mobilities and cavity quality factors of collective modes. Importantly, this
technique enables both the real and imaginary part of the conductivity to be
extracted [63]. This is vital for accurate measurements of parameters such as
the superfluid stiffness, particularly in systems where electronic
correlations shift spectral weight to high frequencies. In such cases, the
applicability of using Kramers-Kronig transformations of the absorption to
extract complex quantities would be limited [64, 65]. The spectroscopically
detected electrodynamics of a vdW heterostructure can additionally be compared
to DC transport measurements on the same device [42]. Furthermore, the
metallic transmission line strips provide sharp electrodynamic boundary
conditions, breaking momentum conservation in a controllable manner that
mediates the coupling between the graphite cavity and collective modes in 2D
materials, as elaborated below.

##  3 Spectrally resolving gate-tunable graphene plasmonic modes

We investigated the behavior of the plasmonic collective modes of hBN-
encapsulated graphene with an electrostatic graphite gate. We designed the
device such that the graphite cavity resonance frequency was intentionally
centered at $\sim 1.13\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$, outside of the
experimental spectral sensitivity range, to facilitate the spectral
characterization of individual graphene plasmonic modes (Fig. 3a). On-chip THz
data taken at $20\text{\,}\mathrm{K}$ show that the transmission decreases as
the sample reflectivity increases with increasing carrier density (Fig. 3b).
Additionally, time-domain oscillations appear, indicating a non-trivial
frequency-dependence of the optoelectronic response.

Figure 3: Probing the spectra of gate-tunable graphene plasmonic modes

Figure 3: a Cavity design for characterizing the plasmonic modes of a graphene
layer ($s=0,1,2$) with a resonance frequency of the bare graphite cavity
($m=0$) mode, $\nu_{0}$, of $\sim
1.13\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$, outside the experimental
spectral sensitivity range (for cavity geometry details see [42]). b
Experimentally measured THz transients of the reference (black line) and
increasing graphene carrier densities (colored lines). c Experimentally
obtained real part of the cavity conductivity spectra. The bare graphite
cavity conductivity (dark blue trace) increases slowly in conductivity with
increasing frequency, and is expected to peak at $\sim$1.13 THz. The resonance
frequencies as a function of graphene carrier density of these three
resonances are shown in the right panel. d Analytical model capturing the
cavity conductivity for the device shown in (a), identifying the three
experimentally detected resonances as the $s=0,1,2$ graphene modes. The
extended frequency range of the theory captures the avoided crossings between
the $s=3,4$ modes to the bare cavity mode (dashed line).

The real (dissipative) part of the cavity conductivity of gate-tunable
monolayer graphene embedded in the microcavity was calculated from the time-
domain data in Fig. 3b [42]. The results are plotted in Fig. 3c. At the lowest
carrier density, close to charge neutrality, the cavity conductivity data grow
slowly with increasing frequency. This is consistent with the expectation that
charge neutral graphene is largely transparent to THz fields. Only the rising
edge of the bare cavity resonance is captured, which peaks outside the
experimentally accessible spectral range.

As the graphene becomes more conductive, three resonances appear on top of
this background, which grow in amplitude and blueshift as the carrier density
increases. This can be explained as the formation of confined plasmonic
standing waves in the graphene ($j_{\text{graphene}}$) [16, 17, 18, 19, 20,
66, 67]. Higher frequency resonances, represented schematically in Fig. 3a,
originate from higher order plasmonic modes. The labelling $s=0,1,2$
correspond to the number of nodes in the current density beneath each coplanar
strip, as determined from finite element electromagnetic field simulations
[42]. The spectra in Fig. 3c were fit with a sum of Lorentzians to extract the
linewidth and resonance frequency of each mode as a function of carrier
density. A representative linewidth for the $s=2$ mode at
$n=$2.85\text{\times}{10}^{12}\text{\,}{\mathrm{cm}}^{-2}$$ corresponds to an
electronic mean free path of $\sim 9.1(0.2)\text{\,}\mathrm{\SIUnitSymbolMicro
m}$, roughly the width of the device [42]. The corresponding quality factor is
$Q\sim{5.7\pm 0.1}$, defined as
$Q\sim\nu_{\textit{s}=2}/\gamma_{\textit{s}=2}$, where $\nu_{\textit{s}=2}$ is
the resonance frequency and $\gamma_{\textit{s}=2}$ is the linewidth. The
scattering rate and quality factor are limited by device dimensions,
illustrating that the graphene in the microstructured cavity exhibited
ballistic transport, thereby achieving the highest possible quality factor for
a device of this geometry in this frequency range.

In an uncoupled system and in the absence of electronic correlations, the
carrier-density dependence of the resonance frequency is determined by the
electronic band dispersion of the sample [68]. This can be monitored by
tracking the fitted resonance frequency as a function of carrier density $n$.
We find that the $s=0$ mode scales as $\nu\propto n^{0.23\pm 0.03}$, in good
agreement with the expectation of a Dirac band structure, where 2D plasmons
scale as $\nu\propto n^{0.25}$ [69, 42]. However, the higher order plasmon
modes, $s=1,2$ scale as $0.12\pm 0.01$ and $0.11\pm 0.01$ respectively.

To understand the suppressed power laws of the $s=1,2$ modes, we extended our
analytical model to describe the light-matter interaction of a cavity
containing a 2D material [42]. Fig. 3d plots the results of the real part of
the simulated cavity conductivity up to
$1.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$, using similar linewidths as
extracted in experiment. Three modes appear below
$1\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ with excellent qualitative
agreement to the experimental data. Examining the simulated conductivity
outside of the experimental range shows that when the graphene $s=3,4$ modes
cross with the bare cavity mode (dashed line, Fig. 3d), they hybridize and
form avoided crossings, splitting into upper and lower polariton branches. The
simulations reveal that at high carrier densities, the observed $s=1,2$ modes
also form lower branches of avoided crossings occurring at carrier densities
outside of the experimentally accessible range [42]. The experimentally
detected, suppressed power-law behaviors of the graphene $s=1,2$ modes are
thus a consequence of the hybridization of graphene and cavity plasmons.

##  4 Ultrastrong cavity-coupling of plasmonic collective modes

Having demonstrated the role of hybridization in a device where resonances
were well-separated in frequency, here we report on the observation of
ultrastrong coupling between graphene and cavity plasmonic modes within the
experimental spectral sensitivity range. To enter this regime, we
intentionally designed the cavity such that four modes, the bare cavity mode
centered at $0.43\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ in addition to three
graphene modes, $s=0,1,2$, overlap in frequency within the experimental
measurement range of the THz circuit, achieved by making both the graphite and
the hBN between graphene and graphite thinner (Fig. 4a) as compared to the
device in Fig. 3.

Figure 4: Ultrastrong coupling between plasmonic modes of graphene and cavity

Figure 4: a Cavity design with a thin hBN layer between graphene and graphite
to center the bare cavity mode to
$0.43\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ (for cavity geometry details see
[42]). The graphene $s=2$ mode ultrastrongly couples to the cavity $m=0$ mode
inside the experimental spectral sensitivity range. b Experimentally-obtained
real part of the near-field THz conductivity spectra of a selection of
measured graphene carrier densities. The data were fit to a sum of four
Lorentzians as indicated by the shaded peaks. c Left panel: Detected resonance
frequencies as a function of carrier density for the four Lorentzians fitted
to the experimental data. The experimentally extracted normalized coupling
strength $\eta_{\text{exp}}=0.12\pm 0.01$, achieving the ultrastrong coupling
regime. The analytical model for this device (center panel) follows similar
trends, with $\eta_{\text{theory}}=0.11\pm 0.01$. The right panel shows the
simulated conductivity spectra at a carrier density of $n=0$, for the
uncoupled $m=0$ mode and for a carrier density of
$n=$1.24\text{\times}{10}^{12}\text{\,}{\mathrm{cm}}^{-2}$$ at which the
avoided crossing appears. d Quantification of experimental spectral weight
transfer between the two highest frequency modes.

The ultrastrong coupling of graphene plasmonic modes to the cavity mode is
evidenced by the near-field carrier-density dependent conductivity spectra,
shown in Fig. 4b. We fit the spectra with a sum of four Lorentzians to capture
the four expected modes of the cavity [42]. As the carrier density increases,
all modes shift to higher frequencies. At a graphene carrier density of
$n=0.5\times 10^{12}$ cm-2, the $s=2$ mode becomes resonant with the $m=0$
mode (dashed line in Fig. 4b). Rather than merely crossing in frequency, as
would be expected given the different momenta of the $m=0$ and $s=2$ modes
[42], they hybridize and split. The spectral separation of these two modes at
this carrier density corresponds to twice the coupling strength $g$, found to
be (50 $\pm$ 5) GHz, corresponding to $\sim
0.2\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. Comparing this to the frequency
of (0.43 $\pm$ 0.01) THz, $\sim 1.8\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$,
at which this coupling occurs, we obtained a normalized coupling strength of
$\eta_{\text{exp}}=0.12\pm 0.01$ characterizing this cavity as in the
ultrastrong coupling regime [42].

The ultrastrong coupling between the graphene and cavity plasmonic modes can
be seen to significantly modify the electrodynamic properties of the system.
Tracking the resonance frequencies of the four modes as a function of carrier
density (seen in Fig. 4c) reveals that the lowest order graphene mode $s=0$ is
again primarily uncoupled from the cavity as it scales with
$\nu_{\text{{s}=0}}\propto n^{0.24\pm 0.01}$. However, the $s=1$ and $s=2$
modes are significantly suppressed when compared to their uncoupled nature
[42].

The spectral weight of each mode was calculated as the area of each Lorentzian
shown in Fig. 4b and plotted as a function of carrier density in Fig. 4d. As
the graphene becomes metallic, the spectral weight of the upper branch
(purple) decreases, as the lower branch (magenta) rapidly increases in
spectral weight. This stands in sharp contrast to the expected spectral weight
behavior of graphene and graphite cavity modes in the absence of coupling,
where the spectral weight of graphene modes should increase in amplitude and
the cavity mode should only be weakly affected when the graphene carrier
density $n$ is increased [42]. The detected transfer in spectral weight across
the avoided crossing is thus a consequence of the hybridization and a
signature of ultrastrong coupling.

The analytical model captures the carrier-density dependence and spectral
weight trends of this device, as shown in the center panel of Fig. 4c. This
theory predicts $\eta_{\text{theory}}=0.11\pm 0.05$ (right panel), consistent
with the experimental observation. The unexpected carrier-density scaling of
the plasmon frequency and spectral weight transfer can be attributed to
tunable, multimodal couplings between the bare cavity mode and graphene modes,
which reach the ultrastrong coupling regime for the $s=2$ mode [42, 70].

##  5 Tunable cavity coupling mechanism

While the devices highlighted in this work exhibited coupling between sample
and cavity modes, it would be useful to also be able to build cavities with
negligible coupling such that the unperturbed electrodynamics of the sample
can be interrogated. The analytical model developed for this work provides
design principles to construct cavities which either minimize or maximize
coupling strengths to either sense or control the THz electrodynamics of
collective modes in vdW materials (Fig. 5).

Figure 5: Coupling strength tunability through geometry: From cavity sensing
to cavity control

Figure 5: a A cavity with a small-gap coplanar strip and no extended sample
regions behaves as a sensing tool. b A cavity with extended 2D material region
$W_{0}$ and gap $W_{2}$ can be used to control collective modes of 2D
materials. $\alpha$ represents a tuning factor for the cavity parameters.
Coplanar strip geometry parameters of $W_{0}=W_{1}=$
$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $W_{2}=$
$8\text{\,}\mathrm{\SIUnitSymbolMicro m}$ were used (for all other device
parameters see [42]). c Simulation of a device as in (a) with a graphene
monolayer and $\alpha=0.001$. Multiple undisturbed graphene plasmonic modes
can be obtained that all scale as $\nu\propto n^{0.25}$. d Simulation of a
device as in (b) with a graphene monolayer and $\alpha=1$. An avoided crossing
with $\eta=0.18$ is observed. e $\eta$ dependence on the scaling factor
$\alpha$ for the geometry defined in (b), with cavity sensing in the small
$\alpha$ limit and cavity control in the large $\alpha$ limit.

The theory predicts that a cavity built with a thick hBN, thick graphite gate,
and with a coplanar strip patterned to cover all but a small gap between the
traces ($\alpha W_{2}\ll W_{1}$, see Fig. 5a) minimizes the coupling between
the graphene and cavity collective modes (Fig. 5c). In such a device, momentum
conservation is preserved and thus the cavity and 2D material plasmon modes
are orthogonal to each other. Our theory and numerical simulations show that
both the current density $j$ and the spatial derivative, $dj/dx$, are
approximately continuous throughout the device [42]. As the THz spectroscopic
platform introduced here is capable of independently capturing the real and
imaginary parts of the uncoupled 2D material conductivity, which are necessary
for extracting spectral weight shifts, measurements of a “sensing cavity” will
provide deeper insight into the dynamics of correlated materials [71]. For
example, properties such as quasiparticle scattering rates, eigenfrequencies
of collective modes, magnitudes of electronic energy gaps, and the density of
superconducting carriers can be captured as a function of temperature,
displacement and magnetic field, and carrier density.

By contrast, the ultrastrong coupling can be enhanced by patterning the
dielectric environment (Fig. 5b). A cavity built with increased conductor
separation and additional vdW heterostructure widths outside of the coplanar
strip maximizes the ultrastrong coupling between modes of different momenta,
with predicted strengths of $\eta>0.18$ (see Fig. 5b,d). This is due to broken
translational symmetry of this cavity design, and thus broken momentum
conservation. While discrete plasmonic and cavity modes are typically
orthogonal to each other, the inhomogenous dielectric environment allows
graphene plasmonic modes to spatially overlap with the cavity mode, through
which they become entangled with each other [70]. Unscreened graphene plasmons
are excited in the gap and mediate the coupling of the screened cavity and
graphene plasmons. The relative widths of the metal traces to the widths of
the unscreened regions provide a tuning knob for the coupling strength. Thus,
the coupling can be controlled through the cavity geometry and carrier
density, allowing for either probing ($\alpha\ll 1$) or controlling
($\alpha\approx 1$) quantum materials (see Fig. 5e). The details of this
coupling are captured by our analytical model [42].

##  6 Cavity control of vdW quantum materials

This work demonstrates the first experimental evidence of the multimodal,
ultrastrong interaction between gate-tunable graphene plasmons and a plasmonic
graphite cavity in the deep sub-wavelength regime. Our results show that
cavity electrodynamics are relevant for vdW heterostructures with graphite
gates of typical dimensions. This raises the possibility that cavity effects
due to vacuum fluctuations may be relevant in vdW heterostructures exhibiting
strongly correlated quantum phases. In common vdW device architectures,
microstructured electrostatic gates and contacts may break momentum
conservation in a similar way to the metallic coplanar strips employed in this
work.

The platform developed here delivers a route to deterministic cavity control
of many-body phenomena in vdW heterostructures. Our advances in on-chip THz
circuitry allow both the light-matter coupling strength and spectral weight
transfer to be probed under the same experimental conditions in which DC
transport is measured. With this, the effects of cavity-mediated vacuum
fluctuations, as probed by transport, could be experimentally related to the
cavity hybridization, as probed by THz measurements. The near-field nature and
broken symmetries in these vdW heterostructure microcavities will enable the
future investigation of cavity control. This could be performed on not only
infrared active, but also optically silent modes [42], such as acoustic and
Raman-active phonons, superconducting collective modes [72], and topological
edge modes. While the cavity demonstrated here achieves ultrastrong coupling
in the linear regime, the possibility to drive the cavity to the nonlinear
regime using Floquet engineering protocols may enable additional control in
the non-equilibrium domain. With such a strategy, these devices provide a
novel and versatile platform for delivering spectroscopic information and
transformative cavity-control of emergent properties in vdW heterostructures.

## Supplementary information

### Acknowledgements

We thank A. Cavalleri, D.N. Basov, J.D. Adelinia, E. Wang and M. Chavez-
Cervantes for their helpful insight throughout the duration of the project. We
also thank J.M. Pizarro, T. Wehling, M. Buzzi, M. Fechner, J. Schmalian and E.
Demler for helpful conversations.

### Funding

We acknowledge support from the Max Planck-New York City Center for Non-
Equilibrium Quantum Phenomena. G.K. and A.R. acknowledges support by the
German Research Foundation through the Cluster of Excellence CUI: Advanced
Imaging of Matter (EXC 2056, project ID 390715994). H.M.B., M.H.M. and M.W.D.
acknowledge support from the Alexander von Humboldt Foundation. H.M.B.
acknowledges financial support from the European Union under the Marie
Sklodowska-Curie Grant Agreement no. 101062921 (Twist-TOC). The fabrication of
THz circuitry optimization devices was supported by the U.S. Department of
Energy, Office of Science, Basic Energy Sciences, under Early Career Award DE-
SC0024334. Development of the cavity coupling theory was supported by the
Center on Programmable Quantum Materials, an Energy Frontier Research Center
funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy
Sciences (BES), under award DE-SC0019443. K.W. and T.T. acknowledge support
from the JSPS KAKENHI (Grant Numbers 20H00354 and 23H02052) and World Premier
International Research Center Initiative (WPI), MEXT, Japan. A.R. acknowledges
support from the Grupos Consolidados (IT1453-22). The Flatiron Institute is a
division of the Simons Foundation. At the early stages of the project, G.K.
and K.K. were supported by SFB 925 - project 170620586/Deutsche
Forschungsgemeinschaft (German Research Foundation). D.M.K acknowledges
funding by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) within the Priority Program SPP 2244 “2DMP” - 443274199 and under
Germany’s Excellence Strategy - Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1 - 390534769. M.A.S. was funded by the
European Union (ERC, CAVMAT, project no. 101124492).

### Duration and volume of study

This project started in January 2020 and ended in March 2024. During this
time, vdW heterostructure microcavity fabrication protocols were developed,
the ultrafast cryogenic measurement setup was built, THz circuitry for probing
the cavity conductivity was developed and optimized, an analysis framework for
extracting the cavity conductivity was constructed, finite element simulations
to simulate the cavity modes were performed, and the analytical theory
describing the cavity modes and their coupling was developed. For this report,
42 THz devices were measured in total: 15 THz circuitry optimization devices,
16 bare graphite cavities, 3 graphene devices without graphite cavity, 2
cavities based on transition metal dichalcogenide gates and 6 gate-tunable
graphene cavity devices.

### Author contributions

J.W.M. conceived the experiment and supervised the project. G.K., B.S., D.H.,
K.K. and S.K. fabricated the THz cavity devices with help from S.M.L., X.L.,
J.H., Y.H. and F.S. B.S. designed and built the measurement setup. H.M.B. and
G.K. performed measurements and analyzed the data, with assistance from
M.W.D., K.K. and B.S. H.M.B. developed the data analysis framework. M.H.M.
developed the theory with support from A.R. and G.K., and G.K. performed the
corresponding numerical simulations. Custom measurement electronics and
circuit simulations were provided by T.M. and G.M.. T.T. and K.W. provided
high-quality hBN single crystals. E.B., D.M.K., M.A.S., A.R., C.J.E., A.M.P.
and G.M. provided valuable conceptual insight. G.K., H.M.B. and J.W.M. wrote
the manuscript with contributions from all authors.

### Competing interests

The authors declare no competing interests.

## References

  *   * Andrei and MacDonald [2020] Andrei, E.Y., MacDonald, A.H.: Graphene bilayers with a twist. Nature materials 19(12), 1265–1275 (2020)
  * Balents et al. [2020] Balents, L., Dean, C.R., Efetov, D.K., Young, A.F.: Superconductivity and strong correlations in moiré flat bands. Nature Physics 16(7), 725–733 (2020)
  * Mak and Shan [2022] Mak, K.F., Shan, J.: Semiconductor moiré materials. Nature Nanotechnology 17(7), 686–695 (2022)
  * Koppens et al. [2011] Koppens, F.H., Chang, D.E., Abajo, F.J.: Graphene plasmonics: a platform for strong light–matter interactions. Nano letters 11(8), 3370–3377 (2011)
  * Grigorenko et al. [2012] Grigorenko, A.N., Polini, M., Novoselov, K.: Graphene plasmonics. Nature photonics 6(11), 749–758 (2012)
  * Basov et al. [2020] Basov, D.N., Asenjo-Garcia, A., Schuck, P.J., Zhu, X., Rubio, A.: Polariton panorama. Nanophotonics 10(1), 549–577 (2020)
  * Low et al. [2017] Low, T., Chaves, A., Caldwell, J.D., Kumar, A., Fang, N.X., Avouris, P., Heinz, T.F., Guinea, F., Martin-Moreno, L., Koppens, F.: Polaritons in layered two-dimensional materials. Nature materials 16(2), 182–194 (2017)
  * Du et al. [2023] Du, L., Molas, M.R., Huang, Z., Zhang, G., Wang, F., Sun, Z.: Moiré photonics and optoelectronics. Science 379(6639), 0014 (2023)
  * Huang et al. [2020] Huang, B., McGuire, M.A., May, A.F., Xiao, D., Jarillo-Herrero, P., Xu, X.: Emergent phenomena and proximity effects in two-dimensional magnets and heterostructures. Nature Materials 19(12), 1276–1289 (2020)
  * Li et al. [2021] Li, T., Jiang, S., Shen, B., Zhang, Y., Li, L., Tao, Z., Devakul, T., Watanabe, K., Taniguchi, T., Fu, L., et al.: Quantum anomalous hall effect from intertwined moiré bands. Nature 600(7890), 641–646 (2021)
  * Park et al. [2023] Park, H., Cai, J., Anderson, E., Zhang, Y., Zhu, J., Liu, X., Wang, C., Holtzmann, W., Hu, C., Liu, Z., et al.: Observation of fractionally quantized anomalous hall effect. Nature 622(7981), 74–79 (2023)
  * Zeng et al. [2023] Zeng, Y., Xia, Z., Kang, K., Zhu, J., Knüppel, P., Vaswani, C., Watanabe, K., Taniguchi, T., Mak, K.F., Shan, J.: Thermodynamic evidence of fractional chern insulator in moiré MoTe2. Nature 622(7981), 69–73 (2023)
  * Cai et al. [2023] Cai, J., Anderson, E., Wang, C., Zhang, X., Liu, X., Holtzmann, W., Zhang, Y., Fan, F., Taniguchi, T., Watanabe, K., et al.: Signatures of fractional quantum anomalous hall states in twisted MoTe2. Nature 622(7981), 63–68 (2023)
  * Song and Gabor [2018] Song, J.C., Gabor, N.M.: Electron quantum metamaterials in van der waals heterostructures. Nature nanotechnology 13(11), 986–993 (2018)
  * Zibrov et al. [2017] Zibrov, A.A., Kometter, C., Zhou, H., Spanton, E., Taniguchi, T., Watanabe, K., Zaletel, M., Young, A.: Tunable interacting composite fermion phases in a half-filled bilayer-graphene landau level. Nature 549(7672), 360–364 (2017)
  * Graef et al. [2018] Graef, H., Mele, D., Rosticher, M., Banszerus, L., Stampfer, C., Taniguchi, T., Watanabe, K., Bocquillon, E., Fève, G., Berroir, J., et al.: Ultra-long wavelength dirac plasmons in graphene capacitors. Journal of Physics: Materials 1(1), 01–02 (2018)
  * Yoshioka et al. [2023] Yoshioka, K., Bernard, G., Wakamura, T., Hashisaka, M., Sasaki, K.-i., Sasaki, S., Watanabe, K., Taniguchi, T., Kumada, N.: On-chip transfer of ultrashort graphene plasmon wavepackets using terahertz electronics. arXiv preprint arXiv:2311.02821 (2023)
  * Zhao et al. [2023] Zhao, W., Wang, S., Chen, S., Zhang, Z., Watanabe, K., Taniguchi, T., Zettl, A., Wang, F.: Observation of hydrodynamic plasmons and energy waves in graphene. Nature 614(7949), 688–693 (2023)
  * Fei et al. [2015] Fei, Z., Goldflam, M., Wu, J.-S., Dai, S., Wagner, M., McLeod, A., Liu, M., Post, K., Zhu, S., Janssen, G., et al.: Edge and surface plasmons in graphene nanoribbons. Nano letters 15(12), 8271–8276 (2015)
  * Lee et al. [2019] Lee, I.-H., Yoo, D., Avouris, P., Low, T., Oh, S.-H.: Graphene acoustic plasmon resonator for ultrasensitive infrared spectroscopy. Nature nanotechnology 14(4), 313–319 (2019)
  * Schlawin et al. [2022] Schlawin, F., Kennes, D.M., Sentef, M.A.: Cavity quantum materials. Applied Physics Reviews 9(1) (2022)
  * Ruggenthaler et al. [2018] Ruggenthaler, M., Tancogne-Dejean, N., Flick, J., Appel, H., Rubio, A.: From a quantum-electrodynamical light–matter description to novel spectroscopies. Nature Reviews Chemistry 2(3), 1–16 (2018)
  * Hübener et al. [2021] Hübener, H., De Giovannini, U., Schäfer, C., Andberger, J., Ruggenthaler, M., Faist, J., Rubio, A.: Engineering quantum materials with chiral optical cavities. Nature materials 20(4), 438–442 (2021)
  * Eckhardt et al. [2024] Eckhardt, C.J., Chattopadhyay, S., Kennes, D.M., Demler, E.A., Sentef, M.A., Michael, M.H.: Theory of resonantly enhanced photo-induced superconductivity. Nature Communications 15(1), 2300 (2024)
  * Sentef et al. [2018] Sentef, M.A., Ruggenthaler, M., Rubio, A.: Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity. Science advances 4(11), 6969 (2018)
  * Curtis et al. [2019] Curtis, J.B., Raines, Z.M., Allocca, A.A., Hafezi, M., Galitski, V.M.: Cavity quantum eliashberg enhancement of superconductivity. Physical review letters 122(16), 167002 (2019)
  * Curtis et al. [2023] Curtis, J.B., Michael, M.H., Demler, E.: Local fluctuations in cavity control of ferroelectricity. arXiv preprint arXiv:2301.01884 (2023)
  * Masuki and Ashida [2023] Masuki, K., Ashida, Y.: Cavity moir$\backslash$’e materials: Controlling magnetic frustration with quantum light-matter interaction. arXiv preprint arXiv:2302.11582 (2023)
  * Bloch et al. [2022] Bloch, J., Cavalleri, A., Galitski, V., Hafezi, M., Rubio, A.: Strongly correlated electron–photon systems. Nature 606(7912), 41–48 (2022)
  * Allocca et al. [2019] Allocca, A.A., Raines, Z.M., Curtis, J.B., Galitski, V.M.: Cavity superconductor-polaritons. Physical Review B 99(2), 020504 (2019)
  * Raines et al. [2020] Raines, Z.M., Allocca, A.A., Hafezi, M., Galitski, V.M.: Cavity higgs polaritons. Physical Review Research 2(1), 013143 (2020)
  * Viñas Boström et al. [2023] Viñas Boström, E., Sriram, A., Claassen, M., Rubio, A.: Controlling the magnetic state of the proximate quantum spin liquid $\alpha$-RuCl3 with an optical cavity. npj Computational Materials 9(1), 202 (2023)
  * Wang et al. [2019] Wang, X., Ronca, E., Sentef, M.A.: Cavity quantum electrodynamical chern insulator: Towards light-induced quantized anomalous hall effect in graphene. Physical Review B 99(23), 235156 (2019)
  * Masuki and Ashida [2023] Masuki, K., Ashida, Y.: Berry phase and topology in ultrastrongly coupled quantum light-matter systems. Physical Review B 107(19), 195104 (2023)
  * Chakraborty and Piazza [2021] Chakraborty, A., Piazza, F.: Long-range photon fluctuations enhance photon-mediated electron pairing and superconductivity. Physical Review Letters 127(17), 177002 (2021)
  * Andolina et al. [2024] Andolina, G.M., De Pasquale, A., Pellegrino, F.M.D., Torre, I., Koppens, F.H., Polini, M.: Amperean superconductivity cannot be induced by deep subwavelength cavities in a two-dimensional material. Physical Review B 109(10), 104513 (2024)
  * Forn-Díaz et al. [2019] Forn-Díaz, P., Lamata, L., Rico, E., Kono, J., Solano, E.: Ultrastrong coupling regimes of light-matter interaction. Reviews of Modern Physics 91(2), 025005 (2019)
  * Frisk Kockum et al. [2019] Frisk Kockum, A., Miranowicz, A., De Liberato, S., Savasta, S., Nori, F.: Ultrastrong coupling between light and matter. Nature Reviews Physics 1(1), 19–40 (2019)
  * Thomas et al. [2021] Thomas, A., Devaux, E., Nagarajan, K., Rogez, G., Seidel, M., Richard, F., Genet, C., Drillon, M., Ebbesen, T.W.: Large enhancement of ferromagnetism under a collective strong coupling of ybco nanoparticles. Nano letters 21(10), 4365–4370 (2021)
  * Jarc et al. [2023] Jarc, G., Mathengattil, S.Y., Montanaro, A., Giusti, F., Rigoni, E.M., Sergo, R., Fassioli, F., Winnerl, S., Dal Zilio, S., Mihailovic, D., et al.: Cavity-mediated thermal control of metal-to-insulator transition in 1T-TaS2. Nature 622(7983), 487–492 (2023)
  * Appugliese et al. [2022] Appugliese, F., Enkner, J., Paravicini-Bagliani, G.L., Beck, M., Reichl, C., Wegscheider, W., Scalari, G., Ciuti, C., Faist, J.: Breakdown of topological protection by cavity vacuum fields in the integer quantum hall effect. Science 375(6584), 1030–1034 (2022)
  * [42] Supplementary information
  * Ashida et al. [2023] Ashida, Y., İmamoğlu, A., Demler, E.: Cavity quantum electrodynamics with hyperbolic van der waals materials. Physical Review Letters 130(21), 216901 (2023)
  * Herzig Sheinfux et al. [2024] Herzig Sheinfux, H., Orsini, L., Jung, M., Torre, I., Ceccanti, M., Marconi, S., Maniyara, R., Barcons Ruiz, D., Hötger, A., Bertini, R., et al.: High-quality nanocavities through multimodal confinement of hyperbolic polaritons in hexagonal boron nitride. Nature Materials, 1–7 (2024)
  * Dirnberger et al. [2023] Dirnberger, F., Quan, J., Bushati, R., Diederich, G.M., Florian, M., Klein, J., Mosina, K., Sofer, Z., Xu, X., Kamra, A., et al.: Magneto-optics in a van der waals magnet tuned by self-hybridized polaritons. Nature 620(7974), 533–537 (2023)
  * Auerbach [2012] Auerbach, A.: Interacting Electrons and Quantum Magnetism. Springer, ??? (2012)
  * Novoselov et al. [2016] Novoselov, K.S., Mishchenko, A., Carvalho, A., Castro Neto, A.: 2d materials and van der waals heterostructures. Science 353(6298), 9439 (2016)
  * Feurer et al. [2007] Feurer, T., Stoyanov, N.S., Ward, D.W., Vaughan, J.C., Statz, E.R., Nelson, K.A.: Terahertz polaritonics. Annu. Rev. Mater. Res. 37, 317–350 (2007)
  * Berman et al. [2010] Berman, O.L., Kezerashvili, R.Y., Lozovik, Y.E., Snoke, D.W.: Bose–einstein condensation and superfluidity of trapped polaritons in graphene and quantum wells embedded in a microcavity. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368(1932), 5459–5482 (2010)
  * Hakala et al. [2018] Hakala, T.K., Moilanen, A.J., Väkeväinen, A.I., Guo, R., Martikainen, J.-P., Daskalakis, K.S., Rekola, H.T., Julku, A., Törmä, P.: Bose–einstein condensation in a plasmonic lattice. Nature Physics 14(7), 739–744 (2018)
  * Di Battista et al. [2024] Di Battista, G., Fong, K., Diez-Carlon, A., Watanabe, K., Taniguchi, T., Efetov, D.: Ultra-low carrier density superconducting bolometers with single photon sensitivity based on magic-angle twisted bilayer graphene. arXiv preprint arXiv:2403.02049 (2024)
  * Bogdanov et al. [2019] Bogdanov, S.I., Boltasseva, A., Shalaev, V.M.: Overcoming quantum decoherence with plasmonics. Science 364(6440), 532–533 (2019)
  * McIver et al. [2020] McIver, J.W., Schulte, B., Stein, F.-U., Matsuyama, T., Jotzu, G., Meier, G., Cavalleri, A.: Light-induced anomalous hall effect in graphene. Nature physics 16(1), 38–41 (2020)
  * Potts et al. [2023] Potts, A.M., Nayak, A.K., Nagel, M., Kaj, K., Stamenic, B., John, D.D., Averitt, R.D., Young, A.F.: On-chip time-domain terahertz spectroscopy of superconducting films below the diffraction limit. Nano Letters 23(9), 3835–3841 (2023)
  * Gallagher et al. [2019] Gallagher, P., Yang, C.-S., Lyu, T., Tian, F., Kou, R., Zhang, H., Watanabe, K., Taniguchi, T., Wang, F.: Quantum-critical conductivity of the dirac fluid in graphene. Science 364(6436), 158–162 (2019)
  * Yoshioka et al. [2022] Yoshioka, K., Wakamura, T., Hashisaka, M., Watanabe, K., Taniguchi, T., Kumada, N.: Ultrafast intrinsic optical-to-electrical conversion dynamics in a graphene photodetector. Nature Photonics 16(10), 718–723 (2022)
  * Wang et al. [2023] Wang, E., Adelinia, J., Chavez-Cervantes, M., Matsuyama, T., Fechner, M., Buzzi, M., Meier, G., Cavalleri, A.: Superconducting nonlinear transport in optically driven high-temperature k3c60. Nature Communications 14(1), 7233 (2023)
  * Parker-Jervis et al. [2020] Parker-Jervis, R., Park, S., Zonetti, S., Sydoruk, O., Wu, J., Wood, C., Li, L., Davies, A., Linfield, E., Cunningham, J.: On-chip terahertz spectroscopy of magnetoplasmons in a two-dimensional electron gas. In: 2020 45th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz), pp. 1–1 (2020). IEEE
  * Zhong et al. [2008] Zhong, Z., Gabor, N.M., Sharping, J.E., Gaeta, A.L., McEuen, P.L.: Terahertz time-domain measurement of ballistic electron resonance in a single-walled carbon nanotube. Nature nanotechnology 3(4), 201–205 (2008)
  * Karadi et al. [1994] Karadi, C., Jauhar, S., Kouwenhoven, L., Wald, K., Orenstein, J., McEuen, P., Nagamune, Y., Sakaki, H.: Dynamic response of a quantum point contact. JOSA B 11(12), 2566–2571 (1994)
  * Dean et al. [2010] Dean, C.R., Young, A.F., Meric, I., Lee, C., Wang, L., Sorgenfrei, S., Watanabe, K., Taniguchi, T., Kim, P., Shepard, K.L., et al.: Boron nitride substrates for high-quality graphene electronics. Nature nanotechnology 5(10), 722–726 (2010)
  * Svintsov and Alymov [2023] Svintsov, D.A., Alymov, G.V.: Refraction laws for two-dimensional plasmons. Physical Review B 108(12), 121410 (2023)
  * Nuss and Orenstein [2007] Nuss, M.C., Orenstein, J.: Terahertz time-domain spectroscopy. Millimeter and submillimeter wave spectroscopy of solids, 7–50 (2007)
  * Armitage [2009] Armitage, N.: Electrodynamics of correlated electron systems. arXiv preprint arXiv:0908.1126 (2009)
  * Basov et al. [2011] Basov, D.N., Averitt, R.D., Van Der Marel, D., Dressel, M., Haule, K.: Electrodynamics of correlated electron materials. Reviews of Modern Physics 83(2), 471 (2011)
  * Alonso-González et al. [2017] Alonso-González, P., Nikitin, A.Y., Gao, Y., Woessner, A., Lundeberg, M.B., Principi, A., Forcellini, N., Yan, W., Vélez, S., Huber, A.J., et al.: Acoustic terahertz graphene plasmons revealed by photocurrent nanoscopy. Nature nanotechnology 12(1), 31–35 (2017)
  * Xu et al. [2023] Xu, S., Li, Y., Vitalone, R.A., Jing, R., Sternbach, A.J., Zhang, S., Ingham, J., Delor, M., McIver, J., Yankowitz, M., et al.: Electronic interactions in dirac fluids visualized by nano-terahertz spacetime mapping. arXiv preprint arXiv:2311.11502 (2023)
  * Hwang and Sarma [2009] Hwang, E., Sarma, S.D.: Plasmon modes of spatially separated double-layer graphene. Physical Review B 80(20), 205405 (2009)
  * Hwang and Sarma [2007] Hwang, E., Sarma, S.D.: Dielectric function, screening, and plasmons in two-dimensional graphene. Physical Review B 75(20), 205418 (2007)
  * Tay et al. [2023] Tay, F., Mojibpour, A., Sanders, S., Liang, S., Xu, H., Gardner, G.C., Baydin, A., Manfra, M.J., Alabastri, A., Hagenmüller, D., et al.: Ultrastrong photon-photon coupling. arXiv preprint arXiv:2308.12427 (2023)
  * Van Loon et al. [2014] Van Loon, E., Hafermann, H., Lichtenstein, A., Rubtsov, A., Katsnelson, M.: Plasmons in strongly correlated systems: spectral weight transfer and renormalized dispersion. Physical Review Letters 113(24), 246407 (2014)
  * Costa et al. [2021] Costa, A., Gonçalves, P., Basov, D., Koppens, F.H., Mortensen, N.A., Peres, N.: Harnessing ultraconfined graphene plasmons to probe the electrodynamics of superconductors. Proceedings of the National Academy of Sciences 118(4), 2012847118 (2021)