* Bittau et al. (2017) Andrea Bittau, Úlfar Erlingsson, Petros Maniatis, Ilya Mironov, Ananth Raghunathan, David Lie, Mitch Rudominer, Ushasree Kode, Julien Tinnes, and Bernhard Seefeld. 2017\.  Prochlo: Strong Privacy for Analytics in the Crowd. In _Proceedings of the 26th Symposium on Operating Systems Principles_ (Shanghai, China) _(SOSP ’17)_. Association for Computing Machinery, New York, NY, USA, 441–459.  https://doi.org/10.1145/3132747.3132769
  * Commission (2022) NYC Taxi & Limousine Commission. 2022\.  _TLC Trip Record Data_.  https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page
  * Ding et al. (2017) Bolin Ding, Janardhan Kulkarni, and Sergey Yekhanin. 2017\.  Collecting Telemetry Data Privately. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ _(NIPS’17)_. 3574–3583.  https://doi.org/10.5555/3294996.3295115
  * Dong and Yi (2021) Wei Dong and Ke Yi. 2021\.  Residual Sensitivity for Differentially Private Multi-Way Joins. In _Proceedings of the 2021 International Conference on Management of Data_ (Virtual Event, China) _(SIGMOD ’21)_. Association for Computing Machinery, New York, NY, USA, 432–444.  https://doi.org/10.1145/3448016.3452813
  * Dwork et al. (2006) Cynthia Dwork, Frank McSherry, Kobbi Nissim, and Adam Smith. 2006.  Calibrating Noise to Sensitivity in Private Data Analysis. In _Theory of Cryptography_ , Shai Halevi and Tal Rabin (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 265–284.
  * Dwork and Roth (2014) Cynthia Dwork and Aaron Roth. 2014.  The Algorithmic Foundations of Differential Privacy.  _Foundations and Trends in Theoretical Computer Science_ 9, 3–4 (2014), 211–407.  https://doi.org/10.1561/0400000042
  * Eichmann et al. (2020) Philipp Eichmann, Emanuel Zgraggen, Carsten Binnig, and Tim Kraska. 2020.  IDEBench: A Benchmark for Interactive Data Exploration. In _Proceedings of the 2020 ACM SIGMOD International Conference on Management of Data_. 1555–1569.  https://doi.org/10.1145/3318464.3380574
  * Gaboardi et al. (2019) Marco Gaboardi, Michael Hay, and Salil Vadhan. 2019.  _A Programming Framework for OpenDP_.  https://projects.iq.harvard.edu/opendp
  * Gaboardi et al. (2016) Marco Gaboardi, James Honaker, Gary King, Jack Murtagh, Kobbi Nissim, Jonathan Ullman, and Salil Vadhan. 2016.  Psi: a private data sharing interface.  (2016).  http://arxiv.org/abs/1609.04340
  * Ge et al. (2019) Chang Ge, Xi He, Ihab F. Ilyas, and Ashwin Machanavajjhala. 2019.  APEx: Accuracy-Aware Differentially Private Data Exploration. In _Proceedings of the 2019 International Conference on Management of Data_ (Amsterdam, Netherlands). Association for Computing Machinery, New York, NY, USA, 177–194.  https://doi.org/10.1145/3299869.3300092
  * Ge et al. (2021) Chang Ge, Shubhankar Mohapatra, Xi He, and Ihab F. Ilyas. 2021\.  Kamino: Constraint-Aware Differentially Private Data Synthesis.  _VLDB_ (2021).
  * Hay et al. (2010) Michael Hay, Vibhor Rastogi, Gerome Miklau, and Dan Suciu. 2010.  Boosting the Accuracy of Differentially Private Histograms through Consistency.  _Proceedings of the VLDB Endowment_ 3, 1–2 (Sept. 2010), 1021–1032.  https://doi.org/10.14778/1920841.1920970
  * Johnson et al. (2018) Noah Johnson, Joseph P. Near, and Dawn Song. 2018.  Towards Practical Differential Privacy for SQL Queries.  _Proceedings of the VLDB Endowment_ 11, 5 (Jan. 2018), 526–539.  https://doi.org/10.1145/3177732.3177733
  * Kohavi (1996) Ron Kohavi. 1996\.  Scaling up the accuracy of naive-bayes classifiers: A decision-tree hybrid. https://archive.ics.uci.edu/ml/datasets/adult. In _KDD_ , Vol. 96. 202–207.
  * Kotsogiannis et al. (2019) Ios Kotsogiannis, Yuchao Tao, Xi He, Maryam Fanaeepour, Ashwin Machanavajjhala, Michael Hay, and Gerome Miklau. 2019.  PrivateSQL: A Differentially Private SQL Query Engine.  _Proceedings of the VLDB Endowment_ 12, 11 (July 2019), 1371–1384.  https://doi.org/10.14778/3342263.3342274
  * Koufogiannis et al. (2016) Fragkiskos Koufogiannis, Shuo Han, and George J Pappas. 2016\.  Gradual Release of Sensitive Data under Differential Privacy.  _Journal of Privacy and Confidentiality_ 7, 2 (2016), 23–52.  https://doi.org/10.29012/jpc.v7i2.649
  * Li et al. (2015) Chao Li, Gerome Miklau, Michael Hay, Andrew Mcgregor, and Vibhor Rastogi. 2015.  The Matrix Mechanism: Optimizing Linear Counting Queries under Differential Privacy.  _The VLDB Journal_ 24, 6 (Dec. 2015), 757–781.  https://doi.org/10.1007/s00778-015-0398-x
  * Machanavajjhala et al. (2008) Ashwin Machanavajjhala, Daniel Kifer, John Abowd, Johannes Gehrke, and Lars Vilhuber. 2008\.  Privacy: Theory meets practice on the map. In _2008 IEEE 24th international conference on data engineering_. IEEE, IEEE, Cancun, Mexico, 277–286.  https://doi.org/10.1109/ICDE.2008.4497436
  * Mazmudar et al. (2022) Miti Mazmudar, Thomas Humphries, Jiaxiang Liu, Matthew Rafuse, and Xi He. 2022.  _CacheDP artifact_.  https://gitlab.uwaterloo.ca/m2mazmud/cachedp-public
  * Mazmudar et al. (2023) Miti Mazmudar, Thomas Humphries, Jiaxiang Liu, Matthew Rafuse, and Xi He. 2023.  Cache Me If You Can: Accuracy-Aware Inference Engine for Differentially Private Data Exploration.  _Proceedings of the VLDB Endowment_.
  * McKenna et al. ([n.d.]) Ryan McKenna, Daniel Sheldon, and Gerome Miklau. [n.d.].  Graphical-model based estimation and inference for differential privacy.  arXiv:1901.09136
  * McSherry (2010) Frank McSherry. 2010\.  Privacy Integrated Queries: An Extensible Platform for Privacy-Preserving Data Analysis.  _Commun. ACM_ 53, 9 (2010), 89–97.  https://doi.org/10.1145/1810891.1810916
  * Mohan et al. (2012) Prashanth Mohan, Abhradeep Thakurta, Elaine Shi, Dawn Song, and David Culler. 2012.  GUPT: Privacy Preserving Data Analysis Made Easy. In _Proceedings of the 2012 ACM SIGMOD International Conference on Management of Data_. Association for Computing Machinery, 349–360.  https://doi.org/10.1145/2213836.2213876
  * Nanayakkara et al. (2022) Priyanka Nanayakkara, Johes Bater, Xi He, Jessica Hullman, and Jennie Rogers. 2022\.  Visualizing Privacy-Utility Trade-Offs in Differentially Private Data Releases.  _CoRR_ (2022).  arXiv:2201.05964
  * of Transportation (2022) United States Department of Transportation. 2022.  Bureau of Transportation Statistics.  https://transtats.bts.gov
  * Peng et al. (2013) S. Peng, Y. Yang, Z. Zhang, M. Winslett, and Y. Yu. 2013.  Query optimization for differentially private data management systems. In _2013 IEEE 29th International Conference on Data Engineering (ICDE)_. IEEE, Brisbane, QLD, Australia, 1093–1104.  https://doi.org/10.1109/ICDE.2013.6544900
  * Qardaji et al. (2013) Wahbeh Qardaji, Weining Yang, and Ninghui Li. 2013.  Understanding hierarchical methods for differentially private histograms. In _Proceedings of the VLDB Endowment_ , Vol. 6. 1954–1965.  http://www.vldb.org/pvldb/vol6/p1954-qardaji.pdf
  * Tantipongpipat et al. (2021) Uthaipon Tantipongpipat, Chris Waites, Digvijay Boob, Amaresh Siva, and Rachel Cummings. 2021\.  Differentially private synthetic mixed-type data generation for unsupervised learning.  _Intelligent Decision Technologies_ (2021).
  * Tao et al. (2021) Yuchao Tao, Ryan McKenna, Michael Hay, Ashwin Machanavajjhala, and Gerome Miklau. 2021.  Benchmarking Differentially Private Synthetic Data Generation Algorithms.  _TPDP_ (2021).  https://tpdp.journalprivacyconfidentiality.org/2021/papers/NingUQKH21.pdf
  * Vesga et al. (2019) Elisabet Lobo Vesga, Alejandro Russo, and Marco Gaboardi. 2019\.  A Programming Framework for Differential Privacy with Accuracy Concentration Bounds.  _CoRR_ (2019).  arXiv:1909.07918
  * Wilson et al. (2020) Royce J Wilson, Celia Yuxin Zhang, William Lam, Damien Desfontaines, Daniel Simmons-Marengo, and Bryant Gipson. 2020\.  Differentially Private SQL with Bounded User Contribution.  _Proceedings on Privacy Enhancing Technologies_ 2020 (2020), 230–250.  Issue 2. https://doi.org/10.2478/popets-2020-0025
  * Xiao et al. (2021) Yingtai Xiao, Zeyu Ding, Yuxin Wang, Danfeng Zhang, and Daniel Kifer. 2021.  Optimizing Fitness-for-Use of Differentially Private Linear Queries.  _Proceedings of the VLDB Endowment_ 14, 10 (2021), 1730–1742.  https://doi.org/10.14778/3467861.3467864
  * Zhang et al. (2016) Jun Zhang, Xiaokui Xiao, and Xing Xie. 2016.  PrivTree: A Differentially Private Algorithm for Hierarchical Decompositions. In _Proceedings of the 2016 International Conference on Management of Data_ _(SIGMOD ’16)_. ACM, 155–170.  https://doi.org/10.1145/2882903.2882928

##  Appendix A Proofs

###  A.1. End-to-end Privacy Proof

####  A.1.1. Proof of Theorem 3.2

Recall the theorem states that CacheDP, as defined in Algorithm 1, satisfies
$\mathcal{B}$-DP.

###### Proof.

(sketch) We begin by addressing the cost estimation phase (lines 4-10). The
cost estimation phases are independent of the data. In addition, each DP
mechanism (MMM or RP) with its corresponding chosen strategy ensures
$\epsilon_{i}$-DP (Proposition 4.1, (Koufogiannis et al., 2016, Theorem 1A)).
At line 11, we check if running the chosen DP mechanism (MMM or RP with the
corresponding chosen strategy) will exceed the total privacy budget
$\mathcal{B}$ by sequential composition (Dwork and Roth, 2014). We only run
the DP mechanism if the total budget is sufficient. This ensures the overall
Algorithm 1 satisfies $\mathcal{B}$-DP.

∎

###  A.2. MMM Module Proofs

####  A.2.1. Proof of Proposition 4.1

Recall the proposition states that the AnswerWorkload interface of MMM
(Algorithm 2) satisfies $\epsilon$-DP, where $\epsilon$ is the output of this
interface.

###### Proof.

(sketch) When the cache is empty, the privacy of MMM (with optional SE)
follows from the matrix mechanisms privacy guarantee in Proposition 2.1. When
there are entries in the cache, we split the strategy into the free and paid
matrix. The paid matrix is private by the same reasoning as above. The privacy
of the free matrix responses and the final concatenation of free and paid
responses follow by the post processing lemma of DP (Dwork and Roth, 2014,
Proposition 2.1). ∎

####  A.2.2. Proof of Proposition 4.2

Given an instant strategy $\mathbf{A}=(\mathbf{F}||\mathbf{P})$ with a vector
of $k$ noise parameters
$\mathbf{b}=\mathbf{b}_{\mathbf{F}}||\mathbf{b}_{\mathbf{P}}$, the error to a
workload $\mathbf{W}$ using the AnswerWorkload interface of MMM (Algorithm 2)
is $\|\mathbf{W}\mathbf{A}^{+}Lap(\mathbf{b})\|$, where $Lap(\mathbf{b})$
draws independent noise from $Lap(\mathbf{b}[1])$, $\ldots,Lap(\mathbf{b}[k])$
respectively. We can simplify its expected total square error as
$\|\mathbf{W}\mathbf{A}^{+}diag(\mathbf{b})\|_{F}^{2}$ where
$diag(\mathbf{b})$ is a diagonal matrix with
$diag(\mathbf{b})[i,i]=\mathbf{b}[i]$.

###### Proof.

The error to the MMM is
$\|\mathbf{W}\mathbf{A}^{+}(\mathbf{A}\mathbf{x}+Lap(\mathbf{b}))-\mathbf{W}\mathbf{x}\|=\|\mathbf{W}\mathbf{A}^{+}Lap(\mathbf{b})\|$.
The expected total square error
$\mathbb{E}[\|\mathbf{W}\mathbf{A}^{+}Lap(\mathbf{b})\|^{2}_{2}]$ can be
expanded to
$\sum_{i=1}^{l}\mathbb{E}[\sum_{j=1}^{k}(\mathbf{W}\mathbf{A}^{+}[i,j]Lap(\mathbf{b}[j]))^{2}]$.
As the $k$ noise variables are independent and has a zero mean, then we have
the error expression equals to

$\sum_{i=1}^{l}\sum_{j=1}^{k}(\mathbf{W}\mathbf{A}^{+}[i,j])^{2}\mathbb{E}[Lap(\mathbf{b}[j]))^{2}]=\sum_{i=1}^{l}\sum_{j=1}^{k}(\mathbf{W}\mathbf{A}^{+}[i,j])^{2}\mathbf{b}[j]^{2}$

which is equivalent to $\|\mathbf{W}\mathbf{A}^{+}diag(\mathbf{b})\|_{F}^{2}$.
∎

####  A.2.3. Proof of Theorem 4.1

The optimal solution to simplified CE problem incurs a smaller privacy cost
$\epsilon$ than the privacy cost $\epsilon_{\mathbf{F}=\emptyset}$ of the
matrix mechanism without cache, i.e., MMM with $\mathbf{F}=\emptyset$.

###### Proof.

(sketch) Let $b^{*}$ be the noise parameter for $\mathbf{A}$ in the matrix
mechanism without cache to achieve the desired accuracy requirement. We can
show that $b^{*}$ is a valid solution to the simplified CE problem: setting
$\mathbf{b}_{\mathbf{F}}=[c.b\in\mathcal{C}~{}|~{}c.\mathbf{a}\in\mathcal{C}\cap\mathbf{A},c.b\leq
b^{*}]$ and $\mathbf{b}_{\mathbf{P}}=[b^{*}~{}|~{}\mathbf{a}\in\mathbf{P}]$
satisfies the accuracy requirement. As
$\|\mathbf{P}\|_{1}\leq\|\mathbf{A}\|_{1}$, the privacy cost
$\epsilon=\frac{\|P\|_{1}}{b^{*}}$ for $b_{\mathbf{P}}=b^{*}$ is smaller than
$\epsilon_{\mathbf{F}=\emptyset}=\frac{\|\mathbf{A}\|_{1}}{b^{*}}$. The
optimal solution to the simplified CE problem has a smaller or the same
privacy cost than a valid solution $b_{\mathbf{P}}=b^{*}$. ∎

###  A.3. Full-rank Transformer (FRT) Proof

Recall Theorem 5.1 states that Given a global strategy $\mathbb{A}^{*}$ in a
$k$-ary tree structure, and an instant strategy
$\mathbb{A}\subseteq\mathbb{A}^{*}$, transformStrategy outputs a strategy
$\mathbf{A}$ that is full rank and supports $\mathbb{A}$. We begin by proving
the following lemma.

###### Lemma A.0.

Given a global strategy $\mathbb{A}^{*}$ in a $k$-ary tree structure, and an
instant strategy $\mathbb{A}\subseteq\mathbb{A}^{*}$, running
getTransformationMatrix($\mathbb{A}$) in Algorithm 4 increases the number of
non-empty buckets in $\mathbb{T}$ by at most 1, for each row
$\mathbb{A}[i]\in\mathbb{A}$.

###### Proof.

If the condition in line 4 is met the result follows trivially. We consider
the remaining case where $\mathbb{A}[i]$ intersects with at least one bucket.
We recall that all entries in $\mathbb{A}$ represent nodes in a $k$-ary tree.
Since adding and subtracting nodes on a $k$-ary tree always results in a
combination of one or more nodes on a $k$-ary tree, we can conclude that at
the end of each loop, all of the buckets are disjoint. For each new row
$\mathbb{A}[i]$, if this row intersects with the buckets in $\mathbb{T}$, then
$\mathbb{A}[i]$ is either (i) a descendant node of one bucket in $\mathbb{T}$,
or (ii) an ancestor node of one or more buckets in $\mathbb{T}$. For the first
case, assume $\mathbb{A}[i]$ is a descendant node of
$\mathbb{t}\in\mathbb{T}$. It cannot be the descendant of other buckets, as
the other buckets are disjoint with $\mathbb{t}$. Then $\mathbb{t}$ will be
replaced by $\mathbb{t}\cap\mathbb{A}[i]$ and $\mathbb{t}-\mathbb{A}[i]$ (line
9), and hence the size of $\mathbb{T}$ increases by 1. For the second case,
assume $\mathbb{A}[i]$ is the ancestor node of multiple buckets in
$\mathbb{T}$. We denote these buckets as
$\\{\mathbb{T}[j_{1}],\ldots,\mathbb{T}[j_{k}]\\}$, then all these buckets
will remain the same (line 9 just adds $\mathbb{t}$), and at most one
additional bucket
$\mathbb{A}[i]-\sum_{\mathbb{t}\in\mathbb{T}^{\prime}\wedge\mathbb{t}\cdot\mathbb{A}[i]\neq
0}\mathbb{t}$ is added in line 10. ∎

We now prove the main result (Theorem 5.1) that transformStrategy ensures full
rank matrices.

###### Proof.

We begin by discussing how $\mathbf{A}$ represents the queries in
$\mathbb{A}$. First, we note that $\mathbf{A}$ and $\mathbb{A}$ have the same
number of rows. The only difference is that $\mathbf{A}$ represents the
queries on the data vector $\mathbf{x}$ where as $\mathbb{A}$ uses
$\mathbb{x}$. By construction, we have
$\mathbb{A}\mathbb{x}=\mathbf{A}\mathbb{T}\mathbb{x}=\mathbf{A}\mathbf{x}$.

Since $\mathbf{A}$ is a representation of $\mathbb{A}$ using the only the
buckets generated in getTransformationMatrix, the number of columns in
$\mathbf{A}$ is the same as the number of buckets, $|\mathbb{T}|$. To show
that $\mathbf{A}$ is always full rank, we will first show that
$|row(\mathbf{A})|\geq|col(\mathbf{A})|$, where $|row(\mathbf{A})|$ and
$|col(\mathbf{A})|$ represent the number of rows and the number of columns in
$\mathbf{A}$. Or equivalently, $|row(\mathbf{A})|\geq|\mathbb{T}|$. This
follows by inductively applying Lemma A.1.

Next we show that $rank(\mathbf{A})=|col(\mathbf{A})|$. We once again show
this by induction on the number of rows in $\mathbb{A}$. In the base case
($|\mathbb{A}|=1$) applying transformStrategy, we get $\mathbf{A}=[1]$ and
thus $rank(\mathbf{A})=1$. Now we assume that
$rank(\mathbf{A})=|col(\mathbf{A})|$ for some $\mathbf{A}$ obtained from
$\mathbb{A}$. We consider adding a new row $r_{new}$ to $\mathbb{A}$ and
define $\bar{\mathbb{A}}=\mathbb{A}\cup r_{new}$. Let $\bar{\mathbf{A}}$
represent the result of applying transformStrategy to $\bar{\mathbb{A}}$.

When adding this row there are two possible cases: First, consider the case
where the newly added row, $r_{new}$, can be represented as a linear
combination of $\mathbb{T}$. That is the buckets created in
getTransformationMatrix are the same for $\mathbb{A}$ and $\bar{\mathbb{A}}$
In this case, $\bar{\mathbf{A}}=\mathbf{A}\cup r^{\prime}_{new}$. Thus
$|col(\bar{\mathbf{A}})|=rank(\bar{\mathbf{A}})$ since the number of linearly
independent columns (or buckets) did not increase or decrease by adding a row.
On the other hand, if the newly added row, $r_{new}$, cannot be represented as
a linear combination of $\mathbb{T}$, then $r_{new}$ must be linearly
independent of $\mathbb{A}$. Thus, $\bar{\mathbf{A}}$ has one additional
linearly independent row and thus $rank(\bar{\mathbf{A}})=rank(\mathbf{A})+1$.
Furthermore by Lemma A.1, we know that adding a row can add at most one new
bucket. Since we assume $r_{new}$ cannot be represented as a linear
combination of $\mathbb{T}$, this means $|col(\bar{A})|=|col(A)|+1$. Thus we
have that $|col(\bar{A})|=rank(\bar{\mathbf{A}})$. Combining the fact that
$|row(\mathbf{A})|\geq|col(\mathbf{A})|$ and
$rank(\mathbf{A})=|col(\mathbf{A})|$, it follows that $\mathbf{A}$ is full
rank. ∎

###  A.4. PQ Module Proofs

####  A.4.1. Proof of Lemma 5.7

Recall the lemma states that the $L_{1}$ norm of the $\mathbb{P}$ matrix is
equal to the subtree norm of the root of the tree with marked nodes
corresponding to $\mathbb{P}$:

(12) $\mathcal{S}_{\mathbb{P}}(\mathcal{T}\text{.root})=\|\mathbb{P}\|_{1}$

###### Proof.

Given that we form $\mathbb{P}$ as shown in Example 5.2, $\|\mathbb{P}\|_{1}$
simply represents the maximum number of overlapping RCQs in $\mathbb{P}$.
Overlapping RCQs on the tree $\mathcal{T}$ must occur on the same path, that
is, they form an ancestor-descendant relationship, since the children of each
node $n$ have non-overlapping ranges. Thus, the maximum number of overlapping
RCQs across all tree paths, is equal to the $\|\mathbb{P}\|_{1}$. ∎

####  A.4.2. Proof of Lemma 5.8

The lemma states that the proactive strategy $\Delta\mathbb{P}$ generated by
generateProactiveStrategy for an input $\mathbb{P}$ satisfies the condition:

(13) $\forall\text{ paths }p\in\mathcal{T},\sum_{v\in
p}\mathcal{M}_{\mathbb{P}\cup\Delta\mathbb{P}}(v)\leq\mathcal{S}_{\mathbb{P}}(\mathcal{T}\text{.root})=\|\mathbb{P}\|_{1}$

###### Proof.

This inequality holds for $\Delta\mathbb{P}=\emptyset$ by applying Lemma 5.7
and Definition 5.6. generateProactiveStrategy only adds a node to
$\Delta\mathbb{P}$ if it satisfies the following condition on line 4, the
remaining path has a length less than $r$. At the root, this condition is set
to be less than $\|\mathbb{P}\|_{1}$. ∎

####  A.4.3. Proof of Theorem 5.3

The theorem states that, given a paid strategy matrix $\mathbb{P}$, Algorithm
6 outputs $\Delta\mathbb{P}$ such that
$\|\mathbb{P}\cup\Delta\mathbb{P}\|_{1}=\|\mathbb{P}\|_{1}$.

###### Proof.

Applying Lemma 5.7 and Definition 5.6 for the matrix
$\mathbb{P}\cup\Delta\mathbb{P}$:

(14)
$\|\mathbb{P}\cup\Delta\mathbb{P}\|_{1}=\mathcal{S}_{\mathbb{P}\cup\Delta\mathbb{P}}(\mathcal{T}\text{.root})=\max_{p\in\text{subtree}(v)}\sum_{v\in
p}\mathcal{M}_{\mathbb{P}\cup\Delta\mathbb{P}}(v)$

Rephrasing Lemma 5.8:

(15) $\max_{p\in\text{subtree}(v)}\sum_{v\in
p}\mathcal{M}_{\mathbb{P}\cup\Delta\mathbb{P}}(v)=\|\mathbb{P}\|_{1}$

Using equations 14 and 15, we get:

(16) $\|\mathbb{P}\cup\Delta\mathbb{P}\|_{1}=\|\mathbb{P}\|_{1}$

∎

##  Appendix B MMM Cache-aware tight bound

Given the cached noise parameters, we propose a theoretical upper bound
$b_{T}$ for the candidate $b_{P0}$ required to satisfy the
$(\alpha,\beta)$-accuracy guarantee.

(17)
$b_{T}\leq\frac{\sqrt{\alpha^{2}\beta/2-\|\mathbf{W}\mathbf{A}^{+}Diag(\vec{b_{\mathbf{F}}})\|_{F}^{2}}}{\|WA^{+}Diag(I_{|\mathbf{P}|})\|_{F}}$

Here, $\vec{b_{\mathbf{F}}}=[b_{1},\ldots,b_{|\mathbf{F}|}]$. In doing so, we
generalize Ge et al.’s tight bound $b_{T}$ for $\mathcal{C}=\phi$ (equation 4)
to consider cached noise parameters.

##  Appendix C Implementation Details

We implement an initial prototype of CacheDP in Python. We use the source code
provided by Ge et al. to evaluate APEx222https://github.com/cgebest/APEx.
Since the authors of Pioneer did not publish their code, we implement it from
scratch in python following the paper. We implement a simple composite plan
using all hyperparamters as described in the paper. We show in Section 8.2.1
that our implementation reproduces similar performance to the results given in
the original paper (Peng et al., 2013)[Section 7.2.1]. We make our full
evaluation including our implementation of related work publicly
available333https://git.uwaterloo.ca/m2mazmud/cachedp-public.git.

We note that pioneer uses variance based accuracy where as APEx uses
$(\alpha,\beta)$ accuracy requirements. This is not a problem for CacheDP as
we accept both types of accuracy requirement. To run Pioneer on workloads with
an $(\alpha,\beta)$ requirement, we use an MC simulation to search for the
variance that satisfies the accuracy requirement. To run APEx on workloads
with a variance based accuracy requirement we utilize a tail bound on the
Laplace distribution (since we find empirically that apex uses the Laplace
Mechanism for all single range queries) to get the alpha beta.

##  Appendix D Evaluation of the SE heuristics

In this section, we reason about the effectiveness of the strategy expander
heuristics, both experimentally and theoretically.

##### Experiments:

We first quantify the probability of SE being selected over MMM using the
results of our frequency analysis, from Table 6. We refer the reader to
Section 8 for all experimental details. Out of all paid workloads for which
MMM and SE were selected (second and fourth rows of the table), we consider
the probability for SE to be selected. SE is overwhelmingly more likely to be
chosen over MMM on both BFS and DFS tasks, with the probability of selection
being $90\%$ for BFS and $95\%$ for DFS. On the other hand, MMM is around
twice as likely to be chosen over SE on both RRQ and IDEBench tasks, as the SE
selection probability is $31\%$ for RRQ and $34\%$ for IDEBench. We conclude
that the likelihood of the SE module being chosen, and thus the success of our
heuristics, depends on the workload sequence since it influences the contents
of our cache. Averaging across all four tasks, the likelihood of the SE module
being chosen over MMM, is around $62\%$.

##### Theoretical Analysis:

We analyze the conditions under which our heuristics result in SE module being
selected. The SE module is only selected if it has a lower cost than the
original strategy in MMM. The privacy cost for our mechanisms is inversely
proportional to the error term, for a given $(\alpha,\beta)$ or
$\alpha^{2}$-expected total square error. Thus, our heuristics are only
successful if they lead to the error term for the SE module to be smaller than
the error term for MMM.

We recall Figure 6 included an example to show that expanding the strategy can
lead to an increased error term. We analyze why such situations can arise and
describe conditions for when our noise parameter-based heuristic can reduce
the error. Our heuristic to choose $b_{\ell+1}<b_{\mathbf{P}}$ strictly
improves the error under the following condition:

###### Theorem D.1.

Given a workload $\mathbf{W}$, a strategy $\mathbf{A}$ of $\ell$ rows, and a
noise vector $\mathbf{b}$, adding a new row to $\mathbf{A}$ to form
$\mathbf{A}_{e}$ reduces the error, i.e.,

(18)
$\|\mathbf{W}\mathbf{A}^{+}Diag(\mathbf{b})\|_{F}^{2}\geq\|\mathbf{W}\mathbf{A}_{e}^{+}Diag(\mathbf{b}||b_{\ell+1})\|_{F}^{2}$

if all entries in $\mathbf{b}$ equal $b^{*}$ and
$b_{\ell+1}<b_{\mathbf{P}}\leq b^{*}$, for any $b^{*}$.

###### Proof.

We recall the following theorem proved by Li et al. (Li et al., 2015) in their
MM paper.

(19)
$\|\mathbf{W}\mathbf{A}^{+}\|_{F}^{2}\geq\|\mathbf{W}\mathbf{A}_{e}^{+}\|_{F}^{2}$

For $\mathbf{b}=[b^{*}\cdots b^{*}]$, we transform Equation 18 to:

(20) $\displaystyle\|\mathbf{W}\mathbf{A}^{+}Diag(\mathbf{b})\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\|\mathbf{W}\mathbf{A}^{+}(b^{*}I)\|_{F}^{2}=(b^{*})^{2}\|\mathbf{W}\mathbf{A}^{+}\|_{F}^{2}$
(21) $\displaystyle\geq$
$\displaystyle(b^{*})^{2}\|\mathbf{W}\mathbf{A}_{e}^{+}\|_{F}^{2}$ (22)
$\displaystyle=$
$\displaystyle\|\mathbf{W}\mathbf{A}_{e}^{+}Diag(\mathbf{b}||b^{*})\|_{F}^{2}$
(23) $\displaystyle\geq$
$\displaystyle\|\mathbf{W}\mathbf{A}_{e}^{+}Diag(\mathbf{b}||b_{\ell+1})\|_{F}^{2}$

where the first inequality comes from applying Li et al.’s result (Equation
19) and the final inequality by applying the condition: $b_{\ell+1}\leq
b^{*}$. ∎

When the noise parameters are not all equal, the sufficiency conditions become
more complicated. For example, if we modify the above counterexample to have a
slightly more accurate expanded row ($b_{\ell+1}=3b<4b$), we get a lower error
than MMM:

$\|\mathbf{W}_{1}\mathbf{A}_{1}^{+}diag(\mathbf{b}_{1})\|=28b^{2},\
\|\mathbf{W}_{1}\mathbf{A}_{1e}^{+}diag(\mathbf{b}_{1e})\|=26.5b^{2}$

Thus our greedy heuristic, which does not simply add entries less than
$b_{p}$, but prioritizes the most accurate rows first, would be effective in
this scenario. We also observe that distribution of the existing noise
parameters in $\mathbf{b}$ is important. If the noise parameters are all close
to each other, strategy expansion is likely to be more beneficial. For example
changing the second entry of $\mathbf{b}_{1}$ from $b$ to $2b$ (i.e.
$\mathbf{b}_{1e}=[b,2b,5b,4b]$), results in a smaller error than MMM:

$\|\mathbf{W}_{1}\mathbf{A}_{1}^{+}diag(\mathbf{b}_{1})\|=31.0b^{2},\
\|\mathbf{W}_{1}\mathbf{A}_{1e}^{+}diag(\mathbf{b}_{1e})\|=30.2b^{2}$

It is evident that both the distribution of the noise parameters in
$\mathbf{b}$ and the new noise parameter $b_{\ell+1}$, influence the success
of our accuracy heuristic. Additionally, we observe that a newly added row
could satisfy our final heuristic by being a parent or a child of _any_ of the
existing rows in $\mathbf{A}_{1}$. We find that the structure of the expanded
strategy $\mathbf{A}_{1e}$ also significantly influences the error term. For
example, changing the final row in $\mathbf{A}_{1e}$ from $(1,1,1)$ to
$(1,0,1)$ or $(0,1,1)$ also reduces the error for $\mathbf{A}_{1e}$ to
$26.5b^{2}$ and $20b^{2}$ respectively, under the exact same noise parameters
as our counterexample. However, changing the final row to $(1,1,0)$ increases
the error to $34.7b^{2}$.

##### Conclusion:

We observe that whether an additional row selected by our noise parameter
heuristics will succeed in decreasing the error for SE over MMM, depends on
how the new row changes elements in $\mathbf{W}\mathbf{A}_{e}^{+}$. However,
we can guarantee that when the noise parameters in $\mathbf{b}$ are
sufficiently similar, strategy expander will reduce the error regardless of
the structure. Our structure-based heuristic only allows strategy queries
related through a parent-child relationship to an existing strategy query in
$\mathbf{A}$. Future research may model the success of this heuristic, by
analyzing the relation between $\mathbf{W}\mathbf{A}^{+}$ and
$\mathbf{W}\mathbf{A}_{e}^{+}$ for $\mathbf{A}$ and $\mathbf{A}_{e}$ that
differ by a parent or child row.


# Van der Waals Heterostructure Polaritons with Moiré -Induced Nonlinearity

Long Zhang1,2, Fengcheng Wu3, Shaocong Hou4, Zhe Zhang1, Yu-Hsun Chou5, Kenji
Watanabe6, Takashi Taniguchi7, Stephen R. Forrest1,4 Hui Deng1<EMAIL_ADDRESS>1 Physics Department, University of Michigan, 450 Church Street, Ann Arbor, MI
48109-2122, USA 2 Department of Physics, Xiamen University, Xiamen 361005,
Fujian, China 3 Condensed Matter Theory Center and Joint Quantum Institute,
Department of Physics, University of Maryland, College Park, Maryland 20742,
USA 4 Department of Electrical Engineering and Computer Science, University of
Michigan, 450 Church Street, Ann Arbor, MI 48109-1040, USA 5 Department of
Photonics, National Cheng Kung University, No. 1, University Road, East
District, Tainan City 70101, Taiwan, ROC 6 Research Center for Functional
Materials, National Institute for Materials Science, Tsukuba, Japan 7
International Center for Materials Nanoarchitectonics, National Institute for
Materials Science, Tsukuba, Japan

Controlling matter-light interactions with cavities is of fundamental
importance in modern science and technology walther_Cavity_2006 . It is
exemplified in the strong-coupling regime, where matter-light hybrid modes
form, with properties controllable via the photon component on the optical-
wavelength scale deng_exciton-polariton_2010 ; schneider_exciton-
polariton_2016 . In contrast, matter excitations on the nanometer scale are
harder to access. In two-dimensional van der Waals heterostructures, a tunable
moiré lattice potential for electronic excitations may form
Bistritzer_Moire_2011 , enabling correlated electron gases in lattice
potentials dean_Hofstadter_2013 ; cao_Unconventional_2018 ;
tang_Simulation_2020 ; regan_Mott_2020a ; shimazaki_Strongly_2020 . Excitons
confined in moiré lattices wu_topological_2017 ; yu_moire_2017 have also been
reported, but cooperative effects have been elusive and interactions with
light have remained perturbative jin_Observation_2019 ; tran_Evidence_2019 ;
seyler_Signatures_2019 ; alexeev_Resonantly_2019a . Here, integrating MoSe2
-WS2 heterobilayers in a microcavity, we establish cooperative coupling
between moiré -lattice excitons and microcavity photons up to liquid-nitrogen
temperature, thereby integrating into one platform versatile control over both
matter and light. The density dependence of the moiré polaritons reveals
strong nonlinearity due to exciton blockade, suppressed exciton energy shift,
and suppressed excitation-induced dephasing, all of which are consistent with
the quantum-confined nature of the moiré excitons. Such a moiré polariton
system combines strong nonlinearity and microscopic-scale tuning of matter
excitations with the power of cavity engineering and long range coherence of
light, providing a new platform for collective phenomena from tunable arrays
of quantum emitters.

##  I Introduction

Controlling matter-light interactions has mostly been implemented with either
microscopic, individual quantum particles, or macroscopic ensembles of free
particles that often can be modelled classically. Bridging the two limits to
realize collective coupling among arrays of quantum particles ushers new
paradigms in quantum many-body physics and quantum simulation research, such
as been pursued using atomic optical lattices and cavities ritsch_cold_2013 ;
gonzalez-tudela_subwavelength_2015 ; leonard_Supersolid_2017 . In solid state
systems, however, such a pursuit has been exceedingly difficult, for
nonlinearity is often too weak, due to screened Coulomb interactions, and
inhomogeneity, too large. In this work, we show that moiré lattices formed in
hetero-bilayer (hBL) transition metal dichalcogenides (TMDCs) crystals may
provide a platform that overcomes these limitations of conventional solids.

When two monolayer (ML) crystals are placed together, a moiré lattice can form
due to a slight mismatch in the lattice constants or crystal orientations of
the MLs Bistritzer_Moire_2011 ; wu_topological_2017 ; yu_moire_2017 . Its
period is tunable by the twist angle between the ML crystals from a few to a
few tens of nanometers – comparable to the range of Coulomb interactions in
semiconductors. The natural formation of the moiré lattice promises the
prospect of uniformity among different moiré cells across the lattice, thereby
potentially enabling a new solid-state platform for cooperative phenomena
between light and arrays of nonlinear quantum particles. While exciton states
induced by the moiré lattice have been reported jin_Observation_2019 ;
tran_Evidence_2019 ; seyler_Signatures_2019 ; alexeev_Resonantly_2019a , the
fundamental question remains whether a sufficiently uniform lattice of
quantized excitations can be formed in a moiré system. Exciton-light
interactions have remained in the perturbative regime, and cooperative
phenomena have not been reported.

##  II Moiré polaritons in WS2 -MoSe2 hetero-bilayers

To enable cooperative coupling between moiré excitons and photons, we use WS2
-MoSe2 hBLs enclosed in a microcavity (Fig. 1a). The hBL is encapsulated in
hexagon boron nitride (hBN) with a twist angle measured to be $56.5^{\circ}\pm
0.9^{\circ}$ (Extended Data Fig.1). The WS2 -MoSe2 hBL is unique among
commonly studied TMDCs hBLs in that the two lowest-energy moiré exciton modes
have large oscillator strengths, inherited from that of the ML MoSe2 A
exciton.alexeev_Resonantly_2019a ; zhang_Twistangle_2020 ; tang_Tuning_2020 .
This not only allows ready identification of moiré excitons via the absorption
spectra, but also suggests both of the moiré excitons may strongly couple with
light and form stable moiré polaritons.

We first identify the existence of moiré lattice and moiré excitons before
enclosing the hBL in a cavity. As shown by comparison of the hBL and ML
reflection spectra in Fig. 1c, the ML MoSe2 A exciton is split into two moiré
excitons in the hBL, X1 and X2, separated by about 40 meV, both with
pronounced absorption, consistent with calculations (Fig. 1b-c) and
temperature dependence of the hBL (Extended Data Fig.2).

When the hBL is placed at the anti-node of a $\lambda$/2 microcavity (Fig. 1d,
see Extended Data Fig.3 for more details of the cavity), in place of the moiré
excitons or the bare cavity, three dispersive modes are observed up to 70 K,
as shown in Fig. 2a-b. The modes anti-cross at the two moiré exciton
resonances, showing clearly strong coupling of both moiré excitons (X1 and X2)
with the cavity photon. Emission, temperature dependence, and time resolved
studies are provided in the Extended Data Fig. 4-6.

The measured dispersions (Fig. 2d-e) are described very well by calculations
based on a three coupled oscillator model with a Hamiltonian:

$H=\begin{bmatrix}E_{X_{1}}&0&\Omega_{1}\\\ 0&E_{X_{2}}&\Omega_{2}\\\
\Omega_{1}&\Omega_{2}&E_{c}\end{bmatrix}.$

Here $E_{c}$ is the energy of the cavity mode, $E_{X_{1},X_{2}}$ are the
energies of the two moiré excitons, and $\Omega_{1,2}$ are their coupling
strengths with the cavity photon. The three distinct dispersions measured
correspond to the three new, light-matter hybrid eigen-modes of the
Hamiltonian: the upper (UP), middle (MP) and lower (LP) polaritons. The fitted
$E_{X_{1},X_{2}}$ agree with the independently measured moiré excitons
energies (Fig. 1c) within 5 meV, where the difference is commonly observed as
a result of strain due to deposition of the top mirror. The fitted
$\Omega_{1,2}=10.1\pm 0.3$ meV and $8.5\pm 0.3$ meV at 4 K and change slightly
to $9.6\pm 0.4$ meV and $8.7\pm 0.4$ meV at 70 K. From independent
measurement, the exciton half linewidths increase from
$\gamma_{X_{1},X_{2}}=7.5$ meV and 8.4 meV at 4 K to 8.8 meV and 8.9 meV at 70
K (Extended Data Fig.5), which are mainly due to inhomogeneous broadening
111The oscillator strengths of the moiré and ML excitons correspond to
radiative linewidths of 100s $\mu$eV and radiative decay rates of 100s fs.
Therefore the measured linewidth of both types of excitons are still dominated
by inhomogeneous broadening. The measured photoluminescence decay times are a
few picoseconds for both moiré and ML excitons (Extended Data Fig.6), which
are consistent with the expected very short radiative lifetime and are
lengthened due to energy relaxation dynamics. The moiré excitons have broader
linewidths (half width of $7-8$ meV) than the ML excitons (about 2 meV),
likely due to inhomogeneity in lattice alignment or strain distribution
introduced during transfer and stacking of the two MLs.. The strong coupling
condition $\Omega_{1,2}>(\gamma_{c}+\gamma_{X_{1,2}})/2$ is satisfied up to 70
K.

In comparison, a similar cavity enclosing a ML MoSe2 also exhibits clearly
strong coupling of the ML MoSe2 A-exciton and the cavity photon (Fig. 2c and
f). Fitting the dispersions of the polaritons yield $\Omega_{XA}$=17.1$\pm$0.1
meV, which is comparable to but slightly larger than
$\sqrt{\Omega_{1}^{2}+\Omega_{2}^{2}}$=13.2 meV. This confirms that most of
the oscillator strength of the ML MoSe2 A-exciton is distributed to the moiré
states X1 and X2. The reduction may be due to additional, higher-energy moiré
states alexeev_Resonantly_2019a ; tang_Simulation_2020 or additional disorder
introduced into the bilayer during stacking of the MLs.

##  III Zero-dimensionality of moiré excitons

With robust polariton modes formed with both moiré hBL and ML excitons, we
study the effect of the underlying moiré lattice via their nonlinear response
to the excitation density $n$. We focus on $n<n_{Mott}=1/a_{B}^{2}\sim
10^{6}~{}\mu$m-2, for $n_{Mott}$ the Mott density and $a_{B}\sim 1$ nm the
Bohr radius chernikov_Population_2015 , and vary $n$ from $10$ $\mu$m-2 up to
$3\times 10^{4}$ $\mu$m-2 using a resonant 150 fs pulsed laser (see Methods
for details) scuri_Large_2018a ; emmanuele_Highly_2020 ; gu_enhanced_2019 ;
kravtsov_Nonlinear_2019 . A few representative reflectance spectra at
different excitation densities are shown in Fig. 3a-b for the hBL- and ML-
cavities, respectively. Fitting the absorption dips in the spectra with
Lorentzian functions, we obtain the polariton energies as plotted in Fig.
3c-d.

With increasing excitation densities, the moiré LP and MP shift symmetrically
toward the moiré exciton resonance $E_{X_{1}}$ while their linewidths remain
the same (Fig. 3a and Fig. 3c). These suggest reduced exciton-photon coupling
strength, yet constant exciton energy or dephasing, which are typical
properties of 0D excitons. In sharp contrast, the ML LP and UP shift together
to higher energies, accompanied by significant linewidth broadening. These
suggest a much weaker saturation of the exciton-photon coupling strength but
pronounced many-body effects, as expected of 2D excitons.

To analyze the density dependence quantitatively, we use the coupled-
oscillator model to extract from the polariton spectra the three basic
properties of excitons: the exciton energy $E_{X}(n)$, linewidth
$\gamma_{X}(n)$, and photon coupling strength $\Omega(n)$ (see Methods for
details). The results can be compared with well-established models for free 2D
excitons:

$\displaystyle E_{X}(n)$ $\displaystyle=E_{X}(0)+\beta_{1}n-\beta_{2}n^{2},$
(1) $\displaystyle\gamma_{X,n}$ $\displaystyle=\gamma_{X}(0)+\alpha n,$ (2)
$\displaystyle\Omega(n)$
$\displaystyle={\Omega(0)}/{\sqrt{1+\frac{n}{n_{s}}}}.$ (3)

These equations describe the energy shift due to two and three-particle
exchange interactions with coefficients $\beta_{1,2}$emmanuele_Highly_2020 ,
the linewidth broadening due to exciton-induced dephasing (EID) with
coefficient $\alpha$ moody_Intrinsic_2015a , and oscillator strength
saturation due to Pauli blocking with a saturation density $n_{s}$
chernikov_Population_2015 ; Huang_Carrier_1990 .

Pronounced differences between moiré and ML-excitons are clearly seen in all
three properties (Fig. 4a-c). For the ML exciton, all three properties are
described very well by Eqs. 1-3 for 2D excitons (blue diamonds and lines in
Fig. 4a-c). The exciton energy blueshifts by 1.5 meV due to exchange
interactions, the linewidth broadens significantly by 3 meV due to EID. The
fitted coefficients, $\beta_{1,ML}=(1.2\pm 0.1)\times
10^{-1}~{}\mu$eV$\cdot\mu$m2, $\beta_{2,ML}=(2.9\pm 0.4)\times
10^{-6}~{}\mu$eV$\cdot\mu$m4, and $\alpha=0.11\pm 0.01~{}\mu$eV$\cdot\mu$m2,
all agree well with reported values emmanuele_Highly_2020 ;
moody_Intrinsic_2015a . The coupling strength decreases slightly by up to
$5\%$; the corresponding $n_{s,ML}=(2.8\pm 0.4)\times 10^{5}~{}\mu$m-2 is
comparable to $1/a_{B}^{2}\sim 10^{6}~{}\mu$m-2. These results confirm that
the ML excitons behave as 2D particles; they also provide a consistency check
of our density calibration.

In contrast to the ML excitons, the moiré excitons show no measurable energy
shift, a much smaller line-broadening of $<1$ meV, and a much stronger
saturation of the coupling strength by up to $20\%$ (Fig. 4a-c, red circles).
These can not be explained with the 2D exciton picture, but are expected of 0D
excitons, as we discuss below.

In a moiré lattice, the exciton center-of-mass wavefunction is no longer a
plane wave, but becomes localized near the potential minimum of each moiré
cell, with a localization length $\ell$ ($a_{B}<\ell<a_{M}/2$, where $a_{M}$
is the moiré period), as illustrated in Fig. 4d. Evaluating $\ell$ for our
device of $a_{M}=4.2$ nm yields $\ell=1.2$ nm. Therefore the confinement leads
to an increase of the effective local density $n_{eff}$ by
$(a_{M}/2\ell)^{2}\sim 3$ at the potential minima, and thus enhanced exchange
and dipole-dipole interactions. If the moiré exciton are 2D-like band
excitons, the enhanced interactions should lead to correspondingly an enhanced
energy shift and enhanced EID. On the contrary, we observe no energy shift and
a much smaller line-broadening. Therefore, our observations show the moiré
excitons in our device are no longer 2D band excitons.

On the other hand, the suppressed energy shift and linewidth broadening are
characteristic of quantum dots with strong exciton blockade. Exciton blockade
takes place when the interaction energy between two excitons in the same cell
becomes greater than the exciton linewidth, so that multiple excitations
correspond to distinct, quantized energy levels. In our hBL, the on-site
exchange interaction $U_{exc,hBL}=\frac{1}{2\pi}(\frac{a_{B}}{w})^{2}E_{b}\sim
40$ meV for binding energy $E_{b}\sim 250$ meV, and exciton wavefunction
extension $w\sim$1 nm. The dipole interaction due to the inter-layer component
is $U_{dd,hBL}\sim\frac{d}{a_{B}}U_{exc,hBL}\sim 30$ meV for dipole length
$d\sim$0.7 nm yu_moire_2017 . These values agree with a more detailed
calculation (Extended Data Fig.9) and are much larger than the exciton full
linewidth of $2\gamma_{X_{1}}\sim 15$ meV. So we indeed expect exciton
blockade in a moiré cell. At the same time, since both the exchange and
dipole-dipole interactions decrease quickly with distance, they are both
suppressed for excitons in different moiré cells (Extended data Fig.9). The
intra-cell exciton blockade, together with suppressed inter-cell interactions,
lead to suppressed many-body effects for the single-exciton resonance, which
manifests as absence of energy shift or EID, in agreement with our
observations.

Consistent with exciton blockade, exciton-photon coupling saturates at one
exciton per moiré cell, or, $n_{s}\sim 1/a_{M}^{2}\sim(\pi
a_{B}^{2}/a_{M}^{2})n_{s,ML}\sim 6\times 10^{4}$, in excellent agreement with
the fitted $n_{s,hBL}=(6.5\pm 0.3)\times 10^{4}~{}\mu$m-2 for $n\geq 1000$
$\mu$m-2 (bottom solid line in Fig. 4c). In contrast, if the moiré excitons
are 2D band excitons, $n_{s}$ would have remained the same as $n_{s,ML}$,
since the total oscillator strength is conserved in the band across the
lattice. We note that the fit does not explain the abnormally strong
saturation at very low densities of $n<1000$ $\mu$m-2, which we will discuss
more later.

##  IV moiré lattice-induced nonlinearity of moiré polaritons

The distinctly different density dependence of the moiré polaritons shown
above suggests the possibility of achieving a much higher polariton
nonlinearity, even for polariton modes that are stable at high temperatures.
Fig. 4e shows the measured $g(n)=\mid{dE(n)/dn}\mid$ for both the moiré and ML
LPs (symbols).

While the nonlinearity increases with decreasing density for both moiré and ML
LPs, the moiré LPs show surprisingly a larger nonlinearity. The nonlinearity
of ML LPs originates primarily from exciton exchange-interactions and $g_{ML-
LP}$ saturates below $n\sim 1000~{}\mu$m-2 to $0.02~{}\mu$eV$\cdot\mu$m2
($0.04~{}\mu$eV$\cdot\mu$m2) at the measured (zero) detuning, in agreement
with reported values emmanuele_Highly_2020 ; kravtsov_Nonlinear_2019 ;
tan_Interacting_2020 . The nonlinearity of moiré LPs results primarily from
exciton blockade. Based on the data at $n>1000~{}\mu$m-2, we obtain $g_{hBL-
LP}$ about four times higher than $g_{ML-LP}$.

At very low densities of $n<1000~{}\mu$m-2, while the ML polaritons or
excitons show no measurable shift, the moiré polaritons show clearly
saturation-induced shifts down to $n\sim 10~{}\mu$m-2 (Fig. 3c), corresponding
to a stronger exciton saturation than expected from exciton blockade, exciton-
interactions, or effects of trions and defect states (Fig. 4c) 222The strong
saturation observed at very low excitation densities cannot be explained by
trions or states in deeper trapping potentials. These states would have a
small initial oscillator strength and lower resonance energy, so the increase
of saturation density would have been accompanied by an increase of the
coupling strength and exciton energy, in contradiction with our observations..
However, this abnormally large nonlinearity is repeatable over multiple
measurements in multiple devices (see Extended Data Fig.7 for another
example), suggesting hidden mechanisms for large polariton nonlinearity in
moiré lattices. Phenomenologically, $\Omega_{X1}(n)$ over the full density
range can be described very well by Eq. 3 if $n_{s}$ is replaced by a density-
dependent effective saturation density $n_{eff,s}=\tilde{n}_{s}(1-Ae^{-n/B})$,
for fitted $\tilde{n}_{s}=(4.4\pm 0.1)\times 10^{4}~{}\mu$m-2, $A=0.98\pm
0.04$ and $B=(7.4\pm 0.5)~{}\mu$m-2 (Fig. 4c).

The polariton nonlinearity is a key figure of merit that distinguishes
polariton systems from pure photon systems. Together with the robust coherence
enforced by the photon component, it gives rise to a wide range of novel
nonlinear many-body and quantum phenomena ballarini_Alloptical_2013 ;
walker_Ultralowpower_2015a ; berloff_realizing_2017 ; delteil_Polariton_2019 ;
munoz-matutano_Emergence_2019 ; kim_Emergence_2020 .

For polaritons made of 2D excitons, however, larger $g$ is obtained only at
the compromise of $\Omega$, or the stability of polariton modes. This is
because the exciton-exchange interaction and the exciton-photon coupling
strength saturation, the two main contributions to $g$, both decrease with the
exciton reduced mass $\mu$: $g_{exc}\sim E_{b}a_{B}^{2}\propto 1/\mu$ and
$g_{sat}\propto a_{B}^{2}\propto 1/\mu^{2}$. Yet the exciton-photon coupling
strength $\Omega\propto\sqrt{f}\propto 1/a_{B}\propto\mu$. The highest
polariton nonlinearity is achieved in single or few quantum-well GaAs
microcavities delteil_Polariton_2019 ; munoz-matutano_Emergence_2019 , with
$g\sim 3~{}\mu$eV$\cdot\mu$m2 and $\Omega$ only $\sim$3 meV. Wide-bandgap
semiconductors, organic crystals and ML TMDCs all feature greater $\Omega$ and
high-temperature polaritons, yet a much weaker polariton nonlinearity
daskalakis_Nonlinear_2014a ; scuri_Large_2018a ; barachati_Interacting_2018 ;
kravtsov_Nonlinear_2019 . Higher order excitations in TMDCs, such as trions
with strong band-filling effect emmanuele_Highly_2020 ; tan_Interacting_2020
and 2s excitons with a larger Bohr radius gu_enhanced_2019 , have shown
enhanced nonlinearity, but again at the compromise of stability.

The moiré polariton system uniquely combines strong nonlinearity, due to
quantum confinement of excitons within each moiré cell, and a large total
photon coupling strength, due to collective coupling among the cells. It
thereby provides a new route to achieve strong nonlinearities simultaneously
with robust exciton-photon coupling.

##  V Conclusion

In summary, we have demonstrated the first polariton system formed via
collective coupling of a 2D lattice of 0D excitons with light in a
microcavity. The system therefore introduces quantum-dot like nonlinearity
into cooperatively coupled solid state system, opening a door to novel quantum
many-body physics and polariton devices ritsch_cold_2013 ; gonzalez-
tudela_subwavelength_2015 ; leonard_Supersolid_2017 ; yu_Electrically_2020 .
Polariton blockade and a strongly-correlated polariton gases may become
possible with reduced inhomogeneous broadening of moiré excitons, improved
cavity quality factors, and a better understanding of the abnormal enhancement
of the polariton nonlinearity at very low excitation densities. Electrical
gating and electrical field tuning of the heterobilayer can be implemented to
further control the nonlinearity and many-body phenomena in the moiré
polariton system.

## Methods

Sample fabrication. The ML MoSe2 , WS2 and few layer hBN flakes were obtained
by mechanical exfoliation from bulk crystals. A polyethylene terephthalate
(PET) stamp was used to pick up the top hBN, WS2 ML, MoSe2 ML, and the bottom
hBN under microscope. The WS2 and MoSe2 monolayer MLs were rotationally
aligned to about $0^{\circ}$ in the heterobilayer.

The bottom half of the cavity consists of 16 pairs of SiO2/TiO2 DBR with a
$\lambda/4$ SiO2 layer. The heterostructure on the PET stamp was stamped onto
the bottom half of the cavity and the PET was dissolved in dichloromethane for
six hours at room temperature. Then 67 nm PMMA film was spin coated on the top
of heterobilayer to form the second half of the $\lambda$/2 cavity. Then a
silver film of 40 nm thick was deposited using electron beam evaporation as
the top mirror of the cavity.

Optical measurements. For low temperature measurements, the sample was kept in
a 4 K cryostat (Montana Instrument). The excitation and collection were
carried out with a home-built confocal microscope with an objective lens with
a numerical aperture (NA) of 0.42. To characterize the dispersion of the
polariton device, we performed angle-resolved reflection measurement using
white light from a tungsten halogen lamp. The white light was focused on the
sample with a beam size of 10 $\mu$m in diameter. To perform power dependent
reflection measurements, we used a 150 fs-pulse laser as the light source with
a repetition rate of 80 MHz and focused beam size on the sample of around 1.5
$\mu$m in diameter. For the hBL, the laser was centered between LP and MP to
simultaneously measure both modes, while the intensity at UP is negligible
(Extended Data Fig.8 for an example spectrum of the laser). For the ML, due to
the larger energy separation of the LP and UP, they were measured separately
with the laser centered at each, respectively. The signal was detected using a
Princeton Instruments spectrometer with a cooled charge-coupled camera.

Theory of moiré excitons The Hamiltonian for the moiré excitons is

$\mathcal{H}=\begin{pmatrix}E_{G}+\frac{\hbar^{2}\bm{k}^{2}}{2M_{\text{X}}}&w(1+e^{-i\bm{b}_{+}\cdot\bm{r}}+e^{-i\bm{b}_{-}\cdot\bm{r}})\\\
w(1+e^{i\bm{b}_{+}\cdot\bm{r}}+e^{i\bm{b}_{-}\cdot\bm{r}})&E_{G}+\delta_{0}+\frac{\hbar^{2}(\bm{k}-\bm{\kappa})^{2}}{2M_{\text{IX}}}\end{pmatrix}$
(4)

where $E_{G}$ is an energy constant, $M_{\text{X}}$ and $M_{\text{IX}}$ are
respectively the effective masses for intralayer and interlayer excitons,
$\bm{\kappa}=(0,4\pi/(3a_{M}))$ accounts for the momentum mismatch between the
two excitons, and $\delta_{0}$ is an energy offset. The off-diagonal terms are
derived from interlayer tunneling, and $\bm{b}_{\pm}=4\pi/(\sqrt{3}a_{M})(\pm
1/2,\sqrt{3}/2)$. The moiré period is
$a_{M}=a_{0}/\sqrt{\delta^{2}+\theta^{2}}$, where
$a_{0}=(a_{\text{MoSe}_{2}}+a_{\text{WS}_{2}})/2=3.26$ Åand
$\delta=|a_{\text{MoSe}_{2}}-a_{\text{WS}_{2}}|/a_{0}\approx 4$%. We use the
following parameter values, $E_{G}=1.61$eV, $\delta_{0}=-20$ meV,$w=14$ meV,
$M_{\text{X}}=0.82m_{0}$, $M_{\text{IX}}=0.71m_{0}$, where $m_{0}$ is the
electron bare mass. The energy spectrum of the moiré excitons is obtained by
diagonalizing the moiré Hamiltonian $\mathcal{H}$ using plane-wave expansion,
and is shown in Fig. 1b. For the lowest-energy exciton X1, its interlayer
component is plotted in Fig. 4d, which shows spatial localization.

Polariton density calibration In this section, we will use the bilayer device
data as an example; the same procedure applies to the monolayer device. To
extract the polariton density: First, we fit the reflection spectra in Fig. 3
using :

$\displaystyle R=1-Absor_{LP}-Absor_{MP},$ (5)

where $Absor_{LP}$ and $Absor_{MP}$ represent absorption by the lower
polariton ($LP$) and middle polariton ($MP$) in the bilayer, respectively, and
they are described by Lorentzian functions:

$\displaystyle
Absor_{LP,MP}(E)=\frac{A_{LP,MP}}{(E-E_{LP,MP})^{2}+\gamma_{LP,MP}^{2}},$ (6)

Resonance energy $E_{LP,MP}$, linewidth (HWHM) $\gamma_{LP,MP}$, and
absorption amplitude $A_{LP,MP}$ are fitting parameters with uncertainties
$\delta_{E_{LP,MP}},\delta_{\gamma_{LP,MP}}$ and $\delta_{A_{LP,MP}}$
corresponding to $95\%$ confidence intervals.

We use pulsed laser with pulse duration of 150 fs and repetition rate $f$ of
80 MHz. We calculate the polariton density $n$ injected per pulse. The average
power of the laser is $P$. The laser profile can be fitted with Gaussian
function:

$\displaystyle
G(E)=\frac{A_{laser}}{\gamma_{laser}\sqrt{\pi/2}}e^{-\frac{2(E-E_{laser})^{2}}{\gamma_{laser}^{2}}}$
(7)

where $A_{laser}$ is the area of the Gaussian function, and $\gamma_{laser}$
is the linewidth. The power absorbed by LP and MP, $P_{LP,MP}$ can be
calculated using the convolution:

$\displaystyle P_{LP,MP}=\frac{P\int G(E)Absor_{LP,MP}}{A_{laser}}=\frac{P\int
G(E)\frac{A_{LP,MP}}{4(E-E_{LP,MP})^{2}+\gamma^{2}_{LP,MP}}}{A_{laser}}$ (8)

The total carrier density including both LP and MP created per pulse and its
error $\delta_{n}$ can be calculated by:

$\displaystyle
n=(n_{LP}+n_{MP})/A_{beam}=(\frac{P_{LP}}{fE_{LP}}+\frac{P_{MP}}{fE_{MP}})/A_{beam}.$
(9)
$\displaystyle\delta_{n}=\sqrt{\delta^{2}_{n_{LP}}+\delta^{2}_{n_{MP}}}/A_{beam}$
(10)

Where, Abeam=1.5 $\mu$m, is the beam area, and
$\delta_{n_{LP,MP}}=n_{LP,MP}\sqrt{(\frac{\delta_{P_{LP,MP}}}{P_{LP,MP}})^{2}+(\frac{\delta_{E_{LP,MP}}}{E_{LP,MP}})^{2}}$

Extraction of exciton energy, linewidth, and coupling strength at different
densities To analyze the results quantitatively, we extract the density
dependence of the exciton properties from the measured polariton spectra. In
the heterobilayer cavity, we focus on the MP and LP modes and neglect changes
caused by the X2, since X2 is at a much higher energy, and its change only
negligibly affects MP and LP (see Extended Fig.10 for details). The cavity
resonance energy $E_{c}$ and linewidth $\gamma_{c}$ are assumed to change
negligibly with exciton density. Using the two coupled oscillator mode, the
energies of the MP and LP of the heterobilayer cavity, as well as of the UP
and LP of the ML-cavity can be extracted. In the following, we will use the
bilayer device data as an example; the same procedure applies to the monolayer
device.

$\displaystyle E_{\textrm{LP,MP}}(n)=$
$\displaystyle\frac{1}{2}\bigg{[}E_{X}(n)+E_{C}+i(\gamma_{C}+\gamma_{X}(n))\bigg{]}$
$\displaystyle\pm$
$\displaystyle\sqrt{\Omega(n)^{2}+\frac{1}{4}\bigg{[}E_{X}(n)-E_{C}+i(\gamma_{C}-\gamma_{X}(n))\bigg{]}^{2}}$
(11)

Here the cavity resonance, $E_{C}$, and cavity half-linewidth $\gamma_{C}$, do
not change with carrier density. Therefore, from the measured polariton
energies and half-linewidth, $\gamma_{LP,MP}$, we obtain the density
dependence of exciton energy, $E_{X}$, half-linewidth, $\gamma_{X}$, and
exciton-photon coupling strength, $\Omega$:

Exciton energy $E_{X}$ and its uncertainty $\delta_{E_{X}}$ :

$\displaystyle E_{X}(n)=E_{LP}(n)+E_{MP}(n)-E_{C}$ (12)
$\displaystyle\delta_{E_{X}}=\sqrt{\delta_{E_{LP}}^{2}+\delta_{E_{MP}}^{2}+\delta_{E_{C}}^{2}}$
(13)

where, $E_{C}$ and $\delta_{E_{C}}$ are the cavity resonance and its
uncertainty obtained by fitting the angle resolved reflection spectrum in Fig.
3, and they do not change with polariton densities.

Exciton linewidth $\gamma_{X}$ and its uncertainty $\delta_{\gamma_{X}}$:

$\displaystyle\gamma_{X}(n)=\gamma_{LP}(n)+\gamma_{MP}(n)-\gamma_{C}$ (14)
$\displaystyle\delta_{\gamma_{X}}=\sqrt{\delta_{\gamma_{LP}}^{2}+\delta_{\gamma_{MP}}^{2}+\delta_{\gamma_{C}}^{2}}$
(15)

where, $\gamma_{C}$ and $\delta_{\gamma_{C}}$ are the cavity linewidth and
uncertainty obtained by fitting the reflection spectrum of the bare cavity in
Fig. 1d, and they do not change with polariton densities.

Coupling strength $\Omega$:

$\displaystyle\Omega(n)=\frac{1}{2}\sqrt{\left[E_{LP}(n)-E_{MP}(n)+i(\gamma_{LP}(n)-\gamma_{MP}(n))\right]^{2}-\left[E_{X}(n)-E_{C}+i(\gamma_{X}(n)-\gamma_{C})\right]^{2}}$
(16)

Dipole-dipole interaction The interlayer component of moiré excitons
contributes to dipole-dipole interactions. The wave function for interlayer
component of a moiré exciton localized near the potential minimum around the
origin can be described as

$W(\bm{r}_{1},\bm{r}_{2})=x_{\text{IX}}\Phi(\frac{\bm{r}_{1}+\bm{r}_{2}}{2})\phi(\bm{r}_{1}-\bm{r}_{2}),$
(17)

where $x_{\text{IX}}$ is the interlayer component weight, $\Phi$ and $\phi$
are respectively the center-of-mass and relative-mass wave function. We
approximate $\Phi$ by a Gaussian
$\Phi(\bm{R})=\frac{1}{\sqrt{\pi}\ell}e^{-R^{2}/(2\ell^{2})}$, and
$\phi(\bm{r})=\sqrt{8/\pi}(1/a_{B})e^{-2r/a_{B}}$, where
$\bm{R}=\frac{\bm{r}_{1}+\bm{r}_{2}}{2}$, $\bm{r}=\bm{r}_{1}-\bm{r}_{2}$,
$\ell$ is the localization length that can be estimated using the exciton
state shown in Fig. 4d, and $a_{B}$ is the Bohr radius.

The dipole-dipole interaction between two excitons localized at different
moiré sites can be approximated as

$\displaystyle U_{dd}(\bm{L})=$ $\displaystyle\int d\bm{r}_{1}\int
d\bm{r}_{2}\int d\bm{r}_{3}\int
d\bm{r}_{4}|W(\bm{r}_{1},\bm{r}_{2})|^{2}|W(\bm{r}_{3}-\bm{L},\bm{r}_{4}-\bm{L})|^{2}$
(18)
$\displaystyle[V_{S}(\bm{r}_{1}-\bm{r}_{3})+V_{S}(\bm{r}_{2}-\bm{r}_{4})-V_{D}(\bm{r}_{1}-\bm{r}_{4})-V_{D}(\bm{r}_{2}-\bm{r}_{3})]$
$\displaystyle=$ $\displaystyle
2x_{\text{IX}}^{2}\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}V(\bm{q})\frac{\exp(-q^{2}\ell^{2}/2)}{[1+\frac{1}{64}a_{B}^{*2}q^{2}]^{3}}e^{-i\bm{q}\cdot\bm{L}},$

where $\bm{L}$ is a moiré lattice vector, $V_{S}(\bm{r})=e^{2}/(\epsilon r)$
is the intralayer Coulomb interaction,
$V_{D}(\bm{r})=e^{2}/(\epsilon\sqrt{r^{2}+d^{2}})$ is the interlayer Coulomb
interaction with $d$ the interlayer distance, and $V(\bm{q})=\frac{2\pi
e^{2}}{\epsilon q}(1-e^{-qd})$. To calculate $U_{dd}(\bm{L})$, we take
$x_{\text{IX}}\approx 1/\sqrt{2}$, $d=0.65$ nm, $\epsilon=5$, and
$a_{B}\approx 1$ nm. Extended Fig.9 shows the onsite and nearest-neighbor
dipole-dipole interactions as a function of twist angle. The onsite repulsion
is sizable ($\sim$ 24 meV for a twist angle of 56.5∘), while the offsite
repulsion is negligible.

Data availability Data are available on request from the authors.

Competing interests The authors declare that they have no competing financial
interests.

Author Contributions H.D., L.Z.conceived the experiment. L.Z. performed the
measurements. F. W. provided theoretical inputs. L.Z. and Z.Z fabricated the
device. L.Z. and H.D. performed data analysis. S.H deposited the silver
mirror. Y.C grew the bottom DBR mirror. K.W. and T.T grew hBN single crystals.
H.D. and S.F. supervised the projects. L.Z and H.D. wrote the paper with
inputs from other authors. All authors discussed the results, data analysis
and the paper.

Acknowledgment We are grateful for helpful discussions with Duncan Steel and
Mack Kira. L. Z., S. H, S. F.and H. D.acknowledge the support by the Army
Research Office under Awards W911NF-17-1-0312. L.Z. and H.D. also acknowledge
the support by the Air Force Office of Scientific Research under Awards
FA2386-18-1-4086 and by the National Science Foundation under Awards
DMR-1838412. S. F. also acknowledges support from the U.S. Department of
Energy, Office of Science, Office of Basic Energy Sciences, under Awards DE-
SC0017971. F. W. is supported by the Laboratory for Physical Sciences. Y. C.
acknowledge the support by the Young Scholar Fellowship Program by the
Ministry of Science and Technology (MOST) in Taiwan, under Grant MOST
108-2636-M-006 -010. K. W. and T. T. acknowledge support from the Elemental
Strategy Initiative conducted by the MEXT, Japan ,Grant Number
JPMXP0112101001, JSPS KAKENHI Grant Numbers JP20H00354 and the
CREST(JPMJCR15F3), JST.

## References

## References

  * (1) Walther, H., Varcoe, B. T. H., Englert, B.-G. & Becker, T.  Cavity quantum electrodynamics.  _Reports on Progress in Physics_ 69, 1325–1382 (2006).
  * (2) Deng, H., Haug, H. & Yamamoto, Y.  Exciton-polariton Bose-Einstein condensation.  _Reviews of Modern Physics_ 82, 1489 (2010).
  * (3) Schneider, C. _et al._ Exciton-polariton trapping and potential landscape engineering.  _Reports on Progress in Physics_ 80, 016503 (2016).
  * (4) Bistritzer, R. & MacDonald, A. H.  Moiré bands in twisted double-layer graphene.  _Proceedings of the National Academy of Sciences_ 108, 12233–12237 (2011).
  * (5) Dean, C. R. _et al._ Hofstadter’s butterfly and the fractal quantum Hall effect in moiré superlattices.  _Nature_ 497, 598–602 (2013).
  * (6) Cao, Y. _et al._ Unconventional superconductivity in magic-angle graphene superlattices.  _Nature_ 556, 43–50 (2018).
  * (7) Tang, Y. _et al._ Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices.  _Nature_ 579, 353–358 (2020).
  * (8) Regan, E. C. _et al._ Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices.  _Nature_ 579, 359–363 (2020).
  * (9) Shimazaki, Y. _et al._ Strongly correlated electrons and hybrid excitons in a moiré heterostructure.  _Nature_ 580, 472–477 (2020).
  * (10) Wu, F., Lovorn, T. & MacDonald, A. H.  Topological exciton bands in moiré heterojunctions.  _Physical Review Letters_ 118, 147401 (2017).
  * (11) Yu, H., Liu, G.-B., Tang, J., Xu, X. & Yao, W.  Moiré excitons: From programmable quantum emitter arrays to spin-orbit–coupled artificial lattices.  _Science Advances_ 3, e1701696 (2017).
  * (12) Jin, C. _et al._ Observation of moiré excitons in WSe2/WS2 heterostructure superlattices.  _Nature_ 567, 76–80 (2019).
  * (13) Tran, K. _et al._ Evidence for moiré excitons in van der Waals heterostructures.  _Nature_ 567, 71–75 (2019).
  * (14) Seyler, K. L. _et al._ Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers.  _Nature_ 567, 66–70 (2019).
  * (15) Alexeev, E. M. _et al._ Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures.  _Nature_ 567, 81–86 (2019).
  * (16) Ritsch, H., Domokos, P., Brennecke, F. & Esslinger, T.  Cold atoms in cavity-generated dynamical optical potentials.  _Reviews of Modern Physics_ 85, 553–601 (2013).
  * (17) González-Tudela, A., Hung, C.-L., Chang, D. E., Cirac, J. I. & Kimble, H. J.  Subwavelength vacuum lattices and atom–atom interactions in two-dimensional photonic crystals.  _Nature Photonics_ 9, 320–325 (2015).
  * (18) Léonard, J., Morales, A., Zupancic, P., Esslinger, T. & Donner, T.  Supersolid formation in a quantum gas breaking a continuous translational symmetry.  _Nature_ 543, 87–90 (2017).
  * (19) Zhang, L. _et al._ Twist-angle dependence of moiré excitons in WS2 /MoSe2 heterobilayers.  _Nature Communications_ 11, 5888 (2020).
  * (20) Tang, Y. _et al._ Tuning layer-hybridized moiré excitons by the quantum-confined Stark effect.  _Nature Nanotechnology_ 1–6 (2020).
  * (21) The oscillator strengths of the moiré+.2222emand ML excitons correspond to radiative linewidths of 100s $\mu$eV and radiative decay rates of 100s fs. Therefore the measured linewidth of both types of excitons are still dominated by inhomogeneous broadening. The measured photoluminescence decay times are a few picoseconds for both moiré+.2222emand ML excitons (Extended Data Fig.6), which are consistent with the expected very short radiative lifetime and are lengthened due to energy relaxation dynamics. The moiré+.2222emexcitons have broader linewidths (half width of $7-8$ meV) than the ML excitons (about 2 meV), likely due to inhomogeneity in lattice alignment or strain distribution introduced during transfer and stacking of the two MLs.
  * (22) Chernikov, A., Ruppert, C., Hill, H. M., Rigosi, A. F. & Heinz, T. F.  Population inversion and giant bandgap renormalization in atomically thin WS2 layers.  _Nature Photonics_ 9, 466–470 (2015).
  * (23) Scuri, G. _et al._ Large Excitonic Reflectivity of Monolayer MoSe2 Encapsulated in Hexagonal Boron Nitride.  _Phys. Rev. Lett._ 120, 037402 (2018).
  * (24) Emmanuele, R. P. A. _et al._ Highly nonlinear trion-polaritons in a monolayer semiconductor.  _Nature Communications_ 11, 3589 (2020).
  * (25) Gu, J. _et al._ Enhanced nonlinear interaction of polaritons via excitonic rydberg states in monolayer WSe2.  _arXiv:1912.12544_ (2019).
  * (26) Kravtsov, V. _et al._ Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum.  _Light: Science & Applications_ 9, 56 (2020).
  * (27) Moody, G. _et al._ Intrinsic homogeneous linewidth and broadening mechanisms of excitons in monolayer transition metal dichalcogenides.  _Nature Communications_ 6, 8315 (2015).
  * (28) Huang, D., Chyi, J.-I. & Morkoç, H.  Carrier effects on the excitonic absorption in GaAs quantum-well structures: Phase-space filling.  _Phys. Rev. B_ 42, 5147–5153 (1990).
  * (29) Tan, L. B. _et al._ Interacting Polaron-Polaritons.  _Physical Review X_ 10, 021011 (2020).
  * (30) The strong saturation observed at very low excitation densities cannot be explained by trions or states in deeper trapping potentials. These states would have a small initial oscillator strength and lower resonance energy, so the increase of saturation density would have been accompanied by an increase of the coupling strength and exciton energy, in contradiction with our observations.
  * (31) Ballarini, D. _et al._ All-optical polariton transistor.  _Nature Communications_ 4, 1778 (2013).
  * (32) Walker, P. M. _et al._ Ultra-low-power hybrid light matter solitons.  _Nature Communications_ 6, 8317 (2015).
  * (33) Berloff, N. G. _et al._ Realizing the classical XY Hamiltonian in polariton simulators.  _Nature Materials_ 16, 1120–1126 (2017).
  * (34) Delteil, A. _et al._ Towards polariton blockade of confined exciton–polaritons.  _Nature Materials_ 18, 219–222 (2019).
  * (35) Muñoz-Matutano, G. _et al._ Emergence of quantum correlations from interacting fibre-cavity polaritons.  _Nature Materials_ 18, 213–218 (2019).
  * (36) Kim, S. _et al._ Emergence of microfrequency comb via limit cycles in dissipatively coupled condensates.  _Physical Review B_ 101, 085302 (2020).
  * (37) Daskalakis, K. S., Maier, S. A., Murray, R. & Kéna-Cohen, S.  Nonlinear interactions in an organic polariton condensate.  _Nature Materials_ 13, 271–278 (2014).
  * (38) Barachati, F. _et al._ Interacting polariton fluids in a monolayer of tungsten disulfide.  _Nature Nanotechnology_ 13, 906–909 (2018).
  * (39) Yu, H. & Yao, W.  Electrically tunable topological transport of moiré polaritons.  _Science Bulletin_ 65, 1555–1562 (2020).

##

Figure 1: Illustration of the moiré polariton system and the constituent moiré
excitons and microcavity. $\mathbb{a}$, Schematic of the moiré polariton
system formed by excitons confined in moiré lattice and coupled with a planar
cavity. $\mathbb{b}$, moiré band structure for hBL excitons, where X1 and X2
are two bright moiré exciton states. See Method (Theory of moiré excitons) for
details of the calculation. $\mathbb{c}$, Top panel: theoretical optical
absorption spectrum calculated from the moiré exciton band structure. Bottom
panel: Reflection contrast spectra near the MoSe2 A exciton resonance, from
the WS2 -MoSe2 hBL (red), MoSe2 ML (blue) and WS2 ML (green). The spectra are
displaced vertically for easier reading. The MoSe2 ML A-exciton
(X${}_{\text{A}}$) splits into two well resolved moiré exciton (X1 and X2) in
the hBN. $\mathbb{d}$ Angle resolved white light reflection spectrum of the
bare cavity, measured in a region where there is no hBL, showing a single
cavity dispersion with a half linewidth $\gamma_{c}=2.7$ meV. Figure 2: Strong
coupling and dispersions of moiré and ML polaritonss. $\mathbb{a-c}$, Angle
resolved white light reflection spectra, demonstrating strong coupling between
moiré exciton and cavity photon at 5K in $\mathbb{a}$ and 70K in $\mathbb{b}$,
in comparison with strong coupling between ML exciton and cavity photon at 5K
in $\mathbb{c}$. The left/right panels show the measured/simulated results,
respectively. $\mathbb{d-f}$, The polariton energies vs. in-plane wavenumber
k// obtained from $\mathbb{a-c}$, respectively. The solid lines are fits to
the polariton dispersions with the coupled harmonic oscillator model. The
dashed lines are the fitted energies of the uncoupled cavity photon and
excitons. The error bars on the energy data correspond to the $95\%$
confidence interval of the Lorentzian fit. Figure 3: Nonlinearity of the moiré
hBL and ML polaritons. $\mathbb{a}$, Reflection spectra of moiré LP and MP at
zero detuning for different pumping densities. From bottom to top, the carrier
density increases from $17~{}\mu$m-2 to $1.8\times 10^{4}~{}\mu$m-2.
$\mathbb{b}$, Reflection spectra of MoSe2 ML LP (left) and UP (right). The
cavity is red detuned from the exciton by 42.5 meV. From bottom to top, the
carrier density increases from $22~{}\mu$m-2 to $1.5\times 10^{4}~{}\mu$m-2.
In $\mathbb{a,b}$, the colored solid lines are examples of Lorentzian fits.
The arrows are guides for the eyes. $\mathbb{c,d}$ Shift of polariton energies
vs. the carrier density (log scale) obtained from $\mathbb{a}$ and
$\mathbb{b}$, respectively. Insets show the density in the linear scale,
plotted over the same ranges of the horizontal and vertical axes as the main
figure. Solid lines are calculations using parameters obtained from fitting
the data in Fig. 4 $\mathbb{a-c}$. The error bars on the energy data
correspond to the $95\%$ confidence interval of the Lorentzian fit.  Figure 4:
Enhanced nonlinearity by moiré lattice confinement. $\mathbb{a-c}$, Shift of
exciton energy $\Delta E$, half-linewidth $\gamma$, and normalized coupling
strength $\Omega/\Omega_{0}$ of the hBL-LP (red) and ML-LP (blue) as a
function of carrier density, obtained from the data in Fig. 3. The hBL $\Delta
E_{hBL-X}$ and $\gamma_{hBL-X}$ (red circles in $\mathbb{a}$ and $\mathbb{b}$)
are approximately constant. The ML $\Delta E_{ML-X}$ and $\gamma_{ML-X}$ (blue
diamonds in $\mathbb{a}$ and $\mathbb{b}$) are fitted by a second order
polynomial and a linear line, respectively (black solid lines in $\mathbb{a}$
and $\mathbb{b}$). In $\mathbb{c}$, the black solid lines are fits with
equation (3) with a constant $n_{s}$ for the moiré excitons at n$>$1000
$\mu$$m^{-2}$ and for ML excitons. The black dashed line is a fit with a
density-dependent effective saturation density $n_{eff,s}$, and is used for
calculating the polariton energies in Fig. 3 $\mathbb{c}$. $\mathbb{d}$, Real-
space distribution of the interlayer component of the X1 state. The white
lines mark a moiré unit cell. The scale bar is 1 nm. $\mathbb{e}$, The
measured (symbols) and fitted (lines) nonlinear coefficient $g$ vs. carrier
density for the moiré hBL-LP (red) and ML-LP (blue). The magenta dashed line
and blue solid line are the calculations using the fitted polariton energies
in Fig. 3 $\mathbb{c}$ and $\mathbb{d}$. The hBL-LP is at zero detuing; the
red solid line and magenta dashed lines correspond to fitted $\Omega_{hBL}$
using a constant and effective saturation density, $n_{s}$ and $n_{eff,s}$, in
Eq (3), respectively. For ML-LP, the blue solid and dashed lines correspond to
the measured detuning and zero detuning, respectively. The error bars in
$\mathbb{a}$-$\mathbb{c}$ and density in $\mathbb{e}$ are explained in
Methods. The error bars of $g$ correspond to the $95\%$ confidence interval of
the fit using $g(n)=\mid{dE(n)/dn}\mid$.

##

Extended Data Fig. 1: Heterobilayer twist angle. Twist angle of hBL in the
main text measured by second harmonic generation (SHG) spectroscopy.
$\mathbb{a}$ The polarization-dependent SHG signal measured on the ML WS2
(green open circles) and MoSe2 (blue filling squares) regions of the hBL, and
the corresponding fits with the sinusoidal functions (green and blue solid
lines). $\mathbb{b}$, the SHG signal from ML WS2 , ML MoSe2 , and hBL regions,
measured with the same experimental configurations. The suppressed SHG signal
from hBL as a result of destructive interference indicates the stacking order
is H-stacking. The twist angle is determined to be $56.5^{\circ}\pm
0.8^{\circ}$ Extended Data Fig. 2: Temperature dependence of moiré exciton.
Temperature dependence of the moiré exciton X1 and X2 measured from a separate
heterobilayer prepared on sapphire substrate. The black dashed lines are
guides for eyes. The two exciton states can be well resolved up to 200K, which
can exclude the possibility of charged exciton or trapped exciton by defect.
Extended Data Fig. 3: Schematic of the device. $\mathbb{a}$ Schematic of the
device shows the different layers of the heterostructure embedded inside a
microcavity that consists of a bottom DBR and a top silver mirror.
$\mathbb{b}$ Microscope image of the hBL on the top of DBR mirror, taken
before depositing the PMMA layer and the silver mirror. $\mathbb{c}$ Thickness
of each layer for the device. Extended Data Fig. 4: Photoluminescence from the
moiré polariton. Left panel: Angle resolved photoluminescence spectrum of hBL
in cavity, excited by a continuous wave laser at the energy of 2.3 eV and
power of 50 $\mu$W. To enhance the visibility of states at higher energy,
emission intensity above 1.607 eV is magnified by 5 times. Right panel:
simulated angle resolved absorption, which agrees well with the measurement.
Extended Data Fig. 5: Transition from strong coupling to weak coupling driven
by thermal broadening. Strong coupling to weak-coupling transition measured by
temperature dependence of $\Omega_{1}$ (red open circles),
$(\gamma_{c}+\gamma_{X1})/2$ (purple left-triangle) in $\mathbb{a}$ and
$\Omega_{2}$ (red square), $(\gamma_{c}+\gamma_{X2})/2$ (purple right-
triangle) in $\mathbb{b}$. $\gamma_{c}=2.7meV$ is constant with temperature.
$\gamma_{X1}$ and $\gamma_{X2}$ are measured independently from bare hBL.
$\Omega_{1}$ and $\Omega_{2}$ drop below the average linewidth at about 100K,
showing the transition to the weak-coupling regime. Extended Data Fig. 6: Time
resolved photoluminescence (TRPL) of ML exciton, ML trion, hBL excitons, and
hBL trion. TRPL spectra measured by streak camera for ML MoSe2 (ML)
$\mathbb{a}$ and hBL WS2 /MoSe2 (hBL) $\mathbb{b}$. (The hBL data is collected
from a different sample from the main text, which is not integrated with
microcavity) The different resonances are labelled with white arrows including
ML exciton (ML-X), ML trion (ML-T), moiré exciton at higher energy (hBL-X2),
moiré excion at lower energy (hBL-X1), and moiré trion (hBL-T). $\mathbb{c}$
and $\mathbb{d}$ shows the time resolved decay of ML-X and hBL-X1 respectively
by integrating the spectrum in the range labeled by the red rectangles in a
and b. The red solid lines in $\mathbb{c}$ and $\mathbb{d}$ are the fits with
single exponential decay function. PL decay time for ML-X and hBL-X1 are 6.7
ps and 8.0 ps respectively. Extended Data Fig. 7: Strong nonlinearity measured
in another device. Measurement taken on a different sample shows
reproducibility of the results. $\mathbb{a}$ Angle resolved white light
reflection spectra taken at 5K on the second sample. White solid lines are the
fits using coupled oscillator model. The dashed white lines are the fitted
energies of the uncoupled cavity photon and excitons. $\mathbb{b}$ Power
dependent reflection spectra for the lower polariton (left panel) and middle
polariton (right panel).$\mathbb{c}$ Shift of polariton energies vs. the
carrier density (log scale) obtained from $\mathbb{b}$. $\mathbb{d}$ Extracted
nonlinear coefficients for lower polariton (red circles) and the calculations
using fitted polariton energies (solid line). Extended Data Fig. 8: Profile of
the laser used for nonlinearity characterization. Profile of pulsed laser (red
dot) used for the nonlinearity measurement of moiré polaritons, compared with
the moiré exciton X1 and X2 (blue dots), and moiré polaritons (black dots).
The laser has a negligibly small tail on X2 and upper polariton.

Extended Data Fig. 9: Dipole-Dipole interaction strength as a function of
twist angle $\theta$.

$U_{dd}^{(0)}$ and $U_{dd}^{(1)}$ are respectively onsite and nearest-neighbor
interaction strength.

Extended Data Fig. 10: Effect on the nonlinearity from the moiré exciton X2.
Energy shifts of $E^{BL}_{LP}$ and $E^{BL}_{MP}$ from the saturation of upper
moiré exciton X2. In the main test, we ignore the effects of moiré exciton X2
on the nonlinearity of lower and middle polaritons. To quantitatively estimate
the effects, we calculate the $\Delta E^{BL}_{LP}$ and $\Delta E^{BL}_{MP}$ as
a function of coupling strength ($\Omega_{2}$) of X2. When the $\Omega_{2}$
changes by $10\%$, the $E_{{LP}}$ and $E_{{MP}}$ will shift within 0.16 meV,
which is less than $8\%$ of the shift observed from the experiments (2meV)
(Fig. 3c in the main text). So the change of $E_{LP}$ and $E_{MP}$ induced by
$\Omega_{2}$ can be safely ignored for simplicity.