# Analysis of $B_{s}\to\phi\nu\bar{\nu}$ at CEPC

Lingfeng Li  Physics Department, Brown University, Providence, RI 02912, USA
Jockey Club Institute for Advanced Study, The Hong Kong University of Science
and Technology, HKSAR, China  Manqi Ruan<EMAIL_ADDRESS>Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
University of Chinese Academy of Sciences, Beijing 100049, China  Yudong Wang
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China
University of Chinese Academy of Sciences, Beijing 100049, China  Yuexin Wang
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China
University of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

The rare $b\to s\nu\bar{\nu}$ decays are sensitive to contributions of new
physics (NP) and helpful to resolve the puzzle of multiple $B$ flavor
anomalies. In this work, we propose to study the $b\to s\nu\bar{\nu}$
transition at a future lepton collider operating at the $Z$ pole through the
$B_{s}\to\phi\nu\bar{\nu}$ decay. Using the $B_{s}\to\phi$ decay form factors
from lattice simulations, we first update the SM prediction of
BR($B_{s}\to\phi\nu\bar{\nu})_{\mathrm{SM}}=(9.93\pm 0.72)\times 10^{-6}$ and
the corresponding $\phi$ longitudinal polarization fraction
$F_{L,{\mathrm{SM}}}=0.53\pm 0.04$. Our analysis uses the full CEPC simulation
samples with a net statistic of $\mathcal{O}(10^{9})$ $Z$ decays. Precise
$\phi$ and $B_{s}$ reconstructions are used to suppress backgrounds. The
results show that BR($B_{s}\to\phi\nu\bar{\nu})$ can be measured with a
statistical uncertainty of $\mathcal{O}(\%)$ and an $S/B$ ratio of
$\mathcal{O}(1)$ at the CEPC. The quality measures for the event
reconstruction are also derived. By combining the measurement of
BR($B_{s}\to\phi\nu\bar{\nu})$ and $F_{L}$, the constraints on the effective
theory couplings at low energy are given.

##  1 Introduction

The rare flavor-changing-neutral-current (FCNC) $b\to s\nu\bar{\nu}$ decays
are widely recognized as important flavor probes. They are suppressed by the
loop factor and the masses of the heavy weak bosons, as shown in Fig. 1. The
inclusive BR($b\to s\nu\bar{\nu}$) is predicted to be $(2.9\pm 0.3)\times
10^{-5}$ according to Standard Model (SM) calculations [1]. The processes of
this mode are one of the most promising probes to test the SM. Even small
contributions from new physics (NP) could significantly alter their branching
fractions. They also offer the possibility to extract the Cabibbo-Kobayashi-
Maskawa (CKM) matrix elements and search for the origin of the $CP$ and $T$
violations. In the absence of non-factorizable corrections and photon-mediated
contributions [2], the theoretical predictions will be much cleaner than $b\to
s\ell\ell$ transitions. Moreover, the differential $b\to s\nu\bar{\nu}$ decay
width becomes smooth without large QCD loop and hadronic resonance
corrections. TABLE 1 summarizes the current experimental constraints and the
corresponding theoretical predictions for various exclusive $b\to
s\nu\bar{\nu}$ decays.

Figure 1:  The penguin and box diagrams of $b\to s\nu\bar{\nu}$ transition at the leading order. | Current Limit | Detector | SM Prediction
---|---|---|---
BR$(B^{0}\to K^{0}\nu\bar{\nu})$ |  $<2.6\times 10^{-5}$ [3] | BELLE |  $(3.69\pm 0.44)\times 10^{-6}$ [1]
BR$(B^{0}\to K^{\ast 0}\nu\bar{\nu})$ |  $<1.8\times 10^{-5}$ [3] | BELLE |  $(9.19\pm 0.99)\times 10^{-6}$ [1]
BR$(B^{\pm}\to K^{\pm}\nu\bar{\nu})$ |  $<1.6\times 10^{-5}$ [4] | BABAR |  $(3.98\pm 0.47)\times 10^{-6}$ [1]
BR$(B^{\pm}\to K^{\ast\pm}\nu\bar{\nu})$ |  $<4.0\times 10^{-5}$ [5] | BELLE |  $(9.83\pm 1.06)\times 10^{-6}$ [1]
BR$(B_{s}\to\phi\nu\bar{\nu})$ |  $<5.4\times 10^{-3}$ [6] | DELPHI | $(9.93\pm 0.72)\times 10^{-6}$
Table 1: Constraints and predictions for various $b\to s\nu\bar{\nu}$ decays.
The updated BR$(B_{s}\to\phi\nu\bar{\nu})$ comes from our calculation, details
in Sec. 2.

Several anomalies are known to be found in other FCNC decays, e.g.,
$R_{K^{(\ast)}}$ anomalies in FCNC $b\to{s}\ell\ell$ transitions [7, 8, 9].
Anomalies also occur in semileptonic $b\to c\tau(\ell)\nu$ decays with flavor-
changed-charged-current (FCCC), such as $R_{D^{(*)}}$ or $R_{J/\psi}$ [10,
11]. See also [12] for an updated calculation of the FCNC $B\to
K^{*}\nu\bar{\nu}$ decay rate by employing the soft-collinear effective theory
(SCET) sum rule predictions of the heavy-to-light $B$-meson decay form
factors. It is natural to look at the relationships between $b\to
c\tau(\ell)\nu$ and $b\to s\ell\ell$ transitions via gauge invariance to check
these anomalies and solve the puzzle. NP Models can be constrained or
investigated by $b\to s\nu\bar{\nu}$, including the supersymmetry [13, 14,
15], leptoquark models [16, 17, 18, 19, 20, 21], compositeness [22, 23, 24,
25], and gauge extensions [26, 27, 18, 28, 29, 30, 31, 32]. Measuring $b\to
s\nu\bar{\nu}$ transitions in multiple exclusive decay channels is therefore
crucial for investigating possible NP models.

The Circular Electron-Positron Collider (CEPC) is a double-ring $e^{+}e^{-}$
collider with a circumference of 100 km and two interaction points (IP),
enabling precise measurements of SM physics and searches for NP effects. It
operates at the $Z$ pole ($\sqrt{s}=$91.2 GeV), at the $W^{+}W^{-}$ threshold
($\sqrt{s}=$161 GeV), and in Higgs factory mode ($\sqrt{s}=$240 GeV) for
electroweak and flavor physics with a nominal integrated luminosity of 16,
2.6, and 5.6 $\mathrm{ab}^{-1}$, respectively. During the $Z$ pole run [33],
about $0.7\times 10^{12}$ on-shell $Z$-bosons will be produced, which could
further increase in the future. This paper focuses on CEPC as a Tera-$Z$
factory ($10^{12}~{}e^{+}e^{-}\to Z$ events). Given the advantages of the high
luminosity and clean collision environment, we expect a significant
improvement in the precision of rare FCNC decays.

Hadrons | Belle II | LHCb (300 fb-1) | CEPC ($10^{12}Z$)
---|---|---|---
$B^{0}$, $\bar{B}^{0}$ | $5.4\times 10^{10}$ | $\sim 3\times 10^{13}$ | $1.2\times 10^{11}$
$B^{\pm}$ | $5.7\times 10^{10}$ | $\sim 3\times 10^{13}$ | $1.2\times 10^{11}$
$B_{s}$, $\bar{B}_{s}$ | $6.0\times 10^{8}$ | $\sim 1\times 10^{13}$ | $3.1\times 10^{10}$
$B_{c}^{\pm}$ | - | $\sim 2\times 10^{11}$ | $1.8\times 10^{8}$
$\Lambda_{b}$, $\bar{\Lambda}_{b}$ | - | $\sim 2\times 10^{13}$ | $2.5\times 10^{10}$
Table 2: The number of $b$-hadrons expected to be produced in Belle II, LHCb,
and CEPC. Here, the Belle II column corresponds to its 50 ab-1 $\Upsilon(4S)$
run and its 5 ab-1 $\Upsilon(5S)$ run. For more details, see [34].

It turns out that the $Z$ factory mode of CEPC is a great new option for
studying flavor physics because of its relatively high production rates and
high efficiency in reconstructing heavy flavor hadrons. First, flavor studies
at the $Z$ pole run benefit from the large $b$ statistics. The abundant energy
at the $Z$ pole allows $b$ quarks to hadronize into different hadrons. As
TABLE 2 shows, the productions of $B^{0}/\bar{B^{0}}$ and $B^{\pm}$ are
comparable to those at Belle II, while $B_{s}$/$\bar{B}_{s}$ is almost two
orders of magnitude more. For even heavier hadrons such as $B_{c}$ and
$\Lambda_{b}$, the advantage of the $Z$ factories is even more pronounced. As
an $e^{+}e^{-}$ collider, CEPC also benefits from negligible pileup, good
geometric coverage of the detector, and a fixed center-of-mass energy that
allows good precision of the missing momentum. The advanced calorimetry [35,
36, 37] and state-of-the-art tracking system [38] proposed for future
detectors further improve the performance in measuring the missing energy.
Given these advantages, accurate measurement of the missing energy of
neutrinos is very likely. The situation is quite different for hadron collider
detectors such as LHCb, where the missing momentum of a given event cannot be
determined directly. In addition, compared to $B$ factories such as Belle II,
the higher $b$ hadron boost from $Z$ decay makes the tracking more accurate.
Therefore, the measurements in terms of energy/momentum [39] and
direction/displacement [33, 40] are more precise and allow better
discrimination of signal and background events.

We focus on the exclusive process $B_{s}(\bar{B_{s}})\to\phi\nu\bar{\nu}$. The
current upper limit of the branching ratio of this channel is about $5.4\times
10^{-3}$, set by the DELPHI detector at LEP [6]. The threshold is much weaker
than other $b\to s\nu\bar{\nu}$ channels listed in TABLE 1. Most $b\to
s\nu\bar{\nu}$ processes are measured by $B$ factories, where $B_{s}$
production is limited. At the $Z$ pole run, extensive statistics of $B_{s}$
and the precise $\phi$ reconstruction [41] are simultaneously fulfilled.
Therefore, we expect that the observation of this channel and the precise
measurements will be realized for the first time in $Z$ factories. The current
projection of BR$(B_{s}\to\phi\nu\bar{\nu})$ at CEPC comes from the luminosity
re-projection of the LEP study [33]. However, the background suppression
$\varepsilon$ at the LEP search is only $\mathcal{O}(10^{-3})$ [6]. For CEPC,
the same strategy leads to a background size of $\gtrsim 10^{7}$, which makes
the analysis vulnerable to background uncertainties. Therefore, we need to
develop a new analysis framework to reduce the SM backgrounds by more than
$\mathcal{O}(10^{-6})$ to provide a healthy signal-to-background ($S/B$) ratio
near $\mathcal{O}(1)$. In such a case, the measurement of the rare
$B_{s}\to\phi\nu\bar{\nu}$ achieves relative precision at the percentage level
and is robust to systematic uncertainties. We have set up another benchmark
for flavor physics at the $Z$ pole with previous phenomenological studies [42,
43, 44, 34, 45, 46, 47, 48]. It is also true that CEPC detector design shares
many commonalities with other proposals for future $Z$ factories, such as the
Tera-$Z$ mode of FCC-$ee$ [49] and the Giga-$Z$ mode of ILC [50]. Therefore,
the methodology and results of this work will also serve as references for
these projects.

This paper is divided into five sections. Section 2 introduces the physical
background and interpretation of the effective theory of
$B_{s}\to\phi\nu\bar{\nu}$ decay. Section 3 describes the detector model,
software framework, and the simulated samples used in this study. Section 4
presents the analysis of $B_{s}\to\phi\nu\bar{\nu}$ at CEPC. Conclusions are
summarized in Section 5.

##  2 Physics of $B_{s}\to\phi\nu\bar{\nu}$

As discussed in the introduction, many NP scenarios could lead to deviations
of $B_{s}\to\phi\nu\bar{\nu}$ from the SM. This section focuses on the model-
independent approach, which describes the contributions of SM and NP as Wilson
coefficients of the low-energy effective theory (LEFT). If there are no BSM
particles lighter than $m_{B_{s}}$, the low-energy effective Hamiltonian fo
$b\to s\nu\bar{\nu}$ could be written as [1, 51]

$\begin{split}\mathcal{H}_{\mathrm{eff}}=-\frac{4G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}(C_{L}\mathcal{O}_{L}+C_{R}\mathcal{O}_{R})+\text{h.c.}\,,\end{split}$
(1)
$\begin{split}\mathcal{O}_{L(R)}=\frac{e^{2}}{8\pi^{2}}(\bar{s}\gamma^{\mu}P_{L(R)}b)(\bar{\nu_{\ell}}\gamma_{\mu}P_{L}\nu_{\ell})\,.\end{split}$
(2)

Only left-handed quarks interact with $W$ bosons and $C_{R}^{\mathrm{SM}}=0$
at leading order in the SM. The SM prediction of
$C_{L}^{\mathrm{SM}}\simeq-6.47$ come primarily from the top-loop diagrams and
preserves QCD and EW corrections [52]. Since the three neutrino flavors are
indistinguishable at the CEPC detector, each contributes equally to the SM
prediction, giving a total of six Wilson coefficients. Here we assume for
simplicity that the lepton flavor violating (LFV) couplings are negligible.
Following the formalism in [1], we denote the dependence of
BR$(B_{s}\to\phi\nu\bar{\nu}$) on the Wilson coefficients as:

$\begin{split}\frac{\mathrm{BR}(B_{s}\to\phi\nu\bar{\nu})}{\mathrm{BR}(B_{s}\to\phi\nu\bar{\nu})_{\mathrm{SM}}}=\frac{1}{3}\sum_{\ell}(1+\kappa_{\eta}\eta_{\ell})\epsilon_{\ell}^{2}~{},(\ell=e,\mu,\tau)\,,\end{split}$
(3)

where $\kappa_{\eta}$ is the coefficient determined by the ratio between
different (axial)vector $B_{s}\to\phi$ form factors [1], and the two real
quantities are

$\begin{split}\epsilon_{\ell}\equiv\frac{\sqrt{|C_{L}^{\ell}|^{2}+|C_{R}^{\ell}|^{2}}}{|C_{L}^{\mathrm{SM}}|}\,,\quad\eta_{\ell}\equiv-\frac{\mathrm{Re}(C_{L}^{\ell}{C_{R}^{\ell}}^{*})}{|C_{L}^{\ell}|^{2}+|C_{R}^{\ell}|^{2}}\,.\end{split}$
(4)

By measuring BR$(B_{s}\to\phi\nu\bar{\nu}$) at the $Z$ pole, we can constrain
the NP-effects in $C_{L}^{\ell,{\mathrm{NP}}}\equiv
C_{L}^{\ell}-C_{L}^{\ell,{\mathrm{SM}}}$ and
$C_{R}^{\ell,{\mathrm{NP}}}=C_{R}^{\ell}$. However, the coefficient
$\kappa_{\eta}$ is not given in the literature. As a theoretical update, both
$\mathrm{BR}(B_{s}\to\phi\nu\bar{\nu})_{\mathrm{SM}}$ and $\kappa_{\eta}$ are
calculated using $B_{s}\to\phi$ form factors from lattice QCD [53, 54],
including their uncertainties and correlations. Finally, we have

$\mathrm{BR}(B_{s}\to\phi\nu\bar{\nu})_{\mathrm{SM}}=(9.93\pm 0.72)\times
10^{-6}~{},$ (5) $\kappa_{\eta}=1.56\pm 0.08\,.$ (6)

The differential decay width $d\Gamma/dq^{2}$ is also calculated using central
values of the form factor, where the quantity
$q^{2}\equiv(p_{B_{s}}-p_{\phi})^{2}=m_{\nu{\bar{\nu}}}^{2}$ is the invariant
mass squared of the neutrino pair. In our prediction, the hadronic
uncertainties dominate both values. Constant factors such as
$C_{L}^{\mathrm{SM}}$ and $|V_{tb}V_{ts}^{\ast}|$ also contribute slightly to
the decay rate uncertainty.

Besides the decay rate, there is also additional information from
$B_{s}\to\phi\nu\bar{\nu}$ decays, such as the longitudinal polarization
fraction of the $\phi$ meson ($F_{L}$). According [1], the $F_{L}$ dependence
of the LEFT Wilson coefficients is as follows.

$F_{L}=F_{L,{\mathrm{SM}}}\frac{\sum_{\ell}(1+2\eta_{\ell})\epsilon_{\ell}^{2}}{\sum_{\ell}(1+\kappa_{\eta}\eta_{\ell})\epsilon_{\ell}^{2}}\,.$
(7)

Using the same form factors and the method used above, we get
$F_{L,{\mathrm{SM}}}=0.53\pm 0.04$. In phenomenology, $F_{L}$ determines the
kinematic distribution of $\phi\to K^{+}K^{-}$ decays as [55]:

$\frac{\,\mathrm{d}\Gamma}{\,\mathrm{d}\cos\theta}=\frac{3}{4}(1-F_{L})\sin^{2}\theta+\frac{3}{2}F_{L}\cos^{2}\theta\,,$
(8)

where $\theta\in[0,\pi)$ is the angle between $K^{+}$ and $B_{s}$ flight
directions in the $\phi$ rest frame. The different dependences of $F_{L}$ and
BR$(B_{s}\to\phi\nu\bar{\nu})$ on $C_{L(R)}^{\ell}$ further constrain NP
effects.

##  3 The CEPC detector and data samples

As shown in TABLE 1, the value of the signal branching ratio, i.e.,
BR$(B_{s}\to\phi\nu\bar{\nu})=9.93\times 10^{-6}$. Considering the $b$-hadron
fragmentation fractions measured in $Z$ decays, $f(b\to B_{s})=0.101$ [56], a
signal of about $3.0\times 10^{5}$ is produced in CEPC. We focus on the
exclusive mode $B_{s}(\bar{B}_{s})\to\phi(\to K^{+}K^{-})\nu\bar{\nu}$, which
accounts for 49.2% of all signal events. Thus $\mathcal{O}(10^{5})$ signal
events are generated by combining Pythia 8 [57] and EvtGen [58] with the
general decay phase space model. All signal events are reweighted according to
the differential decay width ($d\Gamma/dq^{2}$) calculated in Section 2 to
obtain the correct $q^{2}$ distribution.

Only $Z\to q\bar{q}$ ($q=u,\,d,\,c,\,s,\,b$) events are considered, since
leptonic $Z$ decays make a negligible contribution. Moreover, the SM
background is dominated by heavy quarks ($b$ and $c$). All background samples
in this work are from $\mathcal{O}(10^{9})$ inclusive $Z\to q\bar{q}$ events
generated by WHIZARD [59, 60] and Pythia 6 [57]. Because full simulation of
the detector effects is computationally expensive, it is unrealistic to apply
it to all background samples. Instead, only two subsets of the above samples
are run through the full detector simulation to allocate finite resources. In
the first case, the original $\mathcal{O}(10^{9})$ inclusive $Z\to q\bar{q}$
events are refined to truth-level by three cuts: 1) Heavy quarks must be
produced. 2) At least one neutrino must be produced. 3) At least one $\phi\to
K^{+}K^{-}$ decay must occur. The sample size reduces to $\mathcal{O}(10^{7})$
after the above refinement, making the full detector simulation affordable. To
validate the refined samples above, we also apply the full detector simulation
to $\mathcal{O}(10^{7})$ randomly selected inclusive $Z\to q\bar{q}$ events
without any cuts. In practice, the unrefined backgrounds are used in the early
stages of the analysis, where light quarks and random $K^{+}K^{-}$
combinations are still relevant. In later steps (corresponding to those after
the $b$-tag cut in TABLE 3), we turn to refined backgrounds to achieve better
sampling statistics and stability. The background loss from truth-level
refinement is less than $11\%$ for $3\,\sigma$ kaon PID when matching the
yields two methods. This effect is offset by multiplying this factor to the
background yields.

Detector performance for the full simulation follows the CEPC baseline design
[33]. MokkaPlus [61], a GEANT4 [62]-based simulation framework is used. The
track reconstruction is based on Clupatra [63], and the particle flow
reconstruction is based on the Arbor [64, 65] algorithm. Marlin [66] and LCIO
[67] from ilcsoft are used for data management and formatting.

Figure 2:  The topology of FCNC $B_{s}\to\phi\nu\bar{\nu}$ decay at the $Z$
pole.

Realistic particle identifications (PID) are also included. The most important
effect is the large number of charged pions faking charged kaons. Even a low
rate of $K/\pi$ misidentifications can yield many fake $\phi$. Other sources
of fake kaons, such as protons or muons, are neglected because they are much
rarer than pions in our samples. Estimated from Monte Carlo (MC) sampling, the
typical multiplicities for $K^{\pm}$, $\pi^{\pm}$, and $p$ in the event are
about 2.1, 17.2, and 0.9, respectively. Their momentum distributions above
$\sim 15$ GeV range are highly suppressed. The kaon PID is crucial for flavor
physics because it could improve the reconstruction accuracy of hadrons.
According to CEPC CDR [33], the $K/\pi$ separation power [68, 69] can achieve
$3\,\sigma$ or higher if $dE/dx,~{}dN/dx$ and time of flight information are
included. For more details on PID techniques, see also [70]. So a universal
$K/\pi$ separation power $\gtrsim 3\,\sigma$ at CEPC is a reasonable and
conservative assumption. As will be explained in the later section, to ensure
a stable and high accuracy for the reconstruction of hadrons decaying to
kaons, a 3-$\sigma$ $K/\pi$ separation would be necessary. Therefore, we take
the $3\,\sigma$ $K/\pi$ separation power as the benchmark value for the rest
of this paper. However, since an authentic $K/\pi$ PID algorithm is still
under development, the $K/\pi$ separation is simulated using the Gaussian
approximation. Reconstructions of $\phi$ with alternative $K/\pi$ separation
powers are also analyzed. In addition to fake $\phi$, backgrounds from
semileptonic $b$-hadron decays contribute significantly, see discussions in
section 4.2. We adopt the lepton PID algorithm and performance in [71] to
better represent the lepton information.

##  4 Analysis methods

Fig. 2 shows the typical topology of the target process, i.e., the charged
kaon pair produced by the $\phi$ decay and the neutrino-induced missing
energy. The signal identification consists of three steps. First, we
reconstruct $\phi\to K^{+}K^{-}$ decay vertexes. Second, we use various
features such as the $\phi$ kinematics, missing momentum, lepton energy, and
$b$-tagging to separate the signal from backgrounds. Finally, the Boosted
Decision Tree Gradient (BDTG) method is applied to classify the remaining
events and optimize the background reduction.

###  4.1 $\phi$ Reconstruction

As the only visible component in the $B_{s}\to\phi\nu\bar{\nu}$ signal, $\phi$
plays a central role in our analysis. It has a narrow width
($\Gamma_{\phi}\simeq 4.25$ MeV) and a low inclusive production rate $\sim
5\%$ in $Z\to q\bar{q}$ events. The reconstruction chain of the $\phi$
candidate follows the steps listed below:

  * 1)

We reconstruct all charged kaon tracks. With a finite $K/\pi$ separation
power, the reconstructed kaon tracks also contain misidentified pions.

  * 2)

Match all pairs of oppositely charged kaon tracks and use the kinematic
fitting package [72] to reconstruct their vertex.

  * 3)

Choose pairs of kaons with invariant mass $|m_{K^{+}K^{-}}-m_{\phi}|<$8.5 MeV.

  * 4)

The value of the vertex $\chi^{2}$ is calculated by taking the $\chi^{2}$
contribution from each relevant track using the Minuit algorithm [73]:

$\chi^{2}=\sum_{i=1}^{2}\left(\frac{|V_{i}-V_{\mathrm{fit}}|}{\sigma_{i}}\right)^{2}\,,$
(9)

where $V_{\mathrm{fit}}$ is the fitted vertex position, $V_{i}$ is the point
on one track that is closest to the other, and $\sigma_{i}$ is the uncertainty
of the $i$-th track. Only kaon pairs with $\chi^{2}<8$ are selected.

For more details on the algorithm and performance, see [41]. The reconstructed
$\phi$ mass distribution is shown in Fig. 3.

Figure 3:  Fitted candidate $\phi$ mass distributions in $Z\to q\bar{q}$
background samples, scaled to $10^{12}$ $e^{+}e^{-}\to Z$ events. Here
$m_{\phi}=m_{K^{+}K^{-}}$ distribution is fitted by the Crystal Ball function.
The $K/\pi$ separation power used here is $3\,\sigma$. The two vertical dashed
lines represent the optimized window of invariant mass. Figure 4:
Reconstruction performances for the inclusive charged kaon and $\phi$ with
varying $K/\pi$ separation power in $Z\to q\bar{q}$ samples. To avoid
degradation of the performance at low $K/\pi$ separations, the
$m_{K^{+}K^{-}}$ mass window for $\phi$ reconstruction is also optimized
accordingly.

$\begin{split}\varepsilon&=\frac{\mathrm{Number~{}of~{}correctly~{}reconstructed}\mathrm{~{}candidate}~{}\phi}{\mathrm{Number~{}of~{}}\phi\to
K^{+}K^{-}\mathrm{~{}decays}}\,,\\\
p&=\frac{\mathrm{Number~{}of~{}correctly~{}reconstructed}\mathrm{~{}candidate}~{}\phi}{\mathrm{Number~{}of}\mathrm{~{}candidate}~{}\phi}\,.\end{split}$
(10)

The efficiency and purity of candidate $\phi$ are defined in Eq. (10). Similar
definitions apply to reconstructed kaon tracks. The overall efficiency and
purity for candidate $\phi$ are 48% and 76%, respectively. To better
understand the significance of PID, we also plot inclusive kaon and $\phi$
reconstruction performance with varying $K/\pi$ separation power in Fig. 4. We
parameterize the $K$ and $\pi$ PID performance by two Gaussian distributions
with average values $\mu_{K(\pi)}$ and corresponding standard deviations
$\sigma_{K(\pi)}$. The separation power is defined as
$2|\mu_{\pi}-\mu_{K}|/(\sigma_{\pi}+\sigma_{K})$. Without loss of generality,
we set $\sigma_{\pi}=\sigma_{K}$. Compared to the near-perfect PID case with a
$K/\pi$ separation power$>5\,\sigma$, the $\varepsilon\times p$ for the
$3\,\sigma$ benchmark decrease by $\sim 30\%$ for kaon and $\sim 36\%$ for
$\phi$.

###  4.2 Events selection and results

From the kinematics shown in Fig. 2, it is clear that the $\phi$ decay vertex
of the signal shall be in the hemisphere with the lower visible energy
(“signal hemisphere”) and have a distance from the primary vertex (PV)
comparable to the $b$ lifetime. On the other hand, the number of reconstructed
$\phi$ in each event may be zero, one, or even more. It is necessary to
identify these characteristic $\phi$ before applying more sophisticated
selection rules, since those cuts may depend on the choice of $\phi$. We first
divide the space into two hemispheres by the plane perpendicular to the thrust
axis $\hat{n}_{T}$ [74]. Then we define the “signal $\phi$” according to the
following requirements: 1) Its momentum direction must be in the less
energetic hemisphere. 2) The impact parameters of both kaon tracks are larger
than $0.05~{}\mathrm{mm}$. 3) The distance of the $\phi$ vertex to the primary
vertex (PV) should be greater than $0.4~{}\mathrm{mm}$. 4) The $\phi$ has the
highest energy when multiple $\phi$ satisfy the conditions above. The
hemisphere with(without) the signal $\phi$ is then called the signal(tag)
hemisphere for convenience.

Fig. 5 shows the energy distributions of the signal $\phi$ satisfying the
above conditions. Note that both QCD radiation and heavy quark decays
contribute at this stage. In the first case, the $\phi$ is produced at the PV,
typical in light quark events. Therefore, soft $\phi$ with higher impact
parameter uncertainty has a greater chance of getting through all the cuts.
The $\phi$ from the decay of heavy quarks instead carries significant energy
of the parent particle, which is dominant in signal events. The $b\bar{b}$ and
$c\bar{c}$ backgrounds receive both contributions, leading to double-peaked
structures.

Figure 5: The fitted energy distributions of the leading candidate $\phi$ in
the signal hemisphere for different processes. The samples used correspond to
the third row (“Signal $\phi$”) of the TABLE 3 and are scaled to $10^{12}$ $Z$
decays. Figure 6: The nominal $B_{s}$ energy distributions. The samples used
here satisfy all conditions above this cut in TABLE 1. Figure 7: The nominal
$\alpha$ distributions. The samples used here satisfy all conditions above
this cut in TABLE 1. Figure 8: The leading lepton energy distributions in the
signal hemisphere. The samples used here satisfy all conditions above this cut
in TABLE 1.

After selecting a signal $\phi$ for each event, we can further suppress the SM
background using various event features. At this stage, the main SM
backgrounds are semileptonic heavy quark decays with the $\phi$ produced by
$D$ meson decays. Therefore, we choose several variables and corresponding
cuts summarized:

Figure 9: Event distributions in the $\alpha_{1}-\alpha_{2}$ plane. The signal
(left) and $q\bar{q}$ background (right) samples are the same as in Fig. 5 and
are scaled to $10^{12}$ $Z$ decays. Signal features such as significant
missing energy and an energetic $\phi$ in the signal hemisphere are correctly
reflected in the low-$\alpha_{2}$ and high-$\alpha_{1}$ region.

  * •

The energy asymmetry, defined as the total visible energy difference between
the tag and signal hemisphere ($E_{\mathrm{tag}}-E_{\mathrm{sig}}$), should be
larger than 8 GeV.

  * •

The nominal energy of $B_{s}$,
$E_{\mathrm{B_{s}}}^{N}\equiv\sqrt{s}-E_{\mathrm{tot}}+E_{\phi}$ must be
larger than 28 GeV. See as the Fig. 6

  * •

Three parameters, $\alpha_{1}$, $\alpha_{2}$, and $\alpha$ are defined as
follows:

$\begin{split}\alpha_{1}\equiv\frac{E_{\phi}}{{E_{\mathrm{sig}}}}\,,\,\alpha_{2}\equiv\frac{E_{\mathrm{sig}}}{{E_{\mathrm{beam}}}}\,,\,\alpha\equiv\frac{\alpha_{2}}{\alpha_{1}}=\frac{(E_{\mathrm{sig}})^{2}}{E_{\mathrm{beam}}E_{\phi}}\,.\end{split}$
(11)

Considering the topology of the signal decay, most of the energy of the signal
hemisphere should come from the $\phi$, i.e., correspond to a large
$\alpha_{1}$. At the same time, the missing energy from the $B_{s}$ meson
should also be significant, leading to a lower $\alpha_{2}$. We keep only the
events with $\alpha<1.1$, see as the Fig. 7. Meanwhile, Fig. 9 shows the
distributions in the $\alpha_{1}-\alpha_{2}$ plane for the signal and
backgrounds.

  * •

The $b$-tagging score of events (ranging from 0 to 1) must be greater than
0.6, using the same $b$-tagging algorithm described in [75].

  * •

The energy of the leading lepton ($e$ or $\mu$) in the signal hemisphere
should be less than 1.2 GeV. The cut suppresses backgrounds considerably, with
the remaining ones containing leptons softer than 1.2 GeV or hadronic $\tau$.
Fig. 8 shows the energy distribution of the corresponding leading lepton.

  * •

The angle between the missing momentum and the $\phi$ momentum
($\theta^{\mathrm{miss}}_{\phi}$) must be larger than 0.1.

We list the cut flow of the above selection rules corresponding to the second
block in TABLE 3. It is noteworthy that the $b$-tagging score $>$0.6
requirement suppresses the light flavor backgrounds by more than two orders of
magnitude. Even under the conservative assumption that the remaining light
flavor events have similar efficiencies to $b\bar{b}$ in the rest of the
analysis, they contribute at most $\mathcal{O}(10^{-2})$ to the total
background and can be safely ignored.

Cuts | $B_{s}\to\phi\nu\bar{\nu}$ | $u\bar{u}+d\bar{d}+s\bar{s}$ | $c\bar{c}$ | $b\bar{b}$ | total bkg |  $\sqrt{S+B}/S$ (%)
---|---|---|---|---|---|---
CEPC events ($10^{12}Z$) | $3.03\times 10^{5}$ | $4.28\times 10^{11}$ | $1.20\times 10^{11}$ | $1.51\times 10^{11}$ | $6.99\times 10^{11}$ | 276
$N_{\phi(\to K^{+}K^{-})}>0$ | $8.09\times 10^{4}$ | $1.09\times 10^{10}$ | $4.04\times 10^{9}$ | $6.08\times 10^{9}$ | $2.10\times 10^{10}$ | 179
111The candidate $\phi$ here satisfy the following conditions: 1) In the signal hemisphere. 2) The impact parameters of both kaon pair tracks are larger than 0.05 mm. 3) The distance between the decay point of $\phi$ and interact point(IP) is larger than 0.4 mm.“Signal” $\phi$ | $5.38\times 10^{4}$ | $2.52\times 10^{8}$ | $4.09\times 10^{8}$ | $1.69\times 10^{9}$ | $2.35\times 10^{9}$ | 90.9
Energy asymmetry $>8$ GeV | $4.74\times 10^{4}$ | $6.25\times 10^{7}$ | $9.76\times 10^{7}$ | $4.93\times 10^{8}$ | $6.53\times 10^{8}$ | 53.9
$E_{B_{s}}^{N}>28$ GeV | $4.06\times 10^{4}$ | $4.25\times 10^{6}$ | $9.59\times 10^{6}$ | $5.00\times 10^{7}$ | $6.38\times 10^{7}$ | 19.7
$\alpha<1.1$ | $3.03\times 10^{4}$ | $2.41\times 10^{6}$ | $3.10\times 10^{6}$ | $8.47\times 10^{6}$ | $1.40\times 10^{7}$ | 12.4
$b$-tag $>0.6$ | $2.33\times 10^{4}$ | $<2.0\times 10^{4}$ | $2.95\times 10^{5}$ | $5.97\times 10^{6}$ | $6.27\times 10^{6}$ | 10.77
$E_{\mu}$ and $E_{e}<1.2$ GeV | $2.10\times 10^{4}$ | - | $5.85\times 10^{4}$ | $2.10\times 10^{6}$ | $2.16\times 10^{6}$ | 7.03
$\theta_{\phi}^{\mathrm{miss}}>0.1~{}\mathrm{rad}$ | $1.77\times 10^{4}$ | - | $2.75\times 10^{4}$ | $1.38\times 10^{6}$ | $1.41\times 10^{6}$ | 6.75
$q^{2}<14.0~{}\mathrm{GeV}^{2}$ | $1.34\times 10^{4}$ | - | $2.02\times 10^{4}$ | $6.04\times 10^{5}$ | $6.24\times 10^{5}$ | 5.96
BDTG response$>0.89$ | $0.75\times 10^{4}$ | - | $<1\times 10^{2}$ | $1.03\times 10^{4}$ | $1.03\times 10^{4}$ | 1.78
Efficiency | 2.40% | - | - | $6.82\times 10^{-8}$ | $1.47\times 10^{-8}$ | -
Table 3:  The cut chain for the signal and $q\bar{q}$ with full simulation
samples and scaled to the integrated luminosity of the $10^{12}~{}Z$ bosons at
CEPC. The cut chain before the cut of leading lepton energy uses the general
inclusive samples with sizes of $\mathcal{O}(10^{7})$. The light-flavor
contributes less than 3‰ to the total background after the $b$-tagging cut and
is neglected in later steps. Starting from the leading lepton energy cut, the
truth-level refined background samples are used. The background loss due to
the switch to refined samples is $\lesssim$11%. This effect is compensated by
multiplying the subsequent background yields by a factor of 1.11. The kaon PID
are simulated with $3\,\sigma$ $K/\pi$ separation power in the table, see Fig.
4 for the performance of the kaon PID.

After the above cuts, the remaining backgrounds are still an order of
magnitude larger than the signal. It is then necessary to perform a thorough
event reconstruction to better separate them from our signal. The primary goal
is to reproduce the correct $B_{s}$ energy and missing mass squared $q^{2}$
for the signal events using a rational algorithm.

Figure 10: The reconstructed $E_{B_{s}}$ (left) and $q^{2}$ (right)
distributions of the signal and backgrounds before the BDTG cut. For
comparison, we also show the signal distributions at the truth-level.

The reconstruction starts with an updated estimate of $E_{B_{s}}$. In the
previously defined nominal
$E^{N}_{B_{s}}\equiv\sqrt{s}-E_{\mathrm{tot}}+E_{\phi}$, we use global energy
conservation to estimate the missing momentum. However, the calculation
involves the energy measurement errors and neutrino(s) impact in the tag
hemisphere. To reduce the noise in the tag hemisphere, we define a better
approximation of the truth-level $E_{B_{s}}$ as

$E_{B_{s}}^{(0)}=\frac{\sqrt{s}}{2}-E_{\mathrm{sig}}+E_{\phi}\,,$ (12)

where $E_{\mathrm{sig}}$ is the total visible energy in the signal hemisphere.
By this definition, the value $E_{B_{s}}^{(0)}$ is less affected by the tag
hemisphere measurements than $E_{B_{s}}^{N}$. We then assign the $p_{B_{s}}$
direction the same as the displacement of the $\phi$ decay vertex from the PV.
Since $B_{s}$ energy and momentum direction are known, we calculate the four-
momentum $p_{B_{s}}^{(0)}$ after setting the $B_{s}$ on-shell condition. The
value of $q^{2}$ is then calculated by definition as
$(p_{B_{s}}^{(0)}-p_{\phi})^{2}$.

However, the estimate of $E_{B_{s}}$ in Eq. (12) can still be improved. Since
$Z$ hadronic decays are not perfectly symmetric, the total energies at truth-
level in the two hemispheres will not be exactly $\sqrt{s}/2$. An energy
imbalance leads to corrections on top of Eq. (12). Therefore, we introduce the
following relations:

$\begin{split}M_{\mathrm{tag}}&=\sqrt{\left(\sum
p_{\mathrm{tag}}^{\mathrm{vis}}\right)^{2}}\,,\\\
M_{\mathrm{sig}}^{(i)}&=\sqrt{\left(\sum
p_{\mathrm{sig}}^{\mathrm{vis}}+p_{B_{s}}^{(i-1)}-p_{\phi}\right)^{2}}\,,\\\
E_{B_{s}}^{(i)}&=\frac{s+(M_{\mathrm{sig}}^{(i-1)})^{2}-M^{2}_{\mathrm{tag}}}{2\sqrt{s}}-E_{\mathrm{sig}}+E_{\phi}\,,\\\
(q^{2})^{(i)}&=(p_{B_{s}}^{(i-1)}-p_{\phi})^{2}\,,\end{split}$ (13)

where $p^{\mathrm{vis}}_{\mathrm{sig(tag)}}$ are the four-momenta of the
visible particles in the signal (tag) hemisphere. The third equation above
encodes the imbalance of $Z$ decay products in the two hemispheres. Starting
with the initial value $E_{B_{s}}^{(0)}$ in Eq. (12), we solve Eq. (13)
iteratively to obtain a self-consistent signal reconstruction. It turns out
that Eq. (13) converges quickly, leaving little room for improvement after the
first iteration. Therefore, we choose the values of the first iteration
($M_{\mathrm{sig}}^{(1)}$, $E_{B_{s}}^{(1)}$, and $(q^{2})^{(1)}$) as our
event reconstruction results and BDTG inputs.

In Fig. 10 we show the reconstructed $E_{B_{s}}$ and $q^{2}$ distributions for
samples passing all cuts in Section 4.2 to compare the truth-level
distribution. The apparent differences between the signal and the backgrounds
can serve as the input for later analysis. The typical $q^{2}$ and $E_{B_{s}}$
reconstruction errors of signal events, which are defined as the difference
between reconstruction and truth-level values, are $2.5$ GeV2 and 1.7 GeV,
respectively. The complicated and asymmetric response of the detector causes
the overall $E_{B_{s}}$ and $q^{2}$ distributions to deviate slightly from the
truth, which could be partially recovered with a better understanding of the
detector system. For comparison, the error between the nominal $E^{N}_{B_{s}}$
and the truth-level $E_{B_{s}}$ is $5.1~{}\mathrm{GeV}$, three times worse
than the algorithm output. The nominal $q^{2}$ derived from $E_{B_{s}}^{N}$ is
even further from the truth and therefore useless. The accuracy of the
reconstructed $E_{B_{s}}$ and $q^{2}$ thus provides us a simple way to
evaluate the overall CEPC detector performance. In particular, the neutral
hadron/photon momenta suffer larger uncertainties than track momenta. They
contribute significantly to the errors of $E_{B_{s}}$ and $q^{2}$. The
displacement of $\phi$ decay vertex is another source of error since the
reconstruction algorithm relies on the direction of $B_{s}$. Finally, to
further suppress the background, a cut of $q^{2}<14.0~{}\mathrm{GeV}^{2}$ is
imposed based on the above results.

Figure 11: BDTG output distributions for signal and background events, ranging
from -1 to 1. The samples used here passed all the cuts introduced above and
are scaled to $10^{12}~{}Z$ decays.

Besides, the complex relationship among multiple observables is not captured
by simple cuts. As a final step in the analysis, we use the BDTG method of the
TMVA package [76] to train binary event classifiers to optimize measurement
accuracy. The training considers multiple inputs, which are summarized below:

  * •

General event-shape variables: energy asymmetry and $E^{N}_{B_{s}}$.

  * •

The largest impact parameter of all tracks.

  * •

Parameters $\alpha_{1}$ and $\alpha_{2}$ in Eq. (11).

  * •

The angle $\theta^{\mathrm{miss}}_{\phi}$.

  * •

The invariant mass of all visible particles, as well as the visible particle
invariant masses in the tag/signal hemisphere.

  * •

Reconstructed $E_{B_{s}}$ and $q^{2}$.

  * •

The leading electron and muon energies in the signal hemisphere.

  * •

The largest track impact parameter in the signal hemisphere, excluding kaons
from any reconstructed $\phi$.

  * •

Kaon tracks’ impact parameters from the signal $\phi$.

  * •

The signal $\phi$ invariant mass.

Fig. 11 shows the BDTG responses to the test samples, with the signal and
background distributions peaking at $-1.0$ and $1.0$, respectively. With the
optimized cut of the BDTG response at 0.75, we reject over $98\%$ of
$b\bar{b}$ and $c\bar{c}$ backgrounds at the cost of a $44\%$ signal loss. As
summarized in TABLE 3, the $S/B$ ratio reaches 77% after the BDTG cut. The
$1\,\sigma$ Tera-$Z$ sensitivity of the signal strength is estimated by
$\frac{\sqrt{B+S}}{S}$, which corresponds to about $1.78$%. We also evaluate
the sensitivity and $S/B$ ratio with a perfect kaon PID to motivate better
future PID performance. Without any fake kaon tracks and a comparable $S/B\geq
70\%$, the sensitivity of BR($B_{s}\to\phi\nu\bar{\nu}$) is $1.52\%$. The
sensitivity of the branching ratio as a function of the kaon PID is shown in
Fig. 12, which shows stable performance in a wide range of $K/\pi$ separation
power. Besides, taking the benchmark $3\,\sigma$ $K/\pi$ separation power,
Fig. 13 shows the projected sensitivity as a function of
BR($B_{s}\to\phi\nu\bar{\nu}$). Multiple signal features included in the
analysis allow for high sensitivities even in the no kaon PID case.

Figure 12: The sensitivity of BR($B_{s}\to\phi\nu\bar{\nu}$) as a function of
kaon PID, parameterized by the $K/\pi$ separation power.
Figure 13: Projected experimental sensitivity at CEPC ($10^{12}~{}Z$ decays)
as a function of BR($B_{s}\to\phi\nu\bar{\nu}$), shown as the red curve. The
current upper limit from LEP for BR($B_{s}\to\phi\nu\bar{\nu}$) is indicated
by green dashed line. The prediction of SM corresponds to the blue line used
in TABLE 3.

###  4.3 Constraints on Wilson coefficients

The event reconstruction is also effective when measuring the $\phi$
longitudinal polarization fraction $F_{L}$. Fig. 14 shows the distribution of
$\cos\theta$, where $\theta$ is the angle between $B_{s}$ and $K^{+}$(or
$K^{-}$) in the $\phi$ rest frame. Here the truth-level distribution of signal
events is reweighed according to the SM prediction $F_{L,{\mathrm{SM}}}\simeq
0.53$. However, the background statistics after the BDTG cut is insufficient
for a good background fit. Instead, we use the background $\cos\theta$
distribution before the BDTG cut and scale the yields according to the
Tera-$Z$ luminosity. The $p_{B_{s}}$ reconstruction error dominates the
$\sigma_{\theta}$ between reconstructed and the truth values, which is about
$0.047$. Such a $\theta$ reconstruction error corresponds to a difference
$\approx 0.04$ between our $F_{L}$ fit and the truth-value. The estimated
statistical uncertainty of $F_{L}$ is $0.008$ at CEPC, which is subdominant.
Since it is not our goal to thoroughly estimate the differential measurement
of $F_{L}$ in this work, discussions of other systematics are reserved for
future work.

In Fig. 15, we plot the 68% C.L. constraints on the NP contributions to the
LEFT couplings $C_{L}^{\mathrm{NP}}$ and $C_{R}^{\mathrm{NP}}$ at CEPC. We
assume that the Wilson coefficients are a real number and satisfy the LFU,
i.e. $C^{e,\mu,\tau}_{L(R)}=C_{L(R)}$ for simplicity. The
BR$(B_{s}\to\phi\nu\bar{\nu})$ measurement with a relative accuracy of
$1.78\%$ places tight constraints in the
$C_{L}^{\mathrm{NP}}-C_{R}^{\mathrm{NP}}$ plane. We show the regions that
$F_{L}=F_{L,{\mathrm{SM}}}\pm 0.04$ as a suggestive value for the $F_{L}$
measurement. As can be seen in Fig. 15, $|C_{L,R}^{\mathrm{NP}}|$ are limited
to $\lesssim 0.2~{}C_{L}^{\mathrm{SM}}$ after combining the branching ratio
and differential $F_{L}$ measurements. All theoretical uncertainties are
ignored in Fig. 15 to directly illustrate the CEPC flavor physics potential.

##  5 Conclusion

In this paper, we study the phenomenology of the rare FCNC decay
$B_{s}\to\phi\nu\bar{\nu}$ at the $Z$ pole with the full simulated samples of
the CEPC detector profile. The large $B_{s}$ statistic ($\sim 3\times 10^{10}$
from $10^{12}$ $Z$ decays) enables precise measurement of such a rare decay.

We calculated the SM prediction that
BR($B_{s}\to\phi\nu\bar{\nu})_{\mathrm{SM}}=(9.93\pm 0.72)\times 10^{-6}$ with
the lattice $B_{s}\to\phi$ form factors. The hadronic form factors are also
the major contributors to the theoretical uncertainty. The results also
predict the $\phi$ longitudinal polarization fraction $F_{L,{\mathrm{SM}}}$ to
be $0.53\pm 0.04$. In this analysis, $\phi\to K^{+}K^{-}$ vertexes are
primarily reconstructed, with their integrated efficiency and purity reaching
about 48% and 76% under a realistic kaon PID assumption. After a series of
cuts and optimization of the BDTG method, the dominant backgrounds $Z\to
q\bar{q}$ are suppressed by a factor of $\mathcal{O}(10^{-8})$. The remaining
backgrounds are mainly $Z\to b\bar{b}$ events. The final signal efficiency is
$3\%$, resulting in a relative sensitivity of BR($B_{s}\to\phi\nu\bar{\nu}$)
as low as $1.78\%$. The high $S/B$ ratio $\sim 77\%$ makes the measurement
robust against potential systematic uncertainties.

The integrated and differential measurements of this channel are sensitive to
the six dim-6 LEFT operators. The constraints will further contribute to the
global determination of the NP effects behind the $B$ anomalies and allow
discrimination between NP models. We also estimated the $F_{L}$ uncertainty
using the angular distribution of the signal events.

We expect other $b\to s\nu\bar{\nu}$ measurements at CEPC, e.g., $B\to$
pseudoscalar transition $B^{\pm}\to K^{\pm}\nu\bar{\nu}$, can further improve
the NP limit. By studying the specific mode $B_{s}\to\phi\nu\bar{\nu}$, there
is an opportunity to resolve multiple anomalies in the measurements of
$B$-meson decays. The result also allows us to test BSM models and update our
knowledge of QCD.

Figure 14: Fitted $\cos\theta$ distributions for the signal and backgrounds in
total, scaled to $10^{12}$ $Z$ decays. The background samples used here are
before the BDTG cut to ensure sufficient statistics. The error bars shown are
determined by the sample size before scaling.
Figure 15:  Projected $1\,\sigma$ constraints on $C_{L}^{\mathrm{NP}}$ and
$C_{R}^{\mathrm{NP}}$ from $B_{s}\to\phi\nu\bar{\nu}$ measurements at CEPC.
The narrow green band comes from measuring BR($B_{s}\to\phi\nu\bar{\nu}$),
while the orange band corresponds to the suggestive $F_{L}$ uncertainty. Both
values are normalized according to $C_{L}^{\mathrm{SM}}\simeq-6.47$.
Theoretical uncertainties are ignored to better demonstrate experimental
precision.

##  6 acknowledgments

We thank Tao Liu, Wei Wang, Dan Yu and Taifan Zheng for useful discussions.
This study was supported by the International Partnership Program of Chinese
Academy of Sciences (Grant No.113111KYSB20190030). It was also supported by
Innovative Scientific Program of Institute of High Energy Physics, entitled
”New Physics Search at CEPC”. LL was supported by the General Research Fund
(GRF) under Grant No. 16312716, which was issued by the Research Grants
Council of Hong Kong S.A.R..

## References

  * [1] A. J. Buras, J. Girrbach-Noe, C. Niehoff, and D. M. Straub, “$B\to{K}^{\left(\ast\right)}\nu\overline{\nu}$ decays in the Standard Model and beyond,” JHEP, vol. 02, p. 184, 2015.
  * [2] T. Blake, G. Lanfranchi, and D. M. Straub, “Rare $B$ Decays as Tests of the Standard Model,” Prog. Part. Nucl. Phys., vol. 92, pp. 50–91, 2017.
  * [3] J. Grygier et al., “Search for $\boldsymbol{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle,” Phys. Rev. D, vol. 96, no. 9, p. 091101, 2017.  [Addendum: Phys.Rev.D 97, 099902 (2018)].
  * [4] J. P. Lees et al., “Search for $B\to K^{(*)}\nu\overline{\nu}$ and invisible quarkonium decays,” Phys. Rev. D, vol. 87, no. 11, p. 112005, 2013.
  * [5] O. Lutz et al., “Search for $B\to h^{(*)}\nu\bar{\nu}$ with the full Belle $\Upsilon(4S)$ data sample,” Phys. Rev. D, vol. 87, no. 11, p. 111103, 2013.
  * [6] D. Collaboration, “Study of rareb decays with the DELPHI detector at LEP,” Zeitschrift für Physik C: Particles and Fields, vol. 72, no. 2, p. 207, 1996.  Experiment Lastest result.
  * [7] R. Aaij et al., “Test of lepton universality with $B^{0}\rightarrow K^{*0}\ell^{+}\ell^{-}$ decays,” JHEP, vol. 08, p. 055, 2017.
  * [8] R. Aaij et al., “Search for lepton-universality violation in $B^{+}\to K^{+}\ell^{+}\ell^{-}$ decays,” Phys. Rev. Lett., vol. 122, no. 19, p. 191801, 2019.
  * [9] A. Abdesselam et al., “Test of Lepton-Flavor Universality in ${B\to K^{\ast}\ell^{+}\ell^{-}}$ Decays at Belle,” Phys. Rev. Lett., vol. 126, no. 16, p. 161801, 2021.
  * [10] R. Aaij et al., “Measurement of the ratio of branching fractions $\mathcal{B}(B_{c}^{+}\,\to\,J/\psi\tau^{+}\nu_{\tau})$/$\mathcal{B}(B_{c}^{+}\,\to\,J/\psi\mu^{+}\nu_{\mu})$,” Phys. Rev. Lett., vol. 120, no. 12, p. 121801, 2018.
  * [11] A. Abdesselam et al., “Measurement of $\mathcal{R}(D)$ and $\mathcal{R}(D^{\ast})$ with a semileptonic tagging method,” 4 2019.
  * [12] J. Gao, C.-D. Lü, Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, “Precision calculations of $B\to V$ form factors from soft-collinear effective theory sum rules on the light-cone,” Phys. Rev. D, vol. 101, no. 7, p. 074035, 2020.
  * [13] F. Domingo, “Update of the flavour-physics constraints in the NMSSM,” Eur. Phys. J. C, vol. 76, no. 8, p. 452, 2016.
  * [14] D.-Y. Wang, Y.-D. Yang, and X.-B. Yuan, “$b\to c\tau\bar{\nu}$ decays in supersymmetry with $R$-parity violation,” Chin. Phys. C, vol. 43, no. 8, p. 083103, 2019.
  * [15] Q.-Y. Hu, Y.-D. Yang, and M.-D. Zheng, “Revisiting the $B$-physics anomalies in $R$-parity violating MSSM,” Eur. Phys. J. C, vol. 80, no. 5, p. 365, 2020.
  * [16] R. Barbieri, C. W. Murphy, and F. Senia, “B-decay Anomalies in a Composite Leptoquark Model,” Eur. Phys. J., vol. C77, no. 1, p. 8, 2017.
  * [17] R. Barbieri and A. Tesi, “$B$-decay anomalies in Pati-Salam SU(4),” Eur. Phys. J., vol. C78, no. 3, p. 193, 2018.
  * [18] J. Kumar, D. London, and R. Watanabe, “Combined Explanations of the $b\to s\mu^{+}\mu^{-}$ and $b\to c\tau^{-}{\bar{\nu}}$ Anomalies: a General Model Analysis,” Phys. Rev., vol. D99, no. 1, p. 015007, 2019.
  * [19] L. Calibbi, A. Crivellin, and T. Li, “Model of vector leptoquarks in view of the $B$-physics anomalies,” Phys. Rev. D, vol. 98, no. 11, p. 115002, 2018\.
  * [20] M. Blanke and A. Crivellin, “$B$ Meson Anomalies in a Pati-Salam Model within the Randall-Sundrum Background,” Phys. Rev. Lett., vol. 121, no. 1, p. 011801, 2018.
  * [21] A. Crivellin, C. Greub, D. Müller, and F. Saturnino, “Importance of Loop Effects in Explaining the Accumulated Evidence for New Physics in B Decays with a Vector Leptoquark,” Phys. Rev. Lett., vol. 122, no. 1, p. 011805, 2019.
  * [22] D. M. Straub, “Anatomy of flavour-changing $Z$ couplings in models with partial compositeness,” JHEP, vol. 08, p. 108, 2013.
  * [23] C. Niehoff, P. Stangl, and D. M. Straub, “Violation of lepton flavour universality in composite Higgs models,” Phys. Lett. B, vol. 747, pp. 182–186, 2015.
  * [24] F. Sannino, P. Stangl, D. M. Straub, and A. E. Thomsen, “Flavor Physics and Flavor Anomalies in Minimal Fundamental Partial Compositeness,” Phys. Rev. D, vol. 97, no. 11, p. 115046, 2018.
  * [25] P. P. Stangl, Direct Constraints, Flavor Physics, and Flavor Anomalies in Composite Higgs Models.  PhD thesis, Munich, Tech. U., 2018.
  * [26] S. M. Boucenna, A. Celis, J. Fuentes-Martin, A. Vicente, and J. Virto, “Phenomenology of an $SU(2)\times SU(2)\times U(1)$ model with lepton-flavour non-universality,” JHEP, vol. 12, p. 059, 2016.
  * [27] C.-W. Chiang, X.-G. He, J. Tandean, and X.-B. Yuan, “$R_{K^{(*)}}$ and related $b\to s\ell\bar{\ell}$ anomalies in minimal flavor violation framework with $Z^{\prime}$ boson,” Phys. Rev., vol. D96, no. 11, p. 115022, 2017.
  * [28] P. Asadi, M. R. Buckley, and D. Shih, “It’s all right(-handed neutrinos): a new W′ model for the ${R}_{D^{{\left(\ast\right)}}}$ anomaly,” JHEP, vol. 09, p. 010, 2018.
  * [29] A. Greljo, D. J. Robinson, B. Shakya, and J. Zupan, “R(D(∗)) from W′ and right-handed neutrinos,” JHEP, vol. 09, p. 169, 2018\.
  * [30] M. Abdullah, J. Calle, B. Dutta, A. Flórez, and D. Restrepo, “Probing a simplified, $W^{\prime}$ model of $R(D^{(\ast)})$ anomalies using $b$-tags, $\tau$ leptons and missing energy,” Phys. Rev. D, vol. 98, no. 5, p. 055016, 2018.
  * [31] A. Greljo, J. Martin Camalich, and J. D. Ruiz-Álvarez, “Mono-$\tau$ Signatures at the LHC Constrain Explanations of $B$-decay Anomalies,” Phys. Rev. Lett., vol. 122, no. 13, p. 131803, 2019.
  * [32] J. D. Gómez, N. Quintero, and E. Rojas, “Charged current $b\to c\tau\bar{\nu}_{\tau}$ anomalies in a general $W^{\prime}$ boson scenario,” Phys. Rev. D, vol. 100, no. 9, p. 093003, 2019.
  * [33] M. Dong et al., “CEPC Conceptual Design Report: Volume 2 - Physics & Detector,” 2018.
  * [34] L. Li and T. Liu, “$b\to s\tau^{+}\tau^{-}$ physics at future Z factories,” JHEP, vol. 06, p. 064, 2021.
  * [35] M. Dong, “R&D of the CEPC scintillator-tungsten ECAL,” JINST, vol. 13, no. 03, p. C03024, 2018.
  * [36] H. Zhao, C. Fu, D. Yu, Z. Wang, T. Hu, and M. Ruan, “Particle flow oriented electromagnetic calorimeter optimization for the circular electron positron collider,” JINST, vol. 13, no. 03, p. P03010, 2018.
  * [37] J. Jiang, S. Zhao, Y. Niu, Y. Shi, Y. Liu, D. Han, T. Hu, and B. Yu, “Study of SiPM for CEPC-AHCAL,” Nucl. Instrum. Meth. A, vol. 980, p. 164481, 2020\.
  * [38] G. Tang et al., “The circular electron–positron collider beam energy measurement with Compton scattering and beam tracking method,” Rev. Sci. Instrum., vol. 91, no. 3, p. 033109, 2020.
  * [39] N. Berger, M. Kiehn, A. Kozlinskiy, and A. Schöning, “A New Three-Dimensional Track Fit with Multiple Scattering,” Nucl. Instrum. Meth. A, vol. 844, p. 135, 2017.
  * [40] A. Abada et al., “FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2,” Eur. Phys. J. ST, vol. 228, no. 2, pp. 261–623, 2019.
  * [41] T. Zheng, J. Wang, Y. Shen, Y.-K. E. Cheung, and M. Ruan, “Reconstructing $K^{0}_{S}$ and $\Lambda$ in the CEPC baseline detector,” Eur. Phys. J. Plus, vol. 135, no. 3, p. 274, 2020.
  * [42] J. F. Kamenik, S. Monteil, A. Semkiv, and L. V. Silva, “Lepton polarization asymmetries in rare semi-tauonic $b\rightarrow s$ exclusive decays at FCC-$ee$,” Eur. Phys. J., vol. C77, no. 10, p. 701, 2017.
  * [43] M. Dam, “Tau-lepton Physics at the FCC-ee circular e+e- Collider,” SciPost Phys. Proc., vol. 1, p. 041, 2019.
  * [44] T. Zheng, J. Xu, L. Cao, D. Yu, W. Wang, S. Prell, Y.-K. E. Cheung, and M. Ruan, “Analysis of $b_{c}$ $\rightarrow$ $\tau$$\nu_{\tau}$ at cepc,” vol. 45, p. 023001, jan 2021.
  * [45] Y. Amhis, M. Hartmann, C. Helsens, D. Hill, and O. Sumensari, “Prospects for $B_{c}^{+}\to\tau^{+}\nu_{\tau}$ at FCC-ee,” 5 2021.
  * [46] M. Chrzaszcz, R. G. Suarez, and S. Monteil, “Hunt for rare processes and long-lived particles at FCC-ee,” Eur. Phys. J. Plus, vol. 136, no. 10, p. 1056, 2021.
  * [47] R. Aleksan, L. Oliver, and E. Perez, “CP violation and determination of the $bs$ ”flat” unitarity triangle at FCCee,” 7 2021.
  * [48] R. Aleksan, L. Oliver, and E. Perez, “Study of CP violation in $B^{\pm}$ decays to $\overline{D^{0}}(D^{0})K^{\pm}$ at FCCee,” 7 2021.
  * [49] A. Abada et al., “FCC Physics Opportunities,” Eur. Phys. J., vol. C79, no. 6, p. 474, 2019.
  * [50] K. Fujii et al., “Tests of the Standard Model at the International Linear Collider,” 8 2019.
  * [51] W. Altmannshofer et al., “The Belle II Physics Book,” 2018.
  * [52] J. Brod, M. Gorbahn, and E. Stamou, “Two-Loop Electroweak Corrections for the $K\to\pi\nu\bar{\nu}$ Decays,” Phys. Rev. D, vol. 83, p. 034030, 2011\.
  * [53] R. R. Horgan, Z. Liu, S. Meinel, and M. Wingate, “Lattice QCD calculation of form factors describing the rare decays $B\to K^{*}\ell^{+}\ell^{-}$ and $B_{s}\to\phi\ell^{+}\ell^{-}$,” Phys. Rev. D, vol. 89, no. 9, p. 094501, 2014\.
  * [54] R. R. Horgan, Z. Liu, S. Meinel, and M. Wingate, “Rare $B$ decays using lattice QCD form factors,” PoS, vol. LATTICE2014, p. 372, 2015.
  * [55] W. Altmannshofer, A. J. Buras, D. M. Straub, and M. Wick, “New strategies for New Physics search in $B\to K^{*}\nu\bar{\nu}$, $B\to K\nu\bar{\nu}$ and $B\to X_{s}\nu\bar{\nu}$ decays,” JHEP, vol. 04, p. 022, 2009.
  * [56] Y. Amhis et al., “Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2016,” Eur. Phys. J. C, vol. 77, no. 12, p. 895, 2017.
  * [57] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, “An introduction to PYTHIA 8.2,” Comput. Phys. Commun., vol. 191, pp. 159–177, 2015.
  * [58] D. J. Lange, “The EvtGen particle decay simulation package,” Nucl. Instrum. Meth., vol. A462, pp. 152–155, 2001.
  * [59] W. Kilian, T. Ohl, and J. Reuter, “WHIZARD: Simulating Multi-Particle Processes at LHC and ILC,” Eur. Phys. J. C, vol. 71, p. 1742, 2011.
  * [60] M. Moretti, T. Ohl, and J. Reuter, “O’Mega: An Optimizing matrix element generator,” pp. 1981–2009, 2 2001.
  * [61] P. Mora de Freitas and H. Videau, “Detector simulation with MOKKA / GEANT4: Present and future,” in International Workshop on Linear Colliders (LCWS 2002), pp. 623–627, 8 2002.
  * [62] S. Agostinelli et al., “GEANT4–a simulation toolkit,” Nucl. Instrum. Meth. A, vol. 506, pp. 250–303, 2003.
  * [63] F. Gaede, S. Aplin, R. Glattauer, C. Rosemann, and G. Voutsinas, “Track reconstruction at the ILC: the ILD tracking software,” J. Phys. Conf. Ser., vol. 513, p. 022011, 2014.
  * [64] M. Ruan and H. Videau, “Arbor, a new approach of the Particle Flow Algorithm,” in International Conference on Calorimetry for the High Energy Frontier, pp. 316–324, 2013.
  * [65] M. Ruan et al., “Reconstruction of physics objects at the Circular Electron Positron Collider with Arbor,” Eur. Phys. J. C, vol. 78, no. 5, p. 426, 2018.
  * [66] F. Gaede, “Marlin and LCCD: Software tools for the ILC,” Nucl. Instrum. Meth. A, vol. 559, pp. 177–180, 2006.
  * [67] F. Gaede, T. Behnke, N. Graf, and T. Johnson, “LCIO: A Persistency framework for linear collider simulation studies,” eConf, vol. C0303241, p. TUKT001, 2003.
  * [68] C. Lippmann, “Particle identification,” Nucl. Instrum. Meth. A, vol. 666, pp. 148–172, 2012.
  * [69] F. An, S. Prell, C. Chen, J. Cochran, X. Lou, and M. Ruan, “Monte Carlo study of particle identification at the CEPC using TPC dE / dx information,” The European Physical Journal C, vol. 78, no. 6, p. 464, 2018.
  * [70] G. Wilkinson, “Particle identification at FCC-ee,” Eur. Phys. J. Plus, vol. 136, no. 8, p. 835, 2021.
  * [71] D. Yu, M. Ruan, V. Boudry, and H. Videau, “Lepton identification at particle flow oriented detector for the future $e^{+}e^{-}$ Higgs factories,” Eur. Phys. J. C, vol. 77, no. 9, p. 591, 2017.
  * [72] T. Suehara and T. Tanabe, “LCFIPlus: A Framework for Jet Analysis in Linear Collider Studies,” Nucl. Instrum. Meth. A, vol. 808, pp. 109–116, 2016\.
  * [73] F. James and M. Roos, “Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations,” Comput. Phys. Commun., vol. 10, pp. 343–367, 1975.
  * [74] D. Barber et al., “Discovery of Three Jet Events and a Test of Quantum Chromodynamics at PETRA Energies,” Phys. Rev. Lett., vol. 43, p. 830, 1979\.
  * [75] Y. Bai, C. Chen, Y. Fang, G. Li, M. Ruan, J.-Y. Shi, B. Wang, P.-Y. Kong, B.-Y. Lan, and Z.-F. Liu, “Measurements of decay branching fractions of $H\to b\bar{b}/c\bar{c}/gg$ in associated $(e^{+}e^{-}/\mu^{+}\mu^{-})H$ production at the CEPC,” Chin. Phys. C, vol. 44, no. 1, p. 013001, 2020.
  * [76] A. Hocker et al., “TMVA - Toolkit for Multivariate Data Analysis,” 3 2007\.


# High-Throughput Ab Initio Design of Atomic Interfaces using InterMatch

Eli Gerber  School of Applied and Engineering Physics, Cornell University,
Ithaca, NY 14853, USA  Steven B. Torrisi  Department of Physics, Harvard
University, Cambridge, MA 02138, USA  Energy & Materials Division, Toyota
Research Institute, Los Altos, CA 94022, USA  Sara Shabani  Department of
Physics, Columbia University, New York, NY, USA  Eric Seewald  Department of
Physics, Columbia University, New York, NY, USA  Jordan Pack  Department of
Physics, Columbia University, New York, NY, USA  Jennifer E. Hoffman
Department of Physics, Harvard University, Cambridge, MA 02138, USA  John A.
Paulson School of Engineering and Applied Sciences,Harvard University,
Cambridge, MA 02138, USA  Cory R. Dean  Department of Physics, Columbia
University, New York, NY, USA  Abhay N. Pasupathy  Department of Physics,
Columbia University, New York, NY, USA  Eun-Ah Kim  Department of Physics,
Cornell University, Ithaca, NY 14853, USA

( [file: main])

###### Abstract

Forming a hetero-interface is a materials-design strategy that can access an
astronomically large phase space. However, the immense phase space
necessitates a high-throughput approach for an optimal interface design. Here
we introduce a high-throughput computational framework, InterMatch, for
efficiently predicting charge transfer, strain, and superlattice structure of
an interface by leveraging the databases of individual bulk materials.
Specifically, the algorithm reads in the lattice vectors, density of states,
and the stiffness tensors for each material in their isolated form from the
Materials Project. From these bulk properties, InterMatch estimates the
interfacial properties. We benchmark InterMatch predictions for the charge
transfer against experimental measurements and supercell density-functional
theory calculations. We then use InterMatch to predict promising interface
candidates for doping transition metal dichalcogenide MoSe2. Finally, we
explain experimental observation of factor of 10 variation in the supercell
periodicity within a few microns in graphene/$\alpha$-RuCl3 by exploring low
energy superlattice structures as a function of twist angle using InterMatch.
We anticipate our open-source InterMatch algorithm accelerating and guiding
ever-growing interfacial design efforts. Moreover, the interface database
resulting from the InterMatch searches presented in this paper can be readily
accessed through MPContribs.

Figure 1: (a) Role of InterMatch in the materials discovery process. (b) Input
from the bulk database. Lattice vectors $\vec{a}_{\alpha}$,
$\vec{b}_{\alpha}$, density-of-states $g_{\alpha}$ of systems $\alpha=1,2$.
Elastic tensors are additional inputs. (c) $E_{F}^{\alpha}$ are bulk Fermi
levels and $E_{F}^{\prime}$ is new equilibrium Fermi level. (d) $\Delta n$ is
the transferred charge density and $d$ is the interlayer separation between
the two systems, taken to be the sum of the largest van der Waals radii of the
species in each system 1 and 2. (e) Superlattice vectors $\vec{v_{i}}$ (orange
arrows) and their near-equivalent vectors $\vec{u_{i}}$ (magenta arrows).
Candidate supercells, formed in each basis by combining
$\\{(\vec{v}_{i},\vec{u}_{i}),(\vec{v}_{j},\vec{u}_{j})\\}$ pairs specify the
strain $\varepsilon_{i}$ (blue arrows). (f) Optimal supercell minimizes the
elastic energy and the number of atoms in the cell.

With increasing control in interface fabrication, interfacial systems form an
arena of limitless possibilitiesGeim and Grigorieva (2013). Recent
developments with moiré heterostructuresKennes et al. (2021) further enlarged
the phase space to include the twist angle. However, the vast space of
possibilities also implies it is crucial to go beyond serendipitous
discoveries and empirical explorations to effectively harness the intrinsic
potential of interfacial systems. The traditional approach to theoretically
studying interfaces is to carry out density-functional theory (DFT)
calculations on a supercell system consisting of two materialsKomsa and
Krasheninnikov (2013); Terrones and Terrones (2014); Ebnonnasir et al. (2014);
Bokdam et al. (2014); Zhou et al. (2015). While such approaches are rigorous,
computational limitations regularly require imposing unnatural strain to form
a periodic structure. Moreover, the $\mathcal{O}(N^{3})$ scaling of DFT in the
number of electrons $N$ for each such calculation prohibits a comprehensive
exploration. Some of us recently proposed an “intermediate scale” approach
called Mismatched INterface Theory (MINT)Gerber et al. (2020), which can
predict charge transfer and natural strain approximating one layer of the
interface using finite-size scaling of atomic clusters. While MINT
calculations are much more computationally affordable, they are not fast
enough for an exhaustive survey in real-time. Therefore, a comprehensive and
fast approach to scanning the relevant phase space of interfacial combinations
is greatly needed.

With advancements in widely available comprehensive materials databasesHorton
et al. (2021); Levin (2020); Landis et al. (2012); Curtarolo et al. (2012);
Jain et al. (2013); Saal et al. (2013); de Jong et al. (2015); Ashton et al.
(2017); Haastrup et al. (2018); Zhou et al. (2019); Draxl and Scheffler
(2019); Talirz et al. (2020); Choudhary et al. (2020a), it is timely to
establish a high-throughput approach to interface design that can leverage the
information contained in these databases to make predictions of interface
physics. While it is possible to look up a pair of materials in the database
and try to match their properties, the phase space of all possible
combinations amounts easily to $>\mathcal{O}(10^{6})$ possibilities, which is
beyond the scale of what a manual search can reasonably accomplish. At the
same time, it is desirable to benefit from extensive existing databases and
search exhaustively for ideal combinations to reach a given objective. Indeed,
efforts to make interfacial predictions using bulk material databases are
beginning to emergeMathew et al. (2016); Ding et al. (2016); Choudhary et al.
(2020b); Boland and Singh (2022). However, so far, the existing approaches aid
growth and calculation decisions for a specific pair of materials rather than
allowing for a comprehensive query to yield fast, approximate predictions over
a wide range of possible interfaces.

In this paper, we introduce InterMatch, which uses information readily
available from preexisting materials databases such as the Materials Project
and 2DMatPedia databases to predict charge transferMilnes and Feucht (1972),
strainLandau et al. (1986), stabilityMuller (1994), and optimal
superlatticeKaxiras (2003) of an atomic interface. Using these predictions,
InterMatch can narrow the candidate pool from $C>\mathcal{O}(10^{6})$ to
$\mathcal{O}(10)$ that can then be investigated in greater detail using MINT
or supercell DFT (see Fig. 1(a)). We first illustrate how two branches of
InterMatch predict the charge transfer and optimal superlattice after querying
the entries of the Materials Project for each of the constituents of the
interface. We then benchmark InterMatch predictions for the charge transfer
against experimental measurements and supercell DFT predictions. We then
employ InterMatch’s branches to address two bottleneck problems obstructing
design of interfaces towards the goal of discovering new physics: the problem
of doping transition metal dichalcogenides and the problem of predicting
stable interface structure, applied to the graphene/$\alpha$-RuCl3 system. We
comment on many other classes of interfaces that can be optimized using
InterMatch.

Starting with ab initio materials data, InterMatch performs high-throughput
screening of possible heterostructures via pairwise calculation of desired
interface properties including charge transfer $\Delta n$, strain tensor
$\tilde{\varepsilon}$, optimized superlattice vectors $\vec{v_{1}},$
$\vec{v_{2}}$, and number of atoms $N$. Once InterMatch identifies a promising
pool of candidate combinations, one can make a more in-depth analysis of the
smaller pool using MINT or supercell DFT (See Fig 1(a)). The InterMatch
algorithm has two branches to predict two key electronic and mechanical
characteristics of candidate interfaces: charge transfer and optimized
superlattice structure (See Fig. 1). One branch is devoted to calculating
charge transfer (Fig. 1(c)-(d)), and the other is devoted to optimizing
supercell structure by minimizing the number of atoms and elastic energy (Fig.
1(e)-(f)). For the first branch that estimates the direction and magnitude of
charge transfer, we use a simple model to describe the Fermi level shifts
occurring in each material when they are brought together in proximityRuan and
Ching (1987). Fig. 1(c) shows how the Fermi level shifts are determined:
systems 1 and 2 are designated as “donor” or “acceptor” based on their
relative Fermi levels $E_{F}^{1}$ and $E_{F}^{2}$, and the difference of
integrals over $g_{1}(E)$ and $g_{2}(E)$

$\int_{E_{F}^{{}^{\prime}}}^{E_{F}^{1}}dE\;g_{1}(E)=\int_{E_{F}^{2}}^{E_{F}^{{}^{\prime}}}dE\;g_{2}(E)$
(1)

is minimized to determine the equilibrium Fermi level $E_{F}^{{}^{\prime}}$.
We take interaction between the two systems at the interface into account in
the estimation of the charge transfer $\Delta n$ using a simple capacitor
modelLowell and Rose-Innes (1980). Specifically, we model the interface as a
parallel plate capacitor with the separation $d$ given by the sum of the
largest van der Waals radii of the species in each system 1 and 2 (See Fig.
1(d)). The charge transfer depends on the equilibrium Fermi level
$E_{F}^{\prime}$ and the distance $d$ as $e\Delta
n=\varepsilon_{0}E_{F}^{{}^{\prime}}/d$.

Figure 2: Benchmarking InterMatch. (a) Comparison of charge transfer predicted
by InterMatch with measured experimental values for interfaces in RefsAnnadi
et al. (2013); Wang et al. (2020); Sutter et al. (2009); Zheng et al. (2021a).
(b) InterMatch screening of over 10,000 2D materials in heterostructure with
monolayer MoSe2, ranked in descending order of charge transfer
$\absolutevalue{\Delta n}$. (c) Substrate selection based on InterMatch
screening results from red box in (a) according to $\Delta n$, elastic energy
$\mathcal{E}$, and energy above-hull. (d) Comparison of InterMatch predictions
for $\Delta n$ (solid symbols) with supercell DFT calculations (open symbols).

The second branch of the InterMatch algorithm sketched in Fig. 1(e-f)
constructs optimal supercells from a pair of queried systems by calculating
strain and elastic energy at their interface over a series of supercell
configurations. Given the lattice vectors of system 1, $\vec{a}_{1}$,
$\vec{b}_{1}$, and those of system 2, $\vec{a}_{2}$, $\vec{b}_{2}$, the
algorithm searches for pairs of near-equivalent superlattice vectors
$\\{(\vec{u}_{1},\vec{v}_{1}),(\vec{u}_{2},\vec{v}_{2})\\}$:

$\displaystyle\vec{u}_{i}$
$\displaystyle=\mathbf{M}^{i}_{11}\vec{a}_{i}+\mathbf{M}^{i}_{12}\vec{b}_{i}$
(2) $\displaystyle\vec{v}_{i}$
$\displaystyle=\mathbf{M}^{i}_{21}\vec{a}_{i}+\mathbf{M}^{i}_{22}\vec{b}_{i}$

where $\mathbf{M}^{i}$ is a $2\times 2$ matrix of integer coefficients for the
system $i$ and the “near-equivalence” is defined by

$\mathbf{M}^{1}=(\tilde{\varepsilon}^{2}+1)\mathcal{R}_{\theta}\mathbf{M}^{2}.$
(3)

Here $\mathcal{R}_{\theta}$ is an in-plane rotation matrix by angle $\theta$
and $\tilde{\varepsilon}^{2}$ is the strain tensor resulting from straining
the Bravais lattice of system 2 to match that of system 1.

We choose system 2 to be the material with the smallest elements of the
stiffness tensor $\mathbf{C}$ (queried from the Materials Project) in the
strain direction. InterMatch then computes the elastic energy
$\mathcal{E}=\frac{1}{2}C_{ijkl}\varepsilon_{ij}\varepsilon_{kl}$ for the
superlattice candidate according to classical elastic plate theoryLandau et
al. (1986). The optimal supercell is determined by simultaneously minimizing
the elastic energy $\mathcal{E}$ and the cell area
$\absolutevalue{\vec{u}_{i}\crossproduct\vec{v}_{i}}$. The use of elastic
energy goes beyond previous approaches for finding the superlatticeStradi et
al. (2017); Lazić (2015) which only consider geometric strain.

Figure 3: Superlattice structure prediction optimizing elastic energies.
(a)-(b) (Top view) Diagonal stiffness tensor components $C_{11}$ and $C_{22}$
of primitive MoSe2 and ZrTe3 unit cells in Voigt notation. (c)-(d) (Top view)
ZrTe3 layer of two candidate MoSe2/ZrTe3 supercells with the same number of
atoms and average strain $\varepsilon_{av}^{\textnormal{ZrTe}_{3}}$. The solid
black boxes denote the strained ZrTe3 unit cells and the dashed red boxes are
the original unstrained primitive cells. (e) Average strain values
$\varepsilon_{av}$ of the MoSe2/ZrTe3 interfaces in Cells 1 and 2. (f) Elastic
energies $\mathcal{E}$ of the interfaces.

Now we demonstrate how elastic energy considerations can make a difference in
optimization of the superlattice. Consider the MoSe2/ZrTe3 interface. The
primitive unit cells of MoSe2 and ZrTe3 are shown in Fig. 3 (a) and (b)
respectively, along with the diagonal components of their stiffness tensors
$C_{ii}$ in Voigt notation. The anisotropy of the ZrTe3 stiffness tensor is
such that the energetic cost to deforming ZrTe3 along the direction of
$C_{22}$ far exceeds the cost of an equivalent deformation along the direction
of $C_{11}$. Cells 1 and 2 in Fig. 3 (c) and (d) result from an InterMatch
search for low-area, low-strain MoSe2/ZrTe3 supercells. The two cells are
identical in number of atoms and geometric strain $\varepsilon_{av}$ (shown in
Fig. 3 (e)), however, cell 1 is favored energetically due to the different
strains required to make each supercell commensurate with MoSe2 (Fig. 3 (f)).

We now benchmark charge transfer predictions by InterMatch against
experimentally measured charge transfer in known interfaces. Fig. 2 (a) shows
a comparison of InterMatch predictions of charge transfer with experimentally
obtained values for several interfaces: LaAlO3/SrTiO3(1 1 0)Annadi et al.
(2013), GR/$\alpha$-RuCl3Wang et al. (2020), GR/Pt$(1\,1\,1)$Sutter et al.
(2009), and MoS2/MgAl2O4Zheng et al. (2021a). The magnitudes of the $\Delta n$
predicted with InterMatch are at the same order of magnitude as the measured
values, especially given experimental error bars, with the exception of
GR/$\alpha$-RuCl3. However, spin-orbit coupling effects (absent from our
calculations) are known to affect the band structure of $\alpha$-RuCl3Kim et
al. (2015), altering the band alignment with GR and the resulting charge
transfer.

Next we turn to the application of charge transfer prediction to the problem
of doping transition metal dichalcogenides (TMDs). TMDs have emerged as an
exciting van der Waals material platform at the intersection of semiconductor
physics and strong correlation physicsManzeli et al. (2017). Due to the spin-
valley locking Ising spin-orbit coupling, an exotic $p$-wave superconducting
state was proposed for hole-doped TMDsHsu et al. (2017). Recent developments
in TMD moiré systems have further extended the phase space of possibilities.
However, a major bottleneck against testing these proposals is the difficulty
of establishing a good contact. Empirically, it has been established that
doping the contact area can significantly improve the contact resistanceZheng
et al. (2021b). However, gate-based doping does not scale well. While
successful modulation doping using work function difference was established in
graphene/$\alpha$-RuCl3 heterostructuresWang et al. (2020), it is desirable to
perform an exhaustive search of interface possibilities.

We seek 2D substrates for controlling carrier concentration in MoSe2. We
screen all entries of the 2DMatPedia database and 3000 entries from the
Materials Project for stable 2D materials composed of elements making up the
majority of commercially available semiconductors, semimetals, and metals. We
use InterMatch to down-sample from 10,000 candidate 2D substrates based on the
magnitude of the predicted charge transfer $\absolutevalue{\Delta n}$ to MoSe2
in the desired range $\gtrsim\mathcal{O}(10^{13})$ cm-2 (Fig. 2 (b)). We then
select from these the compounds with maximum $\absolutevalue{\Delta n}$,
minimal strain, and minimal above-hull energy (Fig. 2 (c)). Finally, we choose
a small subset of the top interfaces (2H-TaS2, $\beta$-GaSe, and ZrTe3, in the
case of our example) and benchmark InterMatch predictions against supercell
DFT calculations using the optimized supercells generated by InterMatch. Fig.
2 (d) shows a comparison of InterMatch predictions with the results from DFT
for the top three interfaces (for computational details, see the Supplemental
Material). The magnitude of the $\Delta n$ prediction from the two approaches
are within $10^{13}$ cm-2. Moreover, both approaches find consistent relative
magnitude of charge transfer. Given the high-throughput nature of InterMatch,
these agreements encourage using InterMatch as the first pass in searches for
optimal heterostructures.

As an example of the power of Intermatch to understand superlattice structure,
we consider the graphene/$\alpha$-RuCl3 heterostructure (GR/$\alpha$-RuCl3).
This system has attracted great interest due to the presence of strong
modulation dopingWang et al. (2020) and enhancement of $\alpha$-RuCl3’s
proximity to the Kitaev spin liquid phase. However, relatively little
attention has been paid to the atomic scale structure of the heterostructure
and the possible influence on electronic properties. In order to study this
experimentally, we used scanning tunneling microscopy (STM) to investigate the
properties of GR/$\alpha$-RuCl3 heterostructures created by mechanical
exfoliation and colamination, as shown in Fig. 4 (a). The angle between the
$\alpha$-RuCl3 substrate and graphene was not intentionally controlled. Shown
in Fig. 4 (b)-(d) are a set of STM topographs taken at various locations of
the GR/$\alpha$-RuCl3 heterostructure. The locations are within a few microns
of each other on the sample shown in Fig. 4 (a). Intriguingly, all three of
the regions show moiré patterns with large wavelengths - 2.7 nm in Fig. 4 (b),
11.7 nm in Fig. 4 (c) and 25.7 nm in Fig. 4 (d). All three of these
wavelengths are much larger than the wavelength set by the difference in
lattice constants.

Using InterMatch, we perform a comprehensive mapping of the space of
superlattice configurations spanned by $(\theta,L,\mathcal{E})$ where $\theta$
is the twist angle, $L$ is the moiré period, and $\mathcal{E}$ is the elastic
energy of the interface. The resulting spectrum of low-energy superlattice
configurations for $0^{\circ}\leq\theta\leq 30^{\circ}$ and
$0\;\textnormal{nm}\leq L\leq 30\;\textnormal{nm}$ is shown in Fig. 4 (e). We
identify four prominent moiré length scales (blue boxes in Fig. 4 (e))
occurring within a $5^{\circ}$ range between $15^{\circ}$-$20^{\circ}$. Three
of the four length scales coincide with those observed in STM at
$L=2.7,\;11.7,\;25.7$ nm, shown in Fig. 4 (b)-(d). Correctly identifying
energetically favorable GR/$\alpha$-RuCl3 superlattices over a narrow range of
twist angles showcases InterMatch’s capability to predict interfacial
structure of complex (e.g extremely lattice-mismatched) systems.

The presence of an atomic reconstruction at the interface of GR/$\alpha$-RuCl3
can have dramatic consequences for the spectroscopic properties of the
material. Shown in Fig. 4 (f)-(h) are scanning tunneling spectra averaged over
the regions shown in Fig. 4 (b)-(d). These spectra show dramatic differences
from the simple expectation for a doped Dirac spectrum as might be expected
from charge transfer alone. Instead, we observe strong resonances in all three
regions, with the spacing between resonances following the expectation from
Landau levels on a Dirac spectrum. Previously, such spectra have been observed
when graphene has a periodic bucklingMao et al. (2020), where it was ascribed
to periodic strain in the material. In our case, apart from the strain
associated with the moiré latticeShabani et al. (2021), we expect that there
will also be strong periodic variations in the dopingRizzo et al. (2022) that
contribute to the formation of resonances.

Figure 4: (a) Optical image of the measured device of GR/$\alpha$-RuCl3
contacted with bismuth indium tin for the STM measurements. The blue and green
dashed lines show the boundary of $\alpha$-RuCl3 and graphene, respectively.
(b)-(d) STM topographic images (in pm) of GR/$\alpha$-RuCl3 on 2.7 nm (set
points of -100 mV and -100 pA), 11.7 nm and 25.7nm (set points of -1 V and -50
pA) moiré patterns due to atomic reconstruction. (e) InterMatch predictions
for low-energy GR/$\alpha$-RuCl3 superlattice configurations as a function of
period $L$, twist angle $\theta$, and elastic energy $\mathcal{E}$. Left panel
is a projection onto the $L$-axis, bottom panel is a projection onto the
$\theta$-axis. Dashed box indicates interval of $\theta$ containing largest
range of stable superlattice periods, shaded blue boxes indicate regions of
likely superlattice configurations. Purple stars denote the periodicities
extracted from the experiment. (f)-(h) dI/dV measurements corresponding to the
three moiré patterns in (b)-(d) showing strong resonances dependent on moire
wavelengths.

In summary, we introduce and demonstrate InterMatch, a high-throughput
computational framework and database for predicting charge transfer, strain,
and superlattice of an interface between two arbitrary materials. Charge
transfer allows heterostructure-based modulation dopingWang et al. (2020),
which can guide device fabrication and contact design Zheng et al. (2021b).
Efficiently determining the smallest energetically favorable commensurate
supercells from a wide variety of interface configurations is crucial for
accelerating ab initio studies. We showcase the use of InterMatch by
identifying high-charge transfer substrates for doping TMDs, and by predicting
equilibrium moiré superlattice configurations for the lattice-mismatched
GR/$\alpha$-RuCl3 interface that are validated by STM measurements. The
presence of such long-wavelength superlattice modulations at van der Waals
interfaces present new opportunities to tailor bandstructure using materials
that do not have a close match in lattice constants. The evolving interface
database provides open access to InterMatch results which we hope will help
guide future exploration of interfacial systems.

To broadly benefit the community, we made the InterMatch code openly
accessible at https://doi.org/10.5281/zenodo.6823973Gerber and Torrisi (2022).
Moreover, we tabulate InterMatch results in an open-access “interface
database” directly integrated with the Materials Project via the MPContribs
platform. At the time of writing, the database contains $\sim 200,000$
interfaces (and counting) in simple JavaScript Object Notation (JSON) that are
queryable and sortable according to the chemical composition of either
constituent system, charge transfer, strain, and optimized supercell size. In
addition, we generate crystallographic information files (CIF) of interface
supercells with InterMatch which may be readily accessed from the database and
used as inputs for DFT or other first principles studies.

###### Acknowledgements.

The authors thank Kin Fai Mak, Jie Shan, Stephen Carr, Patrick Huck, Matthew
Horton, Jason Munro, and Vidya Madhavan for helpful discussions. EAK and JH
were supported by MURI grant FA9550-21-1-0429. EG was supported by the Cornell
Center for Materials Research with funding from the NSF MRSEC program
(DMR-1719875). SBT was supported by the Department of Energy Computational
Science Graduate Fellowship under grant DE-FG02-97ER25308. The computation was
done using the high powered computing cluster WALLE2 that was established
through the support of Gordon and Betty Moore Foundation’s EPiQS Initiative,
Grant GBMF10436 to EAK and the New Frontier Grant from Cornell University’s
College of Arts and Sciences and hosted and maintained by Cornell Center for
Advanced Computing. STM experiments were supported by NSF DMR-2004691 (ANP)
and AFOSR via grant FA9550-21-1-0378 (SS, ES). Sample synthesis for STM
measurements was supported by the NSF MRSEC program through Columbia in the
Center for Precision-Assembled Quantum Materials (PAQM), grant number
DMR-2011738.

## Supplemental Material

In Section I we present additional InterMatch results identifying optimal
substrate candidates for forming high-charge-transfer interfaces with group-VI
transition metal dichalcogenides and putative spin liquid materials
$\alpha$-RuCl3 and TbInO3. In Section II we provide details of the InterMatch
code. In Section III we provide computational details of the ab initio
density-functional theory (DFT) calculations of MoSe2 supercells in Fig. 2(d)
of the main text.

##  I Charge transfer substrates for WTe2, WSe2, and $\alpha$-RuCl3

We further demonstrate InterMatch by applying it to two types of materials:
group-VI TMDs, and putative spin liquid $\alpha$-RuCl3Plumb et al. (2014).
TMDs’ unique combination of properties makes them highly attractive for
nanoelectronics applications and fundamental studies of novel physical
phenomenaManzeli et al. (2017). However, the realization of many such high-
performance devices and exotic phases is limited by the availability of
systems with high carrier mobility and low contact resistances between metal
contacts and the semiconductor. Quantum spin liquids (QSLs) are interacting
quantum systems in which spins do not order at low temperatures, and have been
theorized to offer insights into high-temperature superconductivity upon
doping. $\alpha$-RuCl3, for example, has been intensively discussed as a
possible candidate for Kitaev physics; however, it orders
antiferromagnetically at low temperatures due to the presence of additional
magnetic couplings extending beyond the pure Kitaev interaction. Doping
$\alpha$-RuCl3 with charge carriers has been predicted to enhance Kitaev
interactions and push $\alpha$-RuCl3 closer to the spin liquid phase. We use
InterMatch to identify stable, high-charge-transfer interfaces for electron-
and hole-doping the TMDs WTe2, WSe2, and MoSe2, and the putative QSL
$\alpha$-RuCl3. Fig. S1 shows sample InterMatch results of substrate
candidates for interfaces with WTe2, WSe2, and $\alpha$-RuCl3, highlighting
those that minimize elastic energy $\mathcal{E}$ and maximize charge transfer
$\Delta n$. In Fig. S1 (a)-(b), we survey all $\sim 70,000$ oxides from the
Materials Project. Oxides have a wide range of charge neutrality levels and
therefore constitute a powerful addition to electrostatic gating or chemical
doping for controlling carrier concentration in heterostructures. In Fig. S1
(c) we survey all entries in the 2DMatpedia database, and determine two-
dimensional (2D) substrate candidates to maximally dope $\alpha$-RuCl3.

Figure S1: Sample InterMatch results for different charge transfer interfaces.
The solid green line indicates the pareto-optimal frontier of substrate
candidates that minimize elastic energy $\mathcal{E}$ and maximize charge
transfer $\absolutevalue{\Delta n}$. Calculated quantities in (a)-(b) use ab
initio data from the Materials Project as inputs, while those in (c) use data
from the 2DMatpedia database.

##  II InterMatch Code

The InterMatch code is written in Python 3.7 and makes extensive use of
pymatgenOng et al. (2013), an open-source Python package of the Materials
Project, for the manipulation and analysis of various structures of interest.
The code is continuously being developed, and the latest version can be
obtained at https://doi.org/10.5281/zenodo.6823973Gerber and Torrisi (2022).
We aim to provide an efficient scheme for computing interface properties
capable of screening a significant fraction of combinations of existing
Materials Project structure entries, returning the results in real-time
(typical run time for the calculation of a single interface is $10\pm 5$
seconds on a 1.7 GHz Intel Core i5 processor at 1333 MHz using 4 GB of RAM,
running macOS Sierra 10.12.6).

##  III Computational Details

All ab initio DFT calculations were carried out within the total-energy plane
wave density-functional pseudopotential approach, using Perdew-Burke-Ernzerhof
generalized gradient approximation functionalsPerdew et al. (1996) and
optimized norm-conserving Vanderbilt pseudopotentials in the SG15
familySchlipf and Gygi (2015). Plane wave basis sets with energy cutoffs of 30
hartree were used to expand the electronic wave functions. We used fully
periodic boundary conditions and a $8\times 8\times 1$ $k$-point mesh to
sample the Brillouin zone. Electronic minimizations were carried out using the
analytically continued functional approach starting with a LCAO initial guess
within the DFT$++$ formalismFreysoldt et al. (2009), as implemented in the
open-source code JDFTxSundararaman et al. (2017) using direct minimization via
the conjugate gradients algorithmPayne et al. (1992). All unit cells were
constructed to be inversion symmetric about $z=0$ with a distance of $\sim 60$
bohr between periodic images of the MoSe2 surface, using coulomb truncation to
prevent image interaction.

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Kennes et al. (2021) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021), ISSN 1745-2481.
  * Komsa and Krasheninnikov (2013) H.-P. Komsa and A. V. Krasheninnikov, Physical Review B 88, 085318 (2013), publisher: American Physical Society.
  * Terrones and Terrones (2014) H. Terrones and M. Terrones, Journal of Materials Research 29, 373 (2014), ISSN 2044-5326.
  * Ebnonnasir et al. (2014) A. Ebnonnasir, B. Narayanan, S. Kodambaka, and C. V. Ciobanu, Applied Physics Letters 105, 031603 (2014), ISSN 0003-6951, publisher: American Institute of Physics.
  * Bokdam et al. (2014) M. Bokdam, T. Amlaki, G. Brocks, and P. J. Kelly, Phys. Rev. B 89, 201404 (2014).
  * Zhou et al. (2015) S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, Physical Review B 92, 155438 (2015).
  * Gerber et al. (2020) E. Gerber, Y. Yao, T. A. Arias, and E.-A. Kim, Physical Review Letters 124, 106804 (2020).
  * Horton et al. (2021) M. K. Horton, S. Dwaraknath, and K. A. Persson, Nature Computational Science 1, 3 (2021), ISSN 2662-8457.
  * Levin (2020) I. Levin, _NIST Inorganic Crystal Structure Database (ICSD)_ (2020).
  * Landis et al. (2012) D. D. Landis, J. S. Hummelshøj, S. Nestorov, J. Greeley, M. Dułak, T. Bligaard, J. K. Nørskov, and K. W. Jacobsen, Computing in Science Engineering 14, 51 (2012).
  * Curtarolo et al. (2012) S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, L. J. Nelson, G. L. Hart, S. Sanvito, M. Buongiorno-Nardelli, et al., Computational Materials Science 58, 227 (2012), ISSN 0927-0256.
  * Jain et al. (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, et al., APL Materials 1, 011002 (2013).
  * Saal et al. (2013) J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, JOM 65, 1501 (2013).
  * de Jong et al. (2015) M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag, J. J. Plata, et al., Scientific Data 2, 150009 (2015).
  * Ashton et al. (2017) M. Ashton, J. Paul, S. B. Sinnott, and R. G. Hennig, Physical Review Letters 118, 106101 (2017).
  * Haastrup et al. (2018) S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, et al., 2D Materials 5, 042002 (2018), ISSN 2053-1583.
  * Zhou et al. (2019) J. Zhou, L. Shen, M. D. Costa, K. A. Persson, S. P. Ong, P. Huck, Y. Lu, X. Ma, Y. Chen, H. Tang, et al., Scientific Data 6, 86 (2019), ISSN 2052-4463.
  * Draxl and Scheffler (2019) C. Draxl and M. Scheffler, Journal of Physics: Materials 2, 036001 (2019).
  * Talirz et al. (2020) L. Talirz, S. Kumbhar, E. Passaro, A. V. Yakutovich, V. Granata, F. Gargiulo, M. Borelli, M. Uhrin, S. P. Huber, S. Zoupanos, et al., Scientific Data 7, 299 (2020).
  * Choudhary et al. (2020a) K. Choudhary, K. F. Garrity, A. C. E. Reid, B. DeCost, A. J. Biacchi, A. R. Hight Walker, Z. Trautt, J. Hattrick-Simpers, A. G. Kusne, A. Centrone, et al., npj Computational Materials 6, 173 (2020a).
  * Mathew et al. (2016) K. Mathew, A. K. Singh, J. J. Gabriel, K. Choudhary, S. B. Sinnott, A. V. Davydov, F. Tavazza, and R. G. Hennig, Computational Materials Science 122, 183 (2016), ISSN 0927-0256.
  * Ding et al. (2016) H. Ding, S. S. Dwaraknath, L. Garten, P. Ndione, D. Ginley, and K. A. Persson, ACS Applied Materials & Interfaces 8, 13086 (2016), ISSN 1944-8244.
  * Choudhary et al. (2020b) K. Choudhary, K. F. Garrity, G. Pilania, and F. Tavazza (2020b).
  * Boland and Singh (2022) T. M. Boland and A. K. Singh, Computational Materials Science 207, 111238 (2022).
  * Milnes and Feucht (1972) A. G. Milnes and D. L. Feucht, _Heterojunctions and metal-semiconductor junctions_ (Academic Press, New York, 1972), ISBN 0-12-498050-3.
  * Landau et al. (1986) L. D. Landau, E. M. Lifshitž, J. B. Sykes, and W. H. Reid, in _Course of Theoretical Physics: Theory of Elasticity_ (Butterworth-Heinemann, 1986), chap. 2, pp. 46–58.
  * Muller (1994) P. Muller, Pure and Applied Chemistry 66, 1077 (1994).
  * Kaxiras (2003) E. Kaxiras, _Atomic and Electronic Structure of Solids_ (Cambridge University Press, 2003).
  * Ruan and Ching (1987) Y.-C. Ruan and W. Y. Ching, Journal of Applied Physics 62, 2885 (1987), ISSN 0021-8979.
  * Lowell and Rose-Innes (1980) J. Lowell and A. C. Rose-Innes, Advances in Physics 29, 947 (1980), ISSN 0001-8732.
  * Annadi et al. (2013) A. Annadi, Q. Zhang, X. Renshaw Wang, N. Tuzla, K. Gopinadhan, W. M. Lü, A. Roy Barman, Z. Q. Liu, A. Srivastava, S. Saha, et al., Nature Communications 4, 1838 (2013).
  * Wang et al. (2020) Y. Wang, J. Balgley, E. Gerber, M. Gray, N. Kumar, X. Lu, J.-Q. Yan, A. Fereidouni, R. Basnet, S. J. Yun, et al., Nano Lett. 20, 8446 (2020), ISSN 1530-6984.
  * Sutter et al. (2009) P. Sutter, J. T. Sadowski, and E. Sutter, Phys. Rev. B 80, 245411 (2009).
  * Zheng et al. (2021a) X. Zheng, E. Gerber, J. Park, D. Werder, O. Kigner, E.-A. Kim, S. Xie, and D. G. Schlom, Applied Physics Letters 118, 093103 (2021a), ISSN 0003-6951.
  * Stradi et al. (2017) D. Stradi, L. Jelver, S. Smidstrup, and K. Stokbro, Journal of Physics: Condensed Matter 29, 185901 (2017), ISSN 0953-8984.
  * Lazić (2015) P. Lazić, Computer Physics Communications 197, 324 (2015), ISSN 0010-4655.
  * Kim et al. (2015) H.-S. Kim, V. S. V., A. Catuneanu, and H.-Y. Kee, Phys. Rev. B 91, 241110 (2015).
  * Manzeli et al. (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nature Reviews Materials 2, 17033 (2017).
  * Hsu et al. (2017) Y.-T. Hsu, A. Vaezi, M. H. Fischer, and E.-A. Kim, Nature Communications 8, 14985 (2017).
  * Zheng et al. (2021b) Y. Zheng, J. Gao, C. Han, and W. Chen, Cell Reports Physical Science 2, 100298 (2021b), ISSN 2666-3864.
  * Mao et al. (2020) J. Mao, S. P. Milovanović, M. Andjelković, X. Lai, Y. Cao, K. Watanabe, T. Taniguchi, L. Covaci, F. M. Peeters, A. K. Geim, et al., Nature 584, 215 (2020).
  * Shabani et al. (2021) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Nature Physics 17, 720 (2021).
  * Rizzo et al. (2022) D. J. Rizzo, S. Shabani, B. S. Jessen, J. Zhang, A. S. McLeod, C. Rubio-Verdú, F. L. Ruta, M. Cothrine, J. Yan, D. G. Mandrus, et al., Nano letters 22, 1946 (2022).
  * Gerber and Torrisi (2022) E. Gerber and S. B. Torrisi, _Basic intermatch scripts_ , Zenodo (2022).
  * Plumb et al. (2014) K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, Physical Review B 90, 041112 (2014).
  * Ong et al. (2013) S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Computational Materials Science 68, 314 (2013), ISSN 0927-0256.
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Schlipf and Gygi (2015) M. Schlipf and F. Gygi, Computer Physics Communications 196, 36 (2015), ISSN 0010-4655.
  * Freysoldt et al. (2009) C. Freysoldt, S. Boeck, and J. Neugebauer, Phys. Rev. B 79, 241103 (2009).
  * Sundararaman et al. (2017) R. Sundararaman, K. Letchworth-Weaver, K. A. Schwarz, D. Gunceler, Y. Ozhabes, and T. Arias, SoftwareX 6, 278 (2017), ISSN 2352-7110.
  * Payne et al. (1992) M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992).