# Teaching Design by Contract using Snap! ††thanks: Identify applicable
funding agency here. If none, delete this.

1st Marieke Huisman  Formal Methods and Tools
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>2nd Raúl E. Monti  Formal Methods and Tools
University of Twente
Enschede, The Netherlands
<EMAIL_ADDRESS>
###### Abstract

With the progress in deductive program verification research, new tools and
techniques have become available to support design-by-contract reasoning about
non-trivial programs written in widely-used programming languages. However,
deductive program verification remains an activity for experts, with ample
experience in programming, specification and verification. We would like to
change this situation, by developing program verification techniques that are
available to a larger audience. In this paper, we present how we developed
prototypal program verification support for Snap!. Snap! is a visual
programming language, aiming in particular at high school students. We added
specification language constructs in a similar visual style, designed to make
the intended semantics clear from the look and feel of the specification
constructs. We provide support both for static and dynamic verification of
Snap! programs. Special attention is given to the error messaging, to make
this as intuitive as possible.

###### Index Terms:

verification, software, education

##  I Introduction

Research in deductive program verification has made substantial progress over
the last years: tools and technique have been developed to reason about non-
trivial programs written in widely-used programming languages, the level of
automation has substantially increased, and bugs in widely-used libraries have
been found [1, 2, 3]. However, the use of deductive verification techniques
remains the field of expert users, and substantial programming knowledge is
necessary to appreciate the benefits of these techniques.

We believe that it is important to make deductive program verification
techniques accessible also to novice programmers. Therefore, we have to teach
the Design-by-Contract [4] (DbC) approach, which requires the programmer to
explicitly specify the assumptions and responsibilities of code in a modular
way, in parallel with actually teaching programming, i.e. DbC should be taught
as an integral part of the process leading from design to implementation. In
this paper, we make the Design-by-Contract idea accessible to high school
students, in combination with appropriate tool support, which is currently
unavailable.

Concretely, this paper presents a Design-by-Contract approach for Snap! [5].
Snap! is a visual programming language targeting high school students. The
design of Snap! is inspired by Scratch, another widely-used visual programming
language. Compared to Scratch, Snap! has some more advanced programming
features. In particular, Snap! provides the possibility to create parametrised
reusable blocks, basically modelling user-defined functions. Also the look and
feel of Snap! targets high school students, whereas Scratch aims at an even
younger age group. Snap! has been successfully integrated in high school
curricula, by its integration in the _Beauty and Joy of Computing_ course [6].
This course combines programming skills with a training in abstract
computational thinking.

The first step to support Design-by-Contract for Snap! is to define a suitable
specification language. The visual specification language that we propose in
this paper is built as a seamless extension of Snap!, i.e. we propose a number
of new specification blocks and natural modifications of existing ones. These
variations capture the main ingredients for the Design-by-Contract approach,
such as pre- and postconditions. Moreover, we also provide blocks to add
assertions and loop invariants in a program and we extend the standard
expression pallets of Snap! with some common expressions to ease
specifications. The choice of specification constructs is inspired by existing
specification languages for Design-by-Contract, such as JML [7], choosing the
most frequently used constructs with a clear and intuitive meaning. Moreover,
all verification blocks are carefully designed to reflect the intended
semantics of the specifications in a visual way.

A main concern for a programmer, after writing the specification of the
intended behaviour of their programs, should be to validate that these
programs behave according to their specification. Therefore, we provide two
kinds of tool support: (i) runtime assertion checking [8], which checks
whether specifications are not violated during a particular program execution,
and (ii) static checking (or deductive verification) [9], which verifies that
all possible program executions respect its specifications. The runtime
assertion checker is built as an extension of the standard Snap! execution
mechanism. The deductive verification support is built by providing a
translation from a Snap! program into Boogie [10].

Another important aspect to take into account for a good learning experience
are the error messages that indicate that a specification is violated. We have
integrated these messages in Snap!’s standard error reporting system, again
sticking to the look and feel of standard Snap!. Moreover, we have put in
effort to make the error messages as clear as possible, so that also a
relatively novice programmer can understand why the implementation deviates
from the specification.

##  II Background

###  II-A Snap!

Snap! is a visual programming language. It has been designed to introduce
children (but also adults) to programming in an intuitive way. At the same
time, it is also a platform for serious study of computer science[11]. Snap!
actually re-implements and extends Scratch [12]. Programming in Snap! is done
by dragging and dropping blocks into the coding area. Blocks represent common
program constructs such as variable declarations, control flow statements
(branching and loops), function calls and assignments. Snapping blocks
together, the user builds a script and visualises its behaviour by means of
turtle graphics visualisations, called sprites. Sprites can change shape,
move, show bubbled text, play music, etc. For all these effects, dedicated
blocks are available.

Figure 1: The Snap! working area.

The Snap! interface divides the working area into three parts: the pallet
area, the scripting area, and the stage area, indicated by labels 1, 2 and 3,
respectively, in Fig. 1. On the left, the various programming blocks are
organised into pallets that describe their natural use. For instance, the
_Variables_ pallet contains blocks for declaring and manipulating variables.
In Snap!, variables are dynamically typed. Blocks are dragged and dropped from
the pallets into the scripting area, located at the centre of the working area
where the Snap! program is constructed. Blocks can be arranged by snapping
them together, or by inserting them as arguments of other blocks. Blocks can
only be used as arguments if their shapes match with the shape of the argument
slots in the target block. These shapes actually provide a hint on the
expected evaluation type of a block, for instance, rounded slots for numbers
and diamond slots for booleans  .

The behaviour of the script is shown with turtle graphics drawings in the
stage area located in the rightmost part of the screen.

In addition, at the bottom of the pallet area, there is a “Make a block”
button. This allows the user to define his or her _Build Your Own Block_
(BYOB) blocks. When pressed, a new floating “Block Editor” window pops out
with a new coding area, in which the behaviour of the personalised block can
be defined (similar to how a script is made in the scripting area). Label $4$
in Fig. 1 shows a BYOB block being edited. Once defined, the BYOB block
becomes available to be used just as any other predefined block.

###  II-B Program Verification

The basis of the Design-by-Contract approach [13] is that the behaviour of all
program components is defined as a contract. For example, a function contract
specifies the conditions under which a function may be called (the function’s
_precondition_), and it specifies the guarantees that the function provides to
its caller (the function’s _postcondition_). There exist several specification
languages that have their roots in this Design-by-Contract approach. For
example the Eiffel programming language has built-in support for pre- and
postconditions [14], and for Java, the behavioural interface language JML [15]
is widely used. As is common for such languages, we use the keyword _requires_
to indicate a precondition, and the keyword _ensures_ to indicate a
postcondition.

If a program behaviour is specified using contracts, various techniques can be
used to validate whether an implementation respects the contract.

Dynamic verification validates an implementation w.r.t. a specification at
runtime. This means that, whenever during program execution a specification is
reached, it will be checked for this particular execution that the property
specified indeed holds. In particular, this means that whenever a function
will be called, its precondition will be checked, and whenever the function
returns, its postcondition will be checked. An advantage of this approach is
that it is easy and fast to use it: one just runs a program and checks if the
execution does not violate the specifications. A disadvantage is that it only
provides guarantees about a concrete execution.

In contrast, static verification aims at verifying that all possible
behaviours of a function respect its contract. This is done by applying Hoare
logic proof rules [16] or using Dijkstra’s predicate transformer semantics
[17]. Applying these rules results in a set of first-order proof obligations;
if these proof obligations can be proven it means that the code satisfies its
specification. Advantage of this approach is that it guarantees correctness of
all possible behaviours. Disadvantage is that it is often labour-intensive,
and often many additional annotations, such as for example loop invariants,
are needed to guide the prover.

##  III Visual Program Specifications

This section discusses how to add visual specification constructs to Snap!.
Our goal was to do this in such a way that (1) the intended semantics of the
specification construct is clear from the way it is visualised, and (2) that
it smoothly integrates with the existing programming constructs in Snap!

Often, Design-by-Contract specifications are added as special comments in the
code. For example, in JML a function contract is written in a special comment,
tagged with an @-symbol, immediately preceding the function declaration. The
tag ensures that the comment can be recognised as part of the specification.
There also exist languages where for example pre- and postconditions are part
of the language (e.g., Eiffel [18], Spec# [19]). We felt that for our goal,
specifications should be integrated in a natural way in the language, rather
than using comments. Therefore, we introduce variations of the existing block
structures, in which we added suitable slots for the specifications. This
section discusses how we added pre- and postconditions, and in-code
specifications such as asserts and loop invariants to Snap!. In addition, to
have a sufficiently expressive property specification language, we also
propose an extension of the expression constructs.

###  III-A Visual Pre- and Postconditions

To specify pre- and postconditions for a BYOB script, we provide a variation
of the initial hat block with a slot for a precondition at the start of the
block, and a slot for a postcondition at the end of the block (Fig. 2).

This shape is inspired by the c-shaped style of other Snap! blocks, such as
blocks for loops. The main advantage is that it visualises at which points in
the execution, the pre- and the postconditions are expected to hold. In
addition, it also graphically identifies which code is actually verified.
Moreover, the shapes are already familiar to the Snap! programmer. If the
slots are not filled, default pre- and postcondition true can be used. Notice
that the pre- and postcondition slots consist of multiple boolean-argument
slots, and we define the property to be the conjunction of the evaluation of
each of these slots. This is similar to how Snap! extends a list or adds
arguments to the header of a BYOB.

Figure 2: Hat block extended with contracts

###  III-B Visual Assertions and Loop Invariants

For static verification, pre- and postconditions are often not sufficient, and
we need additional in-code specifications to guide the prover, such as
assertions, which specify properties that should hold at a particular point in
the program, and loop invariants. Moreover, assertions can also be convenient
for run-time assertion checking to make it explicit that a property holds at a
particular point in the program.

#### Visual Assertions

To specify assertions, both the property specified and the location within the
code are relevant. To allow the specification of assertions at arbitrary
places in a script, we define a special assertion block  similar to all other
control blocks.

#### Visual Loop Invariants

Loop invariants are necessary for static verification [20]. A loop invariant
should hold at the beginning and end of every loop iteration. To account for
this, we provide a (multi-argument boolean) slot to specify the loop invariant
in the traditional Snap! c-shaped loop block. This slot is located just after
the header where the loop conditions are defined. In addition, the c-shaped
loop block repeats the word invariant at the bottom of the block (see Figure
3) to visually indicate that the invariant is checked after each iteration.

Figure 3: Visual loop invariants.

###  III-C Visual Expressions

In addition, we have introduced some specification-only keywords, as commonly
found in Design-by-Contract languages.

  * •

An _old_ expression is used in postconditions to indicate that a
variable/expression should be evaluated in the pre-state of the function. To
support this, we introduced an operator block  with a slot for a variable
name.

  * •

A _result_ expression refers to the return value of a function inside its
postcondition. We support this by introducing a constant  operator, that
allows to specify a property about the result value of a reporter BYOB.

We also introduce syntax to ease the definition of complex Boolean
expressions, by means of the operator blocks  ,  ,  and  , as well as syntax
to write more advanced Boolean expressions, introducing support for quantified
expressions (See Fig.4).

Figure 4: A global quantification expression block

##  IV Graphical approach to verification result reporting

Another important point to consider is how to report on the outcome of the
verification: (1) presenting the verdict of a passed verification, and (2) in
case of failure, giving a concrete and understandable explanation for the
failure. The latter is especially important in our case, as we are using the
technique with inexperienced users.

In order to signal a contract violation, or any assertion invalidated during
dynamic verification, we use Snap!’s pop-up notification windows. These
windows have the advantage that a failing block can be printed inside them
even when the failing script is not currently visible to the user. This allows
to be very precise about the error, even when the BYOB body is not currently
visible.

In order to signal errors while compiling to Boogie, such as making use of
dynamic typing or nested lists in your Snap! BYOB code, we use Snap!’s speech
bubbles that can emerge at specific points in the script while describing the
cause of failure. This has the advantage that the failing block can easily be
singled out by the location of the bubble, while the cause of failure is
described by the text inside the bubble. We find this option less invasive
than a pop-up window but still as precise, and we can be sure that the blocks
involved will be visible since static verification is triggered from the BYOB
editor window (See Fig.5). Notice that currently we do not report the results
of static verification within Snap!, since our extension only returns a
compiled Boogie code which has to be verified with Boogie separately.

Figure 5: Static verification compilation notification.

##  V Tool support

We have developed our ideas into a prototypal extension to Snap! which can be
found at https://git.snt.utwente.nl/montire/verifiedsnap/. This repository
also contains a set of running examples to showcase the new support for
verification. These are available in the _lessons_ folder under the root
directory along with an exercise sheet named _exercises.pdf_. The extension
uses the same technology as the original Snap! and can be run by just opening
the _snap.html_ file in most common web-browsers that support java-script.

Our extension supports both dynamic and static verification of BYOB blocks.
Dynamic verification is automatically triggered when executing BYOB blocks in
the usual way. For static verification, a dedicated button located at the top
right corner of the BYOB editor window allows to trigger the compilation of
the BYOB code into an intended equivalent Boogie code. The compiled code can
be then downloaded and verified with Boogie. Boogie can be run locally or on
the cloud at https://rise4fun.com/Boogie/. _Dynamic verification_ has been
fully integrated into the normal execution flow of a Snap! program, and thus
there is no real restrictions on the characteristics of the BYOB that can be
dynamically verified. For _Static verification_ , we have restricted data
types to be Integers, Booleans and List of Integers. Moreover, we do not
support dynamic typing of variables. Finally, we only focus on compiling an
interesting subset of Snap! blocks for the sake of teaching Design-by-
Contract.

##  VI Conclusions

This paper presented a prototypal program verification extension to Snap!. The
extension is intended to support the teaching of Design-by-Contract in the
later years of high schools. We paid considerable attention to the didactic
aspects of our tool: the looks and feel of the extension should remain
familiar to Snap! users, the syntax and structure of the new blocks should
give a clear intuition about their semantics, and the error reporting should
be precise and expressive.

Our extension allows to analyse BYOB blocks both by runtime assertion checking
and static verification. Runtime assertion checking is fully integrated into
Snap! and there is no limitation on the kind of blocks that can be analysed.
Static verification compiles the Snap! code into a Boogie equivalent code and
the verification needs to be run outside of Snap!. Moreover, we make some
restrictions on the kind of BYOB blocks we can compile, in order to keep the
complexity of the prototype low. As future work we would like to lift these
restrictions as much as possible by integrating the remaining Snap! blocks
into the compilation and by allowing other data types to be used. Also, we
would like to integrate the verification into Snap!, translating Boogie
messages back to the Snap! world, to help student to interpret them.

We would like to carry out an empirical study on our proposed approach. This
will require the development of a concrete study plan and its evaluation in a
Dutch classroom.

Computer science curricula that uses blocks programming is widely and freely
available [21, 22, 23, 24, 25]. Nevertheless, it is hardly spotted that they
include topics around design and verification of code. The words ‘test’ or
‘testing’ are also rare around the curricula and, where mentioned, they are
not sufficiently motivated. The drawbacks of teaching coding with blocks
without paying attention to design nor correctness has already been analysed
[26, 27]. We have not found any work on teaching these concepts in schools,
nor implementations on block programming that support teaching them.

|  Marieke Huisman is a professor in Software Reliability at the University of
Twente. She obtained her PhD in 2001 from the Radboud University Nijmegen.
Afterwards she worked at INRIA Sophia Antipolis, and since 2008 at University
of Twente. Her research interests are in the verification of concurrent
software, as implemented in the VerCors program verifier. She is in particular
interested in making verification usable in a practical setting, and she works
for example on annotation generation, and support for different programming
languages.
---|---
|  Raúl E. Monti received his PhD in 2018 from Universidad Nacional de
Córdoba. He is currently a PostDoc at the Universiteit Twente. His research
interests involve the development and practical application of formal
foundations and tools for analysis and verification of software and hardware
systems by means of model checking and deductive verification. His work
involves interacting with industry to apply his research in the verification
of industrial (embedded) systems and software.
---|---

## References

  * [1] S. De Gouw, J. Rot, F. De Boer, R. Bubel, and R. Hähnle, “OpenJDK’s java.utils.collection.sort() is broken: The good, the bad and the worst case,” in _Proc. 27th Intl. Conf. on Computer Aided Verification (CAV), San Francisco_ , ser. LNCS, D. Kroening and C. Pasareanu, Eds., vol. 9206\. Springer, Jul. 2015, pp. 273–289.
  * [2] W. Oortwijn, M. Huisman, S. J. Joosten, and J. van de Pol, “Automated verification of parallel nested dfs,” in _International Conference on Tools and Algorithms for the Construction and Analysis of Systems_. Springer, 2020, pp. 247–265.
  * [3] M. Safari, W. Oortwijn, S. Joosten, and M. Huisman, “Formal verification of parallel prefix sum,” in _NASA Formal Methods_ , R. Lee, S. Jha, and A. Mavridou, Eds. Cham: Springer International Publishing, 2020, pp. 170–186. [Online]. Available: https://doi.org/10.1007/978-3-030-55754-6_10
  * [4] B. Meyer, “Applying "design by contract",” _Computer_ , vol. 25, no. 10, pp. 40–51, 1992. [Online]. Available: https://doi.org/10.1109/2.161279
  * [5] B. Harvey, D. D. Garcia, T. Barnes, N. Titterton, D. Armendariz, L. Segars, E. Lemon, S. Morris, and J. Paley, “Snap!(build your own blocks),” in _Proceeding of the 44th ACM technical symposium on Computer science education_ , 2013, pp. 759–759.
  * [6] D. Garcia, B. Harvey, and T. Barnes, “The beauty and joy of computing,” _Inroads_ , vol. 6, no. 4, pp. 71–79, 2015. [Online]. Available: https://doi.org/10.1145/2835184
  * [7] G. T. Leavens, A. L. Baker, and C. Ruby, “JML: A notation for detailed design,” in _Behavioral Specifications of Businesses and Systems_ , H. Kilov, B. Rumpe, and I. Simmonds, Eds. Boston, MA: Springer US, 1999, pp. 175–188.
  * [8] Y. Cheon, “A runtime assertion checker for the Java Modeling Language,” Ph.D. dissertation, Department of Computer Science, Iowa State University, Ames, 2003, technical Report 03-09.
  * [9] K. R. M. Leino, “Towards reliable modular programs,” 1995.
  * [10] M. Barnett, B. E. Chang, R. DeLine, B. Jacobs, and K. R. M. Leino, “Boogie: A modular reusable verifier for object-oriented programs,” in _Formal Methods for Components and Objects, 4th International Symposium, FMCO 2005, Amsterdam, The Netherlands, November 1-4, 2005, Revised Lectures_ , ser. Lecture Notes in Computer Science, F. S. de Boer, M. M. Bonsangue, S. Graf, and W. P. de Roever, Eds., vol. 4111. Springer, 2005, pp. 364–387. [Online]. Available: https://doi.org/10.1007/11804192_17
  * [11] B. Harvey and J. Mönig, “Snap! reference manual,” _URL http://snap. berkeley. edu/SnapManual. pdf_ , 2017.
  * [12] M. Resnick, J. Maloney, A. Monroy-Hernández, N. Rusk, E. Eastmond, K. Brennan, A. Millner, E. Rosenbaum, J. S. Silver, B. Silverman, and Y. B. Kafai, “Scratch: programming for all,” _Commun. ACM_ , vol. 52, no. 11, pp. 60–67, 2009. [Online]. Available: https://doi.org/10.1145/1592761.1592779
  * [13] B. Meyer, J.-M. Nerson, and M. Matsuo, “Eiffel: object-oriented design for software engineering,” in _European Software Engineering Conference_. Springer, 1987, pp. 221–229.
  * [14] B. Meyer, “Eiffel: A language and environment for software engineering,” _Journal of Systems and Software_ , vol. 8, no. 3, pp. 199–246, 1988.
  * [15] G. T. Leavens, Y. Cheon, C. Clifton, C. Ruby, and D. R. Cok, “How the design of jml accommodates both runtime assertion checking and formal verification,” _Science of Computer Programming_ , vol. 55, no. 1-3, pp. 185–208, 2005.
  * [16] C. Hoare, “An axiomatic basis for computer programming,” _Communications of the ACM_ , vol. 12, no. 10, pp. 576–580, 1969.
  * [17] E. Dijkstra, _A Discipline of Programming_. Prentice-Hall, 1976.
  * [18] B. Meyer, _Eiffel: The Language_. Prentice-Hall, 1991. [Online]. Available: http://www.eiffel.com/doc/#etl
  * [19] M. Barnett, K. R. M. Leino, and W. Schulte, “The Spec# programming system: An overview,” in _Construction and Analysis of Safe, Secure and Interoperable Smart Devices: Proceedings of the International Workshop CASSIS 2004_ , ser. LNCS, G. Barthe, L. Burdy, M. Huisman, J.-L. Lanet, and T. Muntean, Eds., vol. 3362. Springer, 2005, pp. 151–171.
  * [20] T. Türk, “Local reasoning about while-loops,” in _VSTTE 2010. Workshop Proceedings_ , R. Joshi, T. Margaria, P. Müller, D. Naumann, and H. Yang, Eds. ETH Zürich, 2010, pp. 29–39.
  * [21] “The Beauty and Joy of Computing. An AP CS Principles Course,” https://bjc.edc.org/. Accessed October 2020.
  * [22] “CS First,” https://csfirst.withgoogle.com/s/en/home. Accessed October 2020.
  * [23] “The Creative Computing Curriculum,” http://creativecomputing.gse.harvard.edu/guide/. Accessed October 2020.
  * [24] “An Introduction to Programming. A Pencil Code Teacher’s Manual,” https://manual.pencilcode.net/. Accessed October 2020.
  * [25] P. Factorovich and F. Sawady, _Actividades para aprender a Program. AR: Segundo ciclo de la educación primaria y primero de la secundaria_. Miller Ed Buenos Aires, 2015.
  * [26] O. Meerbaum-Salant, M. Armoni, and M. Ben-Ari, “Habits of programming in scratch,” in _Proceedings of the 16th Annual SIGCSE Conference on Innovation and Technology in Computer Science Education, ITiCSE 2011, Darmstadt, Germany, June 27-29, 2011_ , G. Rößling, T. L. Naps, and C. Spannagel, Eds. ACM, 2011, pp. 168–172. [Online]. Available: https://doi.org/10.1145/1999747.1999796
  * [27] E. Aivaloglou and F. Hermans, “How kids code and how we know: An exploratory study on the scratch repository,” in _Proceedings of the 2016 ACM Conference on International Computing Education Research_ , 2016, pp. 53–61.


# A low-loss ferrite circulator as a tunable chiral quantum system

Ying-Ying Wang  Department of Physics, University of Massachusetts-Amherst,
Amherst, MA, USA  Sean van Geldern  Department of Physics, University of
Massachusetts-Amherst, Amherst, MA, USA  Thomas Connolly  Present address:
Department of Applied Physics, Yale University, New Haven, CT, USA Department
of Physics, University of Massachusetts-Amherst, Amherst, MA, USA  Yu-Xin Wang
Pritzker School of Molecular Engineering, University of Chicago, Chicago, IL,
USA  Alexander Shilcusky  Department of Physics, University of Massachusetts-
Amherst, Amherst, MA, USA  Alexander McDonald  Pritzker School of Molecular
Engineering, University of Chicago, Chicago, IL, USA  Department of Physics,
University of Chicago, Chicago, IL, USA  Aashish A. Clerk  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, IL, USA  Chen Wang
<EMAIL_ADDRESS>Department of Physics, University of Massachusetts-Amherst,
Amherst, MA, USA

###### Abstract

Ferrite microwave circulators allow one to control the directional flow of
microwave signals and noise, and thus play a crucial role in present-day
superconducting quantum technology. They are typically viewed as a black-box,
with their internal structure neither specified nor used as a quantum
resource. In this work, we show a low-loss waveguide circulator constructed
with single-crystalline yttrium iron garnet (YIG) in a 3D cavity, and analyze
it as a multi-mode hybrid quantum system with coupled photonic and magnonic
excitations. We show the coherent coupling of its chiral internal modes with
integrated superconducting niobium cavities, and how this enables tunable non-
reciprocal interactions between the intra-cavity photons. We also probe
experimentally the effective non-Hermitian dynamics of this system and its
effective non-reciprocal eigenmodes. The device platform provides a test bed
for implementing non-reciprocal interactions in open-system circuit QED.

##  I Introduction

Microwave circulators, typically composed of a transmission line Y-junction
with ferrite materials Kord  _et al._ (2020), are ubiquitous in
superconducting circuit QED experiments Devoret and Schoelkopf (2013). They
provide a crucial link in the readout chain of superconducting quantum
processors, by directing the signal traffic while protecting the qubits and
resonators from thermal noise Krantz  _et al._ (2019). They also enable the
interactions between distinct quantum circuit modules to be non-reciprocal
Kurpiers  _et al._ (2018); Axline  _et al._ (2018), a feature which is
important for eliminating long-distance cross-talks in modular quantum
computation architectures. Despite their importance, microwave circulators are
generally treated as broadband black-box devices in experiments. Formulating a
more microscopic quantum description is often challenging, as their internal
modes involving the magnetic spin excitations (magnons) are generally too
lossy and complex to be analyzed using canonical circuit quantization Vool and
Devoret (2017).

On the other hand, there has been growing interest in studying and
manipulating magnon excitations of ferromagnetic/ferrimagnetic materials in
the quantum regime Lachance-Quirion  _et al._ (2019); Awschalom  _et al._
(2021). In particular, the ferromagnetic resonance (FMR) mode of yttrium iron
garnet (YIG), a ferrimagnetic insulator with usage in commercial circulators,
has shown sufficiently high quality factor and coupling cooperativity with
microwave cavities to function as a quantum oscillator mode in strong-coupling
circuit QED Huebl  _et al._ (2013); Tabuchi  _et al._ (2014); Zhang  _et al._
(2014). Notably, coherent coupling of magnons with a superconducting qubit
Tabuchi  _et al._ (2015) and single-shot detection of a single magnon
Lachance-Quirion  _et al._ (2020) have been demonstrated using a millimeter-
sized single-crystalline YIG sphere in a 3D cavity. Furthermore, there is a
plausible pathway towards planar superconducting-magnonic devices Hou and Liu
(2019); Golovchanskiy  _et al._ (2021) to connect circuit QED with spintronics
technologies by advancing fabrication techniques of low-damping YIG films
Heyroth  _et al._ (2019).

It would be interesting to harness these recent advances in the study of
quantum magnonics to revisit the design of microwave circulators, potentially
leading to new kinds of non-reciprocal devices in circuit QED. Our work here
describes a first step in this direction. Here we demonstrate a tunable non-
reciprocal device based on the waveguide circulator loaded with single-
crystalline YIG, which explicitly makes use of well-characterized hybrid
polariton modes. Such modes are the normal modes of coupled magnon-photon
systems Huebl  _et al._ (2013); Tabuchi  _et al._ (2014); Zhang  _et al._
(2014); Boventer  _et al._ (2018); Zhang  _et al._ (2017), and have an
intrinsic chirality that is set by the magnetic field Anderson  _et al._
(2016); Owens  _et al._ (2018); Zhang  _et al._ (2020). While our device
follows the same basic working principles underpinning textbook circulators
Kord  _et al._ (2020); Fay and Comstock (1965), detailed understanding of the
internal modes allows us to incorporate the physical source of non-reciprocity
in the full description of a larger system including two external
superconducting cavities, using a non-Hermitian effective Hamiltonian.

While our device can be configured to operate as a traditional circulator for
its non-reciprocal transmission of travelling waves, the main focus of our
study is to use the device for mediating tunable non-reciprocal interaction
between localized long-lived quantum modes. Such non-reciprocal mode-mode
couplings result in distinct signatures in the eigenvalues and eigenvectors of
the non-Hermitian system Hamiltonian, which is relevant to the more general
study of non-Hermitian dynamics in contexts ranging from classical optics to
quantum condensed matter. Anomalous properties of the eigenvalues and
eigenvectors of a non-Hermitian Hamiltonian have given rise to a number of
striking phenomena such as the existence of exceptional points Heiss (2012);
Özdemir  _et al._ (2019) and the non-Hermitian skin effect Hatano and Nelson
(1997); Yao and Wang (2018); McDonald  _et al._ (2018), but direct
experimental access to the underlying eigenmodes is often difficult. In this
study, we provide comprehensive characterization of the eigenmode structure,
which is a step towards effective Hamiltonian engineering of non-reciprocal
non-Hermitian systems.

The most tantalizing usage of non-reciprocity in quantum systems (such as
entanglement stabilization using directional interactions in chiral quantum
optics setups Stannigel  _et al._ (2012); Lodahl  _et al._ (2017)) require
extremely high quality devices. In particular, they must approach the pristine
limit where undesirable internal loss rates are negligible compared to the
non-reciprocal coupling rates. While many experiments have focused on new
avenues of achieving non-reciprocity Chapman  _et al._ (2017); Lecocq  _et
al._ (2017); Sliwa  _et al._ (2015); Ruesink  _et al._ (2016); Fang  _et al._
(2017); Wang  _et al._ (2019); Xu  _et al._ (2019), this loss-to-coupling
ratio, which can be understood as the quantum efficiency of the non-reciprocal
interactions, has been typically limited to approximately 10% ($\sim$ 0.5 dB)
or more, which is comparable to the linear insertion loss of typical
commercial circulators as measured in modular circuit QED experiments Kurpiers
_et al._ (2018); Axline  _et al._ (2018). This performance lags far behind the
quality of unitary operations between reciprocally coupled quantum components
(i.e. two-qubit gate infidelity $<1$%). Our approach provides a route for
transcending this limitation on the quantum efficiency of non-reciprocal
interactions.

The results of our study have implications in several areas: (1) In the
context of quantum magnonics, we present the first study of polariton modes
with a partially magnetized ferrite material, which features a high quality
factor and low operating field, both of which are crucial for constructing
superconducting-magnonic devices. (2) In the context of modular
superconducting quantum computing, we demonstrate the first circulator with
internal loss well below 1% of the coupling bandwidth, which would enable
high-fidelity directional quantum state transfer. (3) For the general non-
Hermitian physics, we demonstrate an experimental probe of the non-reciprocal
eigenvector composition of a non-Hermitian system. Combining these advances,
we have established an experimental platform that meets the conditions for
future study of nonlinear non-reciprocal interactions with superconducting
qubits.

##  II Experimental Setup

Figure 1: Device and measurement setup. (a) A YIG cylinder (black) is placed
at the center of the intersection of three rounded-rectangular waveguides
placed 120 degrees away from each other. The light grey region is vacuum
inside an oxygen-free copper enclosure. The device can be assembled in two
different configurations: First, a drum-head shaped transition pin can be
attached at the end of each waveguide section to form an impedance matched
waveguide-to-SMA transition (IMT). Alternatively, a short weakly-coupled probe
(WCP) can be attached to each waveguide section to explore the internal modes
of the device. (b) The device is mounted to a mezzanine plate that is
thermalized to the mixing chamber of a dilution refrigerator, and is
positioned at the center of a superconducting solenoid magnet which operates
at 4K. The device is connected to three input cables (with attenuators as
marked) and two output amplifier lines (with directional couplers splitting
signals) for $S$-parameter measurements using a vector network analyzer (VNA).

Our experimental setup is shown in Fig. 1(a). Three rounded-rectangular
waveguides, each with a cross section of 21.0 mm $\times$ 4.0 mm, placed 120
degrees away from each other, intersect to form the body of the circulator. A
$\phi$-5.58 mm $\times$ 5.0 mm single-crystalline YIG cylinder is placed at
the center of the Y-junction, with external magnetic fields applied along its
height (the $z$ axis and the [111] orientation of the YIG crystal). At the end
of the three waveguide sections, we can either attach impedance-matched
waveguide-to-SMA transitions (IMT) to perform standard characterization of the
circulator (as in Section IV), or attach weakly-coupled probe pins (WCP) to
explore the internal modes of this YIG-loaded Y-shaped cavity (as in Section
III). The use of reconfigurable probes in the same waveguide package allows us
to infer the operation condition and the performance of the circulator from
the properties of the internal modes. Furthermore, the copper waveguide
sections can be replaced by superconducting niobium cavities, with details to
be described in Section V and Fig. 5. This modular substitution introduces
additional external high Q modes to the system, and understanding the
resulting Hamiltonian and the hybridized mode structure of the full system
will be a first step towards the study of pristine non-reciprocal interactions
in circuit QED.

The device package is thermalized to the mixing chamber plate ($\sim$20 mK) of
a Bluefors LD-250 dilution refrigerator inside the $\phi$-100 mm bore of a 1 T
superconducting magnet that applies magnetic field along the $z$ axis [Fig.
1(b)]. A vector network analyzer is used to measure the complex microwave
transmission coefficients $S_{ij}$ (from Port $j$ to Port $i$, where
$i,j=1,2,3$) of the device in series with a chain of attenuators, filters and
amplifiers as in typical circuit QED experiments. A magnetic shield made of a
steel sheet is placed outside the bottom half of the refrigerator, and all
data is acquired under the persistent mode of the superconducting magnet to
minimize magnetic-field fluctuations.

##  III Internal mode structure

Figure 2: Internal mode spectrum of the device. (a) VNA transmission
measurement $S_{21}$ of multi-mode photon-magnon hybrid system formed in the
waveguide circulator package with WCP. The blue (red) dashed line plots the
frequency of the clockwise(counterclockwise) mode from a simplified two-mode
model of photon-magnon avoided crossing with $g/2\pi$ = 1.3 GHz (2.1 GHz) to
compare with an observed spectral line. (b) The right panel shows a finer
sweep of $S_{21}$ in the low-field regime. The mode frequencies differ
slightly from (a) since the data was acquired after some modifications to the
device packaging (a piece of Teflon spacer at the top of the YIG cylinder was
removed). The left panel shows the electromagnetic mode structures of the
eigenmode solutions from our HFSS simulation for the WCP with good frequency
agreement to the experimental data (see Fig. 8 in Appendix). The color scale
from red to blue represents electric field strength from high to low in log
scale. The pair of modes around 11 GHz are connected to the circulating modes
of the loaded circulator and (c) shows their linewidths.

We begin by discussing the internal mode structure of the device, as probed by
$S_{21}$ as a function of applied magnetic field $B$ when the device is
installed with WCP [Fig. 2(a)]. A series of electromagnetic modes (relatively
field-independent) are observed to undergo large avoided crossings with a
cluster of magnon modes of the YIG crystal, forming photon-magnon polariton
modes. The magnon mode most strongly coupled to photons is known to correspond
to near-uniform precession of YIG spins, or the Kittel mode of FMR, whose
frequency increases linearly with magnetic field:
$\omega_{m}=\gamma[B+\mu_{0}(N_{x,y}-N_{z})M_{s}]\approx\gamma B$, as marked
by the dashed line in Fig. 2(a). Here $\gamma$ is the gyromagnetic ratio, and
the (volume-averaged) demagnetizing factors $N_{x,y,z}$ in magnetic saturated
state are very close to $1/3$ for the aspect ratio of our YIG cylinder Chen
_et al._ (1991). These avoided crossings are similar to previous experiments
showing strong photon-magnon coupling Tabuchi  _et al._ (2014); Zhang  _et
al._ (2014), but due to the much larger size of the YIG in our experiment, a
large cluster of higher-order magnetostatic modes, most of which have slightly
higher frequency than the Kittel mode Fletcher  _et al._ (1960); Klingler  _et
al._ (2017) also coherently interact with the microwave photons, contributing
to the complex transmission spectra in the vicinity of the crossings.
Nevertheless, to have a coarse estimate of the photon-magnon coupling
strength, it is convenient to model each observed spectral line far away from
the crossing region as a bare electromagnetic mode with frequency
$\omega_{c}/2\pi$ hybridized with a single combined magnon mode. The implied
coupling strengths $g/2\pi$ (in the cavity QED convention) are about 1.2 GHz
and 2.1 GHz for the two modes of particular interest to this study [blue and
red in Fig. 2(a)], placing the mode hybridization in the ultrastrong coupling
regime (see e.g. Marković  _et al._ (2018)) with
$g/(\omega_{c}+\omega_{m})\sim 10\%$. Even at $B=0$, with a photon-magnon
detuning of $\Delta=\omega_{c}-\omega_{m}\approx 2\pi\cdot 10$ GHz, the
participation of magnon excitations in the photon-branch of the polariton
modes remains quite substantial.

Using finite-element simulations (Ansys HFSS, Appendix A), we identify that
the five polariton modes in the frequency range of 8-12 GHz at $B=0$ include
two nearly-degenerate mode pairs with two-fold symmetry and another mode with
three-fold symmetry. Electric field distributions of each of the modes are
illustrated in Fig. 2(b). Each degenerate mode pair can be understood using a
basis of standing-wave modes polarized along the $x$ or $y$ direction. The
application of a magnetic field lifts this $x$-$y$ degeneracy, as the mode
pair forms clockwise and counterclockwise rotating eigenmodes with a frequency
splitting Owens  _et al._ (2018); Zhang  _et al._ (2020); Anderson  _et al._
(2016).

Prior use of these chiral polariton mode pairs have been in the magnetically
saturated regime Anderson  _et al._ (2016); Owens  _et al._ (2018); Zhang  _et
al._ (2020). Here we focus on the low-field regime ($|B|<0.05$ mT) where the
approximately linear increase of frequency splitting between the mode pair
reflects increasing magnetization of YIG under increasing applied magnetic
field. After implementing demagnetization training cycles to suppress a
relatively small hysteretic effect throughout our experiments, we expect an
approximately linear magnetization curve ($M$-$H$) for YIG until it approaches
magnetic saturation. In the limit of high permeability $\mu\gg\mu_{0}$ (with
$\mu_{0}$ being the vacuum permeability), we have $M=B/{\mu_{0}N_{z}}$ (note
that $B$ is the applied magnetic field strength) and $N_{z}\approx 0.285$ is
the $z$-direction demagnetizing factor when the YIG is significantly below
magnetic saturation Chen  _et al._ (1991). Saturation magnetization
$M_{s}=2440$ Oe Solt (1962) of YIG is approached on the scale of
$B\sim\mu_{0}N_{z}M_{s}\approx 70$ mT, which agrees with the changing
curvature of the mode-splitting spectra.

On the other hand, in the completely demagnetized state ($M=0$) at zero field,
the system is expected to satisfy macroscopic time-reversal symmetry. As
supported by numerical simulations, the $x$-$y$ mode pairs should be in
principle exactly degenerate since both the Y-junction geometry and the [111]
YIG crystal has 3-fold rotational symmetry around the $z$ axis. However,
appreciable zero-field splitting is observed experimentally. We attribute this
splitting to some anisotropy in the x-y plane breaking this symmetry and
allowing a preferred magnetization axis of the YIG at 0 field. Some possible
explanations for this anisotropy are a small visible damage to our YIG crystal
on one edge or possible imperfections in eccentricity and alignment. If the
magnetic domains of unsaturated YIG preferentially align with one in-plane
axis compared to its orthogonal axis within the $x$-$y$ plane, this anisotropy
would result in a relative frequency shift between the standing-wave modes
along the in-plane easy and hard axes. This anisotropy-induced frequency shift
$\pm\beta$ for the $x$ and $y$ modes can be modeled in numerical simulations
employing a permeability tensor of unsaturated ferromagnets Schlömann (1970);
Green and Sandy (1974) with certain anisotropic assumption, which can
plausibly explain the data (Appendix A). As $B$ increases, we expect $\beta$
to decay towards 0 when the magnetic domains are increasingly aligned towards
the $z$ direction, thus making any $x$-$y$ plane energetic preference of
negligible effect. We model this decay with a thermodynamic toy model
(Appendix B) whose details do not affect the conclusions of this study.

For the rest of this article, we will focus on the pair of polariton modes
near 11 GHz in Fig. 2(b), and refer to them as “the circulator modes” for
reasons that will become apparent. We can model their frequencies in the
partially magnetized regime ($|B|<50$ mT) using a phenomenological model
accounting for the degeneracy-lifting anisotropy and the field-dependent
magnetization of YIG. Let the zero-field frequencies of the $x$ and $y$ modes
be $\omega_{x}=\omega_{y}$ if the device had perfect 3-fold symmetry, $\beta$
and $\theta/2$ be the magnitude of anisotropy caused degeneracy-lifting and
the direction of the in-plane anisotropy axis (relative to the $x$ axis), and
off-diagonal imaginary coupling term $\pm ikB$ be the magnetic field induced
degeneracy-lifting, linearly increasing with a real coefficient $k$. We use
the following Hamiltonian to characterize the pair of circulator modes in the
basis of $x$ and $y$ mode amplitudes:

$H/\hbar=\begin{pmatrix}\omega_{x}+\beta\cos{\theta}+mB^{2}&\beta\sin{\theta}+ikB\\\
\beta\sin{\theta}–ikB&\omega_{y}-\beta\cos{\theta}+mB^{2}\end{pmatrix}$ (1)

This effective model of the polariton modes has absorbed the magnon
contributions in the regime where they have been adiabatically eliminated. The
formation of clockwise and counterclockwise eigenmodes is due to magnon-
mediated interactions modeled by $\pm ikB$. The level repulsion from the far-
detuned magnon modes is approximated by a small quadratic shift in frequency
$mB^{2}$. The quadratic dependence was empirically chosen because the sum of
the mode frequencies over field displayed a roughly quadratic relationship
with $B$ over the plotted field range. By fitting the mode spectrum in Fig.
2(b), we obtain $\omega_{x}/2\pi=\omega_{y}/2\pi=11.054$ GHz, $k/2\pi=9.82$
GHz/T, $m/2\pi=50$ GHz/T2, $\beta/2\pi=139$ MHz.

Figure 3: Illustration of the circulator working principle and low-temperature
characterization of the non-reciprocity. In the circulator package with IMT,
the frequency splitting of clockwise and counterclockwise rotating modes as
shown in (a) can be tuned such that the phase of the modes are $\pi/6$ and
$-\pi/6$ as shown in (b). This then produces a node at the upper port, thereby
preventing any signal from leaving there at all times where $\omega t=0$ and
$\omega t=\pi/4$ are shown pictorially in (c). (d, e) Measured microwave
transmission (d) $|S_{12}|$ and (e) $|S_{21}|$ spectra as a function of
magnetic field B. (f-i) The isolation performance, (f)
$\mathcal{I}_{12}=|S_{12}/S_{21}|$,(h) $\mathcal{I}_{21}=|S_{21}/S_{12}|$, (g)
$\mathcal{I}_{23}=|S_{23}/S_{32}|$, (i) $\mathcal{I}_{13}=|S_{13}/S_{31}|$.
$S_{21}$ is obtained by measuring the $S_{12}$ at –$B$, which provides a self-
calibrated way to determine the isolation of the circulator.

It is well-known that the FMR modes of partially magnetized ferrimagnetic
insulators, where the magnetic domains are not aligned in equilibrium, have
large damping. Therefore, one may expect broad linewidths for photon-magnon
polariton modes below magnetic saturation. Indeed, we have observed linewidths
exceeding 100 MHz for another polariton mode at 5 GHz at $B<$ 50 mT (not
shown). Surprisingly, the polariton modes at higher frequency display narrow
linewidths, $\kappa_{i}\approx 2$ MHz for the pair of circulator modes [Fig.
2(c)], which corresponds to quality factors on par with some circuit QED
elements such as the readout resonators. The narrow internal linewidth of the
circulator modes is crucial for constructing a low-loss circulator and
eventually achieving high quantum efficiency of non-reciprocal interactions in
circuit QED. It is primarily aided by the use of single crystalline YIG and
the relatively low magnon participation in the circulator modes compared to
commercial circulators. The observed $\kappa_{i}$ may be limited by either the
spin relaxation in YIG or the Ohmic loss in copper. The former remains to be
investigated in this partially magnetized regime, and the latter may be
further reduced through better surface treatment or the use of superconducting
materials in low-field regions of the waveguide package.

##  IV Circulator characterization

The device acts as a circulator when the end of each waveguide section is in
IMT rather than WCP with an applied magnetic field in the $\hat{z}$ direction.
In this configuration, the linewidths of all internal modes are substantially
broadened forming a transmission continuum in the measurement, as shown in
Fig. 3(d,e) for $|S_{12}|$ and $|S_{21}|$. Nevertheless, the operating
condition of the circulator can be conceptually understood as having a pair of
counter-propagating internal modes with their magnetic-field-induced splitting
($\delta$) satisfying the relationship $\delta=2\kappa_{c}/{\sqrt{3}}$ versus
their half linewidths ($\kappa_{c}$) Fay and Comstock (1965). As illustrated
in Fig. 3(a-c), when driven at a frequency in the middle of the two
resonances, the two circulator modes are excited with equal amplitude and a
phase shift of $\pm\frac{\pi}{6}$ relative to the drive. The resultant
standing wave pattern forms a node at the isolation port of the circulator.
This condition can be satisfied by choosing the correct combination of
frequency and magnetic field.

We characterize the non-reciprocity of the circulator by the isolation ratio
$\mathcal{I}_{12}=|S_{12}/S_{21}|$, which may be computed from Fig. 3(d,e).
However, since $S_{12}$ and $S_{21}$ are measured through different cables and
amplifier chains [Fig. 1(b)], it is challenging to calibrate their absolute
values precisely. A much better self-calibrated technique to extract the
isolation ratio in our system is to use the Onsager-Casimir relation Casimir
(1945), $S_{21}(B)=S_{12}(-B)$, resulting from the microscopic time reversal
symmetry. Therefore, we use $\mathcal{I}_{12}=|S_{12}(B)/S_{12}(-B)|$ to
determine the isolation ratio of the circulator as shown in Fig. 3(f), with
the (field-independent) contribution from same transmission chain cancelled
out. The result indicates the circulator working condition is met for the pair
of counter-propagating modes at $\sim$11.2 GHz with external field $\sim$0.022
T. We see $\geq$20 dB of isolation over a bandwidth of about 250 MHz, with
maximum isolation of at least 35 dB.

The same analysis on $S_{21}$ data yields the same isolation property [Fig.
3(h)] as expected. Similarly, $\mathcal{I}_{23}$ and $\mathcal{I}_{13}$ are
measured as in Fig. 3(g) and (i), each showing a slightly different working
field and frequency (possibly due to imperfections of the device geometry) but
similar isolation magnitude and bandwidth. These data are measured at an
estimated circulating photon number on the order of 10’s, but when we lower
the power to below single photon level, the isolation property does not show
notable changes.

An important motivation of our work is to ultimately implement pristine non-
reciprocal interactions between superconducting qubits or cavities. It is
crucial to minimize the ratio between the undesirable internal dissipation
($\kappa_{i}$) and the external bath coupling ($\kappa_{c}$) that enables non-
reciprocity. In the case of a circulator, this ratio sets the limit for the
circulator’s microwave insertion loss $\mathcal{L}_{21}$ Kord  _et al._
(2020); Fay and Comstock (1965):

$\mathcal{L}_{21}=1-|S_{21}|^{2}\geq
1-|S_{21}|^{2}-|S_{11}|^{2}-|S_{31}|^{2}\approx\frac{\kappa_{i}}{\kappa_{c}}$
(2)

Figure 4: Characterization of the internal loss of the circulator at room
temperature. (a) Linewidths of the pair of circulator modes measured at room
temperature. (b) Transmission $S_{21}$ and reflection $S_{11}$ near the
maximum isolation regime of the circulator, measured at $B=24.8$ mT (top
panel) and the internal loss of the circulator calculated from it (bottom
panel) at room temperature.

This lower limit is obtained in principle when the circulator has perfect
impedance matching ($S_{11}=0$) and isolation ratio ($S_{31}=0$). Typical
commercial ferrite circulators used in circuit QED experiments have shown
insertion loss around 10% Kurpiers  _et al._ (2018); Axline  _et al._ (2018),
which is dominated by internal loss. Experimental Josephson circulators so far
have also reported insertion loss of -0.5 dB (11%) or higher Chapman  _et al._
(2017); Lecocq  _et al._ (2017). The lowest quoted insertion loss for a
commercially-listed waveguide circulator is -0.1 dB (or 2.2%) but that is
untested in the quantum regime. In order for the quantum efficiency of a non-
reciprocal two-qubit interaction channel to match the fidelity of state-of-
the-art two-qubit operations, the insertion loss would need to be improved to
the sub-percent level.

To the best of our knowledge, it is an open challenge to calibrate the
insertion loss of a microwave component in a dilution refrigerator with a
precision better than 1%. Even using specialized Thru-Reflect-Line calibration
components and well-characterized cryogenic switches, the resultant precision
would still be limited to about 0.1dB (or 2.3%) Ranzani  _et al._ (2013). In
order to infer the loss of our circulator at 20 mK, we measure its $S$
parameters at room temperature after a careful calibration procedure that uses
attenuators in series to suppress standing waves. We find a conservative upper
bound for room-temperature internal loss of $\leq
1–|S_{21}|^{2}–|S_{11}|^{2}\approx 2\%$, as shown in Fig. 4(b). Assuming
$\kappa_{c}$ does not change as a function of temperature, comparing the
intrinsic linewidth of the circulator mode pair at room temperature versus 20
mK would inform the internal loss at 20 mK. The intrinsic linewidths, measured
in WCP, are 4.1 and 6.3 MHz at room temperature [Fig. 4(a)] and 1.8 and 2.2
MHz at low temperature [Fig. 2(c)], indicating that $\kappa_{c}\gtrsim 260$
MHz and $\kappa_{i}/\kappa_{c}\lesssim 0.8\%$. If we instead use the relation
of $\kappa_{c}=\sqrt{3}\delta/2$, which yields $\kappa_{c}$ in the range of
430 MHz to 550 MHz (and data in Section V would further suggest $\kappa_{c}$
at the high end of this range), or $\kappa_{i}/\kappa_{c}\approx 0.4\%$.
Further improvement of the circulator bandwidth and the coupling ratio can be
achieved by applying impedance transformation techniques to increase
$\kappa_{c}$ Helszajn (2008).

Translating this small internal loss ratio to a sub-percent insertion loss for
a circulator as a peripheral transmission-line device would further require
excellent impedance matching. However, we emphasize that this requirement is
not fundamental if the circulator is modeled as part of the quantum system
itself mediating interactions between other quantum resonance modes. Unlike
most ferrite circulators, our device operates in the regime of partial
magnetization for YIG. It only requires a moderate external magnetic field
that is significantly below the critical field of a variety of superconducting
materials. This allows for 3D integration of superconducting niobium cavities
and shielded transmon qubits for studying circuit QED with non-reciprocal
interactions. In the following section, we demonstrate direct coupling of two
external superconducting cavity modes with the circulator modes and analyze
the resultant non-reciprocal hybrid system as a whole.

##  V Tuning non-reciprocity of eigenmode structure

Figure 5: Waveguide circulator-cavity integration. (a) The photo image, (b) a
schematic top-down view, and (c) a diagrammatic illustration of the effective
Hamiltonian (see Eq. (4), for clarity the $\beta$ and $mB^{2}$ terms have been
neglected in the illustration) of our integrated non-reciprocal device. It is
composed of a Cu waveguide Y-junction loaded with a YIG cylinder, two Nb
cavity segments with weakly-coupled drive ports (Port 1 and 2), and an output
port with IMT (Port 3). For each cavity, the sidewall closest to the copper
Y-junction is formed by a standalone niobium plate in the assembly [enclosed
in blue in (a)], which contains a 5 mm-diameter aperture to create an
evanescent coupling between the superconducting cavity mode and the circulator
modes. One of the cavities is loaded with a transmon qubit [marked as $\times$
in (b)] which stays unused in its ground state in this study.

We integrate superconducting cavities with the ferrite device by replacing the
rectangular waveguide extensions with superconducting 3D cavities made of
niobium [Fig. 5(a)]. Two cavities, attached at Port 1 and 2, are tuned to have
resonance frequencies close to each other, $\omega_{1}\approx\omega_{2}\sim
10.8$ GHz, both of which are within the bandwidth of the circulator. Each
cavity is coupled to the central Y-junction via a coupling aperture. As a
result, the circulator modes will mediate an interaction between these two
external cavities. Crucially, this circulator-mediated interaction can have
both coherent and dissipative aspects, and can be non-reciprocal. The degree
and the direction of non-reciprocity of the coupling can be tuned via the
external magnetic field. Note that Port 3 remains impedance-matched to a
transmission line. This is also essential: it serves as the dominant
dissipative bath that is necessary for achieving non-reciprocal inter-mode
interactions Metelmann and Clerk (2015).

To probe the hybridized mode structure of the composite system, we measure
$S_{31}$ from a weakly-coupled drive port on Cavity 1 to Port 3. The measured
amplitude of $S_{31}$ as a function of magnetic field and frequency is shown
in Fig. 6(a). There are a total of four bare oscillator modes in the vicinity
(within 0.5 GHz) of the frequency range of interest: two superconducting
cavity modes and two internal circulator modes. Since the loaded circulator
modes with very broad linewidths ($>100$ MHz) are difficult to observe in the
presence of the standing-wave background of the coaxial cables, this
spectroscopy measurement primarily reveals the eigenmodes that are localized
in the external superconducting cavities. Indeed, at $|B|>0.03$ T, the
spectrum shows two sharp resonances which we identify as the bare cavity modes
to a good approximation. At lower fields, the cavities appear to more strongly
hybridize with each other and with the lossy circulator modes, but the
spectrum can nontheless be captured relatively well by the sum of two
Lorentzian modes $a$ and $b$:

$S_{31}=\frac{A_{a}e^{i\phi_{a}}}{{-i(\omega-\omega_{a})-\kappa_{a}/2}}+\frac{A_{b}e^{i\phi_{b}}}{{-i(\omega-\omega_{b})-\kappa_{b}/2}}$
(3)

By fitting the spectrum to Eq. (3), we can extract the linewidth
($\kappa_{i}$), frequency ($\omega_{i}$) and amplitude ($A_{i}$) of the two
Lorentzians at each magnetic field, as plotted in Fig. 6(c-e).

The magnetic field dependence of the two prominent Lorentzians can be
connected to the eigenmode solutions of an effective Hamiltonian model of
system. We describe the system using the following 4$\times$4 non-Hermitian
matrix, written in the basis of the amplitudes of the two cavity modes and the
two circulator modes:

$H_{\mathrm{eff}}/\hbar=\\\
\begin{pmatrix}\;\omega_{1}-i\frac{\kappa_{1}}{2}\;&0&g_{y}&g_{x}\\\
0&\;\omega_{2}-i\frac{\kappa_{2}}{2}\;&g_{y}&-g_{x}\\\
g_{y}&g_{y}&\omega_{y}-\beta\cos\theta+mB^{2}-i\frac{\kappa_{3}}{2}&\beta\sin\theta-
ikB\\\ g_{x}\ &-g_{x}\
&\beta\sin\theta+ikB&\omega_{x}+\beta\cos\theta+mB^{2}\\\ \end{pmatrix}$ (4)

The two niobium cavities have bare frequencies $\omega_{1}$, $\omega_{2}$, and
input coupling rates of $\kappa_{1}$ and $\kappa_{2}$. The bottom right block
of Eq. (4) describes the two circulator modes, with their anisotropy
dependence and imaginary coupling due to magnon hybridization following the
same description as in Eq. (1). The zero-field frequencies of the two
circulator modes $\omega_{x}$, $\omega_{y}$ are no longer equal since the
device is no longer 3-fold symmetric. The $y$-mode with frequency $\omega_{y}$
is symmetric with respect to the $y$ axis, and therefore has an equal and in-
phase coupling rate $g_{y}$ with the two cavities. It rapidly leaks to the
waveguide output Port 3, with a decay rate
$\kappa_{3}\gg\kappa_{1},\kappa_{2},g_{x},g_{y}$. $\kappa_{3}$ is related to
$\kappa_{c}$ of the loaded circulator as in Section IV by
$\kappa_{3}=4\kappa_{c}/3$. The $x$-mode is anti-symmetric with respect to the
$y$ axis, preventing it from coupling to the output port. This also leads to a
180∘ phase difference in cavity coupling as accounted for by the negative sign
on two of the $g_{x}$ parameters

Figure 6: Spectroscopy of the hybridized non-reciprocal modes of a circulator-
cavity system. (a) VNA transmission measurement and (b) model prediction of
$|S_{31}|$ frequency spectrum over external magnetic field $B$. Remaining
panels show magnetic field dependence of the system’s eigenmodes and
wavefunctions: (c) eigenmode linewidths $\kappa_{n}/2\pi$, (d) eigenmode
frequencies $\omega_{n}/2\pi$, (e) amplitude parameter $A_{n}$ (c.f. Eq.(3)),
and (f) amplitude ratio (c.f. Eq.(8)) of experimental data (dots) from two-
mode Lorentzian fit (Eq.(3)) and theory predictions (dashed lines). The
symmetry of $\kappa_{n}$ and $\omega_{n}$ ($i.e.~{}$ complex eigen-energy of
the hybrid system) with respect to $B$ exemplifies the microscopic time-
reversal symmetry of the non-Hermitian system. The non-reciprocity is
reflected in the difference in $A_{n}$ at $\pm B$, which reveals the asymmetry
in the left/right eigen-vector structure (c.f. Eq.(8))). The effective
Hamiltonian parameters from the fit are: $\omega_{1}/2\pi=10.8104$ GHz,
$\omega_{2}/2\pi=10.8040$ GHz, $\omega_{x}/2\pi=10.707$ GHz,
$\omega_{y}/2\pi=10.813$ GHz, $\theta=37.7^{\circ}$, $\kappa_{3}/2\pi=730$
MHz, $g_{x}/2\pi=(9.0+0.011\beta)$ MHz, $g_{y}/2\pi=(5.0+0.006\beta)$ MHz with
$\beta/2\pi=139$ MHz at $B=0$ and decays with $|B|$.

This effective non-Hermitian Hamiltonian can be diagonalized as:

$H_{\mathrm{eff}}=\sum_{n}\hbar\omega_{n}\ket{n_{R}}\bra{n_{L}}$ (5)

where $n=a,b,c,d$ are the eigenmode indices of the system, $\omega_{n}$ the
complex eigen-frequencies, and $\ket{n_{R}}$ and $\ket{n_{L}}$ the right and
left eigenvectors of the non-Hermitian Hamiltonian, defined as:
$H_{\mathrm{eff}}\ket{n_{R}}=\hbar\omega_{n}\ket{n_{R}}$ and
$H_{\mathrm{eff}}^{\dagger}\ket{n_{L}}=\hbar\omega_{n}^{*}\ket{n_{L}}$. The
scattering matrix element $S_{ij}$ from Port $j$ to Port $i$ can be generally
drived from the input-output theory relation:
$S_{ij}=\delta_{ij}-i\sqrt{\kappa_{i}\kappa_{j}}G_{ij}^{R}(\omega)$, where the
$4\times 4$ retarded matrix Green’s function is defined as:
$G^{R}(\omega)=(\omega-H_{\mathrm{eff}})^{-1}$, and $\kappa_{i}$ and
$\kappa_{j}$ are the output and input coupling rates, respectively. Applying
this formalism to the $S_{31}$ measurement of our device, we arrive at the
following Lorentzian spectral decomposition to describe the spectrum:

$S_{31}(\omega)=\sum_{n}\frac{-i\sqrt{\kappa_{1}\kappa_{3}}\bra{y}\ket{n_{R}}\bra{n_{L}}\ket{1}}{\omega-\omega_{n}}$
(6)

where the real and imaginary parts of the eigen-frequency $\omega_{n}$
correspond to the observed Lorentzian frequencies and half linewidths,
respectively. The amplitudes of the Lorentzians are proportional to the
product of the left eigenvector overlap with the bare cavity mode $\ket{1}$
and the right eigenvector overlap with the output circulator mode $\ket{y}$.

By fitting the extracted Lorentzian parameters of the two prominent eigenmodes
in Fig. 6(c-e) to the predictions of the $4\times 4$ Hamiltonian model across
all fields (Eq. 4), we can determine all the free Hamiltonian parameters in
this model. This includes $\kappa_{3}=730$ MHz, implying $\kappa_{c}=550$ MHz
for the loaded circulator, consistent with (and at the high end of) the
estimates in Section IV. Somewhat surprisingly, the experimental data strongly
suggests that the cavity-circulator coupling rates $g_{x}$ and $g_{y}$ must be
magnetic field dependent. (For example, it heavily constrains that
$g_{x}/2\pi>16$ MHz near $B$ = 0 and $g_{x}/2\pi<12$ MHz at $|B|>30$ mT.) We
attribute this varying coupling to the change in electromagnetic field
distribution of the $x$\- and $y$-modes around the coupling aperture due to
the $x$-$y$ anisotropy of YIG. Assuming $g_{x}$ and $g_{y}$ contains a
contribution proportional to $\beta(B)$ with the same decay shape over applied
field, the effective Hamiltonian model fits the Lorentzian parameters quite
well and also reproduces the overall transmission spectrum [Fig. 6(b)].

The eigenmode features of the system can be understood intuitively by
considering first the inter-mixing of the $x$, $y$ circulator modes (i.e.
diagonalization of the lower right block of $H_{\mathrm{eff}}$) and then their
mixing with the two cavity modes. At $B=0$, the circulator modes are
relatively close in frequency to the bare cavities, resulting in strong four-
mode hybridization and substantial linewidth-broadening and frequency shift to
Mode $a$. As $B$ increases, the block-diagonalized circulator modes split
further in frequency in response to increasing magnetization of YIG (analogous
to the unloaded internal mode spectrum in Fig. 2b), and become more detuned
from the bare cavities, so the cavity-circulator hybridization is continuously
reduced. This is reflected in the eventual flattening of the frequency and
linewidth of the observed Lorentzians at high fields.

In our device, opposite magnetic fields produce opposite directions of non-
reciprocity, hence the transmission spectra observed at $\pm B$ in Fig. 6(a)
are markedly different. Interestingly, the extracted data in Fig. 6(c,d) shows
that the underlying eigenmode frequency and linewidths at $\pm B$ are equal,
unchanged under the mapping $\mathcal{P}$ of $B:\mapsto-B$. This is no
coincidence, but is rather the direct consequence of microscopic symmetry
requirements. Recall again that the Onsager-Casimir relation Casimir (1945)
requires that the full scattering matrix $S$ satisfy
$S(-B)=S^{\mathrm{T}}(B)$. As $S$ is however directly determined by our non-
Hermitian Hamiltonian, this necessarily implies that
$H^{\phantom{T}}_{\mathrm{eff}}(-B)=H_{\mathrm{eff}}^{\mathrm{T}}(B)$. This in
turn implies that the complex eigenvalues of $H_{\mathrm{eff}}$ are unchanged
under $\mathcal{P}$. Note that the operation $\mathcal{P}$ is not just a
simple time-reversal operation, as it does not involve transforming loss to
gain (and vice versa). This property of $H_{\mathrm{eff}}$ can be easily seen
to hold for our specific model in Eq. (4). Nonetheless, we emphasize that our
experimental observation here of eigenvalue invariance under the mapping
$\mathcal{P}$ is a demonstration of a general physical property; it is by no
means contingent on the specifics of our model.

While the eigenvalues of $H_{\mathrm{eff}}$ do not directly reflect the non-
reciprocal physics of our system, the same is not true of its eigenvectors. As
it involves matrix transposition, the operation $\mathcal{P}$ exchanges the
left and right eigenvectors of the effective Hamiltonian:
$\ket{n_{R}(B)}=\ket{n_{L}(-B)}^{*}$. A defining feature of a non-reciprocal
Hamiltonian is that the left and right eigenvectors generally differ in their
spatial structures (i.e. they look very different when expressed in a basis of
bare modes):

$R_{i,n}=\frac{|\bra{n_{L}}\ket{i}|}{|\bra{i}\ket{n_{R}}|}\neq 1$ (7)

As has been discussed elsewhere McDonald and Clerk (2020); Schomerus (2020),
the $R_{i,n}$ characterizes a fundamental asymmetry in the response of our
system. The numerator characterizes the susceptibility of the eigenmode $n$ to
a perturbation or excitation entering from bare mode $i$. In contrast, the
denominator tells us the amplitude on bare mode $i$ that would result given
that the system eigenmode $n$ is excited. In a Hermitian system these
quantities are necessarily identical, expressing a fundamental kind of
reciprocity between susceptibility and response. In our non-Hermitian system,
the non-unity ratio here reflects the effective non-reciprocity of the inter-
mode interactions.

This non-reciprocal eigenvector structure is experimentally verified by the
asymmetry of the Lorentzian amplitudes with respect to $B$ in Fig. 6(e),

$\frac{A_{n}(B)}{A_{n}(-B)}=\bigg{|}\frac{\bra{n_{L}}\ket{1}}{\bra{1}\ket{n_{R}}}\frac{\bra{y}\ket{n_{R}}}{\bra{n_{L}}\ket{y}}\bigg{|}=\frac{R_{1,n}}{R_{y,n}}$
(8)

We plot this ratio in Fig. 6(f). For our device, a calculation based on Eq.
(4) shows that $R_{y,n}\approx 1$ for most of the field range (near zero field
and $|B|\gtrsim 15$ mT), allowing Fig. 6(f) to be understood as a measurement
of the non-reciprocity ratio $R_{1,n}$, in this field range, showing the role
of Cavity 1 in the two prominent eigenmodes of the system. In particular, the
most pronounced asymmetry is observed near the optimal working point of the
circulator ($B=\pm 28$ mT) for Mode $b$, which can leak through cavity 1 but
cannot be excited from Cavity 1 or vice versa, as expected for a mode
dominated by photons in Cavity 2.

Figure 7: Mediated non-reciprocal coupling rates between the external
superconducting cavities. Red curves show the off-diagonal coupling terms in
the effective two-mode Hamiltonian (c.f. Eq.(9)), $\lvert H_{12}\rvert$ and
$\lvert H_{21}\rvert$, as a function of magnetic field, and blue shows
$r=\sqrt{\lvert H_{21}\rvert/\lvert H_{12}\rvert}$.

It is also interesting to discuss our system in the context of general systems
exhibiting non-reciprocal interactions between constituent parts. Such systems
are commonly described by phenomenological non-Hermitian Hamiltonian matrices
$H$, whose matrix elements in a local basis encode interactions with
directionality: $|H_{ij}|\neq|H_{ji}|$. A prominent example is the Hatano-
Nelson model Hatano and Nelson (1997) of asymmetric tunneling on a lattice. In
our case, we have a microscopically-motivated model that is fully consistent
with the requirements of microscopic reversability, but which encodes non-
reciprocity. As shown in appendix C, one can adiabatically eliminate the
internal circulator modes from our system to obtain an effective two-mode non-
Hermitian Hamiltonian that describes the external cavity modes and their
circulator-mediated interaction:

$\displaystyle
H^{\prime}_{\mathrm{eff}}/\hbar=\begin{pmatrix}\omega_{1,\mathrm{eff}}-i\frac{\kappa_{1,\mathrm{eff}}}{2}&H_{12}\\\
H_{21}&\omega_{2,\mathrm{eff}}-i\frac{\kappa_{2,\mathrm{eff}}}{2}\\\
\end{pmatrix}$ (9)

The field-tunable non-reciprocity can be seen in the asymmetry of the off-
diagonal coupling values $H_{12}$ and $H_{21}$ based on the model, as plotted
in Fig. 7. We note that the scale of $H_{12}$ and $H_{21}$ of a few MHz (which
can be increased by using a larger coupling hole) is much larger than the
achievable internal loss of the superconducting cavities , making the non-
reciprocal coupling the dominant interaction.

Furthermore, this reduced Hamiltonian and its eigenvectors can be mapped to a
reciprocal Hamiltonian and associated eigenvectors using a similarity
transformation $S(r)$, where $\sqrt{\lvert H_{21}\rvert/\lvert
H_{12}\rvert}\equiv r$, as outlined in Appendix C. The similarity transform
effectively localizes the mode participation on the lattice site in the
direction of stronger coupling more than would be expected in the reciprocal
case, causing the amplitude ratios ($R_{i,n}$) to deviate from 1, which can be
viewed as a consequence of the non-Hermitian skin effect on a two site Hatano-
Nelson lattice Yao and Wang (2018). Furthermore,the similarity transformation
can be used to explain the qualitative behavior of the disparate amplitude
ratios seen in Fig. 6(f). In particular, for $B\gtrsim 20$ mT, we have
$R_{1,a}\approx 1$ and $R_{1,b}\approx r^{2}$ (see Appendix C for details).

##  VI Outlook

In this work, we have revisited the working principles of a Y-junction ferrite
circulator Fay and Comstock (1965), a microwave engineering classic from the
1960’s, in an entirely new context of hybrid quantum systems and non-Hermitian
Hamiltonian. The use of reconfigurable probes and single-crystalline YIG in a
low-loss waveguide package allows us to connect the properties of the photon-
magnon polaritons to the circulator performance. We have further leveraged our
direct access to the internal modes and our ability to tune their coupling in
situ to construct a multi-mode chiral system and unambiguously reveal its non-
reciprocal eigenvector structure. An understanding of the circulator modes and
the non-reciprocal eigen-vector structure of the multi-mode chiral system
provides a foundation for future engineering of any target non-Hermitian
Hamiltonian. This is achieved by our creation of a template model that one can
use to couple any circuit QED element to in order to understand how it would
integrate into the non-reciprocal dynamics, as was done here with the
superconducting cavities.

Looking forward, our device architecture provides a versatile testbed for
studying non-reciprocal interactions in circuit QED by integration of
superconducting qubits. This is enabled by two of its highlighted properties:
the low internal loss of the circulator modes ($<$1% of the demonstrated
coupling rates, compatible with potential high-fidelity operations), and the
relatively low-field operation of the circulator ($\sim$25 mT, below
ferrimagnetic saturation). The latter allows niobium waveguides or cavities to
conveniently act as magnetic shields for superconducting qubits. We have
preliminarily tested that the coherence times of a transmon qubit housed in
one of the niobium cavities are unaffected by in-situ application of a global
magnetic field up to at least 0.1 T. We expect a transmon housed in a niobium
waveguide should receive a similar level of protection from magnetic field.

Direct non-reciprocal coupling of superconducting qubits would open a new
frontier in the study of non-reciprocal dynamics currently dominated by linear
systems Sounas and Alù (2017); Fang  _et al._ (2017); Xu  _et al._ (2019);
Ruesink  _et al._ (2016). The physics of a $N$-mode linear non-reciprocal
system can always be described efficiently by a $N\times N$ non-Hermitian
Hamiltonian matrix (exemplified by our application of such a model) and its
dynamics are always in the classical correspondence limit. Direct
participation of multiple nonlinear modes (such as superconducting qubits) in
non-reciprocal coupling, as envisioned in chiral quantum optics Lodahl  _et
al._ (2017), would lead to novel forms of entanglement stabilization and many-
body phases Stannigel  _et al._ (2012); Ramos  _et al._ (2014). Our system
presents another potential platform to implement this regime in circuit QED in
additional to those proposed using dynamic control Guimond  _et al._ (2020);
Gheeraert  _et al._ (2020). Strong coupling of Josephson circuits with low-
loss non-reciprocal elements can even produce degenerate and protected ground
states for robust encoding of qubits Rymarz  _et al._ (2021).

###### Acknowledgements.

We thank Juliang Li and Dario Rosenstock for experimental assistance. This
research was supported by U.S. Army Research Office under grants
W911-NF-17-1-0469 and W911-NF-19-1-0380.

##  Appendix A Numerical simulation of the ferrite device

Finite element analysis software that supports magnetodynamic simulations,
such as Ansys HFSS, can be used to simulate our ciculator system with both
driven mode and eigenmode solutions. Eigenmode analysis can solve for the
frequency and field distributions of our device’s eigenmodes, while driven
mode analysis reports the S-parameters over frequency. Here we discuss
eigenmode simulations, but driven mode analysis can be carried out similarly.

It is well known that when the applied field is large and magnetization is
saturated along $z$ axis, one can generalize to the whole ferrite the
equations of motion derived from the torque experienced by an electron dipole
moment under the presence of an applied field. This approach, augmented by the
small signal approximation of the Landau-Lifshitz equation of motion, yields
the textbook Polder (relative) permeability tensor:

$[\mu]_{z}=\begin{pmatrix}\mu_{r}&i\kappa&0\\\ -i\kappa&\mu_{r}&0\\\ 0&0&1\\\
\end{pmatrix}$ (10)

where $\mu_{r}=1+\frac{\omega_{0}\omega_{m}}{\omega_{0}^{2}-\omega^{2}}$,
$\kappa=\frac{\omega\omega_{m}}{\omega_{0}^{2}-\omega^{2}}$, with
$\omega_{0}=\gamma\mu_{0}H_{0}$ and $\omega_{m}=\gamma\mu_{0}M_{s}$ being the
internal field strength and saturation magnetization converted to frequencies,
respectively.

The Polder permeability tensor is implimented in HFSS by default to solve for
the interaction of a saturated ferrite with an AC microwave field. However, in
our experiment we operate the circulator at a low bias field, where the
ferrite is not fully saturated. We adopt a permeability tensor model proposed
by Sandy and Green Green and Sandy (1974) for a partially-magnetized ferrite:

$\displaystyle[\mu]_{z}=\begin{pmatrix}\mu_{p}&i\kappa_{p}&0\\\
-i\kappa_{p}&\mu_{p}&0\\\ 0&0&\mu_{z}\\\ \end{pmatrix}$ (11)

where

$\displaystyle\mu_{p}$
$\displaystyle=\mu_{d}+(1-\mu_{d})\bigg{(}\frac{|M_{p}|}{M_{s}}\bigg{)}^{3/2}$
(12) $\displaystyle\kappa_{p}$
$\displaystyle=\kappa\bigg{(}\frac{M_{p}}{M_{s}}\bigg{)}$ (13)
$\displaystyle\mu_{z}$
$\displaystyle=\mu_{d}^{\big{(}1-\frac{|M_{p}|}{M_{s}}\big{)}^{5/2}}$ (14)
$\displaystyle\mu_{d}$
$\displaystyle=\frac{1}{3}+\frac{2}{3}\sqrt{1-\bigg{(}\frac{\omega_{m}}{\omega}\bigg{)}^{2}}$
(15)

with $M_{p}$ being the net magnetization of the partially magnetized ferrite.
This model contains functional forms for $\mu_{p}$ and $\mu_{z}$ that are
purely empirical. However, the expressions for $\kappa_{p}$, which dictates
the chiral splitting of the circulator modes, and $\mu_{d}$, which represents
the permeability in fully demagnetized state, are well motivated Schlömann
(1970).

Figure 8: Eigenfrequencies of the device from finite-element simulations. For
the circulator device with WCP, HFSS eigenmode simulation gives the mode
frequency over different magnetic fields(dots) and agree semi-quantitatively
with experimental data.

This model is implemented in simulations by defining materials with the
customized permeability tensor as given above. Whereas HFSS does not by
default support eigenmode simulations for a ferrite under a DC bias field,
manually defining the permeability tensor components allows us to simulate the
circulator’s eigenmode structure at any magnetization (bias field) as shown in
Fig. 8.

The simulation results agree semi-quantitatively with the experimental data in
Fig. 2(a) with a linear relationship between applied magnetic field and
magnetization $M=\mu_{0}B/N_{z}$ mentioned earlier, including a dielectric
resonance mode with steep magnetic field dependence that is visible in the
experimental data.

Relating to the anisotropy mentioned earlier in section III, the simulation is
treating the YIG cylinder as completely isotropic, leading to degenerate Mode
x and y at 0 field. To account for the anisotropy, we introduce a general
energetic preference along the $x$ axis, thus making the domains of the
unsaturated YIG preferentially align along the $x$ axis and breaking the
rotational symmetry.

When all domains are oriented along the $z$ axis with net magnetization of
zero, the permeability is calculated to be

$[\mu]_{z}=\begin{pmatrix}\mu_{\mathrm{eff}}&0&0\\\ 0&\mu_{\mathrm{eff}}&0\\\
0&0&1\\\ \end{pmatrix}$ (16)

where
$\mu_{\mathrm{eff}}=\sqrt{\frac{\omega^{2}-\omega_{m}^{2}}{\omega^{2}}}$. To
get the permeability matrix for domains align along the $x$ and $y$ axes
($[\mu]_{x}$, $[\mu]_{y}$), one can apply a change of coordinates to Eq. (A7).
The matrix for completely random domain orientations would be an equal average
of the three permeability matrices, $[\mu]_{x},[\mu]_{y},[\mu]_{z}$ Schlömann
(1970). Applying a weighted average to the matrices will then allow for
representation of an energetic preference, as shown for a preference along the
$x$ axis:

$[\mu]=(\frac{1}{3}+\delta)[\mu]_{x}+(\frac{1}{3}-\delta)[\mu]_{y}+\frac{1}{3}[\mu]_{z}$
(17)

Using $\delta=0.1$ in Eq. (A8) gives 260 MHz of splitting between Mode x and
Mode y, which is in good agreement with experimental results from Fig. 2(b).

##  Appendix B Modeling YIG anisotropy in system Hamiltonian

As mentioned earlier, there is a clear broken rotational symmetry in the
$x$-$y$ plane apparent from the splitting in Fig. 2(b). Since the exact origin
of the anisotropy is unknown, we will treat it as a general energetic favoring
in the x-y plane. As the $\hat{z}$ bias field is increased, the magnetization
will align more along $\hat{z}$, making $x$-$y$ plane preferences less
impactful. Based off this understanding, we wanted a simple functional form to
describe how the effect of this anisotropy decays with an increase in bias
field strength that we could use to describe the decay of $\beta$. Since we
just want the general form of how the effect of an energetic preference decays
over field, the actual form of the energetic preference in the $x$-$y$ plane
is not important. We chose to use a toy model of a magnetic domain with a
simple Hamiltonian ($H_{\mathrm{an}}$) with a simple energetic preference
given by $K$ along the $x$ axis and a total net magnetic moment $M$:

$H_{\mathrm{an}}=-BM\cos(\theta)-K\sin^{2}(\theta)\cos^{2}(\phi)$ (18)

To see how the effect of this anisotropy changes as we vary the magnetic field
$B$, we utilized classical Boltzmann statistics. We define a partition
function:

$Z=\int_{0}^{\pi}\int_{0}^{2\pi}e^{-H_{\mathrm{an}}(\theta,\phi)/(k_{b}T)}\,\sin(\theta)d\theta
d\phi\ $ (19)

so we can calculate the expectation of the magnetic moment direction using Eq.
(B3) for $A=M_{x},M_{y},M_{z}$;
$M_{x}=M\sin(\theta)\cos(\phi),M_{y}=M\sin(\theta)\sin(\phi),M_{z}=M\cos(\theta)$.

$\langle
A^{2}\rangle=\int_{0}^{\pi}\int_{0}^{2\pi}A^{2}e^{-H_{\mathrm{an}}(\theta,\phi)/(k_{b}T)}\,\sin(\theta)d\theta
d\phi\ /Z$ (20)

We calculate these expectation values numerically, and find that the
difference of $\langle M_{x}^{2}\rangle-\langle M_{y}^{2}\rangle$ follows
approximately a $\sech(B)$ function. This motivates us to use this simple
functional form to model the anisotropy-induced term $\beta$:

$\displaystyle\beta(B)$ $\displaystyle=\beta_{0}\sech(B/B_{0})$ (21)

The scaling factor $B_{0}$ was fit to the $S_{31}$ spectrum giving a value of
18.5 mT. While this is a rather crude phenomenological treatment of the
anisotropy, since the detuning of the circulator modes becomes large enough
that there is little hybridization with the cavities at relatively small
magnetic fields ($\sim$20 mT), the exact dependence on magnetic field becomes
less important to understand the non reciprocal dynamics of the cavities.

##  Appendix C Two-mode Hamiltonian and gauge symmetry

We aim to elucidate the non-reciprocity from the Hamiltonian given in Eq. (4)
by reducing it to the form written in Eq. (9). In order to do this, we
adiabatically integrate out the two circulator modes to reduce the Hamiltonian
to a simple $2\times 2$ matrix ($H^{\prime}_{\mathrm{eff}})$ involving only
the two cavity modes. The adiabatic elimination is justified due to the large
loss rate on the hybridized circulator modes, making their relevant time
scales much faster than the time scale set by the coupling parameters to the
cavities. The form of the of the effective Hamiltonian is written out in Eq.
(9). Due to the complicated dependence on the four mode model parameters, we
have written simple frequency and loss terms on the diagonal entries and
simple non-reciprocal couplings on the off diagonal entries where their
explicit values change as a function of magnetic field. The coupling terms
($H_{12},H_{21}$) along with $r=\sqrt{\lvert H_{21}\rvert/\lvert
H_{12}\rvert}$ are plotted in Fig. 7. The non-reciprocal nature of the system
then becomes immediately apparent as the $H_{21}$ and $H_{12}$ Hamiltonian
terms are different outside of 0 field, showing a clear directionality in the
interaction. We can map this Hamiltonian to a reciprocal one using the
similarity transformation outlined in Eq. (C1) with the transformation matrix
written in Eq. (C2). This new reciprocal Hamiltonain is now symmetric under
flipping the sign of the magnetic field
$H_{\mathrm{rec}}(B)=H_{\mathrm{rec}}(-B)$.

$H^{\prime}_{\mathrm{eff}}\rightarrow SH^{\prime}_{\mathrm{eff}}S^{-1}\equiv
H_{\mathrm{rec}}$ (22) $S=\begin{pmatrix}r^{1/2}&0\\\ 0&r^{-1/2}\\\
\end{pmatrix}$ (23)

This means that plotting the ratio of $R_{i,n,\mathrm{rec}}$ from
$H_{\mathrm{rec}}$ will always yield 1 for all B values. One can also map the
eigenvectors of the original system ($\ket{\psi_{i}}$) to the reciprocal
system ($\ket{\psi_{i,\mathrm{rec}}}$) by $\ket{\psi_{Ri,\mathrm{rec}}}$ =
$S\ket{\psi_{Ri}},\ket{\psi_{Li,\mathrm{rec}}}$ = $S^{-1}\ket{\psi_{Li}}$.
Starting from the ratio $R_{i,n,\mathrm{rec}}=1$ using the eigenvectors of
$H_{\mathrm{rec}}$, it is then apparent that transforming the eigenvectors
back to those of $H^{\prime}_{\mathrm{eff}}$ will allow one to simply
caluclate $R_{i,n}$. To illustrate this, we start with the explicit change in
components from the transformation of the right and left eigenvectors as
written in Eqns. (C3, C4), we can then substitute these in to the earlier
expression for the ratio $R_{i,n}$ and see how the ratio deviates from the
reciprocal case of 1, as done in Eq. (C5) with $i=1$ as an example.

$\displaystyle\ket{\psi_{R,\mathrm{rec}}}=\begin{pmatrix}x\\\ y\\\
\end{pmatrix}\xrightarrow[\text{transform}]{\text{similarity}}\ket{\psi_{R}}=\frac{1}{\sqrt{\frac{|x|^{2}}{r}+|y|^{2}r}}\begin{pmatrix}\frac{1}{\sqrt{r}}x\\\
\sqrt{r}y\\\ \end{pmatrix}$ (24)
$\ket{\psi_{L,\mathrm{rec}}}=\begin{pmatrix}x^{*}\\\ y^{*}\\\
\end{pmatrix}\xrightarrow[\text{transform}]{\text{similarity}}\ket{\psi_{L}}=\frac{1}{\sqrt{|x|^{2}r+\frac{|y|^{2}}{r}}}\begin{pmatrix}\sqrt{r}x^{*}\\\
\frac{1}{\sqrt{r}}y^{*}\\\ \end{pmatrix}$ (25)
$R_{1,n}=\frac{|\bra{n_{L}}\ket{1}|}{|\bra{1}\ket{n_{R}}|}=\frac{|r^{1/2}x\sqrt{|x|^{2}r^{-1}+|y|^{2}r|}}{{|r^{-1/2}x\sqrt{|x|^{2}r+|y|^{2}r^{-1}}|}}$
(26)

It is important to note two simplifying limits for Eq. (C5) that the reader
may verify themselves, for $x/y\gg r$, $R_{1,n}\approx 1$ and for $y/x\gg r$,
$R_{1,n}\approx r^{2}$.

One can use this similarity transformation to understand the qualitative
behavior of the disparate amplitude ratios seen in Fig. 6(f). As mentioned
earlier, the amplitude ratio in this case can be roughly approximated as
$R_{1,n}$ at $|B|\gtrsim$ 15 mT so we can focus primarily on this ratio to
understand the behavior in this field range. At larger fields ($B\gtrsim 20$
mT) Modes a and b are dominated by participation in the bare cavity modes so
we can approximate these modes by using the eigenmode values from
$H^{\prime}_{\mathrm{eff}}$ for the cavity mode components and zeros for the
circulator mode components. Under this approximation we can look at the inner
products in $R_{1,n}$ just from the components in the eigenmodes of
$H^{\prime}_{\mathrm{eff}}$. Mode a is largely dominated by the cavity 1
component with little circulator participation with
$|\bra{a_{\mathrm{rec}}}\ket{1}|/|\bra{a_{\mathrm{rec}}}\ket{2}|\gg r$ for all
values of $r$, thus we can invoke the limit of Eq.(C5) previously mentioned to
find the ratio $R_{1,b}\approx 1$ which is what is seen in Fig. 6(f). The same
argument can be made for mode b, but in the opposite limit of
$|\bra{b_{\mathrm{rec}}}\ket{2}|/|\bra{b_{\mathrm{rec}}}\ket{1}|\gg r$,
leading to the other limit of Eq.(C5), making the ratio $R_{1,b}\approx r^{2}$
which can be seen by comparing Fig. 6(f) with Fig. 7.

## References

  * Kord  _et al._ (2020) A. Kord, D. L. Sounas, and A. Alù, Proceedings of the IEEE 108, 1728 (2020).
  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).
  * Krantz  _et al._ (2019) P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Applied Physics Reviews 6, 021318 (2019).
  * Kurpiers  _et al._ (2018) P. Kurpiers, P. Magnard, T. Walter, B. Royer, M. Pechal, J. Heinsoo, Y. Salathé, A. Akin, S. Storz, J.-C. Besse, S. Gasparinetti, A. Blais, and A. Wallraff, Nature 558, 264 (2018).
  * Axline  _et al._ (2018) C. J. Axline, L. D. Burkhart, W. Pfaff, M. Zhang, K. Chou, P. Campagne-Ibarcq, P. Reinhold, L. Frunzio, S. M. Girvin, L. Jiang, M. H. Devoret, and R. J. Schoelkopf, Nature Physics 14, 705 (2018).
  * Vool and Devoret (2017) U. Vool and M. Devoret, International Journal of Circuit Theory and Applications 45, 897 (2017).
  * Lachance-Quirion  _et al._ (2019) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Applied Physics Express 12, 070101 (2019).
  * Awschalom  _et al._ (2021) D. D. Awschalom, C. H. R. Du, R. He, F. J. Heremans, A. Hoffmann, J. T. Hou, H. Kurebayashi, Y. Li, L. Liu, V. Novosad, J. Sklenar, S. E. Sullivan, D. Sun, H. Tang, V. Tiberkevich, C. Trevillian, A. W. Tsen, L. R. Weiss, W. Zhang, X. Zhang, L. Zhao, and C. W. Zollitsch, arXiv:2102.03222 (2021).
  * Huebl  _et al._ (2013) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, Physical Review Letters 111, 127003 (2013).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Physical Review Letters 113, 083603 (2014).
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Physical Review Letters 113, 156401 (2014).
  * Tabuchi  _et al._ (2015) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Science 349, 405 (2015).
  * Lachance-Quirion  _et al._ (2020) D. Lachance-Quirion, S. P. Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura, Science 367, 425 (2020).
  * Hou and Liu (2019) J. T. Hou and L. Liu, Physical Review Letters 123, 107702 (2019).
  * Golovchanskiy  _et al._ (2021) I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, M. Weides, V. V. Ryazanov, A. A. Golubov, A. V. Ustinov, and M. Y. Kupriyanov, Science Advances 7, eabe8638 (2021).
  * Heyroth  _et al._ (2019) F. Heyroth, C. Hauser, P. Trempler, P. Geyer, F. Syrowatka, R. Dreyer, S. G. Ebbinghaus, G. Woltersdorf, and G. Schmidt, Physical Review Applied 12, 054031 (2019).
  * Boventer  _et al._ (2018) I. Boventer, M. Pfirrmann, J. Krause, Y. Schön, M. Klaui, and M. Weides, Physical Review B 97, 184420 (2018).
  * Zhang  _et al._ (2017) D. Zhang, X.-Q. Luo, Y.-P. Wang, T.-F. Li, and J. Q. You, Nature Communications 8, 1368 (2017).
  * Anderson  _et al._ (2016) B. M. Anderson, R. Ma, C. Owens, D. I. Schuster, and J. Simon, Physical Review X 6, 041043 (2016).
  * Owens  _et al._ (2018) C. Owens, A. LaChapelle, B. Saxberg, B. M. Anderson, R. Ma, J. Simon, and D. I. Schuster, Physical Review A 97, 013818 (2018).
  * Zhang  _et al._ (2020) X. Zhang, A. Galda, X. Han, D. Jin, and V. M. Vinokur, Physical Review Applied 13, 044039 (2020).
  * Fay and Comstock (1965) C. E. Fay and R. L. Comstock, IEEE Transactions on Microwave Theory and Techniques 13, 15 (1965).
  * Heiss (2012) W. D. Heiss,  45, 444016 (2012).
  * Özdemir  _et al._ (2019) S. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nature Materials 18, 783 (2019).
  * Hatano and Nelson (1997) N. Hatano and D. R. Nelson, Physical Review B 56, 8651 (1997).
  * Yao and Wang (2018) S. Yao and Z. Wang, Physical Review Letters 121, 086803 (2018).
  * McDonald  _et al._ (2018) A. McDonald, T. Pereg-Barnea, and A. A. Clerk, Physical Review X 8, 041031 (2018).
  * Stannigel  _et al._ (2012) K. Stannigel, P. Rabl, and P. Zoller, New Journal of Physics 14, 063014 (2012).
  * Lodahl  _et al._ (2017) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, Nature 541, 473 (2017).
  * Chapman  _et al._ (2017) B. J. Chapman, E. I. Rosenthal, J. Kerckhoff, B. A. Moores, L. R. Vale, J. A. B. Mates, G. C. Hilton, K. Lalumière, A. Blais, and K. W. Lehnert, Physical Review X 7, 041043 (2017).
  * Lecocq  _et al._ (2017) F. Lecocq, L. Ranzani, G. A. Peterson, K. Cicak, R. W. Simmonds, J. D. Teufel, and J. Aumentado, Physical Review Applied 7, 024028 (2017).
  * Sliwa  _et al._ (2015) K. M. Sliwa, M. Hatridge, A. Narla, S. Shankar, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Physical Review X 5, 041020 (2015).
  * Ruesink  _et al._ (2016) F. Ruesink, M.-A. Miri, A. Alù, and E. Verhagen, Nature Communications 7, 13662 (2016).
  * Fang  _et al._ (2017) K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Nature Physics 13, 465 (2017).
  * Wang  _et al._ (2019) Y.-P. Wang, J. W. Rao, Y. Yang, P.-C. Xu, Y. S. Gui, B. M. Yao, J. Q. You, and C.-M. Hu, Physical Review Letters 123, 127202 (2019).
  * Xu  _et al._ (2019) H. Xu, L. Jiang, A. A. Clerk, and J. G. E. Harris, Nature 568, 65 (2019).
  * Chen  _et al._ (1991) D.-X. Chen, J. Brug, and R. Goldfarb, IEEE Transactions on Magnetics 27, 3601 (1991).
  * Fletcher  _et al._ (1960) R. C. Fletcher, R. C. LeCraw, and E. G. Spencer, Physical Review 117, 955 (1960).
  * Klingler  _et al._ (2017) S. Klingler, H. Maier-Flaig, C. Dubs, O. Surzhenko, R. Gross, H. Huebl, S. T. B. Goennenwein, and M. Weiler, Applied Physics Letters 110, 092409 (2017).
  * Marković  _et al._ (2018) D. Marković, S. Jezouin, Q. Ficheux, S. Fedortchenko, S. Felicetti, T. Coudreau, P. Milman, Z. Leghtas, and B. Huard, Physical Review Letters 121, 040505 (2018).
  * Solt (1962) I. H. Solt, Journal of Applied Physics 33, 1189 (1962).
  * Schlömann (1970) E. Schlömann, Journal of Applied Physics 41, 204 (1970).
  * Green and Sandy (1974) J. Green and F. Sandy, IEEE Transactions on Microwave Theory and Techniques 22, 641 (1974).
  * Casimir (1945) H. B. G. Casimir, Reviews of Modern Physics 17, 343 (1945).
  * Ranzani  _et al._ (2013) L. Ranzani, L. Spietz, Z. Popovic, and J. Aumentado, Review of Scientific Instruments 84, 034704 (2013).
  * Helszajn (2008) J. Helszajn,  _The Stripline Circulator: Theory and Practice_ , 1st ed. (Wiley-IEEE Press, Hoboken, NJ, 2008).
  * Metelmann and Clerk (2015) A. Metelmann and A. A. Clerk, Physical Review X 5, 021025 (2015).
  * McDonald and Clerk (2020) A. McDonald and A. A. Clerk, Nature Communications 11, 5382 (2020).
  * Schomerus (2020) H. Schomerus, Physical Review Research 2, 013058 (2020).
  * Sounas and Alù (2017) D. L. Sounas and A. Alù, Nature Photonics 11, 774 (2017).
  * Ramos  _et al._ (2014) T. Ramos, H. Pichler, A. J. Daley, and P. Zoller, Physical Review Letters 113, 237203 (2014).
  * Guimond  _et al._ (2020) P.-O. Guimond, B. Vermersch, M. L. Juan, A. Sharafiev, G. Kirchmair, and P. Zoller, npj Quantum Information 6, 1 (2020).
  * Gheeraert  _et al._ (2020) N. Gheeraert, S. Kono, and Y. Nakamura, Physical Review A 102, 053720 (2020).
  * Rymarz  _et al._ (2021) M. Rymarz, S. Bosco, A. Ciani, and D. P. DiVincenzo, Physical Review X 11, 011032 (2021).