consists of a grade algebra $H(\kappa)$, written also $H_{\kappa}$, for every
kind $\kappa\in\mathpzc{S}_{0}$, and a grade algebra homomorphism
$H_{\kappa,\mu}\colon H_{\kappa}\rightarrow H_{\mu}$ for every arrow
$\kappa\sqsubseteq\mu$, respecting composition and identities888 The notation
$H_{\kappa,\mu}$ makes sense, since between $\kappa$ and $\mu$ there is at
most one arrow. , that is, $\kappa\sqsubseteq\mu\sqsubseteq\nu$ implies
$H_{\kappa,\nu}=H_{\mu,\nu}\circ H_{\kappa,\mu}$ and
$H_{\kappa,\kappa}=\mathsf{id}_{H_{\kappa}}$.

Essentially, the homomorphisms $H_{\kappa,\mu}$ realise the refinement
$\kappa\sqsubseteq\mu$, transforming grades of kind $\kappa$ into grades of
kind $\mu$, preserving the grade algebra structure.

Observe that the arrows $I\sqsubseteq\kappa$ and
$\kappa\sqsubseteq\kappa\oplus\mu$ and $\mu\sqsubseteq\kappa\oplus\mu$ in
$\mathpzc{S}$ give rise to the following grade algebra homomorphisms:

$\mathsf{in}_{\kappa}^{H}=H_{I,\kappa}\colon H_{I}\rightarrow
H_{\kappa}\quad\mathsf{lh}_{\kappa,\mu}^{H}=H_{\kappa,\kappa\oplus\mu}\colon
H_{\kappa}\rightarrow
H_{\kappa\oplus\mu}\quad\mathsf{rh}_{\kappa,\mu}^{H}=H_{\mu,\kappa\oplus\mu}\colon
H_{\mu}\rightarrow H_{\kappa\oplus\mu}$

which provide us with a way to map grades of kind $I$ into grades of any other
kind, and grades of kind $\kappa$ and $\mu$ into grades of their composition
$\kappa\oplus\mu$. By functoriality of $H$ and using the commutative
idempotent monoid structure of $\mathpzc{S}$, we get the following equalities
hold in the category $\mathpzc{GrAlg}$, ensuring consistency of such
transformations:

$\displaystyle\mathsf{lh}_{\kappa\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}$ (1)
$\displaystyle\mathsf{rh}_{\kappa\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}\circ\mathsf{lh}_{\mu,\nu}^{H}$
(2) $\displaystyle\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{rh}_{\mu,\kappa}^{H}$ (3)
$\displaystyle\mathsf{lh}_{\kappa,\kappa}^{H}$
$\displaystyle=\mathsf{id}_{H_{\kappa}}$ (4)
$\displaystyle\mathsf{lh}_{\kappa,I}^{H}$
$\displaystyle=\mathsf{id}_{H_{\kappa}}$ (5)
$\displaystyle\mathsf{rh}_{\kappa,I}^{H}$
$\displaystyle=\mathsf{in}_{\kappa}^{H}$ (6)

In the following, we will show how to turn a heterogeneous grade algebra into
a single grade algebra. The procedure we will describe is based on a general
construction due to Grothendieck [16] defined on indexed categories.

Let us assume a grade signature $\mathpzc{S}$ and a heterogeneous grade
algebra $H\colon\mathpzc{S}\rightarrow\mathpzc{GrAlg}$. We consider the
following set:

$|G(H)|=\\{\langle\kappa,r\rangle\mid\kappa\in\mathpzc{S}_{0},\
r\in|H_{\kappa}|\\}$

That is, elements of $G(H)$ will be _kinded grades_ , namely, pairs of a kind
$\kappa$ and a grade of that kind. Note that this is indeed a set because
$\mathpzc{S}$ is small, that is, $\mathpzc{S}_{0}$ is a set. Then, we define a
binary relation $\preceq_{H}$ on $|G(H)|$ as follows:

$\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle\quad\text{iff}\quad\kappa\sqsubseteq\mu\text{
and }H_{\kappa,\mu}(r)\preceq_{\mu}s$

that is, the kind $\kappa$ must be more specific than the kind $\mu$ and,
transforming $r$ by $H_{\kappa,\mu}$, we obtain a grade of kind $\mu$ which is
smaller than $s$. These data define a partially ordered set as the following
proposition shows.

######  Proposition 5.5.

$\langle|G(H)|,\preceq_{H}\rangle$ is a partially ordered set.

The additive structure is given by a binary operation
$+_{H}\colon|G(H)|\times|G(H)|\rightarrow|G(H)|$ and an element
$\mathbf{0}_{H}$ in $|G(H)|$ defined as follows:

$\displaystyle\langle\kappa,r\rangle+_{H}\langle\mu,s\rangle$
$\displaystyle=\langle\kappa\oplus\mu,\mathsf{lh}_{\kappa,\mu}^{H}(r)+_{\kappa\oplus\mu}\mathsf{rh}_{\kappa,\mu}^{H}(s)\rangle$
$\displaystyle\mathbf{0}_{H}$ $\displaystyle=\langle I,\mathbf{0}_{I}\rangle$

This means that, the addition of two elements $\langle\kappa,r\rangle$ and
$\langle\mu,s\rangle$ is performed by first mapping $r$ and $s$ in the most
specific kind generalising both $\kappa$ and $\mu$, namely $\kappa\oplus\mu$,
and then by summing them in the grade algebra over that kind. The zero element
is just the zero of the most specific kind.

######  Proposition 5.6.

$\langle|G(H)|,\preceq_{H},+_{H},\mathbf{0}_{H}\rangle$ is an ordered
commutative monoid.

######  Proposition 5.7.

$\mathbf{0}_{H}\preceq_{H}\langle\kappa,r\rangle$ for every
$\langle\kappa,r\rangle\in|G(H)|$.

Similarly, the multiplicative structure is given by a binary operation
$\cdot_{H}\colon|G(H)|\times|G(H)|\rightarrow|G(H)|$ and an element
$\mathbf{1}_{H}$ in $|G(H)|$ defined as follows:

$\displaystyle\langle\kappa,r\rangle\cdot_{H}\langle\mu,s\rangle$
$\displaystyle=\begin{cases}\langle\kappa\oplus\mu,\mathsf{lh}_{\kappa,\mu}^{H}(r)\cdot_{\kappa\oplus\mu}\mathsf{rh}_{\kappa,\mu}^{H}(s)\rangle&\langle\kappa,r\rangle\neq\mathbf{0}_{H}\text{
and }\langle\mu,s\rangle\neq\mathbf{0}_{H}\\\
\mathbf{0}_{H}&\text{otherwise}\end{cases}$ $\displaystyle\mathbf{1}_{H}$
$\displaystyle=\langle I,\mathbf{1}_{I}\rangle$

Notice that the definitions above follow almost the same pattern as additive
operations, but we force that multiplying by $\mathbf{0}_{H}$ we get again
$\mathbf{0}_{H}$, which is a key property of grade algebras.

######  Proposition 5.8.

$\langle|G(H)|,\preceq_{H},\cdot_{H},\mathbf{1}_{H}\rangle$ is an ordered
monoid.

Altogether, we finally get the following result.

######  Theorem 5.9.

$G(H)=\langle|G(H)|,\preceq_{H},+_{H},\cdot_{H},\mathbf{0}_{H},\mathbf{1}_{H}\rangle$
is a grade algebra.

##  6 Grades as Java expressions

In Section 4 we described how a Java-like language could be equipped with
_grades_ decorating types, taken in an arbitrary grade algebra. Moreover, in
Section 5, we have shown that such grade algebra could have been obtained by
composing, in a specific way determined by providing a minimal collection of
grade homomorphisms, different grade algebras corresponding to different ways
to track usage of resources. In this way, _heterogeneous_ grades can coexist
in the same program, and do not need to be fixed once and for all, but can be
enriched by user-defined grade algebras and homomorphisms. In this section, we
consider the issue of providing a linguistic support to this aim. This could
be done by using an ad-hoc configuration language, however, we believe an
interesting solution is that the grade annotations in types could be written
themselves in Java.

The key idea is that such grade annotations are Java expressions of (classes
implementing) a predefined interface Grade, analogously to Java exceptions
which are expressions of (subclasses of) Exception. Moreover, grade
homomorphisms are user-defined as well. Namely, a user program can include:

  * •

pairs of _grade classes_ and _grade factory classes_ , each one modeling a
grade algebra desired for the specific application, with the factory class
providing its constants

  * •

_grade homomorphism classes_ , each one modeling a homomorphism from a grade
algebra (class) to another.

When typechecking code with grade annotations, the grades internally used by
the type system are those obtained by combining all the declared grade
algebras (classes) by means of the declared grade homomorphism classes, with
the construction described in Section 5. Operations on grades in the same
grade algebra (class) are derived from user-defined methods, as discussed more
in detail below, whereas operations on heterogeneous grades are derived as in
the construction in Section 5.

In the following we describe each point.

###  6.1 Grade and grade factory classes

They are classes implementing the following generic interfaces,
respectively:999We omit access modifiers to make the code lighter.

⬇

interface Grade<T extends Grade<T>> {

boolean leq(T x);

T sum(T x);

T mult(T x);

}

interface GradeFactory<T extends Grade<T>>{

T zero();

T one();

}

For instance, considering Example 5.1, we show below the complete
implementation of affinity grades, and only the skeletons of the other
classes, for brevity:

⬇

abstract class Affinity implements Grade<Affinity>{}

class AffinityFactory implements GradeFactory<Affinity>{

Affinity zero(){return new ZeroA();}

Affinity one(){return new One();}

}

class ZeroA extends Affinity{

boolean leq(Affinity x){return true;}

Affinity sum(Affinity x){return x;}

Affinity mult(Affinity x){return this;}

}

class One extends Affinity{

boolean leq(Affinity x){return !(x instanceof ZeroA);}

Affinity sum(Affinity x){

if (x instanceof ZeroA) return this;

return new Omega();

}

Affinity mult(Affinity x){return x;}

}

class Omega extends Affinity{

boolean leq(Affinity x){return (x instanceof Omega);}

Affinity sum(Affinity x){return this;}

Affinity mult(Affinity x){

if (x instanceof ZeroA) return x;

return this;

}

}

abstract class Privacy implements Grade<Privacy> {...}

class ZeroP extends Privacy{...}

class Private extends Privacy{...}

class Public extends Privacy{...}

class PrivacyFactory implements GradeFactory<Privacy>{...}

abstract class PPrivacy implements Grade<PPrivacy> {...}

class ZeroPP extends PPrivacy{...}

class LevelA extends PPrivacy{...}

...

class LevelD extends PPrivacy{...}

class PPrivacyFactory implements GradeFactory<PPrivacy>{...}

class AffinityPrivacy implements Grade<AffinityPrivacy>{

Affinity first; Privacy second;

...

}

class AffinityPrivacyFactory implements GradeFactory<AffinityPrivacy>{

...

}

Note that for Affinity, Privacy, and PPrivacy we could have used an enum class
with constants as well; here we preferred to use a hierarchy of classes, which
is a core object-oriented feature provided in the FJ calculus. More in
general, note that the implementation, as expected, depends very much on the
features of the target language; for instance, in Haskell we could express
Grade as a typeclass, and grade algebras as its instances, and the same for
GradeHom below.

###  6.2 Grade homomorphism classes

They are classes implementing the following generic interface:

⬇

interface GradeHom<T extends Grade<T>, R extends Grade<R>> {

R apply(T x);

}

For instance, the grade homomorphisms of Example 5.2 can be implemented as
follows:

⬇

class PPtoPGradeHom implements GradeHom<PPrivacy,Privacy>{

Privacy apply (PPrivacy p){

if (p instanceof LevelD) return new Public();

return new Private();

}

}

class APtoAGradeHom implements GradeHom<AffinityPrivacy,Affinity>{

Affinity apply (AffinityPrivacy ap){return ap.first;}

}

class APtoPGradeHom implements GradeHom<AffinityPrivacy,Privacy>{

Privacy apply (AffinityPrivacy ap){return ap.second;}

}

###  6.3 Predefined grade classes

The grade algebra $\mathsf{Nat}$ can be implemented by the following classes,
assumed to be predefined.

⬇

abstract class Nat implements Grade<Nat> {

abstract <R extends Grade<R>> R transform(GradeFactory<R> factory);

}

class Zero extends Nat {

<R extends Grade<R>> R transform(GradeFactory<R> factory){

return factory.zero();

}

boolean leq(Nat x){return true;}

Nat sum(Nat x){return x;}

Nat mult(Nat x){return this;}

}

class Succ extends Nat {

Nat pred;

Succ(Nat pred){this.pred=pred;}

<R extends Grade<R>> R transform(GradeFactory<R> factory) {

return factory.one().sum(pred.transform(factory));

}

boolean leq(Nat x){

if (x instanceof Zero){ return false;}

return pred.leq(((Succ) x).pred);

}

Nat sum(Nat x){return new Succ(pred.sum(x));}

Nat mult(Nat x){return pred.mult(x).sum(x);}

}

class NatFactory implements GradeFactory<Nat>{

Nat zero(){return new Zero();}

Nat one(){return new Succ(zero());}

}

class NtoGradeHom<R extends Grade<R>> implements GradeHom<Nat,R> {

GradeFactory<R> factory;

NtoGradeHom(GradeFactory<R> factory){this.factory=factory;}

R apply(Nat n) {return n.transform(this.factory);}

}

Order, sum, and multiplication on natural numbers are defined in the expected
inductive way. An object of class NtoGradeHom<R> models the unique
homomorphism $\iota_{\mathit{R}}$ from $\mathsf{Nat}$ to the grade algebra
(modeled by) R, which maps a natural number $n$ to the sum in $\mathit{R}$ of
$n$ copies of $\mathbf{1}_{\mathit{R}}$, as defined at page 2.2. To achieve
this, the object is created specifying a constructor argument of the
corresponding factory class.101010 To this end, we need a convention to
construct names of factory classes. When the homomorphism is applied to a
natural number n (invoking apply with argument n), the auxiliary method
transform constructs the object of class R corresponding to n, using the zero
and one method provided by the factory.

Analogously, for the grade algebra $\mathsf{Triv}$ we have:

⬇

class Triv implements Grade<Triv>{

boolean leq(Triv t){return true;}

Triv sum(Triv t){return this;}

Triv mult(Triv t){return this;}

}

class TrivFactory implements GradeFactory<Triv>{

Triv zero(){return new Triv();}

Triv one(){return new Triv();}

}

class toTrivGradeHom<R extends Grade<R>>implements GradeHom<R,Triv>{

Triv apply(R r){return new Triv();}

}

###  6.4 Grade annotations as Java expressions

Assuming to have the classes decribed above, the programmer could write code
with grade annotations being Java expressions. For instance, Example 4.1 could
be written as follows, where annotations are in square brackets:

⬇

class Pair { A[new One()] first; A[new One()] second;

A[new ZeroA()] getFirstZero() [new One()]{this.first}

A[new One()] getFirstAffine() [new One()]{this.first}

A[new Omega()] getFirst() [new One()]{this.first}

}

Typechecking could then be performed in two steps:

  1. 1.

Code defining grades, which is assumed to be standard (that is, non-graded)
Java code, is typechecked by the standard compiler.

  2. 2.

Graded code (containing grade annotations written in Java) is typechecked
accordingly to the graded type system in Figure 4, where the underlying grade
algebra is obtained by composing, by the construction described in Section 5,
the user-defined grade algebras through the user-defined grade homomorphisms.
Each user-defined algebra has as carrier (set of grades) the Java values which
are instances of the corresponding class, and the operations are computed by
executing user-defined methods in such class. For instance, to compute the sum
$\mathit{v}_{1}+\mathit{v}_{2}$ of two grades which are values of a grade
class, we evaluate
$\mathit{v}_{1}{\texttt{.}}\texttt{sum}\texttt{(}\mathit{v}_{2}\texttt{)}$.
Analogously to compute the result of a grade homomorphism.

For the whole process to work correctly, the following are responsabilities of
the programmer:

  * •

Grade classes, grade factory classes, and grade homomorphism classes should
satisfy the axioms required for the structures they model, e.g., that the sum
derived from sum methods is commutative and associative. The same happens, for
instance, in Haskell, when one defines instances of Functor or Monad.

  * •

Code defining grades should be _terminating_ , since, as described above, the
second typechecking step requires to _execute_ code typechecked in the first
step.

  * •

Finally, the relation among grade classes implicitly defined by declaring
grade homomorphism classes should actually be a direct refinement relation,
that is, should satisfy the two requirements: (1) there exists at most one
path between two grade classes, and (2) each two grade classes with a common
ancestor have a _least_ common ancestor. These are requirements easy to check,
similarly to the check that inheritance is acyclic, or that there are no
diamonds in multiple inheritance.

An interesting point is that implementations could use in a parametric way
auxiliary tools, notably a termination checker to prevent divergence in
methods implementing grade operations, and/or a verifer to ensure that they
provide the required properties.

##  7 Related work

The two contributions which have been more inspiring for the work in this
paper are the instrumented semantics proposed in [7] and the Granule language
[22]. In [7], the authors develop GraD, a graded dependent type system that
includes functions, tensor products, additive sums, and a unit type. Moreover,
they define an instrumented operational semantics which tracks usage of
resources, and prove that the graded type system is sound with respect to such
instrumented semantics. In this paper, we take the same approach to define a
resource-aware semantics, parametric on an arbitrary grade algebra. However,
differently from [7], where such semantics is defined on typed terms, with the
only aim to show the role of the type system, the definition of our semantics
is given _independently_ from the type system, as is the standard approach in
calculi. That is, the aim is also to provide a simple purely semantic model
which takes into account usage of resources.

Granule [22] is a functional language equipped with graded modal types, where
different kinds of grades can be used at the same time, including naturals for
exact usage, security levels, intervals, infinity, and products of coeffects.
We owe to Granule the idea of allowing different kinds of grades to coexist,
and the overall objective to exploit graded modal types in a programming
language. Concerning heterogeneous grades, in this paper we push forward the
Granule approach, since we do not want this grade algebra to be fixed, but
extendable by the programmer with user-defined grades. To this aim we define
the construction in Section 5. Concerning the design of a graded programming
language, here we investigate the object-oriented rather than functional
paradigm, taking some solutions which seem more adequate in that context,
e.g., to have once-graded types and no boxing/unboxing. The design and
implementation of a real Java-like language are not objectives of the current
paper; however, we outline in Section 6 a possible interesting solution, where
grade annotations are written in the language itself.

Coming more in general to resource-aware type systems, coeffects were first
introduced by [23] and further analyzed by [24]. In particular, [24] develops
a generic coeffect system which augments the simply-typed $\lambda$-calculus
with context annotations indexed by _coeffect shapes_. The proposed framework
is very abstract, and the authors focus only on two opposite instances:
structural (per-variable) and flat (whole context) coeffects, identified by
specific choices of context shapes.

Most of the subsequent literature on coeffects focuses on structural ones, for
which there is a clear algebraic description in terms of semirings. This was
first noticed by [6], who developed a framework for structural coeffects for a
functional language. This approach is inspired by a generalization of the
exponential modality of linear logic, see, e.g., [5]. That is, the distinction
between linear and unrestricted variables of linear systems is generalized to
have variables decorated by coeffects, or grades, that determine how much they
can be used. In this setting, many advances have been made to combine
coeffects with other programming features, such as computational effects [13,
22, 10], dependent types [2, 7, 21], and polymorphism [1]. Other graded type
systems are explored in [2, 14, 1], also combining effects and coeffects [13,
22]. In all these papers, the process of tracking usage through grades is a
powerful method of instrumenting type systems with analyses of irrelevance and
linearity that have practical benefits like erasure of irrelevant terms
(resulting in speed-up) and compiler optimizations (such as in-place update).

As already mentioned, [21] and [26] observed that contexts in a structural
coeffect system form a module over the semiring of grades, event though they
do not use this structure in its full generality, restricting themselves to
free modules, that is, to structural coeffect systems. Recently, [4] shows a
significant non-structural instance, namely, a coeffect system to track
sharing in the imperative paradigm.

##  8 Conclusion

The contributions of the paper can be summarized as follows:

  * •

Resource-aware extension of FJ reduction, parametric on an arbitrary grade
algebra.

  * •

Resource-aware extension of the type system, proved to ensure soundness-may of
the resource-aware semantics.

  * •

Formal construction which, given grades of different kinds and grade
transformations corresponding to a refinement relation among kinds (formally,
a functor over a grade signature), provides a grade algebra of _heterogeneous
grades_.

  * •

Notion of _direct refinement_ allowing a minimal and easy to check way to
specify the above functor.

  * •

Outline of a Java extension where grades are user-defined, and grade
annotations are written in the language itself.

As already noted, the key novel ideas in the contributions above are mostly
independent from the language. So, a first natural direction for future work
is to explore their incarnation in another paradigm, e.g., the functional one.
That would include the definition of a parametric resource-aware reduction
independent from types, the design of a type system with once-graded types,
and possibly the design of user-defined grades in a functional language, e.g.,
in Haskell by relying on the typeclass feature. Though the overall approach
should still apply, we expect the investigation to be significant due to the
specific features of the paradigm.

The resource-aware operational semantics defined in this paper requires
_annotations_ in subterms, with the only aim to fix their reduction grade in
the reduction of the enclosing term. As mentioned in Section 3, adopting a
big-step style would clearly remove the need of such technical artifice; only
annotations in constructor subterms should be kept, since they express a true
constraint on the semantics. Thus, a very interesting alternative to be
studied is a big-step version of resource-aware semantics, allowing a more
abstract and clean presentation. With this choice, we should employ, to prove
soundness-may, the techniques recently introduced in [9, 8].

Coming back to Java-like languages, the FJ language considered in the paper
does not include imperative features. Adding mutable memory leads to many
significant research directions. First, besides the model presented in this
paper, and in general in literature, where “using a resource” means “replacing
a variable with its value”, another view is possible where the resource is the
memory and “using” means “interacting with the memory”. Moreover, we would
like to investigate more in detail how to express by grade algebras forms of
usages which are typical of the imperative paradigm, such as the readonly
modifier, and, more in general, _capabilities_ [17, 15].

## References

  * [1] Andreas Abel and Jean-Philippe Bernardy.  A unified view of modalities in type systems.  Proceedings of ACM on Programming Languages, 4(ICFP):90:1–90:28, 2020.  doi:10.1145/3408972.
  * [2] Robert Atkey.  Syntax and semantics of quantitative type theory.  In Anuj Dawar and Erich Grädel, editors, IEEE Symposium on Logic in Computer Science, LICS 2018, pages 56–65. ACM Press, 2018.  doi:10.1145/3209108.3209189.
  * [3] Riccardo Bianchini, Francesco Dagnino, Paola Giannini, and Elena Zucca.  A Java-like calculus with user-defined coeffects.  In Ugo Dal Lago and Daniele Gorla, editors, ICTCS’22 - Italian Conference on Theoretical Computer Science, volume 3284 of CEUR Workshop Proceedings, pages 66–78. CEUR-WS.org, 2022.
  * [4] Riccardo Bianchini, Francesco Dagnino, Paola Giannini, Elena Zucca, and Marco Servetto.  Coeffects for sharing and mutation.  Proceedings of ACM on Programming Languages, 6(OOPSLA):870–898, 2022.  doi:10.1145/3563319.
  * [5] Flavien Breuvart and Michele Pagani.  Modelling coeffects in the relational semantics of linear logic.  In Stephan Kreutzer, editor, 24th EACSL Annual Conference on Computer Science Logic, CSL 2015, volume 41 of LIPIcs, pages 567–581. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2015.  doi:10.4230/LIPIcs.CSL.2015.567.
  * [6] Aloïs Brunel, Marco Gaboardi, Damiano Mazza, and Steve Zdancewic.  A core quantitative coeffect calculus.  In Zhong Shao, editor, European Symposium on Programming, ESOP 2013, volume 8410 of Lecture Notes in Computer Science, pages 351–370. Springer, 2014.  doi:10.1007/978-3-642-54833-8\\_19.
  * [7] Pritam Choudhury, Harley Eades III, Richard A. Eisenberg, and Stephanie Weirich.  A graded dependent type system with a usage-aware semantics.  Proceedings of ACM on Programming Languages, 5(POPL):1–32, 2021\.  doi:10.1145/3434331.
  * [8] Francesco Dagnino.  A meta-theory for big-step semantics.  ACM Transactions on Computational Logic, 23(3):20:1–20:50, 2022\.  doi:10.1145/3522729.
  * [9] Francesco Dagnino, Viviana Bono, Elena Zucca, and Mariangiola Dezani-Ciancaglini.  Soundness conditions for big-step semantics.  In Peter Müller, editor, European Symposium on Programming, ESOP 2020, volume 12075 of Lecture Notes in Computer Science, pages 169–196. Springer, 2020.  doi:10.1007/978-3-030-44914-8\\_7.
  * [10] Ugo Dal Lago and Francesco Gavazzo.  A relational theory of effects and coeffects.  Proceedings of ACM on Programming Languages, 6(POPL):1–28, 2022\.  doi:10.1145/3498692.
  * [11] Rocco De Nicola and Matthew Hennessy.  Testing equivalences for processes.  Theoretical Computer Science, 34(1):83 – 133, 1984.  doi:https://doi.org/10.1016/0304-3975(84)90113-0.
  * [12] Werner Dietl, Sophia Drossopoulou, and Peter Müller.  Generic universe types.  In Erik Ernst, editor, European Conference on Object-Oriented Programming, ECOOP 2007, volume 4609 of Lecture Notes in Computer Science, pages 28–53. Springer, 2007.
  * [13] Marco Gaboardi, Shin-ya Katsumata, Dominic A. Orchard, Flavien Breuvart, and Tarmo Uustalu.  Combining effects and coeffects via grading.  In Jacques Garrigue, Gabriele Keller, and Eijiro Sumii, editors, ACM International Conference on Functional Programming, ICFP 2016, pages 476–489. ACM Press, 2016.  doi:10.1145/2951913.2951939.
  * [14] Dan R. Ghica and Alex I. Smith.  Bounded linear types in a resource semiring.  In Zhong Shao, editor, European Symposium on Programming, ESOP 2013, volume 8410 of Lecture Notes in Computer Science, pages 331–350. Springer, 2014.  doi:10.1007/978-3-642-54833-8\\_18.
  * [15] Colin S. Gordon.  Designing with static capabilities and effects: Use, mention, and invariants (pearl).  In Robert Hirschfeld and Tobias Pape, editors, European Conference on Object-Oriented Programming, ECOOP 2020, volume 166 of LIPIcs, pages 10:1–10:25. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020.  doi:10.4230/LIPIcs.ECOOP.2020.10.
  * [16] Alexander Grothendieck.  Catégories fibrées et descente.  In Revêtements étales et groupe fondamental, pages 145–194. Springer, 1971.
  * [17] Philipp Haller and Martin Odersky.  Capabilities for uniqueness and borrowing.  In Theo D’Hondt, editor, European Conference on Object-Oriented Programming, ECOOP 2010, volume 6183 of Lecture Notes in Computer Science, pages 354–378. Springer, 2010.
  * [18] Atsushi Igarashi, Benjamin C. Pierce, and Philip Wadler.  Featherweight Java: A minimal core calculus for Java and GJ.  In ACM Symp. on Object-Oriented Programming: Systems, Languages and Applications 1999, pages 132–146. ACM Press, 1999.  doi:10.1145/320384.320395.
  * [19] Saunders Mac Lane.  Categories for the working mathematician, volume 5.  Springer Science & Business Media, 2013.
  * [20] Daniel Marshall, Michael Vollmer, and Dominic Orchard.  Linearity and uniqueness: An entente cordiale.  In Ilya Sergey, editor, European Symposium on Programming, ESOP 2022, volume 13240 of Lecture Notes in Computer Science, pages 346–375. Springer, 2022.  doi:10.1007/978-3-030-99336-8\\_13.
  * [21] Conor McBride.  I got plenty o’ nuttin’.  In Sam Lindley, Conor McBride, Philip W. Trinder, and Donald Sannella, editors, A List of Successes That Can Change the World - Essays Dedicated to Philip Wadler on the Occasion of His 60th Birthday, volume 9600 of Lecture Notes in Computer Science, pages 207–233. Springer, 2016.  doi:10.1007/978-3-319-30936-1\\_12.
  * [22] Dominic Orchard, Vilem-Benjamin Liepelt, and Harley Eades III.  Quantitative program reasoning with graded modal types.  Proceedings of ACM on Programming Languages, 3(ICFP):110:1–110:30, 2019.  doi:10.1145/3341714.
  * [23] Tomas Petricek, Dominic A. Orchard, and Alan Mycroft.  Coeffects: Unified static analysis of context-dependence.  In Fedor V. Fomin, Rusins Freivalds, Marta Z. Kwiatkowska, and David Peleg, editors, Automata, Languages and Programming, ICALP 2013, volume 7966 of Lecture Notes in Computer Science, pages 385–397. Springer, 2013.  doi:10.1007/978-3-642-39212-2\\_35.
  * [24] Tomas Petricek, Dominic A. Orchard, and Alan Mycroft.  Coeffects: a calculus of context-dependent computation.  In Johan Jeuring and Manuel M. T. Chakravarty, editors, ACM International Conference on Functional Programming, ICFP 2014, pages 123–135. ACM Press, 2014.  doi:10.1145/2628136.2628160.
  * [25] Emily Riehl.  Category theory in context.  Courier Dover Publications, 2017.
  * [26] James Wood and Robert Atkey.  A framework for substructural type systems.  In Ilya Sergey, editor, European Symposium on Programming, ESOP 2022, volume 13240 of Lecture Notes in Computer Science, pages 376–402. Springer, 2022.  doi:10.1007/978-3-030-99336-8\\_14.

##  Appendix A Proofs of Section 5

######  Proof A.1 (Proof of Proposition 2.3).

Item 2 is straightforward as the singleton set is a terminal object in the
category of sets and functions. Towards a proof of Item 1, let
$f\colon\mathsf{Nat}\rightarrow\mathit{R}$ be a grade algebra homomorphism and
note that, since $n=1+\cdots+1$ ($n$ times), for all $n\in\mathbb{N}$, and $f$
preserves sums and the unit, we get
$f(n)=f(1)+_{\mathit{R}}\cdots+_{\mathit{R}}f(1)=\mathbf{1}_{\mathit{R}}+_{\mathit{R}}\cdots+_{\mathit{R}}\mathbf{1}_{\mathit{R}}$
($n$ times). That is, we have $f(n)=\iota_{\mathit{R}}(n)$, for all
$n\in\mathbb{N}$. Therefore, to conclude, we just have to show that the map
$\iota_{\mathit{R}}$ is a grade algebra homomorphism. The fact that
$\iota_{\mathit{R}}(0)=\mathbf{0}_{\mathit{R}}$ and
$\iota_{\mathit{R}}(1)=\mathbf{1}_{\mathit{R}}$ is immediate. The fact that
$\iota_{\mathit{R}}(n+m)=\iota_{\mathit{R}}(n)+_{\mathit{R}}\iota_{\mathit{R}}(m)$
and ${\iota_{\mathit{R}}(n\cdot
m)}=\iota_{\mathit{R}}(n)\cdot_{\mathit{R}}\iota_{\mathit{R}}(m)$ follows from
a straightforward induction on $n$, using distributivity and nullity
properties of the grade algebra $\mathit{R}$. Finally, to prove monotonicity,
consider $n\leq m$ and proceed by induction on $m-n$. If $m-n=0$, then $n=m$
and so the thesis is trivial. If $m-n=k+1$, we have $m-(n+1)=k$, then by
induction hypothesis we get
$\iota_{\mathit{R}}(n+1)\preceq_{\mathit{R}}(\iota_{\mathit{R}}(m)$. Since
$\iota_{\mathit{R}}(n+1)=\iota_{\mathit{R}}(n)+_{\mathit{R}}\iota_{\mathit{R}}(1)$
and $\mathbf{0}_{\mathit{R}}\preceq_{\mathit{R}}\iota_{\mathit{R}}(1)$, we get
$\iota_{\mathit{R}}(n)=\iota_{\mathit{R}}(n)+_{\mathit{R}}\mathbf{0}_{\mathit{R}}\preceq_{\mathit{R}}\iota_{\mathit{R}}(n)+_{\mathit{R}}\iota_{\mathit{R}}(1)\preceq_{\mathit{R}}\iota_{\mathit{R}}(m)$,
as needed.

######  Proof A.2 (Proof of Proposition 5.5).

We have to prove that $\preceq_{H}$ is reflexive, transitive and
antisymmetric.

Given an element $\langle\kappa,r\rangle\in|G(H)|$, since
$\kappa\sqsubseteq\kappa$ and $H_{\kappa,\kappa}=\mathsf{id}_{H_{\kappa}}$, by
functoriality of $H$, we have

$H_{\kappa,\kappa}(r)=\mathsf{id}_{H_{\kappa}}(r)=r$

hence $\langle\kappa,r\rangle\preceq_{H}\langle\kappa,r\rangle$, which proves
reflexivity.

Given
$\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle\preceq_{H}\langle\nu,t\rangle$,
we know that $\kappa\sqsubseteq\mu\sqsubseteq\nu$ and
$H_{\kappa,\mu}(r)\preceq_{\mu}s$ and $H_{\mu,\nu}(s)\preceq_{\nu}t$ and, by
functoriality of $H$, $H_{\kappa,\nu}=H_{\mu,\nu}\circ H_{\kappa,\mu}$.
Therefore, we get

$\displaystyle H_{\kappa,\nu}(r)$
$\displaystyle=H_{\mu,\nu}(H_{\kappa,\mu}(r))\preceq_{\nu}H_{\mu,\nu}(s)\preceq_{\nu}t$

hence $\langle\kappa,r\rangle\preceq_{H}\langle\nu,t\rangle$, which proves
transitivity.

Given
$\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle\preceq_{H}\langle\kappa,r\rangle$,
we know that $\kappa\sqsubseteq\mu\sqsubseteq\kappa$ and
$H_{\kappa,\mu}(r)\preceq_{\mu}s$ and $H_{\mu,\kappa}(s)\preceq_{\kappa}r$.
Since $\sqsubseteq$ is antisymmetric, we get $\kappa=\mu$, hence
$H_{\kappa,\mu}=H_{\mu,\kappa}=H_{\kappa,\kappa}$, which, and, by
functoriality of $H$, is equal to the identity $\mathsf{id}_{H_{\kappa}}$.
This implies that $r\preceq_{\kappa}s$ and $s\preceq_{\kappa}r$, which implies
$r=s$ by antisymmetry of $\preceq_{\kappa}$.

######  Proof A.3 (Proof of Proposition 5.6).

We check the four properties.

Monotonicity of $+_{H}$

Consider $\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle$ and
$\langle\kappa^{\prime},r^{\prime}\rangle\preceq_{H}\langle\mu^{\prime},s^{\prime}\rangle$,
then we know that $\kappa\sqsubseteq\mu$ and
$\kappa^{\prime}\sqsubseteq\mu^{\prime}$ and $H_{\kappa,\mu}(r)\preceq_{\mu}s$
and
$H_{\kappa^{\prime},\mu^{\prime}}(r^{\prime})\preceq_{\mu^{\prime}}s^{\prime}$.
By monotonicity of $\mathsf{lh}_{\mu,\mu^{\prime}}^{H}$ and
$\mathsf{rh}_{\mu,\mu^{\prime}}^{H}$ and $+_{\mu\oplus\mu^{\prime}}$, we get

$\mathsf{lh}_{\mu,\mu^{\prime}}^{H}(H_{\kappa,\mu}(r))+_{\mu\oplus\mu^{\prime}}\mathsf{rh}_{\mu,\mu^{\prime}}^{H}(H_{\kappa^{\prime},\mu^{\prime}}(r^{\prime}))\preceq_{\mu\oplus\mu^{\prime}}\mathsf{lh}_{\mu,\mu^{\prime}}^{H}(s)+_{\mu\oplus\mu^{\prime}}\mathsf{rh}_{\mu,\mu^{\prime}}^{H}(s^{\prime})$

By monotonicity of $\oplus$, we get
$\kappa\oplus\kappa^{\prime}\sqsubseteq\mu\oplus\mu^{\prime}$, then, by
functoriality of $H$ and by definition of $\mathsf{lh}^{H}$ and $injr$, we
have

$\mathsf{lh}_{\mu,\mu^{\prime}}^{H}\circ
H_{\kappa,\mu}=H_{\kappa\oplus\kappa^{\prime},\mu\oplus\mu^{\prime}}\circ\mathsf{lh}_{\kappa,\kappa^{\prime}}^{H}\qquad\mathsf{rh}_{\mu,\mu^{\prime}}^{H}\circ
H_{\kappa^{\prime},\mu^{\prime}}=H_{\kappa\oplus\kappa^{\prime},\mu\oplus\mu^{\prime}}\circ\mathsf{rh}_{\kappa,\kappa^{\prime}}^{H}$

therefore, we get

$\begin{array}[]{l}H_{\kappa\oplus\kappa^{\prime},\mu\oplus\mu^{\prime}}(\mathsf{lh}_{\kappa,\kappa^{\prime}}^{H}(r)+_{\kappa\oplus\kappa^{\prime}}\mathsf{rh}_{\kappa,\kappa^{\prime}}^{H}(r^{\prime}))\preceq_{\mu\oplus\mu^{\prime}}\\\
\mathsf{lh}_{\mu,\mu^{\prime}}^{H}(H_{\kappa,\mu}(r))+_{\mu\oplus\mu^{\prime}}\mathsf{rh}_{\mu,\mu^{\prime}}^{H}(H_{\kappa^{\prime},\mu^{\prime}}(r^{\prime}))\preceq_{\mu\oplus\mu^{\prime}}\\\
\mathsf{lh}_{\mu,\mu^{\prime}}^{H}(s)+_{\mu\oplus\mu^{\prime}}\mathsf{rh}_{\mu,\mu^{\prime}}^{H}(s^{\prime})\end{array}$

This proves
$\langle\kappa,r\rangle+_{H}\langle\kappa^{\prime},r^{\prime}\rangle\preceq_{H}\langle\mu,s\rangle+_{H}\langle\mu^{\prime},s^{\prime}\rangle$,
as needed.

Associativity of $+_{H}$

Consider elements $\langle\kappa,r\rangle$, $\langle\mu,s\rangle$ and
$\langle\nu,t\rangle$ in $|G(H)|$. Using Equations 1 and 2, we have the
following

$\displaystyle(\langle\kappa,r\rangle+_{H}\langle\mu,s\rangle)+_{H}\langle\nu,t\rangle$
$\displaystyle=\langle\kappa\oplus\mu,\mathsf{lh}_{\kappa,\mu}^{H}(r)+\mathsf{rh}_{\kappa,\mu}^{H}(s)\rangle+_{H}\langle\nu,t\rangle$
$\displaystyle=\langle(\kappa\oplus\mu)\oplus\nu,\mathsf{lh}_{\kappa\oplus\mu,\nu}^{H}(\mathsf{lh}_{\kappa,\mu}^{H}(r)+\mathsf{rh}_{\kappa,\mu}^{H}(s))+\mathsf{rh}_{\kappa\oplus\mu,\nu}^{H}(t)\rangle$
$\displaystyle=\langle\kappa\oplus(\mu\oplus\nu),\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}(r)+(\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}(\mathsf{lh}_{\mu,\nu}^{H}(s)+\mathsf{rh}_{\mu,\nu}^{H}(t))\rangle$
$\displaystyle=\langle\kappa,r\rangle+_{H}(\langle\mu,s\rangle+_{H}\langle\nu,t\rangle)$

Commutativity of $+_{H}$

Consider elements $\langle\kappa,r\rangle$ and $\langle\mu,s\rangle$ in
$|G(H)|$. Using Equation 3, we get the following

$\displaystyle\langle\kappa,r\rangle+_{H}\langle\mu,s\rangle$
$\displaystyle=\langle\kappa\oplus\mu,\mathsf{lh}_{\kappa,\mu}^{H}(r)+_{\kappa\oplus\mu}\mathsf{rh}_{\kappa,\mu}^{H}(s)\rangle=\langle\kappa\oplus\mu,\mathsf{rh}_{\kappa,\mu}^{H}(s)+_{\kappa\oplus\mu}\mathsf{lh}_{\kappa,\mu}^{H}(r)\rangle$
$\displaystyle=\langle\mu\oplus\kappa,\mathsf{lh}_{\mu,\kappa}^{H}(s)+_{\mu\oplus\kappa}\mathsf{rh}_{\mu,\kappa}^{H}(r)\rangle=\langle\mu,s\rangle+_{H}\langle\kappa,r\rangle$

Neutrality of $\mathbf{0}_{H}$

Consider $\langle\kappa,r\rangle$ in $|G(H)|$. Using Equations 5, 6 and 3, we
the the following

$\displaystyle\langle\kappa,r\rangle+_{H}\mathbf{0}_{H}$
$\displaystyle=\langle\kappa\oplus
I,\mathsf{lh}_{\kappa,I}^{H}(r)+_{\kappa\oplus
I}\mathsf{rh}_{\kappa,I}^{H}(\mathbf{0}_{I})\rangle=\langle\kappa,r+_{\kappa}\mathsf{in}_{\kappa}^{H}(\mathbf{0}_{I})\rangle$
$\displaystyle=\langle\kappa,r+_{\kappa}\mathbf{0}_{\kappa}\rangle=\langle\kappa,r\rangle$

######  Proof A.4 (Proof of Proposition 5.7).

Since $\mathsf{in}_{\kappa}^{H}$ is a grade algebra homomorphism, we have
$\mathsf{in}_{\kappa}^{H}(\mathbf{0}_{I})=\mathbf{0}_{\kappa}$ and by
definition of grade algebra we have $\mathbf{0}_{\kappa}\preceq_{\kappa}r$.
Therefore, $\mathsf{in}_{\kappa}^{H}(\mathbf{0}_{I}\preceq_{\kappa}r$, which
proves the thesis by definition of $\preceq_{H}$.

######  Proof A.5 (Proof of Proposition 5.8).

To prove the equational axioms of monoid (associativity, and neutrality) the
proof is the same as Proposition 5.6 when all the involved elements are
different from $\mathbf{0}_{I}$, and it is trivial otherwise. Indeed, if one
of such elements is $\mathbf{0}_{I}$, then the whole multiplication gives
$\mathbf{0}_{I}$ by definition.

To prove monotonicity of $\cdot_{H}$, consider
$\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle$ and
$\langle\kappa^{\prime},r^{\prime}\rangle\preceq_{H}\langle\mu^{\prime},s^{\prime}\rangle$
in $|G(H)|$. If they are all different from $\mathbf{0}_{I}$, the proof goes
as in Proposition 5.6. If either $\langle\kappa,r\rangle=\mathbf{0}_{I}$ or
$\langle\kappa^{\prime},r^{\prime}\rangle=\mathbf{0}_{I}$, then
$\langle\kappa,r\rangle\cdot_{H}\langle\kappa^{\prime},r^{\prime}\rangle=\mathbf{0}_{I}$
and so the thesis follows by Proposition 5.7. If
$\langle\mu,s\rangle=\mathbf{0}_{I}$, then
$\langle\kappa,r\rangle\preceq_{H}\langle\mu,s\rangle=\langle
I,\mathbf{0}_{I}\rangle$ implies $\kappa\sqsubseteq I$ and
$H_{\kappa,I}(r)\preceq_{I}\mathbf{0}_{I}$ and, since $I\sqsubseteq\kappa$ by
definition of grade signature, we get $\kappa=I$. Therefore, by functoriality
of $H$, we have $H_{\kappa,I}=\mathsf{id}_{H_{I}}$, hence
$r\preceq_{I}\mathbf{0}_{I}$ which implies $r=\mathbf{0}_{I}$, since $H_{I}$
is a grade algebra. This proves that $\langle\kappa,r\rangle=\langle
I,\mathbf{0}_{I}\rangle=\mathbf{0}_{H}$ and so we get
$\langle\kappa,r\rangle\cdot_{H}\langle\kappa^{\prime},r^{\prime}\rangle=\mathbf{0}_{H}=\langle\mu,s\rangle\cdot_{H}\langle\mu^{\prime},s^{\prime}\rangle$,
as needed. Finally, the case $\langle\mu^{\prime},s^{\prime}\rangle=\langle
I,\mathbf{0}_{I}\rangle$ is analogous, hence we get the thesis.

######  Proof A.6 (Proof of Theorem 5.9).

By Propositions 5.6 and 5.8 we have both the additive and multiplicative
monoid structures. Proposition 5.7 proves that $\mathbf{0}_{H}$ is the least
element of the order $\preceq_{H}$. The fact that multiplying by
$\mathbf{0}_{I}$ we get again $\mathbf{0}_{I}$ holds by definition. Hence, it
remains to prove that $\cdot_{H}$ distributes over $+_{H}$. To this end,
consider $\langle\kappa,r\rangle$, $\langle\mu,s\rangle$ and
$\langle\nu,t\rangle$ in $|G(H)|$ and assume they are all different from
$\mathbf{0}_{I}$. Using Equations 1, 2, 3 and 4, we get the following
equations:

$\displaystyle\mathsf{lh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{(\kappa\oplus\mu)\oplus\kappa,\nu}^{H}\circ\mathsf{lh}_{\kappa\oplus\mu,\kappa}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa\oplus(\kappa\oplus\mu),\nu}^{H}\circ\mathsf{rh}_{\kappa,\kappa\oplus\mu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{(\kappa\oplus\kappa)\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa\oplus\kappa,\mu}^{H}\circ\mathsf{lh}_{\kappa,\kappa}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}=\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}$
$\displaystyle\mathsf{rh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}\circ\mathsf{lh}_{\kappa,\nu}^{H}$
$\displaystyle=\mathsf{lh}_{(\kappa\oplus\mu)\oplus\kappa,\nu}^{H}\circ\mathsf{rh}_{\kappa\oplus\mu,\kappa}^{H}=\mathsf{lh}_{\kappa\oplus(\kappa\oplus\mu),\nu}^{H}\circ\mathsf{lh}_{\kappa,\kappa\oplus\mu}^{H}$
$\displaystyle=\mathsf{lh}_{(\kappa\oplus\kappa)\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa\oplus\kappa,\mu}^{H}\circ\mathsf{lh}_{\kappa,\kappa}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa\oplus\mu,\nu}^{H}\circ\mathsf{lh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}$
$\displaystyle\mathsf{lh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}\circ\mathsf{rh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{(\kappa\oplus\mu)\oplus\kappa,\nu}^{H}\circ\mathsf{lh}_{\kappa\oplus\mu,\kappa}^{H}\circ\mathsf{rh}_{\kappa,\mu}^{H}$
$\displaystyle=\mathsf{lh}_{\kappa\oplus(\kappa\oplus\mu),\nu}^{H}\circ\mathsf{rh}_{\kappa,\kappa\oplus\mu}^{H}\circ\mathsf{rh}_{\kappa,\mu}^{H}=\mathsf{lh}_{(\kappa\oplus\kappa)\oplus\mu,\nu}^{H}\circ\mathsf{rh}_{\kappa\oplus\kappa,\mu}^{H}$
$\displaystyle=\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}\circ\mathsf{lh}_{\mu,\nu}^{H}$
$\displaystyle\mathsf{rh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}\circ\mathsf{rh}_{\kappa,\nu}^{H}$
$\displaystyle=\mathsf{rh}_{(\kappa\oplus\mu)\oplus\kappa,\nu}^{H}=\mathsf{rh}_{\kappa\oplus\mu,\nu}^{H}$
$\displaystyle=\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}\circ\mathsf{rh}_{\mu,\nu}^{H}$

which imply the following

$\displaystyle(\langle\kappa,r\rangle$
$\displaystyle\cdot_{H}\langle\mu,s\rangle)+_{H}(\langle\kappa,r\rangle\cdot_{H}\langle\nu,t\rangle)=$
$\displaystyle=\langle(\kappa\oplus\mu)\oplus(\kappa\oplus\nu),\mathsf{lh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}(\mathsf{lh}_{\kappa,\mu}^{H}(r)+\mathsf{rh}_{\kappa,\mu}^{H}(s))+\mathsf{rh}_{\kappa\oplus\mu,\kappa\oplus\nu}^{H}(\mathsf{lh}_{\kappa,\nu}^{H}(r)+\mathsf{rh}_{\kappa,\nu}^{H}(t))\rangle$
$\displaystyle=\langle\kappa\oplus(\mu\oplus\nu),(\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}(r)\cdot\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}(\mathsf{lh}_{\mu,\nu}^{H}(s)))+(\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}(r)\cdot\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}(\mathsf{rh}_{\mu,\nu}^{H}(t)))\rangle$
$\displaystyle=\langle\kappa\oplus(\mu\oplus\nu),\mathsf{lh}_{\kappa,\mu\oplus\nu}^{H}(r)\cdot\mathsf{rh}_{\kappa,\mu\oplus\nu}^{H}(\mathsf{lh}_{\mu,\nu}^{H}(s)+\mathsf{rh}_{\mu,\nu}^{H}(t))\rangle$
$\displaystyle=\langle\kappa,r\rangle\cdot_{H}(\langle\mu,s\rangle+_{H}\langle\nu,t\rangle)$

which proves distributivity when all the elements are different from
$\mathbf{0}_{H}$.

Now, suppose that $\langle\kappa,r\rangle=\mathbf{0}_{H}$, then we have
$\langle\kappa,r\rangle\cdot_{H}(\langle\mu,s\rangle+_{H}\langle\nu,t\rangle)=\langle\kappa,r\rangle\cdot_{H}\langle\mu,s\rangle=\langle\kappa,r\rangle\cdot_{H}\langle\nu,t\rangle=\mathbf{0}_{H}$,
hence distributivity trivialy holds. Finally, suppose
$\langle\mu,s\rangle=\mathbf{0}_{H}$ (the case
$\langle\nu,t\rangle=\mathbf{0}_{H}$ is similar), then we have
$\langle\mu,s\rangle+_{H}\langle\nu,t\rangle=\langle\nu,t\rangle$ and
$\langle\kappa,r\rangle\cdot_{H}\langle\mu,s\rangle=\mathbf{0}_{H}$.
Therefore, we get

$\displaystyle\langle\kappa,r\rangle\cdot_{H}(\langle\mu,s\rangle+_{H}\langle\nu,t\rangle)$
$\displaystyle=\langle\kappa,r\rangle\cdot_{H}\langle\nu,t\rangle=\mathbf{0}_{H}+_{H}\langle\kappa,r\rangle\cdot_{H}\langle\nu,t\rangle$
$\displaystyle=\langle\kappa,r\rangle\cdot_{H}\langle\mu,s\rangle+_{H}\langle\kappa,r\rangle\cdot_{H}\langle\nu,t\rangle$


# Dynamics of van der Waals Charge Qubit in 2D Bilayers: Ab initio Quantum
Transport and Qubit Measurement

Jiang Cao  Integrated Systems Laboratory, ETH Zürich, 8092 Zürich, Switzerland
Guido Gandus  Integrated Systems Laboratory, ETH Zürich, 8092 Zürich,
Switzerland  Empa, Swiss Federal Laboratories for Materials Science and
Technology, Überlandstrasse 129, 8600 Dübendorf, Switzerland  Tarun Agarwal
Department of Electrical Engineering, IIT Gandhinagar, Gandhinagar, India
Mathieu Luisier  Integrated Systems Laboratory, ETH Zürich, 8092 Zürich,
Switzerland  Youseung Lee  Corresponding author<EMAIL_ADDRESS>Integrated Systems Laboratory, ETH Zürich, 8092 Zürich, Switzerland

###### Abstract

A van der Waals (vdW) charge qubit, electrostatically confined within two-
dimensional (2D) vdW materials, is proposed as building block of future
quantum computers. Its characteristics are systematically evaluated with
respect to its two-level anti-crossing energy difference ($\Delta$). Bilayer
graphene ($\Delta$ $\approx$ 0) and a vdW heterostructure ($\Delta$ $\gg$ 0)
are used as representative examples. Their tunable electronic properties with
an external electric field define the state of the charge qubit. By combining
density functional theory and quantum transport calculations, we highlight the
optimal qubit operation conditions based on charge stability and energy-level
diagrams. Moreover, a single-electron transistor (SET) design based on
trilayer vdW heterostructures capacitively coupled to the charge qubit is
introduced as measurement setup with low decoherence and improved measurement
properties. It is found that a $\Delta$ greater than 20 meV results in a rapid
mixing of the qubit states, which leads to a lower measurement quantity, i.e.
contrast and conductance. With properly optimized designs, qubit architectures
relying on 2D vdW structures could be integrated into an all-electronic
quantum computing platform.

quantum computing, charge qubit, two-dimensional materials, bilayer graphene,
quantum transport

##  I Introduction

Emerging nanoelectronic devices are widely seen as key enablers of future
quantum computers that will be able to solve intractable problems for
classical machines. A wide range of qubit platforms have already been
experimentally demonstrated, e.g. charge, spin, superconducting, and trapped-
ion qubits [1, 2, 3, 4]. In all these approaches, a large number of physical
qubits are required to realize one single logical qubit [5]. The ideal qubit
structure should therefore be scalable and integrable with the existing CMOS
technology to build a large-scale quantum computer.

Quantum dots (QD) offer a great platform to generate CMOS-compatible qubits
based on the spin and/or electronic orbital degrees of freedom of electrons
[6, 7, 8]. For example, silicon spin qubits show very long coherence times up
to milliseconds at low temperature [2]. However, they are rather challenging
to handle due to the rapid flip-flop of single electron spins. Charge qubits,
on the other hand, can simplify the read-out of qubits thanks to well-
developed technologies such as single-electron transistors (SET). Moreover,
these qubits can be directly operated by electrodes up to gigahertz rates.

Figure 1:  (a) Sketch of a charge qubit system implemented on a 2D bilayer
material to which a vertical electric field $\vec{E}$ is applied. The wave
functions localized in the bottom and top layers are labeled $\ket{B}$ and
$\ket{T}$, respectively. (b) Schematic of the proposed qubit measurement setup
where three QDs labeled QDL,C,R are confined inside a vdW trilayer
MoS2-WSe2-MoS2 heterostructure controlled by a top (TG) and bottom gate (BG)
and connected to a source and drain contact forming altogether a single-
electron transistor (SET). The trilayer MoS2-WSe2-MoS2 SET is electrostaticlly
coupled to the charge qubit located on the BLG. The distance between the SET
and qubit is $L$.

Two-dimensional (2D) materials represent a promising alternative to
conventional semiconductors to host qubits owing to the wide tunability of
their electronic properties and the manyfold quantum phenomena that emerge in
the atomically thin limit such as massless Dirac fermions [9], anyons [10], or
the quantum spin Hall effect [11]. In addition, the van der Waals (vdW)
stacking of 2D materials provides versatility in growth, functionalization,
and heterostructure integration, which facilitates the processing of complex
systems. AB-stacked bilayer graphene (BLG), a vdW homojunction, is an
excellent example: its band gap can be gradually modified by applying an
external electric field [12, 13, 14, 15, 16, 17, 18]. Moreover, the electric
field can modify the quantum superposition of wave functions localized in
different layers, affecting both the conduction and valence bands of BLG. The
same effect can be obtained in other vdW structures as well if the band edges
of the individual 2D materials that constitute the stack get nearly aligned on
an absolute energy scale [19]. Based on this observation, vdW structures are
attractive to implement charge qubits capable of switching between two states
with the help of a vertical electric field [20]. However, for vdW charge
qubits to become a viable option, an all-electronic platform that can host
multiple QDs should be available. The coupling between neighboring QDs should
be tunable, while it should be possible to measure the qubit state with low
decoherence.

In this work, we propose a platform for vdW charge qubits that takes advantage
of quantum superposition of the top and bottom charge states formed by vdW
bilayers, as illustrated in Fig. 1(a), and we analyze its potential through ab
initio device simulations. Differently from the conventional charge qubits
that require two QDs, e.g. the left and right QDs, a single QD of bilayer acts
as an individual qubit. As a promising candidate, we showcase the
functionality of a BLG vdW charge qubit system. The role of two-level anti-
crossing energies $\Delta$ is analyzed for two exemplary vdW materials that
illustrate the extreme cases of $\Delta$: BLG ($\Delta\approx 0$) and
ZrS2-HfS2 ($\Delta=0.174$ eV). As a first step, density functional theory
(DFT) and quantum transport calculations are combined to study the charge
stability diagram and energy-level spectrum of multiple vdW charge qubits with
tunable coupling, formed on an electrostatically confined BLG device. It
should be pointed out that directly coupled read-out leads, as originally
proposed in the paper of Lucatto et al. [20], can be prominent sources of
decoherence due to the finite tunnel coupling of the reservoirs to the charge
qubit. As a consequence, it might be experimentally difficult to measure the
qubit state. To overcome this issue, we designed a novel measurement setup
using triple-QD SETs based on a trilayer vdW heterostructure electrostatically
coupled to our vdW charge qubits, as depicted in Fig. 1(b). Here, its
influence will be simulated through the Lindblad master equation with the
QuTiP package [21] using a phenomenological model including the important
physics of decoherence. Our results highlight that 2D vdW structures can
provide highly integrable all-electronic universal qubits. Moreover, our
findings indicate that devices with $\Delta\gg$ 20 meV pose specific design
challenges related to the material combination and geometrical parameters.
Hence, a careful selection of the $\Delta$ value is necessary to design
properly working vdW charge qubits.

##  II vdW charge qubit manipulation

We first examine the effect of a vertical electric field on the selected 2D
vdW bilayer materials, BLG and ZrS2-HfS2, and show that the electric-field
externally applied to the bilayers can be used to control the wave function
localized in the top or bottom layer to construct the qubit states. The BLG
band structure without and with electric field is shown in Figs. 2(a) and
2(b), respectively. The detailed DFT simulation procedure that was used to
compute those quantities can be found in Appendix A. The conduction band
minimum (CBM) and valence band maximum (VBM) are chosen as the two-level qubit
states ($\ket{0}$ and $\ket{1}$) in the case of BLG, while the first and
second lowest conduction bands of the ZrS2-HfS2 heterostructure play the same
role. One importance for vdW bilayers to properly act as qubits is the
existence of two bands lying quite close to each other and that can be easily
inverted by applying a practical electric field.

In the BLG system without electric field $\vec{E}$, the $\ket{0}$ and
$\ket{1}$ states are made of exactly equal contributions from each graphene
layer due to the inversion symmetry of this material. A bandgap can be opened
by applying a symmetry-breaking vertical $\vec{E}$. The vdW charge qubit
states, $\ket{0}$ and $\ket{1}$, are constructed by combining the wave
functions localized in the top and bottom layers, $\ket{T}$ and $\ket{B}$,
respectively.

Figure 2: Electronic bandstructure of BLG (a) without and (b) with an electric
field of 0.5 V/nm. The color of the different states indicates the
$\absolutevalue{b_{0,1}}^{2}$ weights, which result from the projection of the
atomic wave functions onto the bottom graphene layer. The valence band maximum
and conduction band minimum define the $\ket{0}$ and $\ket{1}$ states. (c)
Energy-level vs. electric field for the $\ket{0}$ and $\ket{1}$ states in BLG.
(d) Corresponding $\absolutevalue{b_{0,1}}^{2}$ and
$\absolutevalue{t_{0,1}}^{2}$ weights vs. electric field. (e) Same as (c) but
for ZrS2-HfS2. (f) Same as (d) but for ZrS2-HfS2.

By applying a positive $\vec{E}$ along the $+z$ direction, as in Fig. 1(a),
the bands of the bottom graphene layer are lowered in energy, which increases
the $\ket{B}$ character of $\ket{0}$, while $\ket{1}$ becomes more localized
in the top graphene layer. The opposite behavior happens with $\vec{E}$
pointing along the $-z$ direction. The contribution of each layer to the Bloch
wave function of $\ket{0}$ and $\ket{1}$ strongly depends on the vertical
electric field. An electron wave function $\ket{\psi(\vec{E})}$ can be written
as the superposition of the wave functions $\ket{T}$ and $\ket{B}$ for a given
$\vec{E}$

$\ket{\psi(\vec{E})}=t_{\psi}(\vec{E})\ket{T}+b_{\psi}(\vec{E})\ket{B},$ (1)

where $t_{\psi}(\vec{E})$ and $b_{\psi}(\vec{E})$ are complex coefficients
depending on $\vec{E}$. Due to the vdW gap between the two layers, the
$\ket{T}$ and $\ket{B}$ states have a relatively small overlap and form a
quasi-orthogonal basis. The $t_{\psi}(\vec{E})$ and $b_{\psi}(\vec{E})$
coefficients satisfy the normalization condition
$\absolutevalue{t_{\psi}}^{2}+\absolutevalue{b_{\psi}}^{2}=1$.

Figure 2(c) shows the energy level of the $\ket{0}$ and $\ket{1}$ states and
the $\ket{B}$ contribution to them as a function of $\vec{E}$. In Fig. 2(d)
the $\absolutevalue{t_{0,1}}^{2}$ and $\absolutevalue{b_{0,1}}^{2}$ weights of
the $\ket{0}$ and $\ket{1}$ states are plotted, highlighting a complementary
relation between them. Here, we note that by tuning $\vec{E}$ we can readily
address any coherent superposition of the basic quantum states $\ket{B}$ and
$\ket{T}$. The state $\ket{\psi}$ can be equivalently described by a Bloch
vector in the standard Bloch sphere representation, where the coefficients are
given by a polar $\theta$ and azimuthal angle $\phi$ such that
$t_{\psi}=\cos(\theta/2)$ and $b_{\psi}=\sin(\theta/2)e^{i\phi}$. Moreover,
under a reasonably strong positive/negative electric field, the $\ket{0}$ and
$\ket{1}$ states can be localized in either the top or bottom layer, thus
enabling the “must-have” qubit initialization process.

Figure 3:  (a) Simulation domain corresponding to a bilayer graphene (BLG)
with a total length of 260 nm along the transport direction $x$ and a width of
2.1 nm along $y$ for a total of 42400 atoms, with the periodic boundary
condition along $y$. A vertical electric field $\vec{E}$ acts on the BLG and
is applied through a back gate (BG) contact. The BLG is embedded between two
$h$BN layers. The red regions under L and R indicate the position of the
electrostatically confined quantum dots. The electrostatic potential is
further controlled by three top gates labeled L, M, and R. They are separated
from the BLG by a SiO2 oxide and a $h$BN layer. The bottom inset shows a
typical potential profile obtained by superposing Gaussian functions induced
by each top electrode. (b) Local density-of-states of the electron and hole
populations. Positive (negative) values refer to the top (bottom) layer,
respectively. $\vec{E}=0.5$ V/nm, $V_{\rm M}=\pm 4$ mV, and $V_{\rm L}=V_{\rm
R}=\pm 0.01$ V where + stands for hole while - for electron. The dashed lines
represent the conduction and valence band edges. (c) Schematics of the band
profile along the BLG sheet illustrating the formation of dots: (I) using the
M gate, two n-doped islands form the desired electron dots (ED) with a
tunneling barrier between them; (II) at $V_{\rm M}=0$, one large ED is formed;
(III) at larger $V_{\rm L}$ and $V_{\rm R}$, two additional hole dots (HD) are
created under these gates. (d) Charge stability diagram: logarithmic
conductance at a drain-to-source voltage $V_{\rm DS}$=0.1 mV, Fermi level
$E_{F}$=-2.0525 eV, $\vec{E}=0.5$ V/nm, $V_{\rm M}$=0 V and $k_{B}T$=1 $\mu$eV
as a function of the L and R gate voltages. The equilibrium charge number in
HDL, ED and HDR are depicted as (N0,N1,N2).

A similar behavior is observed for the first and second lowest conduction band
states of ZrS2-HfS2, as depicted in Fig. 2(e) and (f). Due to the energy
offset in the natural band alignment of the vdW heterostructure, a non-zero
$\vec{E}_{ac}$ is required to reach the crossing point where the layer
contributions are equal. At this point, a so-called anti-crossing energy
$\Delta$ can be introduced. The latter is then a material parameter and
determined as the difference between the energies of the $\ket{1}$ and
$\ket{0}$ eigenstates [1]. Table I summaries the values of $\Delta$ for the
considered vdW structures. For BLG, we note that $\Delta\gtrsim 20$ $\mu$eV
was observed in previous experiments [14, 22, 23] due to the presence of
symmetry-breaking factors even in the absence of an electric field, for
example spin-orbit coupling, surrounding substrate atoms, or local strain.

By comparing our results to the ZrSe2-SnSe2 vdW heterostructure investigated
in the original paper of Lucatto et al. [20], we notice that the Bloch vector
of the BLG and ZrS2-HfS2 systems covers a larger $\theta$ interval for
electric fields within $\pm$3 V/nm. For BLG, the full range of $\theta$ can be
achieved with a small electric field of less than 1 V/nm, thus making BLG the
most technologically appealing candidate.

A quantum control technique over the angle $\phi$ (i.e. relative phase of the
wave functions) applicable to the vdW charge qubits has been proposed in Ref.
24. It relies on the Landau-Zener-Stückelberg (LZS) interference when the
system non-adiabatically sweeps through the intersection point. Based on this
technique, general rotations in the Bloch sphere and universal single qubit
operations have already been experimentally demonstrated for Si and GaAs DQD
charge qubits [1, 24].

##  III Quantum transport simulations

To create a clear two-level system with a $\ket{0}$ and $\ket{1}$ state, the
QDs in the vdW bilayers must be formed electrically through gate electrodes
that induce a local confinement, as illustrated in Fig. 3(a). In the proposed
device setup, a bottom gate uniformly acts on the BLG encapsulated within two
hBN dielectric layers to electrostatically control the Fermi level. A
transport channel can be defined through the BLG. Separated by an oxide, three
top gates labeled L, M, and R are placed. By applying different potentials to
the bottom and top gates, the vertical electric field acting on the BLG can be
precisely controlled. In addition, the potential landscape can be modified by
tuning the top L, M, and R gates. The simulated device structure measures 260
nm along the transport direction, possesses zigzag edges, and is made of 42400
atoms in total. The three top gates are 10 nm long and separated by 50 nm.
Charge transport is solved with a quantum transport simulator [25] based on
the ab initio Non-equilibrium Green’s function (NEGF) formalism, as described
in Appendix A. In this study, the potential profiles are approximated as
Gaussian functions that are plotted in the bottom inset of Fig. 3(a),
Poisson’s equation being not self-consistently solved due to the large
computational burden associated with the consideration of thousands of gate
voltage configurations.

Figure 3(b) shows that by applying $\vec{E}=0.5$ V/nm and by properly
adjusting the applied top gate voltages both electron and hole QDs can be
generated in BLG, without spurious bandgap states and with a strong
localization of the $\ket{0}$ and $\ket{1}$ states in the top and bottom
graphene layer. The projected local density-of-states (LDOS) confirm that the
qubit state defined in Fig. 2 is not altered by the confinement potential. By
applying an appropriate bias to the top M gate (or plunger gate), a gradual
transition from two QDs hosting electrons or holes to a single larger QD can
be achieved, thus providing a way to tune the interdot tunnel coupling. This
transition occurs between regime I and II in Fig. 3(c). The charge stability
is demonstrated by the distinct resonances in the conductance map as a
function of the gate voltages applied to the L and R electrodes, $V_{\rm
L,R}$, setting $V_{\rm M}=0$. $V_{\rm L,R}$ creates two potential barriers
that confine an electron dot (ED) in the middle and that tune its energy
levels. Each time an energy level crosses the energy window defined by the
Fermi energies of the source and drain contacts, a resonance peak can be seen
in the conductance map (see Fig. 3(d)). In addition to creating n and p-type
QDs in the central region, the electrostatic potential can also be shaped to
form ambipolar QDs through the application of large enough L and R gate
voltages, as highlighted by regime III in Fig. 3(c). The L/R gates closely
control the energy levels in the hole dot (HDL/HDR) underneath it. Whenever
the energy levels of the ED and HDL/HDR cross the energy window of interest, a
crossing of two lines of maximum conductance can be observed in Fig. 3(d). We
note that the observed high LDOS on the QD boundaries, which originated from
the coupling to the continuum of states in the left and right contacts due to
the NEGF open boundary condition, is expected to be reduced in a fully self-
consistent Schrödinger-Poisson solution. However, we do not expect this to
alter qualitatively our result because the conductance is limited by the
electron states inside the quantum dot rather than in the boundaries.

Figure 4: Local density-of-states of the electron and hole population of BLG
at a positive electric field of 0.5 V/nm. Positive (negative) values refer to
the top (bottom) layer. (a) At gate voltages $V_{M}$=0 V, $V_{R}$=-0.06 V, and
$V_{L}$=-0.01 V. (b) Same as (a), but at $V_{L}$=-0.06 V. The dashed lines
represent the conduction and valence band edges. The corresponding charge
density distributions in the top and bottom layers with $E_{F}$=-2.048 eV are
plotted in the upper panels.

In order to investigate the system scalability towards a highly integrable
qubit platform, we also show the posibility of creating multiple qubits within
one BLG device in Fig. 4. Since the electron and hole dots possess well-
separated LDOS and charge densities localized in the top and bottom layers,
they can be exploited for hosting $\ket{1}$ and $\ket{0}$ qubit states. It is
worth noting that between the QDs, the Fermi energy lies inside the band gap,
which leads to a finite region of zero charge density and a tunneling barrier.
Figure 4(a) shows an n-p double-QD, while Fig. 4(b) presents a p-n-p triple-QD
formed in the channel. A chain of QDs of $\ket{1}$ and $\ket{0}$ can be
constructed via $N$ top gates, thus promising easily scalable multi-qubit
systems [8]. Contrary to the traditional charge qubits where two QDs are
required to form a single charge qubit [1, 26], a single vdW-bilayer QD is
sufficient to form one qubit because the wave function in each dot is a
quantum superposition of $\ket{T}$ and $\ket{B}$. Moreover, a single global
bottom gate electrode reduces the number of gate electrodes, which greatly
simplifies the geometrical complexity and increases the integration density of
qubits. Coupling remote charge qubits has been experimentally demonstrated for
graphene [27, 28] using microwave resonators, offering the possibility to
entangle more-than-nearest-neighbor QDs. While the previous experiments show a
viable path to implement a multi-qubit coupling scheme for the vdW-bilayer
QDs, this remains a challenging task that will need detailed investigations in
the future. They go far beyond the scope of our present work.

##  IV Time-dependent qubit measurement

After demonstrating the qubit operation conditions and obtaining charge
stability for multiple QDs, the procedure to read out the qubit state is now
discussed. For time-dependent simulations, the qubit is modeled by a two-level
Hamiltonian constructed from the DFT bandstructure of 2D bilayers under
different vertical electric fields

$\begin{split}\hat{H}_{\rm
Q}(\vec{E})=&\epsilon_{T}(\vec{E})\ket{T}\bra{T}+\epsilon_{B}(\vec{E})\ket{B}\bra{B}\\\
&-\frac{\Delta}{2}(\ket{T}\bra{B}+\ket{B}\bra{T}),\end{split}$ (2)

where $\epsilon_{T}(\vec{E})$ and $\epsilon_{B}(\vec{E})$ are the on-site
energies of the top and bottom layer that depend linearly on $\vec{E}$, and
$\Delta$ is the anti-crossing energy. These values are obtained by fitting the
VBM and CBM energy levels for varying $\vec{E}$ (see Fig. 2) to the
eigenvalues of the quadratic Hamiltonian.

Figure 5: (a) Electronic bandstructure of the trilayer vdW heterostructure
MoS2-WSe2-MoS2 for the SET measuring setup. The color of the states indicates
the contribution from each individual MoS2 and WSe2 monolayer to the overall
bandstructure. The inset represents the simulated vdW structure. (b)
Conceptual representation of the whole system consisting of one vdW qubit (red
circled) and a SET (blue circled) viewed from the y-z cross-section in Fig.
1(b). The relevant energy levels of QDL,C,R are $E_{L,C,R}$, respectively.
They are shown in the top sub-plot. The tunneling probabilities between the
source, QDL,C,R and the drain are $\Gamma_{L}$, $t_{L,R}$, and $\Gamma_{R}$,
respectively. The SET is coupled through Coulomb interaction to the charge
qubit located within a bilayer system at a distance $L$.

A SET formed by a linear array of three QDs capacitively coupled to the qubit
is chosen for the measurement, as suggested in Ref. [29] and shown in Fig.
5(b). Such a configuration has an enhanced readout sensitivity as compared to
the more common single QD SET. The designed SET is described by the following
Hamiltonian operator,

$\begin{split}\hat{H}_{\rm
SET}=&E_{L}\ket{L}\bra{L}+E_{R}\ket{R}\bra{R}+E_{C}\ket{C}\bra{C}\\\
&-t_{L}(\ket{C}\bra{L}+\ket{L}\bra{C})\\\
&-t_{R}(\ket{C}\bra{R}+\ket{R}\bra{C})\\\
&+\sum_{i,j=L,C,R}W_{i,j}\ket{i,j}\bra{i,j},\end{split}$ (3)

where $\ket{L}$, $\ket{R}$, and $\ket{C}$ denote the single-electron states
when the electron is located in QDL, QDR, and QDC, respectively, and $E_{\rm
L}$, $E_{\rm R}$, and $E_{\rm C}$ for the on-site energies; $t_{L}$ and
$t_{R}$ indicate the tunneling probabilities between the QDs; $\ket{L,R}$,
$\ket{L,C}$ and $\ket{C,R}$ are the two-electron states, while $W_{i,j}$
represents the Coulomb interaction between dots $i$ and $j$. The electron
density in the QD is assumed as a point charge. In this case the energies of
the electrostatic interaction have the following values $W_{i,j}=2{\rm
Ry}\times a_{B}/d_{i,j}$, where $d_{i,j}$ is the distance between the centers
of the QDs $i$ and $j$, $a_{B}=0.52\times 10^{-10}$ m is the Bohr radius, Ry =
13.6 eV is the Rydberg energy. This approximation works well for small-sized
charge dots as the ones considered in this work [29]. The SET is only coupled
to the qubit $via$ the Coulomb interactions without any direct tunneling,

$\hat{H}_{\rm Q-SET}=\sum_{i=T,B;j=L,C,R}W_{i,j}\ket{i,j}\bra{i,j}.$ (4)

The time evolution of the density matrix operator $\hat{\rho}$ corresponding
to the qubit-SET system was simulated by solving the Lindblad equation [29,
30]. More details are given in Appendix B.

The optimal condition for the aforementioned triple-QD SET requires first that
the energetic detunings between the external (QDL, QDR) and the central (QDC)
dots, $\delta_{\rm L,R}=E_{\rm L,R}-E_{\rm C}$, are equivalent. Second, the
electrostatic couplings $W$ between the SET QDs and the qubit should be
symmetric (asymmetric) for the electron in the top (bottom) layer,
respectively. More explicitly, $W_{T,L}=W_{T,R}$ and $W_{B,L}\neq W_{B,R}$. We
therefore propose a vdW heterostructure made of three 2D materials
(MoS2-WSe2-MoS2) as SET design, as illustrated in Fig. 1(b).

The projected bandstructure in Fig. 5(a) reveals that the CBM and VBM are
located within one of the 2D layers. The degenerate CBM of the top and bottom
MoS2 layers define the energies of the external QDs L and R, and the VBM of
WSe2 realizes the central QD C. The direct bandgap at K in the Brillouin Zone
corresponds to the symmetric detuning $\delta_{\rm L,R}$. The vertical
alignment of the top qubit layer with WSe2 guarantees that the presence of an
electron in this layer causes a symmetric renormalization of the energies in
the SET. In contrast, an electron in the qubit bottom layer creates an
imbalance of detunings and causes the current to cease. The triple-QD in the
vdW heterostructure is formed electrically by a top gate electrode that
induces a local confinement. To asses the impact of electron-phonon and
electron-electron interactions on the measurement sensitivity of the coupled
qubit and SET systems, we included the relaxation and dephasing rates
$\gamma_{rel}$ and $\gamma_{dep}$, accounting for those scattering processes,
respectively. The corresponding collapse operator in the Lindblad equation is
constructed as

$\begin{split}\mathcal{C}=&\gamma_{rel}[\mathcal{D}(\ket{C}\bra{L})+\mathcal{D}(\ket{C}\bra{R})],\\\
+&\gamma_{dep}[\mathcal{D}(\ket{L}\bra{L}-\ket{C}\bra{C})+\mathcal{D}(\ket{R}\bra{R}-\ket{C}\bra{C})\\\
+&\mathcal{D}(\ket{T}\bra{T}-\ket{B}\bra{B})],\end{split}$ (5)

where
$\mathcal{D}(\hat{O})\equiv\hat{O}\hat{\rho}\hat{O}^{\dagger}-\frac{1}{2}[\hat{O}^{\dagger}\hat{O},\hat{\rho}]$.

Figure 6: (a) Contrast $C$ of the designed SET as a function of the anti-
crossing energy $\Delta$ and the distance $L$ between the qubit and SET. (b)
Same as (a), but for the differential conductance $g_{m}$ in unit of
($e^{2}/\hbar$). (c) Contrast $C$ and (d) differential conductance $g_{m}$ as
a function of $\Delta$ and of the dephasing rate $\gamma_{dep}$ with $L$=12.5
nm. $\Delta$ relates to the thermodynamic temperature $T=\Delta/2k_{B}$
limiting qubit operation.

The qubit is initialized by the application of an external electric field into
one of the two states $\ket{0}$ or $\ket{1}$. The system then evolves until it
reaches a steady state at time $\sim 1/\Gamma_{R}$ where $\Gamma_{R}$ is the
escaping rate of the electrons to the drain (see Fig. 5(b)). The output
current $I=2\pi e\Gamma_{R}\rho_{R}/\hbar$ is proportional to the population
in the SET $\ket{R}$ state $\rho_{R}$. In what follows, we set $\Gamma_{R}$ =
$\Gamma_{L}$ = 0.12 THz in agreement with experimental studies [1]. To assess
the measurement capabilities of the proposed SET we record its contrast $C$
and conductance $g_{m}$ during the qubit dynamics. Here, the contrast
determines how well both qubit states can be distinguished. It is defined as

$C=\absolutevalue{I(1)-I(0)}/\absolutevalue{I(1)+I(0)},$ (6)

where $I(0)$ and $I(1)$ denote the output current when measuring a qubit in
the $\ket{0}$ and $\ket{1}$ states, respectively. The differential conductance
quantifying the SET current changes caused by variations of the qubit states
is defined as [29]

$g_{m}=Ce[I(1)+I(0)]/\absolutevalue{W_{T,R}-W_{B,R}}.$ (7)

The resulting contrast and differential conductance are plotted in Fig. 6 as a
function of the qubit anti-crossing energy $\Delta$ and distance $L$ between
the qubit and SET with a characteristic time of $20/\Gamma_{R}$, including the
aforementioned dissipative processes ($\gamma_{rel}=0.05$ THz based on low-
temperature transport measurement [31], while $\gamma_{dep}=2\gamma_{rel}$ to
account for the fact that dephasing happens faster than relaxation). The
bottom left corner of Fig. 6(a) shows the limiting range of $L$ (up to 13 nm)
and $\Delta$ (up to 25 meV) to obtain high contrasts. Figure 6(b) also reports
that $g_{m}$ decreases as a function of $\Delta$, but in a limited range of
$L$ values between $12$ and $25$ nm. Combining both sub-plots, it clearly
appears that the optimal conditions for qubit measurement are $L$=12.5 nm and
$\Delta$ between 0 and 25 meV. The maximum $g_{m}$ value is higher than the
experimental value reported in Ref. [32] using a different SET setup. A widely
accepted idea is that a large $\Delta$ is desired between the qubit states, in
order to make them energetically distinguishable from each other. A necessary,
but insufficient condition is that $\Delta$ must be larger than two times the
thermal energy ($2k_{B}T$) to suppress thermalization effects [33]. Therefore,
a smaller $\Delta$ value necessarily entails a lower operation temperature of
the qubit, as indicated by the upper abscissa axis in Fig 6. Here, we show
that a normally functional qubit also imposes an upper-limit on $\Delta$, thus
reducing the qubit design space.

Figure 7: (a) Time evolution of the drain current $I_{D}(0)$ and $I_{D}(1)$
flowing through a SET electrostatically coupled to a vdW qubit system with two
logical states ($\ket{0}$ and $\ket{1}$) and of the top layer ($\rho_{T}(0)$
and $\rho_{T}(1)$) and bottom layer ($\rho_{B}(0)$ and $\rho_{B}(1)$)
population for the BLG system. (b) Same as (a), but for the ZrS2-HfS2 system.
The inset shows the short-time evolution after the beginning of measurement.

To examine the intrinsic feasibility of vdW charge qubits, we study $C$ and
$g_{m}$ as a function of the dephasing rate and anti-crossing energy in Fig.
6(c) and (d). The distance between the SET and qubit is fixed to $12.5$ nm to
ensure that both measurement quantities have high resolutions. It can be seen
that the dephasing rate has a small impact on the degradation of $C$ and
$g_{m}$ as long as $\Delta$ remains below 20 meV. However, when $\Delta$
exceeds 20 meV, $g_{m}$ and $C$ start to significantly decrease. This
parameter region should therefore be avoided for the qubit design. This means
that the anti-crossing energy is an intrinsically limiting factor, which
suggests that the charge qubit should be defined with $\Delta\leq 20$ meV and
the operating temperature lower than 120 K.

Table 1:  Summary of the material properties and qubit measurement results ($C$ and $g_{m}$) for the BLG and ZrS2-HfS2 systems with $\vec{E}=-3~{}V/nm$, $\gamma_{dep}=2\gamma_{rel}=0.1$ THz. The $\Delta$ is the anti-crossing energy, whereas $\vec{E}_{ac}$ denotes the electric field at the anti-crossing point when $t_{0,1}=b_{0,1}$. | BLG | ZrS2 (T) HfS2 (B)
---|---|---
Contrast [%] | 71 | 26
$g_{m}$ [$\mu$S] | 102 | 35.7
$\Delta$ [meV] | 0.02 | 174
$\vec{E}_{ac}$ [V/nm] | 0 | 1.6

A typical time evolution of the drain current in this SET and the
corresponding charge population in the qubit $\ket{T}$ ($\rho_{T}$) and
$\ket{B}$ ($\rho_{B}$) states are presented in Fig. 7 for BLG and for
ZrS2-HfS2. The BLG exhibits a much higher contrast of saturation current and
almost no mixing. For ZrS2-HfS2, the strong $\Delta$ results in a fast qubit
state mixing and consequently requires very rapid measurements within the
picosecond time scale. By reducing the strength of $\Delta$, for example by
inserting an $h$BN monolayer between the vdW bilayers, the qubit measurement
quality can be recovered. These results underline the fact that high contrast
qubit measurements are intrinsically not possible if the charge qubit is
characterized by a large $\Delta$. Table I summaries the calculated qubit
measurement results.

##  V Conclusion and Discussion

In this work, we proposed a vdW charge qubit architecture based on 2D vdW
materials. Our results indicate that the qubit states can be readily
controlled by external electric fields. The two-level anti-crossing energy,
which is related to the vdW interlayer coupling strength, can significantly
vary between different material systems. Our quantum transport simulations
suggest that quantum dots can be confined by electrically gating the vdW
bilayer and that the coupling between neighboring quantum dots can be tuned by
the gate potential. Furthermore, we also proposed a single-electron transistor
(SET) design based on a trilayer vdW heterostructure and capacitively coupled
to the quantum dot to measure the charge qubit state. Our time evolution
simulations of the coupled qubit and SET systems reveal that the measurement
sensitivity and contrast are affected by the intrinsic anti-crossing energy
and by the dephasing processes in the vdW heterostructure. Our study
identifies the optimal design parameter space to simultaneously achieve a low
decoherence of the qubit states and enhanced measurement contrast and
conductance.

In conclusion, our theory explains that upon the SET measurement the vdW-
heterostructure qubit states rapidly degrade due to the large coupling
intrinsically induced by anti-crossing energies higher than 20 meV. This
phenomenon was not captured by previous theoretical studies. Amoung the
systems we studied, bilayer graphene is the one maintaining the highest
contrast between qubit states. Our findings show that 2D vdW structures can
provide highly integrable all-electronic universal qubits. Our design of the
coupled vdW qubit and SET system allows for an enhanced sensitivity of the
charge qubit, thus paving the way for a scalable and CMOS-compatible quantum
computing platform.

##  Appendix A Ab initio quantum transport simulation

Following Ref. [20], a charge qubit can be defined in 2D bilayer materials
through the superposition of a top ($\ket{T}$) and bottom ($\ket{B}$) orbital
state, whose properties depend on a vertically applied electric field
($\vec{E}$). Since each local orbital occupies a different electrically-
induced energy state, two distinct energy levels with different weights
($t(\vec{E})$ and $b(\vec{E})$) can be obtained (Fig. 1(a)). They are used as
building blocks for the Bloch sphere and its $\ket{0}$ and $\ket{1}$ states,
which form the two-level quantum system shown in Fig. 1(b). The latter is
simulated at the ab initio level.

Through the QUANTUM ESPRESSO DFT package [34], the electronic structure of 2D
bilayer materials is first obtained. The generalized gradient approximation of
Perdew, Burke, and Ernzerhof (PBE) [35] is used as exchange and correlation
functional, while the vdW interactions are taken into account with the DFT-D2
scheme of Grimme [36]. The plane-wave DFT Hamiltonian is then converted into a
maximally-localised Wannier function (MLWF) basis with the Wannier90 code
[37]. These results are scaled up [38] to construct the real-space device
Hamiltonian and to perform realistic qubit device simulations based on the
Non-equilibrium Green’s Function (NEGF) formalism. The following equation is
solved for the retarded Green’s function $\mathbf{G^{R}}$,

$\left[E\mathbf{I}-\mathbf{H}_{\rm MLWF}-\mathbf{\Sigma}_{\rm S/D}^{\rm
R}\right]\mathbf{G}^{\rm R}=\mathbf{I},$ (8)

where $E$ is the electron energy, $\mathbf{I}$ the identity matrix,
$\mathbf{H}_{\rm MLWF}$ the Wannier Hamiltonian matrix of the 2D bilayer
graphene in the MLWF basis, and $\mathbf{\Sigma}_{\rm S/D}^{\rm R}$ the
retarded boundary self-energy accounting for the source and drain semi-
infinite contacts. The observables of the system can be obtained from the
retarded Green’s function $\mathbf{G}^{\rm R}$. The periodic boundary
condition is imposed along the $y$ direction as shown in Fig. 3(a). By
including an electrostatically-induced confinement potential into
$\mathbf{H}_{\rm MLWF}$, a wide range of QD systems can be realized and their
operation verified by inspecting the local density-of-states (LDOS). The
electrical current of the device is calculated with the Landauer-Büttiker
formula from the transmission probability $T(E)$ between two contacts [39, 40]

$I=-\frac{e}{\hbar}\int\frac{dE}{2\pi}T(E)[f(E,E_{fL})-f(E,E_{fR})],$ (9)

with the elementary charge $e$. The two contacts are characterized by their
Fermi distribution function $f$ and Fermi levels $E_{fL}$ and $E_{fR}$, whose
difference depends on the applied source-to-drain voltage
$E_{fL}-E_{fR}=eV_{DS}$. The conductance $G$ is calculated as $dI/dV_{DS}$.

##  Appendix B Evaluation of time-dependent Lindblad equation for qubit
readout

The dynamics of the charge qubits is simulated with the QuTiP [21] package,
which solves the Lindblad equation [29],

$\begin{split}\frac{d\hat{\rho}}{dt}=&-i\hbar\left[\hat{H}_{\rm
Q}+\hat{H}_{\rm SET}+\hat{H}_{\rm Q-SET},\hat{\rho}\right]\\\
&+\mathcal{C}+\Gamma_{R}\ket{vac_{R}}\bra{R}+\Gamma_{L}\ket{L}\bra{vac_{L}},\\\
\end{split}$ (10)

where the qubit Hamiltonian ($H_{\rm Q}$) is a two-level quantum system
constructed from the electronic structure of a 2D bilayer material, while a
three-level QD is used to represent the SET Hamiltonian ($H_{\rm SET}$). The
qubit and SET parts are only electrostatically coupled through the electron-
electron Coulomb interaction $H_{\rm Q-SET}$, assuming no tunneling between
them. The dissipative processes including relaxation and dephasing are
described by the collapse operators $\mathcal{C}$ in Eq. (5). The meaning and
definition of each operator are given in Section IV. An arbitrary state as in
Eq. (1) is evolved by integrating the set of ordinary differential equations
that define Eq. (10). The expectation value of the populations in the SET are
recorded at each time step by means of the operators
$\rho_{L}=\ket{L}\bra{L}$, $\rho_{R}=\ket{R}\bra{R}$, and
$\rho_{C}=\ket{C}\bra{C}$. The output at the last time step is used in Eqs.
(6) and (7) to asses the contrast and conductance of the proposed SET.

###### Acknowledgements.

This work was supported by the NCCR MARVEL and the NCCR SPIN of the Swiss
National Science Foundation (SNSF), by SNSF under Grant No. 175479 (ABIME), by
the Marie Skłodowska-Curie Grant No. 885893, and by a grant from the Swiss
National Supercomputing Centre (CSCS) under Project s1119.

## References

  * Gorman  _et al._ [2005] J. Gorman, D. G. Hasko, and D. A. Williams, Phys. Rev. Lett. 95, 090502 (2005).
  * Maurand  _et al._ [2016] R. Maurand, X. Jehl, D. Kotekar-Patil, A. Corna, H. Bohuslavskyi, R. Laviéville, L. Hutin, S. Barraud, M. Vinet, M. Sanquer, and S. De Franceschi, Nat. Commun. 7, 13575 (2016).
  * Gambetta  _et al._ [2017] J. M. Gambetta, J. M. Chow, and M. Steffen, Npj Quantum Inf. 3, 2 (2017).
  * Bruzewicz  _et al._ [2019] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, Appl. Phys. Rev. 6, 021314 (2019).
  * Van Meter and Horsman [2013] R. Van Meter and C. Horsman, Commun. ACM 56, 84–93 (2013).
  * Eich  _et al._ [2018] M. Eich, R. Pisoni, A. Pally, H. Overweg, A. Kurzmann, Y. Lee, P. Rickhaus, K. Watanabe, T. Taniguchi, K. Ensslin, and T. Ihn, Nano Letters 18, 5042 (2018), publisher: American Chemical Society.
  * Kurzmann  _et al._ [2019] A. Kurzmann, H. Overweg, M. Eich, A. Pally, P. Rickhaus, R. Pisoni, Y. Lee, K. Watanabe, T. Taniguchi, T. Ihn, and K. Ensslin, Nano Letters 19, 5216 (2019), publisher: American Chemical Society.
  * Banszerus  _et al._ [2020] L. Banszerus, S. Möller, E. Icking, K. Watanabe, T. Taniguchi, C. Volk, and C. Stampfer, Nano Lett. 20, 2005 (2020).
  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Xiaofeng  _et al._ [2014] Q. Xiaofeng, L. Junwei, F. Liang, and L. Ju, Science 346, 1344 (2014).
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006).
  * McCann  _et al._ [2007] E. McCann, D. S. Abergel, and V. I. Fal’ko, Solid State Communications 143, 110 (2007), exploring graphene.
  * Aoki and Amawashi [2007] M. Aoki and H. Amawashi, Solid State Communications 142, 123 (2007).
  * Castro  _et al._ [2007] E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. C. Neto, Phys. Rev. Lett. 99, 216802 (2007).
  * Oostinga  _et al._ [2008] J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo, and L. M. K. Vandersypen, Nature Materials 7, 151 (2008).
  * Mak  _et al._ [2009] K. F. Mak, C. H. Lui, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 102, 256405 (2009).
  * Kumar and Guo [2011] S. B. Kumar and J. Guo, Appl. Phys. Lett. 98, 222101 (2011).
  * Koda  _et al._ [2018] D. S. Koda, F. Bechstedt, M. Marques, and L. K. Teles, Phys. Rev. B 97, 165402 (2018).
  * Lucatto  _et al._ [2019] B. Lucatto, D. S. Koda, F. Bechstedt, M. Marques, and L. K. Teles, Phys. Rev. B 100, 121406 (2019).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Comput. Phys. Commun. 183, 1760 (2012).
  * Konschuh  _et al._ [2012] S. Konschuh, M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. B 85, 115423 (2012).
  * Banszerus  _et al._ [2021] L. Banszerus, S. Möller, C. Steiner, E. Icking, S. Trellenkamp, F. Lentz, K. Watanabe, T. Taniguchi, C. Volk, and C. Stampfer, Nature Communications 12, 5250 (2021).
  * Cao  _et al._ [2013] G. Cao, H.-O. Li, T. Tu, L. Wang, C. Zhou, M. Xiao, G.-C. Guo, H.-W. Jiang, and G.-P. Guo, Nat. Commun. 4, 1401 (2013).
  * Luisier [2014] M. Luisier, Chem. Soc. Rev. 43, 4357 (2014).
  * Kim  _et al._ [2015] D. Kim, D. R. Ward, C. B. Simmons, J. K. Gamble, R. Blume-Kohout, E. Nielsen, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, and M. A. Eriksson, Nat. Nanotechnol. 10, 243 (2015).
  * Deng  _et al._ [2015a] G.-W. Deng, D. Wei, S.-X. Li, J. R. Johansson, W.-C. Kong, H.-O. Li, G. Cao, M. Xiao, G.-C. Guo, F. Nori, H.-W. Jiang, and G.-P. Guo, Nano Lett. 15, 6620 (2015a).
  * Deng  _et al._ [2015b] G.-W. Deng, D. Wei, J. R. Johansson, M.-L. Zhang, S.-X. Li, H.-O. Li, G. Cao, M. Xiao, T. Tu, G.-C. Guo, H.-W. Jiang, F. Nori, and G.-P. Guo, Phys. Rev. Lett. 115, 126804 (2015b).
  * Tsukanov [2019] A. V. Tsukanov, Phys. Rev. A 100, 062305 (2019).
  * Manzano [2020] D. Manzano, AIP Advances 10, 025106 (2020).
  * Engels  _et al._ [2014] S. Engels, B. Terrés, A. Epping, T. Khodkov, K. Watanabe, T. Taniguchi, B. Beschoten, and C. Stampfer, Phys. Rev. Lett. 113, 126801 (2014).
  * Barthel  _et al._ [2010] C. Barthel, M. Kjærgaard, J. Medford, M. Stopa, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Phys. Rev. B 81, 161308 (2010).
  * Makhlin  _et al._ [2001] Y. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).
  * Giannozzi  _et al._ [2017] P. Giannozzi _et al._ , J. Phys.: Condens. Matter 29, 465901 (2017).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme [2006] S. Grimme, Journal of computational chemistry 27, 1787 (2006).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Szabó  _et al._ [2019] A. Szabó, A. Jain, M. Parzefall, L. Novotny, and M. Luisier, Nano Letters 19, 3641 (2019), pMID: 31079463.
  * Landauer [1988] R. Landauer, IBM Journal of Research and Development 32, 306 (1988).
  * Büttiker [1990] M. Büttiker, Phys. Rev. Lett. 65, 2901 (1990).