Further author information:
Vinicius G. Goecks, Graduate Research Assistant, Vehicle Systems $\&$ Control
Laboratory<EMAIL_ADDRESS>
Grayson Woods, Graduate Research Assistant, Unmanned Systems Laboratory,
<EMAIL_ADDRESS>
Niladri Das, Graduate Research Assistant, Intelligent Systems Research
Laboratory<EMAIL_ADDRESS>
John Valasek, Professor and Director, Vehicle Systems $\&$ Control Laboratory,
<EMAIL_ADDRESS>

# Combining Visible and Infrared Spectrum Imagery using Machine Learning for
Small Unmanned Aerial System Detection

Vinicius G. Goecks  Department of Aerospace Engineering, Texas A&M University
Grayson Woods  Department of Mechanical Engineering, Texas A&M University
John Valasek  Department of Aerospace Engineering, Texas A&M University

###### Abstract

There is an increasing demand for technology and solutions to counter
commercial, off-the-shelf small unmanned aerial systems (sUAS). Advances in
machine learning and deep neural networks for object detection, coupled with
lower cost and power requirements of cameras, led to promising vision-based
solutions for sUAS detection. However, solely relying on the visible spectrum
has previously led to reliability issues in low contrast scenarios such as
sUAS flying below the treeline and against bright sources of light.
Alternatively, due to the relatively high heat signatures emitted from sUAS
during flight, a long-wave infrared (LWIR) sensor is able to produce images
that clearly contrast the sUAS from its background. However, compared to
widely available visible spectrum sensors, LWIR sensors have lower resolution
and may produce more false positives when exposed to birds or other heat
sources. This research work proposes combining the advantages of the LWIR and
visible spectrum sensors using machine learning for vision-based detection of
sUAS. Utilizing the heightened background contrast from the LWIR sensor
combined and synchronized with the relatively increased resolution of the
visible spectrum sensor, a deep learning model was trained to detect the sUAS
through previously difficult environments. More specifically, the approach
demonstrated effective detection of multiple sUAS flying above and below the
treeline, in the presence of heat sources, and glare from the sun. Our
approach achieved a detection rate of 71.2 $\pm$ 8.3%, improving by 69% when
compared to LWIR and by 30.4% when visible spectrum alone, and achieved false
alarm rate of 2.7 $\pm$ 2.6%, decreasing by 74.1% and by 47.1% when compared
to LWIR and visible spectrum alone, respectively, on average, for single and
multiple drone scenarios, controlled for the same confidence metric of the
machine learning object detector of at least 50%. With a network of these
small and affordable sensors, one can accurately estimate the 3D position of
the sUAS, which could then be used for elimination or further localization
from more narrow sensors, like a fire-control radar (FCR). Videos of the
solution’s performance can be seen at https://sites.google.com/view/tamudrone-
spie2020/.

###### keywords:

Automatic Target Recognition, ATR, Small Unmanned Aerial Systems, LWIR, Deep
Learning, Machine learning.

##  1 INTRODUCTION

Small unmanned aircraft systems (sUAS), commonly called drones, were once only
thought of as military aircraft. We have seen images from the government or
the news media showing comparatively larger aircraft such as the General
Atomics MQ-9 Reaper, which are being controlled by an operator sitting
hundreds of miles away. They are effective in performing surveillance but can
also deploy weapons without risking its operator. Increasingly, the technology
to control and operate such unmanned aircraft is becoming affordable and
easily available. We are seeing a rise in the manufacturing of consumer-grade
and commercial-grade sUAS. Companies such as DJI[1], Parrot[2], and Autel
Robotics[3] have drones equipped with high definition RGB and infrared
cameras, sophisticated mission planners, and are capable of completing its
mission autonomously, that can be bought for as little as $\$$1000 USD. This
capability is posing crucial privacy and security threats to people and
organizations around the world, specially when these drones are retrofitted
with sophisticated cameras and possibly weapons, which sparked a need for
counter-sUAS, or counter-drone, solutions.

One of the components of a counter-drone solution is the detection of the
drone. Due to the wide availability of cameras and advancements in deep
learning and computer vision techniques, vision-based solutions became more
appealing to address the detection problem. Vision based-object detection
techniques were mainly developed to detect objects on the visible spectrum
while using RGB cameras. However, accuracy strongly decreases when detecting
objects under fog, rain, or in poorly illuminated environments. As an
alternative, long-wave infrared (LWIR) cameras are more robust to these
conditions but are less accurate due to lower resolution and false positives
when detecting objects with similar heat signatures such as differentiating at
long distance drones from birds or light bulbs.

To leverage the advantages of each sensor, our proposed solution combines
images captured from RGB and LWIR cameras to train a machine learning model
for the detection of sUAS. Our main contribution is the insight of utilizing
the heightened background contrast from the LWIR sensor combined and
synchronized with the relatively increased resolution of the RGB sensor to
train a deep learning model to detect sUAS through previously difficult
environments. More specifically, the approach demonstrated effective detection
of multiple sUAS flying above and below the treeline, in the presence of heat
sources, and glare from the sun.

##  2 Related Work

For large unmanned aerial systems (UAS) detection, the majority of techniques
use radar imaging [4, 5, 6, 7, 8, 9]. However, the size and speed of small
unmanned aerial systems (sUAS) present challenges for radar detection, which
must be able to scan a large area with discrimination to prevent false alarms
from birds and moving ground objects [10]. Deep learning techniques have been
developed to help detect and classify sUAS by recognizing patterns in the
Doppler signature of micro UAS [11]. However, the confidence with which these
systems can be practically employed for counter-sUAS is still up for debate.

Image-based sUAS detection provides another solution for the detection of
sUAS. Visual sensors have been used to track targets moving over complex
backgrounds [12]. In addition, temporal information from subsequent images can
be supplemented to make the tracking robust to when the background or target
object is changing in appearance [13]. Detection of sUAS with visual sensors
can even be employed when the sensors themselves are moving (for example, when
onboard another sUAS [14]). Image processing techniques such as feature
detection and corner detection can also be used to analyze an object or scene
from two different positions and find features that would correspond to a
real-world 3D object or scene [15], offering an increased ability to classify
the target object.

Fusing multi-modal sensor outputs such as RGB and infrared video streams are
gaining prominence to detect objects under changing environmental conditions.
If one sensing modality is impaired in a particular condition, other sensing
modalities might not be affected. The fusion of sensor data for object
detection is typically done in any of these three phases: at the pixel level,
at the feature level, or at the decision level. RGB and infrared sensor data
fused at a pixel level [16] can then be segmented using either parametric
[17]techniques such as Gaussian Mixture Model or non-parametric techniques
such as kernel-based density estimators. A new non-parametric technique based
on codebook model [18] has shown enhanced object detection performance [16]
when applied to pixel-level fused data. The image segmentation accuracy as
measured using the detection rate metric is considerably higher for the fused
data compared to RGB and infrared data separately[16]. Moreover the processing
time for fused data was found to be less compared to the summation of the
processing time required for RGB and infrared data individually. The detection
accuracy is shown to be further improved if we allow a context-aware, quality-
based fusion [19]. Context sensitive indicators are introduced [19] to weigh
in the soft segmentation performed over individual RGB and infrared data.

More recently, machine learning approaches were developed to remove the
necessity of hand-designing image features and corners to, instead, learn
these patterns automatically from large datasets. Simple Online and Realtime
Tracking with a Deep Association Metric (Deep SORT) [20, 21] trained a
convolutional neural network on a large-scale person identification dataset to
learn an association metric for each target identified. Paired with Kalman
filtering in the image-space and frame-by-frame data association, this
approach achieves favorable performance when tracking objects for longer
periods of time under occlusion. One of the most common approaches, You Only
Look Once (YOLO) [22, 23, 24], frames object detection as a regression problem
where a neural network predicts, in real time, bounding boxes and class
probabilities directly from full images in only one evaluation step. Improving
upon YOLO, EfficientNet [25] proposes the compound scaling method that
uniformly scales depth, width, and resolution of convolution neural networks
following a fixed ratio to achieve state-of-the-art object detection while
using smaller and faster networks. Combined with machine learning, data from
visual sensors can be used to not only track sUAS [26] but also distinguish
them from birds [13] and other background objects. Our approach, detailed
below, utilizes many of these advancements to provide a robust framework for
sUAS detection.

##  3 Methods

This research work proposes combining the advantages of the LWIR and visible
spectrum sensors using machine learning for vision-based detection of sUAS.
Due to the power and battery discharge rates required to fly small unmanned
quadrotors, these vehicles hold a characteristic heat signature when imaged
with LWIR sensors that contrasts with common cluttered backgrounds suck as
buildings, trees, sky, and the terrain. The same does not happen in the
visible spectrum where changes in lighting and cluttered background can
completely camouflage a flying small unmanned quadrotor. Conversely, visible
spectrum sensors are widely available in the market and have higher
resolution, which aids detection at long-ranges assuming that the vehicle
contrasts with the background. Our prototype, as seen in Figure 1a consists of
two calibrated cameras aligned based on the position of their image capture
sensor and attached to a custom mount system that allows them to be attached
to a camera tripod. The cameras, from left to right, are a FLIR infrared
camera BOSON 640 coupled with 4.9mm lens and a FLIR RGB camera CM3-U3-50S5C-CS
Chameleon3 USB3 coupled with a 3.5mm C Series fixed focal length lens. The
camera specifications are shown in Table 1. All camera data was streamed and
recorded using ROS (Robot Operating System) [27], synchronized to run at $30$
Hz. A sample RGB and respective LWIR frame demonstrating the data capture
system can be seen in Figure 1b.

(a)

(b)

Figure 1: (a) Picture of the mounted RGB and LWIR cameras used to capture data for this research work and (b) sample RGB and its respective LWIR frame after synchronization and before fusion. Table 1: Camera Specifications. | Camera Type
---|---
Specifications | RGB | Infrared
Model | FLIR CM3-U3-50S5C-CS | FLIR BOSON 640
Lens focal length (mm) | 3.5 | 4.9
Frame Rate (Hz) | $35$ | $60$
Resolution (pixels) | $2448\times 2048$ | $640\times 512$
Sensor Type | CMOS | N/A
Sensor Format | $2/3"$ | N/A
Thermal Sensitivity (mK) | N/A | $<60$

The data collection process consisted of imaging two different sUAS, DJI Mavic
Pro and DJI Phantom 4, while flying at different times of the day and weather
conditions, such as during sunset and sunrise in a clear sky or during the
afternoon of a cloudy day, and under different lighting conditions in order to
create a diverse dataset to be used as training data for the machine learning
model. A total of 30 minutes of data at 30 frames per second were recorded for
a total of 52800 frames for each camera, LWIR and RGB. From that, 1275 frames
in total were manually labeled by drawing a bounding box surrounding the
vehicle. The main assumption during this process was that the person should
label LWIR and RGB data independently so there is no information leakage
between one spectrum to another.

To combine the calibrated and synchronized LWIR and RGB images, we first
computed the homography transformation from the LWIR to the RGB image frame so
each pixel in the infrared spectrum could be associated with one pixel in the
visible spectrum. The pixels were combined by computing the equally weighted
sum between them so the resulting pixel carries information from both spectra.
Figure 2 illustrates our proposed approach to combined LWIR+RGB data in
comparison to LWIR and RGB alone. The green box indicates that the machine
learning object detector has reached the desired level of confidence to
identify a drone in that location, to be explained in more detail in the
following paragraphs.

Figure 2: Sample of frames illustrating our proposed approach (LWIR+RGB,
leftmost image) compared to LWIR and RGB alone to detect sUAS. The proposed
method synchronizes and combines both LWIR and RGB images leveraging the
advantages of both sensors for object detection.

In possession of a synchronized and annotated dataset of LWIR+RGB, LWIR, and
RGB images, we trained a machine learning model for drone detection for each
separate dataset using YOLOv3 [22, 23, 24] official implementation111YOLOv3 on
Github: https://github.com/pjreddie/darknet.. The LWIR and RGB models used
their respective images and labels to train the YOLO object detector while the
LWIR+RGB model used the combined frames, as seen in the leftmost image of
Figure 2, and all LWIR and RGB labels. As illustrated in Figure 3, YOLO splits
the input image in $S\times S$ grid cells, predicting one object and a fixed
number of bounding boxes per cell using dimension clusters as anchor boxes
[28]. Each prediction of a bounding box is accompanied by a box confidence
score and each grid cell by a conditional class probability, which, in our
case, represents the probability of containing or not a drone. These scores
are condensed with a dot product in a final class confidence score that
measures the confidence in both the classification (correct class of
containing or not a drone) and the localization of the object (coordinates of
the corners of the bounding box surrounding the drone). YOLO was initially
designed to be used with RGB images but no modification in the architecture
was required since the LWIR+RGB frames still have 3 image channels, similar to
an RGB image. The same applies to the LWIR frames, displayed in colors, which
also have 3 image channels. Due to the modularity of this approach, it can be
integrated into existing surveillance and camera systems capable of recording
LWIR and RGB data.

Figure 3: Illustration of the YOLO method for object detection. Adapted from
original YOLO diagram [23].

##  4 Results and Discussion

The proposed approach is evaluated in four drone detection scenarios: drones
crossing above and below the horizon and treeline, drones flying in front of
cameras pointed directly the sun during sunrise, drones flying together with
other heat sources such as public street lighting, and drone detection at
approximately 300 meters away, here called long-range since the sUAS is barely
visible on the RGB spectrum and also to differentiate of the previous cases
where the vehicles were flying at less than 300 meters away. We compared our
approach of combining LWIR and RGB images and detecting drones using a custom
trained YOLO model against the same YOLO architecture but trained only using
LWIR and RGB images alone. All YOLO models for the baselines and the proposed
method were trained on the same scenarios and using the same labels, only
changing the input video appropriately. Performance is computed in terms of
Detection Rate (DR) and False Alarm Rate (FAR) [18]. DR is computed as the
ratio between true positives and the sum of true positives (TP) and false
negatives (FN), as seen in Equation 1, and FAR is computed as the ratio
between false positives (FP) and the sum of true and false positives, as seen
in Equation 2. Since YOLO is a probabilistic model, it outputs a confidence
metric for each object detected in the image frame. In this work, we evaluate
the performance metrics DR and FAR for all baselines and the proposed method
for four confidence thresholds: $25$, $50$, $75$, and $90\%$. Each specific
condition is measured, explained, and illustrated bellow and qualitative
results can be seen in the complete prediction videos for each condition shown
in the project webpage222Project page:
https://sites.google.com/view/tamudrone-spie2020/..

$DR=\frac{TP}{TP+FN}$ (1)

$FAR=\frac{FP}{TP+FP}$ (2)

Our first experiment, shown in Figure 4, has a drone crossing above and below
the treeline and horizon. This scenario is problematic for solutions that
solely rely on the visible spectrum because of the drone’s low contrast
relative to the background. The infrared sensor, shown in the middle of Figure
4, has a significantly higher contrast compared to the RGB sensor, due to
salient heat emission from the drone batteries and motors. This allows the
sensor to keep track of the drone even when it crosses above and below the
skyline. In comparison, our proposed method, LWIR+RGB, presents a better
horizon contrast when compared to LWIR and RGB alone.

Figure 4: Frames illustrating our proposed approach (LWIR+RGB, leftmost image)
compared to LWIR and RGB alone to detect sUAS when it is not visible on the
RGB spectrum as, for example, when flying below the horizon in front of trees.

Figure 5 illustrates our second experiment with the drone flying in front of
cameras pointed directly at the sun during sunrise. This is another scenario
where the RGB sensor does not have enough contrast for the proper detection of
the sUAS. Using both the visible and infrared spectrum, our solution can
detect the drone with much less difficulty, while also rejecting the false
positives that arise from purely relying on the infrared spectrum.

Figure 5: Frames illustrating our proposed approach (LWIR+RGB, leftmost image)
compared to LWIR and RGB alone to detect sUAS when it is not visible on the
RGB spectrum as, for example, during sunrise when cameras are pointed to the
sun.

Our third experiment, shown in Figure 6, illustrates the solution’s robustness
to other heat sources such as public street lighting. This scenario is
problematic for solutions that solely rely on the infrared spectrum because of
the similarities of drone’s and lamp’s heat signatures at low resolution. By
combining both LWIR and RGB images, our method is still capable to leverage
the heat signature of the drone to improve detection but also filters out
light bulbs because they are clearly distinct from drones in the visible
spectrum. In this scenario, RGB alone still fails when the drones are below
the horizon and in front of the trees.

Figure 6: Frames illustrating our proposed approach (LWIR+RGB, leftmost image)
compared to LWIR and RGB alone to detect sUAS when objects with similar heat
signatures to drones, such as lamps, are present in the image.

In our final experiment, shown in Figure 7, we demonstrate our proposed method
capabilities detecting drones when the sUAS is barely visible on the RGB
spectrum, here called long-range to differentiate of the previous cases where
the vehicles were flying at less than 300 meters away. This scenario is
challenging to LWIR alone due to the reduced resolution and still challenging
for RGB alone depending on how the drone contrasts with the background. As
before, our proposed method still presents higher performance metrics for
being able to combine both sensors, as explained in the following paragraph.

Figure 7: Frames illustrating our proposed approach (LWIR+RGB, leftmost image)
compared to LWIR and RGB alone to detect sUAS flying beyond 300 meters from
the cameras, here defined as long-range since the sUAS is barely visible on
the RGB spectrum.

Table 2 and 3 show the numerical results in terms of mean and standard
deviation of Detection Rate (DR) and False Alarm Rate (FAR), respectively for
each table, averaging the metrics for all scenarios and also showing them
separately for cases where a single or multiple drones were flying at the same
time, for all confidence threshold (CT) conditions of the machine learning
model output. In bold we highlighted the best performing method for each
condition, showing that our method outperforms the baselines 87.5% of the
time. Improvements in the DR ranges from 12.9 to 880.4% while the reduction in
FAR ranges from -15.4 to -100.0%, depending on the scenario. Overall, setting
CT to $50\%$ presents a good balance of DR and FAR for all conditions
investigated. Qualitative results can be seen in the complete prediction
videos for each condition shown in the project webpage333Project page:
https://sites.google.com/view/tamudrone-spie2020/..

Table 2: Detection Rate (DR) of our proposed method combining LWIR and visible spectrum imagery (LWIR+RGB) compared to using LWIR and the visible spectrum (RGB) alone. Results, in terms of the mean and standard deviation (SD), were computed for different scenarios containing single and multiple drones in the field of view of the cameras. Results are also controlled for the confidence threshold (CT) used to filter the outputs of the machine learning model. The “LWIR+RGB Variation” columns quantifies the difference between the baseline and our proposed method. Condition |  CT (%) |  | LWIR+RGB
---
(Ours)
LWIR | RGB
Mean $\pm$ SD | Mean $\pm$ SD |  | LWIR+RGB
---
Variation
Mean $\pm$ SD |  | LWIR+RGB
---
Variation
Combined | 25 | 0.821$\pm$ 0.067 | 0.617 $\pm$ 0.489 | 33.0% | 0.633 $\pm$ 0.144 | 29.6%
50 | 0.712$\pm$ 0.083 | 0.421 $\pm$ 0.374 | 69.0% | 0.546 $\pm$ 0.151 | 30.4%
75 | 0.522$\pm$ 0.1 | 0.279 $\pm$ 0.249 | 86.7% | 0.402 $\pm$ 0.136 | 29.7%
90 | 0.279$\pm$ 0.161 | 0.137 $\pm$ 0.137 | 103.0% | 0.2 $\pm$ 0.131 | 39.8%
Single Drone | 25 | 0.849$\pm$ 0.076 | 0.973 $\pm$ 0.024 | -12.8% | 0.644 $\pm$ 0.2 | 31.8%
50 | 0.742$\pm$ 0.1 | 0.657 $\pm$ 0.266 | 12.9% | 0.578 $\pm$ 0.2 | 28.3%
75 | 0.57$\pm$ 0.093 | 0.43 $\pm$ 0.196 | 32.4% | 0.436 $\pm$ 0.159 | 30.7%
90 | 0.385$\pm$ 0.066 | 0.203 $\pm$ 0.146 | 89.9% | 0.276 $\pm$ 0.109 | 39.5%
Multiple Drones | 25 | 0.778$\pm$ 0.025 | 0.082 $\pm$ 0.05 | 842.6% | 0.616 $\pm$ 0.042 | 26.2%
50 | 0.667$\pm$ 0.029 | 0.068 $\pm$ 0.044 | 880.4% | 0.498 $\pm$ 0.061 | 34.0%
75 | 0.449$\pm$ 0.069 | 0.052 $\pm$ 0.036 | 751.2% | 0.351 $\pm$ 0.124 | 27.7%
90 | 0.121$\pm$ 0.104 | 0.04 $\pm$ 0.026 | 202.9% | 0.085 $\pm$ 0.035 | 41.3%
Table 3: False Alarm Rate (FAR) of our proposed method combining LWIR and visible spectrum imagery (LWIR+RGB) compared to using LWIR and the visible spectrum (RGB) alone. Results, in terms of mean and standard deviation (SD), were computed for different scenarios containing single and multiple drones in the field of view of the cameras. Results are also controlled for the confidence threshold (CT) used to filter the outputs of the machine learning model. The “LWIR+RGB Variation” columns quantifies the difference between the baseline and our proposed method. Condition |  CT (%) |  | LWIR+RGB
---
(Ours)
LWIR | RGB
Mean $\pm$ SD | Mean $\pm$ SD |  | LWIR+RGB
---
Variation
Mean $\pm$ SD |  | LWIR+RGB
---
Variation
Combined | 25 | 0.059 $\pm$ 0.043 | 0.18 $\pm$ 0.124 | -66.7% | 0.08 $\pm$ 0.047 | -25.4%
50 | 0.027 $\pm$ 0.026 | 0.107 $\pm$ 0.075 | -74.1% | 0.052 $\pm$ 0.041 | -47.1%
75 | 0.009 $\pm$ 0.009 | 0.062 $\pm$ 0.065 | -85.1% | 0.026 $\pm$ 0.03 | -64.4%
90 | 0 $\pm$ 0.001 | 0.031 $\pm$ 0.031 | -97.2% | 0.018 $\pm$ 0.025 | -95.2%
Single Drone | 25 | 0.041 $\pm$ 0.04 | 0.248 $\pm$ 0.106 | -83.2% | 0.049 $\pm$ 0.028 | -15.4%
50 | 0.018 $\pm$ 0.025 | 0.135 $\pm$ 0.077 | -86.2% | 0.023 $\pm$ 0.014 | -19.4%
75 | 0.006 $\pm$ 0.009 | 0.073 $\pm$ 0.083 | -90.9% | 0.004 $\pm$ 0.002 | 46.3%
90 | 0.001 $\pm$ 0.002 | 0.032 $\pm$ 0.04 | -95.6% | 0.001 $\pm$ 0.002 | -25.2%
Multiple Drones | 25 | 0.087 $\pm$ 0.043 | 0.077 $\pm$ 0.063 | 12.3% | 0.126 $\pm$ 0.019 | -31.2%
50 | 0.041 $\pm$ 0.029 | 0.065 $\pm$ 0.068 | -36.5% | 0.096 $\pm$ 0.013 | -57.1%
75 | 0.013 $\pm$ 0.01 | 0.047 $\pm$ 0.05 | -71.6% | 0.058 $\pm$ 0.007 | -77.3%
90 | 0 $\pm$ 0 | 0.029 $\pm$ 0.027 | -100.0% | 0.042 $\pm$ 0.026 | -100.0%

##  5 Conclusion

In this paper, we have demonstrated that using state-of-the-art deep learning
algorithms for object detection with the insight of combining long-wave
infrared (LWIR) and visible spectrum (RGB) sensors can significantly improve
the detection of sUAS in previously difficult situations as, for example, when
sUAS flies crossing above and below the treeline and horizon line, in the
presence of other heat sources, with cameras pointed directly to the sun, and
in long-range scenarios. We compared our approach of combining LWIR and RGB
images and detecting drones using a custom trained YOLO model against the same
YOLO architecture but trained only using LWIR and RGB images alone. Our
approach achieved a detection rate of 71.2 $\pm$ 8.3%, improving by 69% when
compared to LWIR and by 30.4% when visible spectrum alone, and achieved false
alarm rate of 2.7 $\pm$ 2.6%, decreasing by 74.1% and by 47.1% when compared
to LWIR and visible spectrum alone, respectively, on average, for single and
multiple drone scenarios, controlled for the same confidence metric of the
machine learning object detector of at least 50%. We believe that the further
development of this technology has serious benefits to the area of counter-
sUAS and can supplement existing technology in practical applications and
enable counter-sUAS systems to operate in a 24/7 continuous regime and prevent
malicious use of the sUAS.

##  6 Limitations and Future Work

There are still some inherent limitations to the vision-based approaches for
object detection. One of those limitations stems from scenarios where the
drone is indistinguishable from objects with similar visual characteristics
and heat signatures. Although the multi-modal approach mitigated this issue,
the system had difficulty distinguishing the drone from the birds during long-
range detection since both objects have a similar heat signatures at low
resolution. Besides using higher-resolution sensors and collecting more
training data, one avenue we are exploring to further mitigate this limitation
uses a filtering algorithm, like Deep SORT [20], to track drones and other
objects over the individual frames. From there, each object’s flight
characteristics are used as another input for detection and provide a more
educated estimate for indiscernible heat signatures. In combination with other
sensor modalities like a radar and spectrum analyzer, the solution presented
could offer an accurate and robust solution to sUAS detection.

###### Acknowledgements.

The authors would like to thank every other participant of the “A-Team”,
winning team of the A-Hack-of-the-Drones 2018, where this project started:
Edan Coben, Emily Fojtik, Garrett Jares, Humberto Ramos, Nicholas Waytowich,
Niladri Das, Sunsoo Kim, Vedang Deshpande, Venkata Tadiparthi, and Vernon
Lawhern; Texas A&M University Professors: Dr. Daniel Ragsdale, Dr. John
Hurtado, and Dr. Raktim Bhattacharya, for supporting the students in the
event; Rodney Boehm from Texas A&M University Engineering Entrepreneurship
program for supporting the work after the event; and MD5 and Army Future
Commands for sponsoring the event and following up with the team to continue
the work after the hackathon.

## References

  * [1] “DJI drones.” https://www.dji.com/ (2006).  Accessed: 2020-03-20.
  * [2] “Parrot Drones.” https://www.parrot.com/us/drones (1994).  Accessed: 2020-03-20.
  * [3] “Autel Robotics.” https://auteldrones.com/collections/drones (2015).  Accessed: 2020-03-20.
  * [4] Moffatt, D. and Mains, R., “Detection and discrimination of radar targets,” IEEE Transactions on Antennas and Propagation 23(3), 358–367 (1975).
  * [5] Bullard, B. D. and Dowdy, P. C., “Pulse doppler signature of a rotary-wing aircraft,” IEEE Aerospace and Electronic Systems Magazine 6(5), 28–30 (1991).
  * [6] Khan, R., Gamberg, B., Power, D., Walsh, J., Dawe, B., Pearson, W., and Millan, D., “Target detection and tracking with a high frequency ground wave radar,” IEEE Journal of Oceanic Engineering 19(4), 540–548 (1994).
  * [7] Khan, R. H. and Power, D., “Aircraft detection and tracking with high frequency radar,” in [Proceedings International Radar Conference ], 44–48, IEEE (1995).
  * [8] Headrick, J. and Thomason, J., “Applications of high-frequency radar,” Radio Science 33(4), 1045–1054 (1998).
  * [9] Brown, J., Woodbridge, K., Stove, A., and Watts, S., “Air target detection using airborne passive bistatic radar,” Electronics letters 46(20), 1396–1397 (2010).
  * [10] Poitevin, P., Pelletier, M., and Lamontagne, P., “Challenges in detecting uas with radar,” in [2017 International Carnahan Conference on Security Technology (ICCST) ], 1–6, IEEE (2017).
  * [11] Mendis, G. J., Randeny, T., Wei, J., and Madanayake, A., “Deep learning based doppler radar for micro uas detection and classification,” in [MILCOM 2016-2016 IEEE Military Communications Conference ], 924–929, IEEE (2016).
  * [12] Rozantsev, A., Lepetit, V., and Fua, P., “Detecting flying objects using a single moving camera,” IEEE transactions on pattern analysis and machine intelligence 39(5), 879–892 (2016).
  * [13] Aker, C. and Kalkan, S., “Using deep networks for drone detection,” in [2017 14th IEEE International Conference on Advanced Video and Signal Based Surveillance (AVSS) ], 1–6, IEEE (2017).
  * [14] Opromolla, R., Fasano, G., and Accardo, D., “A vision-based approach to uav detection and tracking in cooperative applications,” Sensors 18(10), 3391 (2018).
  * [15] Rosten, E. and Drummond, T., “Machine learning for high-speed corner detection,” in [European conference on computer vision ], 430–443, Springer (2006).
  * [16] St-Laurent, L., Maldague, X., and Prevost, D., “Combination of colour and thermal sensors for enhanced object detection,” in [2007 10th International Conference on Information Fusion ], IEEE (jul 2007).
  * [17] Stauffer, C. and Grimson, W., “Adaptive background mixture models for real-time tracking,” in [Proceedings. 1999 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No PR00149) ], IEEE Comput. Soc (1999).
  * [18] Kim, K., Chalidabhongse, T. H., Harwood, D., and Davis, L., “Real-time foreground–background segmentation using codebook model,” Real-Time Imaging 11, 172–185 (jun 2005).
  * [19] Alldieck, T., Bahnsen, C., and Moeslund, T., “Context-aware fusion of RGB and thermal imagery for traffic monitoring,” Sensors 16, 1947 (nov 2016).
  * [20] Bewley, A., Ge, Z., Ott, L., Ramos, F., and Upcroft, B., “Simple online and realtime tracking,” in [2016 IEEE International Conference on Image Processing (ICIP) ], 3464–3468 (2016).
  * [21] Wojke, N., Bewley, A., and Paulus, D., “Simple online and realtime tracking with a deep association metric,” CoRR abs/1703.07402 (2017).
  * [22] Redmon, J., “Darknet: Open source neural networks in c.” http://pjreddie.com/darknet/ (2013–2016).
  * [23] Redmon, J., Divvala, S., Girshick, R., and Farhadi, A., “You only look once: Unified, real-time object detection,” in [Proceedings of the IEEE conference on computer vision and pattern recognition ], 779–788 (2016).
  * [24] Redmon, J. and Farhadi, A., “Yolov3: An incremental improvement,” (2018).
  * [25] Tan, M. and Le, Q. V., “Efficientnet: Rethinking model scaling for convolutional neural networks,” CoRR abs/1905.11946 (2019).
  * [26] Wang, N. and Yeung, D.-Y., “Learning a deep compact image representation for visual tracking,” in [Advances in neural information processing systems ], 809–817 (2013).
  * [27] Quigley, M., Conley, K., Gerkey, B., Faust, J., Foote, T., Leibs, J., Wheeler, R., and Ng, A. Y., “Ros: an open-source robot operating system,” in [ICRA workshop on open source software ], 3(3.2), 5, Kobe, Japan (2009).
  * [28] Redmon, J. and Farhadi, A., “Yolo9000: better, faster, stronger,” in [Proceedings of the IEEE conference on computer vision and pattern recognition ], 7263–7271 (2017).


# Ultrastrong capacitive coupling of flux qubits

Mar\́mathrm{i}a Hita-Pérez  Institute of Fundamental Physics IFF-CSIC, Calle
Serrano 113b, 28006 Madrid, Spain  Gabriel Jaumà  Institute of Fundamental
Physics IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain  Manuel Pino
Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, 28006 Madrid,
Spain  Juan José Garc\́mathrm{i}a-Ripoll  Institute of Fundamental Physics
IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain

###### Abstract

A flux qubit can interact strongly when it is capacitively coupled to other
circuit elements. This interaction can be separated in two parts, one acting
on the qubit subspaces and one in which excited states mediate the
interaction. The first term dominates the interaction between the flux qubit
and an LC-resonator, leading to ultrastrong couplings of the form
$\sigma^{y}(a+a^{\dagger}),$ which complement the inductive
$\sigma^{x}i(a^{\dagger}-a)$ coupling. However, when coupling two flux qubits
capacitively, all terms need to be taken into account, leading to complex non-
stoquastic ultrastrong interaction of the $\sigma^{y}\sigma^{y}$,
$\sigma^{z}\sigma^{z}$ and $\sigma^{x}\sigma^{x}$ type. Our theory explains
all these interactions, describing them in terms of general circuit
properties—coupling capacitances, qubit gaps, inductive, Josephson and
capactive energies—, that apply to a wide variety of circuits and flux qubit
designs.

## Introduction.–

Flux qubits are superconducting loops with one or several Josephson junctions
that, when the qubit is threaded by a magnetic flux, create a frustrated
inductive energy landscape. The qubit’s low-energy space is built from quantum
superpositions of persistent current states with opposite directions [1, 2, 3,
4]. Flux qubits exhibit strong magnetic interactions and large anharmonicities
while retaining good coherence times [5, 6, 7]. These are useful properties to
implement fast qubit gates [8, 9], perform quantum annealing [10, 11, 12, 13,
14] or to simulate strongly coupled quantum systems [15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25].

A flux qubit is usually described by two observables: the tunneling between
persistent current states $\sigma^{z}$, and the qubit’s magnetic dipole moment
$\sigma^{x}.$ The former accounts for the energy splitting $\Delta$ between
current superpositions due to tunneling at the symmetry point
$H=\Delta\sigma^{z}/2$ [2]. The dipole moment appears in the inductive
coupling to microwave photons $g\sigma^{x}i(a^{\dagger}-a)$ and to other flux
qubits $J_{xx}\sigma^{x}_{1}\sigma^{x}_{2}$. An open question is how to escape
the narrow framework provided by these interactions, allowing flux qubits to
simulate non-stoquastic Hamiltonians, the ones that exhibit a sign problem
when solving them via Quantum Monte-Carlo [26, 27, 28, 29, 30] and that enable
universal adiabatic quantum computation [31, 32, 33, 34, 35].

There is still no complete satisfactory framework to implement different flux-
qubit couplings. One may obtain a ZZ interaction by converting the main qubit
junction into a SQUID, and coupling those SQUIDs to each other inductively
[1]. However, the resulting interaction is weak and decreases monotonically
with the gap, making it unsuitable for general quantum simulation and quantum
annealing [36]. A similar coupling for qubits displaying non-trivial
topological effects [37, 38, 39] could give the desired interaction [36], but
those qubits may suffer from an enhanced susceptibility to charge noise [1,
40]. The most promising approach so far is the capacitive coupling between
flux qubits. Experiments with flux qubits [41, 30] have demonstrated capacitve
interactions along more than one direction [42], but the coupling strength
seems to be limited and there is no analytical framework to understand the
range of available interactions.

Figure 1:  (a) Two superconducting circuits, described by flux variables
$\phi_{i}$, coupled by a capacitance $C_{g}.$ (b) The first circuit is a flux
qubit and the second is either an identical qubit or an LC-resonator. The flux
qubit, operated at full frustration, is described by an inductive potential
$\mathcal{U}_{q}$ with two identical wells (solid line). The relevant energy
scales are the qubit $\Delta$, the energy differences between states in the
same well $\hbar\omega_{q}$ and, for deep potentials, the splitting
$\Delta_{e}$ between excited states. The resonator has a quadratic potential
$\mathcal{U}_{r}$ (dashed) and a harmonic spectrum with equispaced energies,
separated by $\hbar\omega_{r}$. (c) Our simulations use three Josephson
junction flux qubits, with a central junction $\alpha$ times smaller, a
possible shunting capacitance $\beta$, and a magnetic flux
$\Phi\simeq\Phi_{0}/2$. The charging and Josephson energies of the large
junctions are denoted by $E_{J}$ and $E_{C}=e^{2}/(2C)$

.

In this work we introduce an analytical model for the capacitive coupling of
flux qubits to other superconducting circuit elements [cf. Fig. 1],
complemented by a non-perturbative numerical treatment. We obtain the shape
and scaling of the interaction between a qubit and a microwave resonator, and
also between two flux qubits. We provide evidence of ultrastrong qubit-
resonator interaction mediated by the electric dipole moment $\sigma^{y}$,
thereby extending the family of ultrastrong inductive couplings [15, 16, 21],
which are mediated by $\sigma^{x}$ terms. For two capacitively coupled qubits,
we explain the appearance of complex interactions along multiple directions,
YY, ZZ and XX. This is a rich landscape of spin Hamiltonians for quantum
simulation, quantum annealing and quantum computation, which exceeds the
simple picture from spectroscopic characterizations [30] and complements
earlier numerical studies [42] with scalings based on qubit’s design
parameters and a modelization of interactions mediated by excited states.

## Model.–

We study on equal footing both the qubit-qubit and qubit-resonator capacitive
couplings, as shown in Fig. 1. Each element, qubit or resonator, is
represented by a flux $\phi$ and a charge $q$ operator, and the type of
circuit is determined by the inductive potential $\mathcal{U}$. The
Hamiltonian $H=H^{(0)}+\epsilon V$ splits into bare circuits and interaction
[43]

$\displaystyle H^{(0)}$
$\displaystyle=\sum_{i=1,2}\frac{q_{i}^{2}}{2C_{i}}+\mathcal{U}_{i}(\phi_{i}),$
(1) $\displaystyle V$
$\displaystyle=-\frac{q_{1}q_{2}}{\sqrt{C_{1}C_{2}}}-\sqrt{\frac{C_{2}}{C_{1}}}\
\frac{q_{1}^{2}}{C_{1}}-\sqrt{\frac{C_{1}}{C_{2}}}\ \frac{q_{2}^{2}}{C_{2}}.$
(2)

$C_{1,2}$ are the capacitances of bare circuits, $C_{g}$ is the coupling
capacitance and
$\overline{C}_{\text{od}}=\left[C_{1}C_{2}+(C_{1}+C_{2})C_{g}\right]/C_{g}$ is
the off-diagonal of the inverse capacitance matrix. The strength of the
coupling is controlled by
$\epsilon=\sqrt{C_{1}C_{2}}/\overline{C}_{\text{od}},$ with
$\epsilon\sim\mathcal{O}(C_{g})$ for $C_{g}/C_{1,2}\ll 1.$

The first circuit is always a flux qubit $C_{1}=C_{q}$. The second circuit
will be either an identical qubit $C_{2}=C_{1}$, or a microwave resonator
$C_{2}=C_{r}$. Without loss of generality, we study three Josephson Junction
flux qubits (3JJQ) [cf. Fig. 1(c)], operating at full frustration, with half a
flux quantum $\Phi_{q}=\Phi_{0}/2$ in the loop. In this situation, this or any
other similar qubit will exhibit an inductive potential $\mathcal{U}_{1,2}$
with local minima at $\phi=\pm\varphi_{*}(\Phi_{0}/2\pi)$. Each minima is
associated to one persistent current state and a local excitation energy
$\hbar\omega_{q}$. We present results in terms of the 3JJQ’s relative coupling
$\gamma=C_{g}/C$. For the resonator we use a quadratic inductive potential
$\mathcal{U}_{2}=\phi_{2}^{2}/(2L_{r})$ with resonator frequency
$\omega_{r}=1/\sqrt{L_{r}C_{r}}.$

## Methods.–

We model the whole system as a an effective qubit-resonator or qubit-qubit
Hamiltonian, using a Schrieffer-Wolff transformation $U(\epsilon)$ [44] that
maps the eigenspaces of $H$ to the eigenspaces of the bare model $H^{(0)}$. In
our analytical treatment, we start from a projector $P_{0}$ onto a low-energy
subspace—e.g. the 4-dimensional space of two qubits, or a tensor product of a
qubit and resonator spaces—and develop the effective Hamiltonian as a
perturbative series

$H_{\text{eff}}=P_{0}H^{(0)}P_{0}+\sum_{n=1}\epsilon^{n}\mathcal{M}_{n}.$ (3)

The first order term is the projection of the interaction onto the qubit
subspace $\mathcal{M}_{1}=P_{0}^{\dagger}VP_{0}$, while $\mathcal{M}_{i\geq
2}$ describe interactions mediated by virtual transitions [43].

Numerically, we could imitate this procedure [42], but instead we sum the
series to all orders [45, 44], as $H_{\text{eff}}=P_{0}UHU^{\dagger}P_{0}$.
The unitary transformation $UP=P_{0}U$ is derived from the projector $P$ and
$P_{0}$, onto the numerically exact eigenstates of $H$ and $H^{(0)}$. We then
expand the effective model $H_{\text{eff}}$ using Pauli and Fock operators.
This allows us to compare the effective Hamiltonian to the predictions from
perturbation theory, validating the type and scaling of the coupling terms.

## First order capacitive interaction.–

The interaction $V$ has two terms that renormalize the bare circuits, and only
one that entangles their dynamics $V_{c}=-q_{1}q_{2}/\sqrt{C_{1}C_{2}}.$ To
develop the first order correction $\mathcal{M}_{1}$ we must express the
charges $q_{i}$ in the qubit and resonator basis. For the resonator
$q_{2}=\sqrt{\hbar/2Z}(a^{\dagger}+a)$ exactly, in terms of Fock operators
$\\{a,a^{\dagger}\\}$ and the resonator impedance $Z$. For the flux qubit
$q_{1}$ we assume that the renormalized Hamiltonian is anharmonic and
$H_{1}\simeq\hbar\Delta\sigma_{1}^{z}/2$. We derive the voltage operator
projected onto the qubit subspace $\mathcal{V}_{1}$ as the derivative of the
flux $\mathcal{V}_{1}\simeq
i[H_{1},\phi_{1}]/\hbar=(\Phi_{0}/2\pi)\frac{\varphi_{\star}\Delta}{\hbar}\sigma^{y}_{1},$
approximating the flux operator as the flux jump between qubit states
$\phi_{1}=(\Phi_{0}/2\pi)\varphi_{\star}\sigma^{x}.$ Since
$\mathcal{V}_{1}=q_{1}/\overline{C}_{1}$ for the bare renormalized qubit, the
projected charge operator is
$q_{1}\simeq\Phi_{0}\overline{C}_{1}\Delta\varphi_{\star}\sigma_{1}^{y}/h.$

Using this method, we obtain the first order effective interaction between two
qubits $H^{(1)}_{qq}=g^{(1)}_{qq}\sigma_{1}^{y}\sigma^{y}_{2}$ with

$\displaystyle\frac{g^{(1)}_{qq}}{\Delta}$
$\displaystyle=\frac{\overline{C}_{q}\varphi_{\star}^{2}}{\overline{C}_{\text{od}}}\frac{\Delta}{8E^{q}_{C}}.$
(4)

and the interaction between a qubit and a resonator
$H^{(1)}_{qr}=ig^{(1)}_{qr}\sigma_{1}^{y}(b^{\dagger}-b)$ with

$\displaystyle\frac{g^{(1)}_{qr}}{\Delta}$
$\displaystyle=\frac{\overline{C}_{q}}{\overline{C}_{\text{od}}}\frac{\varphi^{\star}}{2}\sqrt{\frac{1}{2\pi
G_{0}\mathcal{Z}}}.$ (5)

Everything depends on the qubit’s renormalized gap $\Delta$, the resonator
impedance $Z$, the conductance quantum $G_{0}$ and the qubit’s charging energy
$E_{C}^{q}=e^{2}/(2\overline{C}_{1})$.

This treatment neglects higher order terms in the perturbation series (3),
generated by matrix elements of the qubit’s charge operator $(\openone-
P_{0})q_{1}P_{0}$ connecting qubit states with excited states delocalized
among the inductive potential wells. Those elements grow as
$q_{1}\sim\sqrt{\omega_{q}}$, requiring us to analyze their effect in a case-
by-case basis.

Figure 2:  A 3JJQ capacitively couple to a LC-resonator. (a) Coupling divided
by qubit gap $\Delta$ as a function of the first order corrections (the dashed
line is the theoretical prediction for those corrections). The design
parameters for all the panels are qubit large junction energy ratio
$r_{q}=E_{J}/E_{C}$, energy ratio for the resonator
$r_{r}=E^{r}_{J}/E^{r}_{C}$ and ratio of resonator and qubit Josephson
energies $r_{qr}=E_{J}/E_{J}^{r}.$ The Josephson energy for the resonator as a
function of its inductance $L_{r}$ is
$E_{J}^{r}=\frac{1}{L_{r}}\left(\frac{\hbar}{2e}\right)^{2}.$ (b) Coupling
divided by qubit gap (solid line) and resonator energy (dashed) as a function
of the ratio $\gamma=C_{g}/C$ between shared and qubit’s large junction
capacitances.

## Strong qubit-resonator coupling.

Let us discuss the capacitive coupling of a 3JJQ [c.f. Fig. 1(c)] with an LC-
resonator. The second order corrections to the capacitive coupling involve a
simultaneous excitation of the qubit and the resonator, which acquire an
energy $\hbar\omega_{q}$ due to leaving the qubit space and $\hbar\omega_{r}$
due to the creation or annihilation of a photon. While the amplitude of these
processes in $\epsilon V_{c}$ grows as $\hbar\sqrt{\omega_{q}\omega_{r}}$, the
resonator cannot easily absorb the energy $\hbar\omega_{q}.$ Thus, higher-
order corrections $\mathcal{M}_{i\geq 2}$ only renormalize the qubit’s self-
energy and the capacitive coupling is fully captured by Eq. (5).

We confirm this hypothesis with the full Schrieffer-Wolff transformation of
the qubit-resonator circuit model (2). In terms of qubit and photon operators,
it takes the form $H_{\text{eff}}\simeq\hbar\Delta\sigma^{z}/2+\hbar\omega
a^{\dagger}a+g\sigma^{y}(a^{\dagger}+a)$ for moderate numbers of photons. As
shown in Fig. 2(a), all the coupling constants from different designs of qubit
and resonator collapse into the first order correction derived analytically
(5), up to coupling strengths $g/\Delta\approx 1$ beyond the perturbative
limit. Since our model only demands a qubit anharmonic spectrum, we conclude
that (5) is a general theory for the capacitive interaction between flux
qubits and microwave resonators.

We also have evidence that the capacitive coupling can enter the ultrastrong
coupling regime, where $g/\omega_{r}\sim g/\Delta\geq 12\%$. Unlike the
inductive case [15, 16, 21], exploring the actual designs where this happens
is complicated, because one has to consider the renormalization of the qubit’s
gap, while the resonator remains more or less unperturbed. Fig. 2(c) shows
evidence of this regime using a qubit with an intermediate ratio of qubit
Josephson and charging energies $r_{q}=E_{J}/E_{C}=20$ which slows down the
gap renormalization. Doing so, we achieve the ultra-strong coupling regime
$g\approx 0.15\Delta$ at zero detuning $\Delta=\omega_{r}$. We expect future
simulations will reveal more favorable situations, by tuning both the qubit’s
and the resonator’s impedances.

Figure 3:  (a) Coupling strength in gap units and (b) low energy spectrum for
two equal 3JJQ capacitively coupled, both having $\alpha=0.65$ and ratio
$E_{J}/E_{C}=50,$ as a function of the ratio $\gamma=C_{g}/C$ between shared
and qubit’s large junction capacitances. The $E_{J},E_{C}$ are Josephson and
charging of the large junction qubit. The inset of (b) shows the relative
anharmonicity of each qubit approximated from the full spectrum as
$\alpha_{r}=(E_{3}-E_{0}-\Delta)/\Delta.$

## Qubit-Qubit coupling.–

Let us now discuss the capacitive interaction between two flux qubits. We
begin by presenting the numerically exact Schrieffer-Wolff transformation for
two 3JJQ’s. In Fig. 3(a) we plot the interaction coefficients that result from
expanding $H_{\text{eff}}$ in the basis of Pauli matrices. The numerical model
clearly shows the first order terms associated to the $\sigma^{y}$ charge
dipole operator in Eq. (4). However, the flux qubits acquire also a comparable
ZZ interaction that enables the tunneling of current states, and we also find
a residual inductive XX coupling that explodes once $\gamma=C_{g}/C\gg 1$.
Note while the capacitive term induces a renormalization of the qubit’s gap
$\Delta$, the qubit nature is preserved by an improvement in the qubit’s
relative anharmonicity $\alpha_{r}$ [c.f. Fig. 3(b)]. Only at very large
$\gamma,$ the form of the interactions approaches
$-g(\tilde{\sigma}^{+}_{1}\tilde{\sigma}^{+}_{2}+\tilde{\sigma}^{-}_{1}\tilde{\sigma}^{-}_{2}),$
with ladder operators in the persistent current base
$\tilde{\sigma}^{\pm}_{i}=\tilde{\sigma}_{i}^{z}\pm i\tilde{\sigma}^{y}_{i}$.
In this limit, the coupling produces a large mass in the direction
$\phi_{1}-\phi_{2}$ of the two qubit system and strongly suppress transitions
of the form
$\tilde{\sigma}^{+}_{1}\tilde{\sigma}^{-}_{2}+\tilde{\sigma}^{-}_{1}\tilde{\sigma}^{+}_{2}$
[46].

Figure 4:  Diagrams for first (a), second (b) and third (c) order corrections
to the capacitive coupling of two 3JJQ. Qubit and excited states are
represented by solid and dashed lines, respectively. Wiggle lines are qubit
interactions, which can be diagonal (shorter blue line) and off-diagonal
(longer red) in the qubit subspace. First order corrections are the projection
of charge-charge operators to the qubit subspace. It gives a $YY$ type of
coupling. Second order produces a $ZZ$ type of coupling. We only show one of
the possible diagrams for third order corrections. It can be seen to give $XX$
coupling using the lower order diagrams. Data collapse of the (e) first and
(d) second order corrections for several qubit design parameters.

We interpret these results using the perturbation theory (3) processes
$\mathcal{M}_{1}$, $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$, sketched in Figs.
4(a-c). The horizontal lines denote qubit (solid) and excited states (dashed).
These are connected by interactions (wiggly lines), which can be qubit terms
$P_{0}VP_{0}$ (blue) or connect to excitations $P_{0}V(\openone-P_{0})$ and
$(\openone-P_{0})VP_{0}$ (red). To first order, the capacitor produces YY
terms in the qubit space. To second order, the operator $V_{c}$ enables
virtual transitions where both circuits momentarily excite, acquiring an
energy $\sim 2\times\hbar\omega_{q}$. Despite the large energy difference,
these processes are assisted by matrix elements in $V_{c}$ that grow as
$\omega_{q}$, and cannot be neglected. We estimate the second order term
$H_{qq}^{(2)}=g^{(2)}_{qq}\sigma_{1}^{z}\sigma_{2}^{z}$ with

$\displaystyle
g^{(2)}_{qq}=\frac{\overline{C}_{q}}{\overline{C}_{\text{od}}}\frac{(\Delta_{e}-\Delta)^{2}}{\hbar\omega_{q}},$
(6)

where $\Delta_{e}$ is the approximate splitting between excited states due to
tunneling [See [43] and Fig. 1a]. Finally, we can also deduce that third order
terms create an XX coupling. This coupling becomes dominant when the
renormalization of the qubit’s capacitance enables phase-slip transitions
between the qubit’s unit cells—as opposed to the tunneling between current
states enabled by $\Delta$—. However, this mechanism leads to an enhancement
of charge noise [1, 40]. We can stay far from this regime by limiting
$\gamma\sim O(1)$ and choosing $\alpha\ll 0.9$, away from the single qubit
phase slip regime [1]. This way we can still recover an XX interaction by
coupling the qubits inductively, which in combination with YY an ZZ
interactions gives a fully non-stoquastic spin model.

As in the resonator case, we compare our predictions (4) and (6) to the
numerically exact couplings, for different qubit designs. Fig. 3(d) shows how
the $g^{yy}$ coupling collapses to the theoretical prediction (4) only for
small enough couplings $g/\Delta<0.05$. Similarly, Fig. 3(e) shows that
$g^{zz}$ follows the expected scaling, up to factors $\mathcal{O}(1)$.
Perturbation theory thus captures the overall dependency of the couplings on
the qubit’s parameters, but fails to estimate the non-perturbative corrections
that account for the full interaction. This contrasts with the qubit-resonator
model, and highlights the relevance of qubit-qubit interactions mediated by
excited states.

## Conclusions.–

We have presented a non-perturbative study of capacitive interactions between
a flux qubit and other circuits. This study reveals that the flux qubit charge
operator is the sum of two dipole moments: one connecting qubit states and one
enabling transitions to higher energy excitations. The first term describes
the coupling between a qubit and a resonator and supports ultrastrong qubit-
photon interactions along directions orthogonal to the ones created by
inductive terms. The second term combines with the first one to enable a rich
family of non-stoquastic qubit-qubit couplings, including YY, ZZ and XX
interactions.

The capacitive interactions combined with inductive ultrastrong qubit-photon
interactions [15, 16, 20, 21] open new regimes in light-matter and light-
mediated interactions, such as the ultrastrong coupling limit of the Jaynes-
Cummings model [47], new regimes of the spin-boson model [20] with two
transverse couplings, and new models of coherent and dissipative interactions
mediated by photon exchange, beyond those in Ref. [48].

Regarding quantum simulation, our study confirms the idea that flux qubits
exhibit rich families of non-stoquastic interactions. These may appear
combined, as in the two-qubit model, or they may be pure YY interactions, if
we use resonators to mediate the coupling [49, 24]. From a fundamental point
of view, it would be interesting to explore the long-range interactions at the
hardware level, without embeddings [50]. Indeed, capacitively coupled flux
qubits in 2D geometries support YY interactions that only decay
logarithmically up to a length $\xi=\sqrt{C_{g}/C_{q}}$ and then exponentially
fast [51]. As in classical spin-glasses, such long-range interactions will
produce hard to solve quantum models [52, 53, 54, 55, 56].

Finally, recent works [57, 58] suggest that the ground state properties of
many superconducting circuits, including capacitively coupled flux qubits,
could be efficiently simulable in the semiclassical charge-flux
representation. Our work opens a rigorous avenue to study these circuits and
their effective models in a non-perturbative fashion. This will help us
understand whether the ground states of circuits are in some sense
trivial—e.g. low-energy states are classical spin configurations—, or whether
the energy scales and types of interactions reveal other kinds of
obstructions, different from the sign problem, that prevent the classical
simulation of the device.

###### Acknowledgements.

This work has been supported by European Commission FET-Open project AVaQus GA
899561 and CSIC Quantum Technologies Platform PTI-001. Financial support by
Fundación General CISC (Programa Comfuturo) is acknowledged. The numerical
computations have been performed in the cluster Trueno of the CSIC.

## References

  * Orlando  _et al._ [1999] T. Orlando, J. Mooij, L. Tian, C. H. Van Der Wal, L. Levitov, S. Lloyd, and J. Mazo, Superconducting persistent-current qubit, Physical Review B 60, 15398 (1999).
  * Mooij  _et al._ [1999] J. Mooij, T. Orlando, L. Levitov, L. Tian, C. H. Van der Wal, and S. Lloyd, Josephson persistent-current qubit, Science 285, 1036 (1999).
  * Chiorescu  _et al._ [2003] I. Chiorescu, Y. Nakamura, C. M. Harmans, and J. Mooij, Coherent quantum dynamics of a superconducting flux qubit, Science 299, 1869 (2003).
  * Peltonen  _et al._ [2018] J. Peltonen, P. Coumou, Z. Peng, T. Klapwijk, J. Tsai, and O. Astafiev, Hybrid rf squid qubit based on high kinetic inductance, Scientific reports 8, 1 (2018).
  * You  _et al._ [2007] J. Q. You, X. Hu, S. Ashhab, and F. Nori, Low-decoherence flux qubit, Phys. Rev. B 75, 140515 (2007).
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nature Physics 7, 565 (2011).
  * Braumüller  _et al._ [2020] J. Braumüller, L. Ding, A. P. Vepsäläinen, Y. Sung, M. Kjaergaard, T. Menke, R. Winik, D. Kim, B. M. Niedzielski, A. Melville, _et al._ , Characterizing and optimizing qubit coherence based on squid geometry, Physical Review Applied 13, 054079 (2020).
  * Liu  _et al._ [2006] Y.-x. Liu, L. F. Wei, J. S. Tsai, and F. Nori, Controllable coupling between flux qubits, Phys. Rev. Lett. 96, 067003 (2006).
  * Saito  _et al._ [2009] S. Saito, T. Tilma, S. J. Devitt, K. Nemoto, and K. Semba, Experimentally realizable controlled not gate in a flux qubit/resonator system, Physical Review B 80, 224509 (2009).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse ising model, Phys. Rev. E 58, 5355 (1998).
  * Bunyk  _et al._ [2014] P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, and et al., Architectural considerations in the design of a superconducting quantum annealing processor, IEEE Transactions on Applied Superconductivity 24, 1–10 (2014).
  * Boixo  _et al._ [2014] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nature Physics 10, 218–224 (2014).
  * Weber  _et al._ [2017] S. J. Weber, G. O. Samach, D. Hover, S. Gustavsson, D. K. Kim, A. Melville, D. Rosenberg, A. P. Sears, F. Yan, J. L. Yoder, W. D. Oliver, and A. J. Kerman, Coherent coupled qubits for quantum annealing, Phys. Rev. Applied 8, 014004 (2017).
  * Hauke  _et al._ [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Perspectives of quantum annealing: Methods and implementations, Reports on Progress in Physics 10.1088/1361-6633/ab85b8 (2020).
  * Niemczyk  _et al._ [2010] T. Niemczyk, F. Deppe, H. Huebl, E. Menzel, F. Hocke, M. Schwarz, J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, _et al._ , Circuit quantum electrodynamics in the ultrastrong-coupling regime, Nature Physics 6, 772 (2010).
  * Forn-D\́mathrm{missing}iaz  _et al._ [2010] P. Forn-D\́mathrm{i}az, J. Lisenfeld, D. Marcos, J. J. Garcia-Ripoll, E. Solano, C. Harmans, and J. Mooij, Observation of the bloch-siegert shift in a qubit-oscillator system in the ultrastrong coupling regime, Physical review letters 105, 237001 (2010).
  * Peropadre  _et al._ [2010] B. Peropadre, P. Forn-D\́mathrm{i}az, E. Solano, and J. J. Garc\́mathrm{i}a-Ripoll, Switchable ultrastrong coupling in circuit qed, Physical review letters 105, 023601 (2010).
  * Kakuyanagi  _et al._ [2016] K. Kakuyanagi, Y. Matsuzaki, C. Déprez, H. Toida, K. Semba, H. Yamaguchi, W. J. Munro, and S. Saito, Observation of collective coupling between an engineered ensemble of macroscopic artificial atoms and a superconducting resonator, Physical review letters 117, 210503 (2016).
  * Pino  _et al._ [2015] M. Pino, A. Tsvelik, and L. Ioffe, Unpaired majorana modes in josephson-junction arrays with gapless bulk excitations, Physical review letters 115, 197001 (2015).
  * Forn-D\́mathrm{missing}iaz  _et al._ [2017] P. Forn-D\́mathrm{i}az, J. J. Garc\́mathrm{i}a-Ripoll, B. Peropadre, J.-L. Orgiazzi, M. Yurtalan, R. Belyansky, C. M. Wilson, and A. Lupascu, Ultrastrong coupling of a single artificial atom to an electromagnetic continuum in the nonperturbative regime, Nature Physics 13, 39 (2017).
  * Yoshihara  _et al._ [2017] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Superconducting qubit–oscillator circuit beyond the ultrastrong-coupling regime, Nature Physics 13, 44 (2017).
  * De Bernardis  _et al._ [2018] D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Breakdown of gauge invariance in ultrastrong-coupling cavity qed, Physical Review A 98, 053819 (2018).
  * Zhu  _et al._ [2011] X. Zhu, S. Saito, A. Kemp, K. Kakuyanagi, S.-i. Karimoto, H. Nakano, W. J. Munro, Y. Tokura, M. S. Everitt, K. Nemoto, _et al._ , Coherent coupling of a superconducting flux qubit to an electron spin ensemble in diamond, Nature 478, 221 (2011).
  * Pino and Garc\́mathrm{missing}ia-Ripoll [2018] M. Pino and J. J. Garc\́mathrm{i}a-Ripoll, Quantum annealing in spin-boson model: from a perturbative to an ultrastrong mediated coupling, New Journal of Physics 20, 113027 (2018).
  * Pino and García-Ripoll [2020] M. Pino and J. J. García-Ripoll, Mediator-assisted cooling in quantum annealing, Physical Review A 101, 032324 (2020).
  * Bravyi  _et al._ [2006] S. Bravyi, D. P. Divincenzo, R. I. Oliveira, and B. M. Terhal, The complexity of stoquastic local hamiltonian problems, arXiv preprint quant-ph/0606140 (2006).
  * Susa  _et al._ [2017] Y. Susa, J. F. Jadebeck, and H. Nishimori, Relation between quantum fluctuations and the performance enhancement of quantum annealing in a nonstoquastic hamiltonian, Physical Review A 95, 042321 (2017).
  * Albash [2019] T. Albash, Role of nonstoquastic catalysts in quantum adiabatic optimization, Physical Review A 99, 042334 (2019).
  * Hormozi  _et al._ [2017] L. Hormozi, E. W. Brown, G. Carleo, and M. Troyer, Nonstoquastic hamiltonians and quantum annealing of an ising spin glass, Physical review B 95, 184416 (2017).
  * Ozfidan  _et al._ [2020] I. Ozfidan, C. Deng, A. Smirnov, T. Lanting, R. Harris, L. Swenson, J. Whittaker, F. Altomare, M. Babcock, C. Baron, _et al._ , Demonstration of a nonstoquastic hamiltonian in coupled superconducting flux qubits, Physical Review Applied 13, 034037 (2020).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Reviews of Modern Physics 90, 015002 (2018).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. Shen, M. N. Vyalyi, and M. N. Vyalyi,  _Classical and quantum computation_ , 47 (American Mathematical Soc., 2002).
  * Kempe and Regev [2003] J. Kempe and O. Regev, 3-local hamiltonian is qma-complete, arXiv preprint quant-ph/0302079 (2003).
  * Kempe  _et al._ [2006] J. Kempe, A. Kitaev, and O. Regev, The complexity of the local hamiltonian problem, SIAM Journal on Computing 35, 1070 (2006).
  * Oliveira and Terhal [2005] R. Oliveira and B. M. Terhal, The complexity of quantum spin systems on a two-dimensional square lattice, arXiv preprint quant-ph/0504050 (2005).
  * Kerman [2019] A. J. Kerman, Superconducting qubit circuit emulation of a vector spin-1/2, New Journal of Physics 21, 073030 (2019).
  * Friedman and Averin [2002] J. R. Friedman and D. V. Averin, Aharonov-casher-effect suppression of macroscopic tunneling of magnetic flux, Physical review letters 88, 050403 (2002).
  * Bell  _et al._ [2016] M. Bell, W. Zhang, L. Ioffe, and M. Gershenson, Spectroscopic evidence of the aharonov-casher effect in a cooper pair box, Physical review letters 116, 107002 (2016).
  * Kalashnikov  _et al._ [2020] K. Kalashnikov, W. T. Hsieh, W. Zhang, W.-S. Lu, P. Kamenov, A. Di Paolo, A. Blais, M. E. Gershenson, and M. Bell, Bifluxon: Fluxon-parity-protected superconducting qubit, PRX Quantum 1, 010307 (2020).
  * Chirolli and Burkard [2006] L. Chirolli and G. Burkard, Full control of qubit rotations in a voltage-biased superconducting flux qubit, Phys. Rev. B 74, 174510 (2006).
  * Satoh  _et al._ [2015] T. Satoh, Y. Matsuzaki, K. Kakuyanagi, K. Semba, H. Yamaguchi, and S. Saito, Ising interaction between capacitively-coupled superconducting flux qubits, arXiv preprint arXiv:1501.07739 (2015).
  * Consani and Warburton [2020] G. Consani and P. A. Warburton, Effective hamiltonians for interacting superconducting qubits: local basis reduction and the schrieffer–wolff transformation, New Journal of Physics 22, 053040 (2020).
  * sup [2021] Supplementary material for “ultrastrong capacitive coupling of flux qubits” (2021).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–wolff transformation for quantum many-body systems, Annals of physics 326, 2793 (2011).
  * [45] M. Hita-Pérez, G. Jaumà, M. Pino, and J. J. Garc\́mathrm{i}a-Ripoll, 3-josephson junctions flux qubit couplings, in preparation.
  * Levitov  _et al._ [2001] L. Levitov, T. Orlando, J. Majer, and J. Mooij, Quantum spin chains and majorana states in arrays of coupled qubits, arXiv preprint cond-mat/0108266 (2001).
  * Huang  _et al._ [2020] J.-F. Huang, J.-Q. Liao, and L.-M. Kuang, Ultrastrong jaynes-cummings model, Phys. Rev. A 101, 043835 (2020).
  * Gonzalez-Tudela  _et al._ [2011] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Entanglement of two qubits mediated by one-dimensional plasmonic waveguides, Phys. Rev. Lett. 106, 020501 (2011).
  * Kurcz  _et al._ [2014] A. Kurcz, A. Bermudez, and J. J. Garcia-Ripoll, Hybrid quantum magnetism in circuit qed: from spin-photon waves to many-body spectroscopy, Physical review letters 112, 180405 (2014).
  * Venturelli  _et al._ [2015] D. Venturelli, S. Mandrà, S. Knysh, B. O’Gorman, R. Biswas, and V. Smelyanskiy, Quantum optimization of fully connected spin glasses, Physical Review X 5, 10.1103/physrevx.5.031040 (2015).
  * Ortuño  _et al._ [2015] M. Ortuño, A. Somoza, V. Vinokur, and T. Baturina, Electronic transport in two-dimensional high dielectric constant nanosystems, Scientific reports 5, 1 (2015).
  * Fernandez  _et al._ [2018] L. A. Fernandez, E. Marinari, V. Martin-Mayor, G. Parisi, and J. J. Ruiz-Lorenzo, Out-of-equilibrium 2d ising spin glass: Almost, but not quite, a free-field theory, Journal of Statistical Mechanics: Theory and Experiment 2018, 103301 (2018).
  * Fernandez  _et al._ [2019] L. A. Fernandez, E. Marinari, V. Martin-Mayor, G. Parisi, and J. J. Ruiz-Lorenzo, An experiment-oriented analysis of 2d spin-glass dynamics: a twelve time-decades scaling study, Journal of Physics A: Mathematical and Theoretical 52, 224002 (2019).
  * Marshall  _et al._ [2016] J. Marshall, V. Martin-Mayor, and I. Hen, Practical engineering of hard spin-glass instances, Physical Review A 94, 10.1103/physreva.94.012320 (2016).
  * Pino  _et al._ [2012] M. Pino, A. M. Somoza, and M. Ortuño, Quantum coulomb gap in low dimensions, Physical Review B 86, 094202 (2012).
  * Pino [2020] M. Pino, Scaling up the anderson transition in random-regular graphs, Phys. Rev. Research 2, 042031 (2020).
  * Ciani and Terhal [2020] A. Ciani and B. M. Terhal, Stoquasticity in circuit-qed, arXiv preprint arXiv:2011.01109 (2020).
  * Halverson  _et al._ [2020] T. Halverson, L. Gupta, M. Goldstein, and I. Hen, Efficient simulation of so-called non-stoquastic superconducting flux circuits, arXiv preprint arXiv:2011.03831 (2020).

##  Appendix A APPENDIX 1: Hamiltonian of two flux qubits coupled via a
capacitor

We derive the Hamiltonian for a system of two flux qubits coupled via a
capacitor. We assume two identical qubits, with qubit capacitance $C_{q}$,
which are coupled by a capacitor $C_{g}$ as in Fig. 1 of the main text. The
Lagrangian in term of the capacitance matrix:

$\displaystyle\mathcal{L}=\frac{1}{2}\vec{\dot{\phi}}\ \mathcal{C}\
\vec{\dot{\phi}}-\sum_{i=1,2}\mathcal{U}_{i}$ (7)

with the vector of fluxes of each qubit $\vec{\phi}=(\phi_{1},\phi_{2})$ and
capacitance matrix:

$\displaystyle\mathcal{C}=\begin{bmatrix}C_{q}+C_{g}&-C_{g}\\\
-C_{g}&C_{q}+C_{g}\end{bmatrix}.$ (8)

The Hamiltonian $H=\sum_{i=1,2}q_{i}\dot{\phi}_{i}-\mathcal{L}$ can be written
using flux and its conjugate charge:

$\displaystyle
H=\sum_{i=1,2}\frac{q_{i}^{2}}{2\overline{C}_{q}}+\mathcal{U}(\phi_{i})+\frac{q_{1}q_{2}}{\overline{C}_{\text{od}}}.$
(9)

The inverse of capacitance matrix $\mathcal{C}^{-1}$ is used to obtain
$1/\overline{C}_{q}=(\mathcal{C}^{-1})_{11}=(\mathcal{C}^{-1})_{22}$ and
$1/\overline{C}_{\text{od}}=(\mathcal{C}^{-1})_{12}=(\mathcal{C}^{-1})_{21}$.
We can obtain the following expression for the capacitances involved in Eq.
(9) as:

$\displaystyle\frac{1}{\overline{C}_{q}}$
$\displaystyle=\frac{1}{C_{q}}\frac{C_{q}+C_{g}}{2C_{g}+C_{q}}$ (10)
$\displaystyle\frac{1}{\overline{C}_{\text{od}}}$
$\displaystyle=\frac{1}{C_{q}}\frac{C_{g}}{2C_{g}+C_{q}}$ (11)

Taking into account that the renormalized qubit matrix is
$\frac{1}{\overline{C}_{q}}=\frac{1}{C_{q}}-\frac{1}{\overline{C}_{\text{od}}},$
we express the Hamiltonian as a sum of the non-coupled Hamiltonian $H_{0}$
plus a perturbation $V$, as in Eq. (2) of the manuscript:

$\displaystyle H$ $\displaystyle=H^{(0)}+\epsilon V$ (12) $\displaystyle
H^{(0)}$
$\displaystyle=\sum_{i=1,2}\frac{q_{i}^{2}}{2C_{q}}+\mathcal{U}(\phi_{i})$
(13) $\displaystyle V$
$\displaystyle=\frac{q_{1}q_{2}}{C_{q}}-\sum_{i=1,2}\frac{q_{i}^{2}}{2C_{q}}$
(14)

where $\epsilon=C_{q}/\overline{C}_{\text{od}}$. The second term in the
previous Hamilton is a perturbation to the non-coupled system when
$\epsilon\ll 1$. The Hamiltonian for the general case Eq. (2) with asymmetric
capacitances $C_{1},C_{2}$ can be obtained from the inverse of the capacitance
matrix $1/\overline{C}_{1}=1/C_{1}-C_{2}/(\overline{C}_{\text{od}}C_{1})$ and
$1/\overline{C}_{2}=1/C_{2}-C_{1}/(\overline{C}_{\text{od}}C_{2})$.

##  Appendix B APPENDIX 2: Perturbation theory for the capacitive coupling of
two flux qubits

We use the Schrieffer–Wolff transformation [44] and expand the effective
Hamiltonian in series of the small parameter $\epsilon$, as in Eq. (3) of the
main part of the manuscript. The perturbation series up to third order can be
expressed as:

$\displaystyle
H_{\text{eff}}=H^{(0)}+\epsilon\mathcal{M}_{1}+\frac{\epsilon^{2}}{2}\mathcal{M}_{2}+\frac{\epsilon^{3}}{2}\mathcal{M}_{3}$
(15)

with matrix:

$\displaystyle\mathcal{M}_{1}$ $\displaystyle=P_{0}VP_{0}$ (16)
$\displaystyle\mathcal{M}_{2}$
$\displaystyle=P_{0}\hat{S}(V_{\text{od}})P_{0}$ (17)
$\displaystyle\mathcal{M}_{3}$
$\displaystyle=P_{0}\hat{V}_{\text{od}}\mathcal{L}\hat{V}_{d}(S)P_{0}$ (18)

where the adjoint representation is $\hat{Y}(x)=[Y,X]$ and the operators
$P_{0},Q_{0}$ project onto the qubit and excited subspaces of the unperturbed
system. The notation $O_{\text{od}}$ is used for the non-diagonal part of an
operator $O_{\text{od}}=P_{0}OQ_{0}+Q_{0}OP_{0}$ and
$S=\mathcal{L}(V_{\text{od}}),$ being
$\mathcal{L}(O)=\sum_{ij}\left(O_{\text{od}}\right)_{ij}/E_{ij}\ket{j}\bra{i}$
(matrix elements are denoted by $O_{ij}=\braket{i}{O}{j}$ and energy
differences $E_{ij}=E_{i}-E_{j}$).

We now write the explicit forms for the matrix involved in the computation of
the effective Hamiltonian up to third order. We do so by employing Latin and
Greek letters to denote unperturbed qubit and excited states, respectively:

$\displaystyle\mathcal{M}_{1}$ $\displaystyle=\sum_{i,j}V_{ij}\ket{j}\bra{i}$
(19) $\displaystyle\mathcal{M}_{2}$
$\displaystyle=\sum_{i,j,\alpha}\frac{V_{i\alpha}V_{\alpha
j}}{E_{i\alpha}}\left(\ket{j}\bra{i}-\ket{i}\bra{j}\right).$ (20)
$\displaystyle\mathcal{M}_{3}$
$\displaystyle=\left(\sum_{\alpha,j,\beta,k}\frac{V_{j\alpha}V_{\alpha\beta}V_{\beta
k}}{E_{\alpha j}E_{\beta
j}}-\sum_{\alpha,i,j,k}\frac{V_{ji}V_{i\alpha}V_{\alpha k}}{E_{\alpha
i}E_{\alpha j}}\right)\left(\ket{k}\bra{j}-\ket{j}\bra{k}\right)$ (21)

the braket of the interactions in the qubit subspace is
$V_{ij}=\braket{i}{V}{j},$ on the excited one
$V_{\alpha\beta}=\braket{\alpha}{V}{\beta}$ and the off diagonal between qubit
and excited states $V_{\alpha i}=\braket{\alpha}{V}{i}.$

As discussed in the main body of the work, the first order corrections Eq.19
are given by the projection of the perturbation to the qubit subspace. They
produce a reonormalization of the gap mass and a coupling of the $YY$ type.
The higher order corrections that couple the qubits involve operators that
move the state of the superconducting circuit from qubit to excited states. We
analyze in the following those corrections.

##  Appendix C APPENDIX 3: Second and third order corrections to qubit-qubit
coupling.

We focus now on the second order correction for two three Josephson junctions
qubits (3JJQ) coupled via a capacitor. We do not analyze contributions that
renormalizes each qubit gap, only contributions that couple the two qubit
system. The relevant diagram is shown in panel in Fig. 3(b) of the main
manuscript, where the qubits visit high energy states due to the off-diagonal
part of the interaction. We need to make strong simplifications on the qubit
spectrum, as we explain in what follows, in order to treat analytically second
order corrections. Each qubit is approximated by the first fourth eigenstate
of the uncoupled superconducting circuit. We use the notation
$\ket{\pm_{g}}=(\ket{L_{e}}\pm\ket{R_{e}})$ for the qubit eigenstates and
$\ket{\pm_{e}}=(\ket{L_{e}}\pm\ket{R_{e}})$ for second and third excited
states, being $\ket{L_{g}}$ and $\ket{L_{e}}$ the ground and first excited
state inside the left well of the potential and similar notation for the right
well $\ket{R_{g}},\ \ket{R_{e}}$, see Fig. 1(b) of the main text for a picture
of the spectrum.

The unperturbed projected Hamiltonian onto the subspace expanded by the first
four states of each qubit is:

$\displaystyle
H^{(0)}=\sum_{i=1,2}\frac{\Delta}{2}P_{0i}^{\dagger}\sigma_{i}^{z}P_{0i}+\sum_{i=1,2}Q_{0i}^{\dagger}\left(\hbar\omega+\frac{\Delta_{e}}{2}\sigma_{i}^{z}\right)Q_{0i}$
(22)

where $P_{0i}$ and $Q_{0i}$ projects onto the unperturbed low and excited
subspace of qubit $i$, respectively. We have extended the domain of Pauli
matrices so that they act on ground
$\sigma_{i}^{z}\ket{\pm_{g}}=\pm\ket{\pm_{g}}$ and excited states
$\sigma_{i}^{z}\ket{\pm_{e}}=\pm\ket{\pm_{e}}.$ In order to obtain the second
order term, we need to compute the non-diagonal energy between the low-energy
and high energy sector. We recall that, in the case treated here,
$V_{\text{od}}=P_{0}VQ_{0}+Q_{0}VP_{0}$, so that we can approximate the off-
diagonal elements of the interaction as:

$\displaystyle
V_{i_{g},\alpha_{e}}=\braket{i_{g}}{q_{1}q_{2}}{\alpha_{e}}=\hbar\omega_{q}\delta_{i\alpha}$
(23)

where $i,\alpha=1,2,3,4.$ The low energy sector of the coupled system is
$\ket{1_{g}}=\ket{+_{g},+_{g}},$ $\ket{2_{g}}=\ket{+_{g},-_{g}},$
$\ket{3_{g}}=\ket{-_{g},+_{g}}$ and $\ket{4_{g}}=\ket{-_{g},-_{g}}$.
Similarly, we have for the excited states $\ket{1_{e}}=\ket{+_{e},+_{e}},$
$\ket{2_{e}}=\ket{+_{e},-_{e}},\dots$. The effective harmonic frequency of
each of the single-qubit potential wells is $\omega_{q},$ see Fig. 1(b) of the
main text. We do not consider matrix process that only take one qubit outside
of the low-energy subspace and bring it back. Those processes have an
amplitude that scale as $\sqrt{\hbar\omega_{q}\Delta},$ which is much smaller
than the one which takes the two qubits outside of the qubit subspace in Eq.
(23) (we are in the anharmonic limit $\Delta\ll\hbar\omega_{q}$). We then
approximate the second order corrections as:

$\displaystyle\mathcal{M}_{2}=2(\hbar\omega_{q})^{2}\left(\frac{\ket{+_{g}+_{g}}\bra{+_{g}+_{g}}}{\hbar\omega_{q}+2(\Delta_{e}-\Delta)}+\frac{\ket{-_{g}-_{g}}\bra{-_{g}-_{g}}}{\hbar\omega_{q}-2(\Delta_{e}-\Delta)}+\frac{\ket{-_{g}+_{g}}\bra{-_{g}+_{g}}+\ket{+_{g}-_{g}}\bra{+_{g}-_{g}}}{\hbar\omega_{q}}\right)$
(24)

Expanding the denominator up to second order in $\Delta_{e}-\Delta$ and using
the identity $\ket{\pm}\bra{\pm}=(1\pm\sigma^{z})/2,$ we found that:

$\displaystyle H^{(2)}$
$\displaystyle=g^{(2)}_{qq}\sigma^{z}_{1}\sigma_{2}^{z}.$ (25) $\displaystyle
g^{(2)}_{qq}$
$\displaystyle=\frac{\epsilon^{2}}{2}\frac{(\Delta_{e}-\Delta)^{2}}{\hbar\omega_{q}}$
(26)

In the case of two 3JJQs, we can express previous formula using the ratio
between large Josephson junction and coupling capacitances $\gamma=C_{g}/C$,
instead of the parameter in the perturbative expansion via
$\epsilon=C_{q}/\overline{C}_{\text{od}}=\gamma/(1+2(\alpha+\gamma+\beta))$.
The second order corrections are then:

$\displaystyle
g^{(2)}_{qq}=\frac{\gamma^{2}}{8\hbar\omega_{q}}\left(\frac{\Delta_{e}-\Delta}{1+2\alpha+2\beta}\right)^{2},$
(27)

This is the formula that is employed to plot results in Fig. 4 of the main
text.

Using similar approximations as before, we can find that third order
corrections induces a coupling of the type $\sigma_{1}^{x}\sigma_{2}^{x}$. To
do so, we analyze the third order correction correspond to the second term
inside the parenthesis at the right hand side of Eq. (21). These corrections
can be approximated as:

$\displaystyle
g^{(3)}_{qq}\approx\sum_{\alpha,i,j,k}\frac{V_{ji}V_{i\alpha}V_{\alpha
k}}{E_{\alpha i}E_{\alpha
j}}=\frac{\mathcal{M}_{1}}{\hbar\omega_{q}}\mathcal{M}_{3}$ (28)

Taking into account the shape of the first two corrections, this terms would
give a contribution to coupling $\sigma_{1}^{x}\sigma_{2}^{x}$ as explained in
Fig. 3. Although there are other third order corrections, the analysis
performed here shows that the coupling $XX$ dominates the third order. This is
exactly the case in our numeric in Fig. 3 of the main part of the manuscript,
where $g_{xx}$ depends on $\epsilon^{3}$ at small $\epsilon$.