(eccv) Package eccv Warning: Package ‘hyperref’ is loaded with option
‘pagebackref’, which is *not* recommended for camera-ready version

11institutetext: 1National Supercomputing Center in Shenzhen 2Tsinghua
University
3The University of Hong Kong 4Sun Yat-Sen University 5SenseTime Group
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# A2-MAE: A Spatial-temporal-spectral Unified Remote Sensing Pre-training
Method Based on Anchor-aware Masked Autoencoder

Lixian Zhang 1122Those authors contribute equally to this work.Those authors
contribute equally to this work. Yi Zhao 2211footnotemark: 111footnotemark: 1
Runmin Dong 2211footnotemark: 111footnotemark: 1Corresponding
authors.Corresponding authors. Jinxiao Zhang 22 Shuai Yuan 33 Shilei Cao 44
Mengxuan Chen 22 Juepeng Zheng 4422footnotemark: 222footnotemark: 2 Weijia Li
44 Wayne Zhang 55 Wei Liu 55 Litong Feng 55 Haohuan Fu 2222footnotemark:
222footnotemark: 2

###### Abstract

Vast amounts of remote sensing (RS) data provide Earth observations across
multiple dimensions, encompassing critical spatial, temporal, and spectral
information which is essential for addressing global-scale challenges such as
land use monitoring, disaster prevention, and environmental change mitigation.
Despite various pre-training methods tailored to the characteristics of RS
data, a key limitation persists: the inability to effectively integrate
spatial, temporal, and spectral information within a single unified model. To
unlock the potential of RS data, we construct a Spatial-Temporal-Spectral
Structured Dataset (STSSD) characterized by the incorporation of multiple RS
sources, diverse coverage, unified locations within image sets, and
heterogeneity within images. Building upon this structured dataset, we propose
an Anchor-Aware Masked AutoEncoder method (A2-MAE), leveraging intrinsic
complementary information from the different kinds of images (featuring
different resolutions, spectral compositions, and acquisition times) and geo-
information to reconstruct the masked patches during the pre-training phase.
A2-MAE integrates an anchor-aware masking strategy and a geographic encoding
module to comprehensively exploit the properties of RS images. Specifically,
the proposed anchor-aware masking strategy dynamically adapts the masking
process based on the meta-information of a pre-selected anchor image, thereby
facilitating the training on images captured by diverse types of RS sources
within one model. Furthermore, we propose a geographic encoding method to
leverage accurate spatial patterns, enhancing the model generalization
capabilities for downstream applications that are generally location-related.
Extensive experiments demonstrate our method achieves comprehensive
improvements across various downstream tasks compared with existing RS pre-
training methods, including image classification, semantic segmentation, and
change detection tasks. The dataset and pre-training model will be released.

Figure 1: An illustrative overlook of the proposed global-scale dataset STSSD
and pre-training method A2-MAE. STSSD is a comprehensive remote sensing
dataset structured and characterized by the inclusion of diverse spatial,
temporal, and spectral coverage. A2-MAE facilitates the efficient utilization
of the intrinsic complementarity information in STSSD within one unified
spatial-temporal-spectral model. Lon., Lat., and GSD indicate longitude,
latitude, and ground spatial distance information, respectively.

##  1 Introduction

Earth observations through remote sensing (RS) constitute a fundamental tool
for monitoring the evolution of global-scale phenomena, including the
urbanization process [38, 53], land-use change [25, 26], and biodiversity loss
[43, 45]. Over the past half-century, there has been a substantial increase in
the volume of RS data, resulting in spatial, temporal, and spectral
diversities within extensive RS image archives. The spatial-temporal-spectral
diversities inherent in RS images offer critical and complementary information
for comprehensive analysis and recognition of objects and scenes.
Consequently, RS plays an pivotal role in various operational and complex
research domains within the field of geoscience.

In response to the challenges encountered by extensive RS downstream
applications arising from the scarcity of expensive annotations [46], self-
supervised learning (SSL) [21, 20] has emerged as a promising technique for
deriving robust feature representations from a vast repository of satellite
images. The acquired feature representations can subsequently be fine-tuned
with limited labeled data for specific downstream applications. Despite
considerable efforts in constructing large pre-trained models through methods
like masked autoencoders (MAE) or contrastive learning (CL), most of the
existing RS SSL methodologies have been custom-tailored for specific
scenarios, such as temporal SatMAE [13], multi-spectral SatMAE [13], multi-
resolution ScaleMAE [34], and spatiotemporal foundation models [50]. These
methods enhance performance in specific downstream tasks but fall short of
achieving comprehensive improvements across various downstream tasks. Besides,
the SSL methods underutilize geographical information, a powerful prior for
leveraging spatial patterns. In the work, we address the pivotal question: How
can we present a single spatial-temporal-spectral unified RS pre-training
method to effectively leverage a diverse collection of RS images?

The key to this question lies in two aspects. The first aspect involves the
construction of a location-unified and extensive RS dataset encompassing
images with varying temporal, spatial resolutions, and spectral compositions.
Presently available RS datasets are typically derived from one or two
satellite sources, offering limited spatial-temporal-spectral coverage. For
instance, the Million-AID dataset exclusively covers optical RS images with
RGB bands [30, 31]. SEN12MS [37] and SSL4EO-S12 [47] exhibit constrained
temporal coverage. The SeCo [33] and CACo [32] datasets are pre-trained
exclusively on Sentinel-2 images. However, real-world RS data exhibits
significant variations in spatial resolution, temporal coverage, and spectral
composition. The SSL models trained on homogeneous RS data struggle to provide
effective representation for fine-tuning of downstream tasks involving
different RS sources.

The second aspect entails mining the intrinsic relevance from the images with
different spatial, temporal, and spectral characteristics through SSL
techniques. A straightforward approach is to design separate backbones for
different types of sources and align the representations of different types
[19]. However, this method leads to a linear escalation in model parameters
and computational overhead with the expansion of source type count. As there
are a large number of types of RS sources, such as Landsat-8 with 7 bands,
Sentinel-2 with 13 bands, Gaofen-2 with 4 bands, and WorldView-2 with 8 bands,
it is difficult to simultaneously model the relationship across different
types of sources.

To address these challenges, we introduce STSSD (Figure 1), a global-scale RS
dataset containing half a million sampling locations with 2.5 million spatial-
temporal-spectral structured image sets collected from multiple multi-spectral
sources. Each image set is meticulously crafted to exhibit different spatial
resolutions, temporal and spectral compositions for the same location. Our
data processing method preserves heterogeneity within images and diversity
across images for SSL. To harness the rich and varied representation features
within STSSD effectively, we propose an Anchor-Aware Masked AutoEncoder method
(A2-MAE), including an anchor-aware masking strategy and geographic encoding
module (Figure 1). The proposed anchor-aware masking strategy enables training
on images captured by diverse sources within one unified spatial-temporal-
spectral model. Besides, the proposed geographic encoding method allows the
model to leverage accurate spatial patterns, unleashing the potential of geo-
location priors for downstream tasks. Experiments verify that our method
achieves comprehensive improvements across various downstream tasks compared
with state-of-the-art RS SSL methods. Taking DynamicEarthNet as an example,
the performance can be further enhanced by over 8.4% on mIoU through the
introduction of geographic information during the fine-tuning process (refer
to Sec. 5.3).

In summary, our contributions are as follows:

  * •

We build the STSSD, a globally spatial-temporal-spectral structured RS dataset
featuring high diversity, unification, and heterogeneity. STSSD is
meticulously curated to encompass diverse land-use types spatially, capture
landscape changes temporally, and incorporate various band compositions
spectrally.

  * •

We propose a pre-training method, A2-MAE, designed to accommodate various
types of RS sources within a unified backbone architecture. A2-MAE leverages
spatial-temporal-spectral relationships and geographical information to
improve model representation and generalization capabilities.

  * •

Experiments verify the effectiveness and advantages of A2-MAE compared to
existing RS pre-training models with similar complexities across image
classification, semantic segmentation, and change detection tasks.

##  2 Related Work

###  2.1 Large-scale datasets for remote sensing imagery pre-training

Inspired by the achievements of computer vision (CV) datasets [15, 29, 16, 1],
researchers have introduced several large-scale RS datasets [11, 39, 40, 41,
47, 5].

These datasets exhibit a gradual expansion in the volume of data, starting the
fMoW [11] encompassing 1 million images, progressing to BigEarthNet-MM [40]
with 1.2 million images, and further expanding to SSL4EO-S12 [47] comprising 3
million images. Additionally, there has been a progression in the diversity of
spectral sources in datasets, transitioning from datasets like BigEarthNet
[39] solely from Sentinel-1, to BigEarthNet-MM [40] combining Sentinel-1/2
pairs, to SatlasPretrain [5], which incorporates data from Sentinel-1/2 and
NAIP, and then to DynamicEarthNet [41] containing diverse spatial-temporal-
spectral images with constrained sampling locations. Therefore, there is an
urgent need to construct a large-scale spatial-temporal-spectral structured RS
dataset encompassing more multi-spectral sources and diverse coverage. In this
work, we introduce STSSD for spatial-temporal-spectral unified learning,
surpassing the DynamicEarthNet dataset by 10 times and incorporating data from
4 multi-spectral sources.

###  2.2 Self-supervised learning for satellite imagery

SSL primarily focuses on generating supervisory signals from unlabeled data,
through the design of various pretext tasks such as masked patches
reconstruction [20, 42, 51, 3, 48] and contrasting semantically similar inputs
[9, 21, 18, 10, 7]. Furthermore, SSL enables the acquisition of semantic
information without human annotation. Therefore, SSL plays a vital role in the
RS domain [46], where annotation demands specialized expertise and incurs high
costs. Existing RS pre-training methods leverage different properties of RS
images or specific RS tasks [46]. For instance, Ayush et al.[2] leverage
spatially aligned but temporally separated images as positive pairs to learn
feature representation for 10m multi-spectral images. Similarly, Mall et
al.[32] propose a new SSL loss for CL to distinguish between short-term and
long-term changes in multi-spectral images. Cong et al.[13] introduce SatMAE
to leverage temporal or multi-spectral information in data through positional
encoding. Reed et al.[34] present Scale-MAE to reconstruct both low and high-
frequency images to learn robust multi-scale representations for RS imagery.
Nevertheless, these studies are customized for specific types of RS images and
cannot simultaneously utilize RS images from different kinds of multi-spectral
sources in one unified model. To fill this gap, we propose an anchor-aware
masking strategy to leverage intrinsic complementarity information from an
image set, which can be easily extended to various multi-spectral sources.

###  2.3 Geography-aware learning

RS images offer essential metadata records containing geographic information,
such as geographic location and ground sample distance (GSD) [49, 28]. This
prior information enables the capture of geographic patterns and fosters a
robust linkage between fine-tuning data and models pre-trained globally [52,
12]. Consequently, it is anticipated to bolster the representational capacity
of the pre-trained model [27, 4]. While a few studies have leveraged recorded
geographic data [2, 36], they are constrained in efficiently utilizing such
information on a large scale [35]. one-hot geo-encoding [11] offers limited
encoding outcomes, while GSD scaling encoding [34] cannot be jointly
integrated with geo-location data. Another alternative (_i.e_., geo-context
prototype learning [19]) demands additional computational resources while
yielding encoding outcomes unsuitable for varying spatial resolutions. To
bridge this gap, we introduce a Geographic Encoding Module in A2-MAE,
providing more accurate geographical priors without additional computation
overhead (_i.e_., latitude, longitude, and GSD), thereby improving the
generalization of applications on a global scale.

##  3 Data

###  3.1 Overview

Figure 2: The compositions of image sets and the global sampling location
distribution. (a) shows the S2-L8 image sets. "ng" denotes the non-growth
period, and "g" denotes the growth period within one year. (b) shows the GF-S2
image sets. (c) is the sampling location distribution (i.e., purple circles
for S2-L8 in urban areas, green circles for S2-L8 in nature reserves, and pink
circles for GF-S2).

We introduce STSSD, a large-scale RS dataset designed for spatial-temporal-
spectral unified SSL. This dataset is meticulously curated through data
pruning from an initial pool of 4.2 million original images collected at
1,045K sampling locations. The resulting STSSD comprises 510K image sets, each
containing up to six images collected from two sources with different
resolutions and spectral compositions. As shown in Figure 2, STSSD consists of
4 kinds of image sets, featuring diverse sources, spatial resolutions,
coverage, and acquisition times. It is characterized by the following four key
attributes:

Diversity. STSSD exhibits comprehensive source diversity, containing 4
different satellite sources with 3 different band compositions. It also boasts
spectral diversity, manifested in different band compositions with resolutions
ranging from 0.8m/pixel to 30m/pixel, derived from various and distinct data
sources.

Coverage. STSSD owns dynamic and diverse geographical features, involving more
than 12,000 urban centers and 10,000 nature reserves around the world. This
global coverage enhances the model’s abilities in downstream tasks with
diverse coverage and geographical characteristics.

Unification. STSSD integrates spatial, temporal, and spectral contexts to
generate image sets sourced from different origins and acquisition times.
These image sets are spatial-temporal-spectral unified, providing more
dimensions of information compared to single-temporal or single-source data,
thereby significantly contributing to the unification and robustness of the RS
foundation model.

Heterogeneity. Employing a clustering-based data pruning strategy to eliminate
redundant types (_e.g_., desert) and low-quality images, STSSD balances the
diversity across images and heterogeneity within an image. The heterogeneity
can increase the difficulty of image reconstruction for SSL, thereby
empowering the model with enhanced feature representation capabilities.

###  3.2 STSSD Construction

In the pursuit of constructing a unified RS dataset characterized by diverse
sources, we strategically opt for Gaofen-1 (4 bands with 1 m/pixel), Gaofen-2
(4 bands with 0.8 m/pixel), Sentinel-2 (13 bands with 10 m/pixel), and
Landsat-8 (7 bands with 30 m/pixel) as our sources to maximize coverage while
considering sources’ accessibility, each contributing diverse resolutions,
with variations up to 37.5$\times$, and multiple band compositions.

Building upon the texture-rich images of urban areas proved by previous work
[32], we further expand to include nature reserves [43, 45], so as to enhance
the model’s understanding and capabilities of a more diversified and dynamic
planet. Consequently, we meticulously select over the original 1,045K sampling
locations, spanning nature reserves (depicted in green) and main cities
(depicted in purple), to collect Sentinel-2 and Landsat-8 image sets (S2-L8),
as illustrated in Figure 2 (c). To capture the dynamic nature of geographical
features, a time series of images are provided for each sampling location,
ranging from the year 2020 to 2023, with periodic seasonal revisits.
Furthermore, we utilize the locations of the available Gaofen images to gather
the corresponding Sentinel-2 images, subsequently forming Sentinel-2 and
Gaofen image sets (GF-S2) to enhance the representation ability for higher-
resolution data (depicted in pink in Figure 2).

The structuring of these image sets is designed to ensure optimal resolution
and band gaps for effective model learning. Specifically, there are 2 kinds of
image sets: S2-L8 image sets and GF-S2 image sets. For S2-L8 (Figure 2 (a)),
the image sets collected from main cities comprise 6 images, involving 3
Sentinel-2 and 3 Landsat-8 images annually from 2021 to 2023, to capture the
temporal changes in land cover. As for nature reserves, we conduct the image
sets comprising 4 images, including 2 Sentinel-2 and 2 Landsat-8 images during
both the growth and non-growth periods in 2020, to showcase the phenological
characteristics. For GF-S2 (Figure 2 (b)), each image set integrates 3 images,
including a Gaofen-1 or Gaofen-2 image and 2 Sentinel-2 images captured at
different time points. Note that each image set comprises two distinct data
sources with different sources and temporal snapshots. The integration of
diverse image sources, characterized by varying spatial resolutions, within a
single image set enables multi-scale observation of the same geographical
areas. This simultaneous consideration of fine-grained details and broader
contextual views facilitates a more comprehensive feature representation,
consequently enhancing performance across diverse downstream tasks such as
building extraction and land cover mapping. Moreover, the incorporation of
multi-temporal information ensures accessibility to temporal dynamics, thereby
fortifying the robustness of temporal variations for downstream tasks such as
change detection.

Since the original STSSD contains observations from areas with high
homogeneities, such as deserts which do not contribute significantly to the
diversity and complexity of the dataset due to their uniform nature, we employ
a data pruning strategy to remove redundant contents and filter out low-
quality images, resulting in a more refined and curated collection of data.
This process ensures that the images in STSSD are high-quality and
heterogeneous. After pruning, the final STSSD owns over 510K sampling
locations with 2.5 million curated images. Refer to supplementary material for
more details about STSSD, such as the data retrieving and pre-processing.

##  4 Methodology

###  4.1 Overall Architecture

Figure 3: The overall framework of A2-MAE. A2-MAE incorporates an anchor-aware
masking strategy and a geographic encoding module, allowing for the efficient
utilization of spatial-, temporal-, and spectral-variant information in large-
scale RS imagery.

As illustrated in Figure 3, A2-MAE is a self-supervised pre-training method
based on the MAE [20], which makes two key contributions to the MAE framework
to unlock representative potentials in STSSD. First, A2-MAE presents an
anchor-aware masking (AAM) strategy to utilize the image sets collected from
different sources. The AAM dynamically adjusts the masking strategy according
to the meta-information of a pre-selected anchor image for each training
iteration. This adaptation learning allows the model to leverage the intrinsic
complementarity of spatial-temporal-spectral information to reconstruct the
masked patches, thereby improving the model representation ability. In
addition, A2-MAE introduces a geographic encoding module (GEM) to obtain the
geo-embedding of the given image set to provide accurate geographical priors
for A2-MAE, improving the model generalization ability.

###  4.2 Setup

Since the image sets in the STSSD are gathered by geo-locations, we denote
$P_{i}=\\{I_{i}^{1,1},...,I_{i}^{t,s}\\}$ represents the image set at the
$i^{th}$ location. $I_{i}^{t,s}\in\mathbb{R}^{H\times W\times C}$ represents
an RS image with $H$ height, $W$ width, and $C$ channels, which is captured bu
source $s$ at time $t$. Note that images from different sources $s$ have
different representative features, including different spectral compositions
and GSDs. Three images of $P_{i}$ are randomly selected as input $I_{in}$ for
A2-MAE, ensuring a minimum of 2 different $s$ and 2 different $t$ to capture
sufficient diversity in spatio-temporal-spectral relationships while balance
computational cost. The A2-MAE then patchifies the selected $I_{in}$ into
three sets of sequence $Seq$ of independent patches. After randomly removing a
fraction of the obtained patches, the A2-MAE reconstructs the removed patches
by leveraging the complementarity information within the remaining patches
from $Seq$. Unlike the traditional MAE architecture, the A2-MAE includes an
AAM to encourage the A2-MAE to implicitly leverage the intrinsic
complementarity information within $I_{in}$ and a GEM to introduce the
geographic-related information.

###  4.3 Anchor-Aware Masking Strategy

Existing RS SSL methods are often tailored for specific scenarios, limiting
the capabilities to leverage symbiotic features among images and increasing
the costs for transferring to other scenarios. In contrast, our method works
towards a unified pre-training method that potentially benefits various
representative features between images with different spatial resolutions,
temporal, and spectral compositions in $P_{i}$. However, this symbiotic and
diverse complementarity information within the image sets also poses
challenges for obtaining robust and generalized RS representative features due
to the complexity of spatial-temporal-spectral relationships.

To jointly utilize the spatial-temporal-spectral information within the image
sets, a straightforward method is to apply the random masking strategy to
different spectral combinations of the input image set. However, when input
image sets have diverse combinations of spatial resolutions and temporal
compositions, the random masking strategy may lead to feature leakage from the
remaining high-resolution patches during the reconstruction of the removed
low-resolution patches at the same position, resulting in shortcut learning
during model pre-training. To this end, we propose the AAM to dynamically
adjust the masking strategy of images for each input $I_{in}$, enabling
training with images from diverse sources while preventing feature leakage.
Specifically, we adopt a consistent masking strategy for images from different
sources $s$ at the same retrieval time $t$, a mutually-exclusive masking
strategy for images from the same $s$ at different $t$, and a random masking
strategy for the other circumstances. For a quantitative ablation study on
AAM, please refer to Section 5.3

In the example depicted in Figure 3, three images is randomly sampled from an
image set $P_{i}$, specifically, $I_{in}=\\{I_{i}^{2020,Sen2}$,
$I_{i}^{2020,Lan8}$, $I_{i}^{2022,Sen2}\\}$. Taking the middle image
$I_{i}^{2020,Sen2}$ as the referenced anchor, A2-MAE explicitly _aware_ the
meta-information (_i.e_., source $s$ and time $t$) to opt specific masking
strategy for removing patches from the other two images (_i.e_.,
$I_{i}^{2020,Lan8}$ and $I_{i}^{2022,Sen2}$). We first randomly select three
bands of the anchor image $I_{i}^{2020,Sen2}$ to encompass a substantial
diversity of band compositions while balancing the computational costs. If an
image in $I_{in}$ differs in $s$ but the same in $t$ as the anchor image
(_i.e_., the bottom image $I_{i}^{2020,Lan8}$), a consistent masking strategy
is employed, obtaining a patch sequence $Seq_{i}^{2020,Lan8}$ where the
patches are removed from the same position as those in $Seq_{i}^{2020,Sen2}$.
This ensures that remaining patches maintain their coarsest version for the
same source $s$, preventing the A2-MAE from feature leakage during pre-
training. To address temporal disparities, if an image in $I_{in}$ has a
different time $t$ from the same $s$ (_i.e_., the upper image
$I_{i}^{2022,Sen2}$), a mutually-exclusive masking strategy is adopted to
ensure the removed patch retains positional uniqueness with
$Seq_{i}^{2020,Sen2}$, enhancing A2-MAE’s capacity to leverage multi-temporal
symbiotic features. Additionally, the incorporation of $I_{in}$ offers
sufficient diversity in spatio-temporal-spectral relationships, encouraging
A2-MAE to effectively leverage multi-scale symbiotic features for patch
reconstruction.

###  4.4 Geographic Encoding

Figure 4: The proposed geographic encoding module in A2-MAE. By efficiently
encoding the geographic metadata (_i.e_., a $Geo_{GSD}$ and four sets of
$(Geo_{Lat}^{c},Geo_{Lon}^{c})$ in RS images, the geographic encoding module
encourages A2-MAE to be explicitly aware of this crucial geographic prior for
varying geographic information in images.

The metadata stored with the RS images contains geographic information,
including the latitude, longitude, and GSD. The latitude and longitude
information indicates the absolute location of the retrieved image on the
Earth, which is of significance in leveraging the geographic pattern in the
pre-training on worldwide RS datasets. The GSD indicates the ground scale of
the RS image, which is critical to understanding the spatial ranges and
frequency specificity of the image. For example, an image with a low GSD has
more details in high frequency than a high GSD image does. Therefore, we
propose the GEM to explicitly incorporate essential geographic priors into the
MAE model, thereby enhancing its generalization capabilities for downstream
applications.

As illustrated in Figure 4, given an RS image, the corresponding metadata
records one GSD $Geo_{GSD}$ and four sets of latitude and longitude (_i.e_.,
the four corners $(Geo_{Lat}^{c},Geo_{Lon}^{c}),c\in\\{TL,TR,BL,BR\\}$).
Instead of directly utilizing the decimal geographic information, the GEM
views the RS image as a group of squared grids and encodes it to achieve
better representative geo-encoding features. Let $\mathbb{G}$ be the set of
grids formed with latitudes and longitudes. The $\mathbb{G}$ contains several
levels of mesh in a pyramid way, which are composed of equally subdivided
grids with integer coding. Let _Level 0_ mesh $G_{0}=\\{g_{0}\\}\in\mathbb{G}$
be a $512^{\circ}\times 512^{\circ}$ grid which covers the whole global area.
_Level 1_ mesh $G_{1}=\\{g_{1},...,g_{n1}\\}\in\mathbb{G}$ is defined as
equally subdivided four grids, each of which has $256^{\circ}\times
256^{\circ}$ in the height and width. Following this logic, _Level k_ mesh
$G_{k}$ is obtained by a quadtree division of _Level k-1_ mesh $G_{k-1}$. The
higher level mesh represents finer resolution with low GSD. In all, given the
metadata of an RS image, we first query the closest _Level k_ according to its
$Geo_{GSD}$, then four sets of latitude and longitude are embedded as
sequences of binary array. For example, The $Geo_{GSD}$ of Landsat-8 images is
30 m, which can be referred to grids in _Level 21_ (32 m of the size of each
grid in the equator). Considering that the $Geo_{GSD}$ is utilized as an
approximate version in this encoding strategy, we further encode the precise
$Geo_{GSD}$ by replacing the positional embedded vector with the ground scaled
positional encoding vector (inspired by [34]), which can be embedded as
follows:

$v_{gsd,x}(pos,2i)=sin\frac{\widehat{Geo_{GSD}}}{Geo_{GSD}}\frac{pos}{10000^{\frac{2i}{D}}},$
(1)

$v_{gsd,y}(pos,2i+1)=cos\frac{\widehat{Geo_{GSD}}}{Geo_{GSD}}\frac{pos}{10000^{\frac{2i}{D}}},$
(2)

where $pos$ is the position of the embedded patch along the given axis, $i$ is
the patch index, and $D$ is the number of embedded dimensions, exactly as
introduced in [44]. $\widehat{Geo_{GSD}}$ is the reference GSD (nominally set
to 1 m).

As a result, the proposed GEM efficiently embeds the geographic metadata,
providing unique embedded features for images with specific locations and
GSDs.

##  5 Experiments

###  5.1 Implementation Details and Baselines

We adopt ViT-Large architecture as the backbone for the proposed A2-MAE, and
pre-train the A2-MAE using the constructed STSSD. We employ a progressive
training strategy [24] by starting with S2-L8 data and then progressively
transitioning to GF-S2 in STSSD. The patch size is fixed to $16\times 16$
pixels. A2-MAE is pre-trained for 130 epochs with a batch size of 1,024 on 8
NVIDIA A800 GPUs. AdamW optimizer [22] is utilized with an initial learning
rate of 0.0001, coupled with a half-cycle cosine decay schedule. According to
existing works [20, 13], we adopt a masking ratio of 75%, balancing training
efficiency and pretext task difficulty. Five self-supervised learning methods
with the officially released pre-trained weights are selected as competing
methods in this study, including 2 ResNet-50 based methods (SeCo [33], CACo
[32]) and 3 ViT-Large based methods (vanilla MAE pre-trained on ImageNet-1k
[20], SatMAE (the version for spectral data) [13], and ScaleMAE [34]). Full
fine-tuning is employed in the downstream tasks for all methods, of which the
first layer is modified to fit the data structure of specific datasets.

###  5.2 Comparison Results

We conduct experiments on 7 datasets with diverse distributions in spatial,
temporal, and spectral coverage, encompassing various data sources. This
ensures a comprehensive assessment of the capabilities in efficiently
utilizing spatial-, temporal-, and spectral-variant features, involving
different downstream tasks, including classification, segmentation, and change
detection. We employ the encoder of the competing pre-trained models for all
downstream tasks, and details of the training setups for fine-tuning of each
task are included in the supplementary materials. As shown in Table 1, A2-MAE
achieves comprehensive improvements across downstream tasks, indicating the
effectiveness of A2-MAE in exploiting the properties of RS images.

Land Cover Classification: We perform the scene classification task on EuroSAT
[23] and the multi-label classification task on BigEarthNet [39]. EuroSAT
comprises 27K Sentinel-2 images with 13 bands collected from 34 European
countries. BigEarthNet encompasses 590K Sentinel-2 images with 13 bands
collected from 10 countries. As presented in Table 1, despite that several
competing methods (_i.e_., SeCo, CACo, and SatMAE) are custom-tailored and
pre-trained on Sentinel-2 image dataset, A2-MAE still outperforms all
competing methods in both EuroSAT and BigEarthNet datasets, highlighting the
proposed A2-MAE’s effectiveness in leveraging diverse RS information within
one unified model.

Table 1: Comparison results (%) of classification, segmentation, and change
detection. The vanilla MAE model was pre-trained on ImageNet-1k.

Methods | Classification | Segmentation | Change detection
---|---|---|---
EuroSAT | BigEarthNet | Sen1Floods11 | CropSeg | LEVIR-CD | OSCD | DynamicEarthNet
Accuracy | mAP | mIoU | mIoU | mIoU | F1 | mIoU
SeCo[33] | 95.63 | 82.6 | 79.10 | 32.78 | 82.37 | 46.78 | 43.1
CACo[32] | 95.90 | 82.1 | 84.62 | 32.83 | 82.90 | 51.74 | 41.2
Vanilla MAE[20] | 98.78 | 80.1 | 85.42 | 44.51 | 80.41 | 34.74 | 39.8
SatMAE[13] | 98.98 | 82.1 | 88.77 | 44.21 | 82.11 | 47.63 | 44.7
ScaleMAE[34] | 98.31 | 82.2 | 85.05 | 44.04 | 83.07 | 48.70 | 44.3
A2-MAE | 99.09 | 83.0 | 88.87 | 44.81 | 84.32 | 53.97 | 46.0

Semantic Segmentation: We perform experiments on Sen1Floods11 [6] and
CropSeg[17]. Sen1Floods11 is a surface water segmentation dataset including
4,831 Sentinel-2 imagery with 13 bands covering 120,406 km2 and spans 14
biomes and 6 continents of the world across 11 flood events. CropSeg is a
cropland segmentation dataset containing 3,854 Harmonized Landsat-Sentinel
imagery with 7 bands at 30 m resolution across the Contiguous United States.
Given that SeCo and CACo are pre-trained on datasets covering only urban
regions, A2-MAE, which is pre-trained on STSSD covering diverse land cover
types, achieves significant improvements by 9.77%/12.03% against SeCo and
4.25%/11.98% against CACo, highlighting its superior generalization ability
when pre-trained on STSSD with diverse coverage and geographical
characteristics.

Change Detection: We conduct experiments on the LEVIR-CD [8], OSCD [14], and
DynamicEarthNet datasets [41]. LEVIR-CD comprises 637 image pairs with a
resolution of 0.5 m and a time span ranging from 5 to 14 years. The OSCD
dataset comprises Sentinel-2 images with 13 bands collected from 24 urbanized
regions worldwide. DynamicEarthNet provides daily images from Planet with 4
bands at 3 m resolution and monthly images from Sentinel-2 with 13 bands
across approximately 75 areas of interest worldwide. Table 1 presents the
quantitative evaluation results of baselines and A2-MAE. A2-MAE outperforms
competing methods across various backbones and self-supervised architectures
by 1.25% on mIoU / 2.23% on F1 / 1.3% on mIoU against the second-best results
in LEVIR-CD, OSCD, and DynamicEarthNet, respectively. These observed
improvements underscore the effectiveness of the proposed A2-MAE in harnessing
multi-temporal RS images within a single unified model.

###  5.3 Ablation Study

To efficiently and fairly investigate the key contributions of A2-MAE, _i.e_.,
AAM and GEM, we conduct ablation experiments by pre-training and fine-tuning
A2-MAE on DynamicEarthNet, which is split into 55 locations for training and
10 locations for testing. The pre-training and the training phase of fine-
tuning are conducted on the split training locations. A scratch version of
SatMAE, which is pre-trained and further fine-tuned on DynamicEarthNet
training sites, is viewed as the baseline. We also conduct comparisons in
terms of masking strategy (_i.e_., random masking and tube masking strategies
[13]) and geographic encoding method (_i.e_., one-hot geographic encoding [11]
and scale encoding [34]). We further encode the geographic information using
the proposed GEM during fine-tuning on DynamicEarthNet. This configuration
(denoted as A2-MAE+) gives us a promising glance at the potential of fully
utilizing the GEM when the geographic metadata of RS images is available in
downstream RS tasks.

Table 2: Performance comparison (%) between different sub-components of A2-MAE
on DynamicEarthNet change detection. SatMAE* denotes a scratch version of
SatMAE architecture pre-trained on DynamicEarthNet and further fine-tuned on
DynamicEarthNet. A2-MAE+ denotes introducing the geographic encoding during
fine-tuning.

Pre-training | Pix. Acc. | mIoU
---|---|---
SatMAE* [13] | 72.2 | 45.9
Van. MAE (random masking) | 69.2 | 40.2
\+ tube masking [13] | 67.4 | 38.6
\+ AAM | 72.7 | 46.8
\+ One-hot geo-encoding[11] | 71.9 | 47.1
\+ scale encoding [34] | 72.8 | 47.3
\+ full GEM (A2-MAE) | 73.8 | 47.7
A2-MAE+ | 76.4 | 56.1

As shown in Table 2, A2-MAE outperforms the scratched SatMAE by 1.6/0.8 of
Pix. Acc./mIoU, decoupling and showcasing the contributions of the proposed
pre-training method. Besides, both the AAM (+ 6.6% on mIoU) and GEM (+ 0.9% on
mIoU) contribute to the significant performance improvements of A2-MAE.
Specifically, for masking strategies, the proposed AAM outperforms random
masking and tube masking strategies by 6.6% on mIoU, indicating the
effectiveness of AAM. For geographic encoding methods, comparisons against
one-hot geographic encoding [11] show improvements by 1.9%/0.6% on Acc./mIoU.
Further ablation studies reveal enhancements by 0.1%/0.5% on Acc./mIoU for GSD
embedding and 1.0%/0.4% on Acc./mIoU for Lat/Lon embedding.

Furthermore, by introducing the geographic information via the GEM in fine-
tuning phase, A2-MAE+ achieves a notable improvement of 8.4% on mIoU against
A2-MAE, indicating the promising potential of the GEM. It reveals that for
downstream tasks that provide raw geographic metadata, introducing the GEM
during fine-tuning can improve the results by a large margin.

###  5.4 Discussion on Model Efficiency in Exploiting RS Images with Diverse
Characteristics

Given the vast quantity and varied characteristics of RS images, it is crucial
to efficiently exploit the intrinsic relevance of images with diverse spatial,
temporal, and spectral attributes. Previous studies have predominantly focused
on addressing specific facets of this diversity, such as multi-spectral [24]
and multi-resolution images [34], or by achieving separate pre-training models
[13]. Consequently, achieving comprehensive and generalizable improvements
across downstream tasks that span spatial, temporal, and spectral dimensions
remains challenging. In response, a contemporaneous study, Skysense [19],
designs separate backbones for three types of sources in a larger model with
2.06 billion parameters, which is trained on 80 A100 GPUs. However, this
approach faces difficulties in scaling model parameters and computational
overhead to accommodate expanding RS sources, which is evidently
unsustainable. Different from this method, A2-MAE explores the utilization of
various multi-spectral sources in a unified backbone. Benefiting from the
proposed anchor-aware masking strategy, A2-MAE enables the efficient
exploitation of the intrinsic complementarity information within RS images
from different multi-spectral sources within one unified spatial-temporal-
spectral model. Moreover, it requires 6$\times$ fewer model parameters than
[19], thus facilitating efficient pre-training utilizing only 8 A800 GPUs,
significantly saving computational costs while achieving comprehensive
improvements across various downstream tasks.

##  6 Conclusion

In this study, we introduce STSSD, a spatial-temporal-spectral structured
dataset comprising 510K sampling locations with 2.5 million structured images
collected from multiple RS sources. To exploit different kinds of multi-
spectral sources in one unified backbone, we propose an anchor-aware masking
strategy to harness the intrinsic complementary information from different
kinds of images, thus achieving more powerful feature representations.
Furthermore, we propose the geographic encoding module to leverage geographic
information, thereby improving the model generalization ability. Experiments
verify the effectiveness and advantages of our method compared to existing RS
pre-training models with the same parameter amount across image
classification, semantic segmentation, and change detection tasks. In future
work, we will expand the diversity of modalities such as Synthetic Aperture
Radar and hyperspectral images in STSSD and A2-MAE.

## References

  * [1] Abu-El-Haija, S., Kothari, N., Lee, J., Natsev, P., Toderici, G., Varadarajan, B., Vijayanarasimhan, S.: Youtube-8m: A large-scale video classification benchmark. arXiv preprint arXiv:1609.08675 (2016)
  * [2] Ayush, K., Uzkent, B., Meng, C., Tanmay, K., Burke, M., Lobell, D., Ermon, S.: Geography-aware self-supervised learning. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 10181–10190 (2021)
  * [3] Bachmann, R., Mizrahi, D., Atanov, A., Zamir, A.: Multimae: Multi-modal multi-task masked autoencoders. In: European Conference on Computer Vision. pp. 348–367. Springer (2022)
  * [4] Bai, L., Huang, W., Zhang, X., Du, S., Cong, G., Wang, H., Liu, B.: Geographic mapping with unsupervised multi-modal representation learning from vhr images and pois. ISPRS Journal of Photogrammetry and Remote Sensing 201, 193–208 (2023)
  * [5] Bastani, F., Wolters, P., Gupta, R., Ferdinando, J., Kembhavi, A.: Satlaspretrain: A large-scale dataset for remote sensing image understanding. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 16772–16782 (2023)
  * [6] Bonafilia, D., Tellman, B., Anderson, T., Issenberg, E.: Sen1floods11: A georeferenced dataset to train and test deep learning flood algorithms for sentinel-1. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops. pp. 210–211 (2020)
  * [7] Caron, M., Touvron, H., Misra, I., Jégou, H., Mairal, J., Bojanowski, P., Joulin, A.: Emerging properties in self-supervised vision transformers. In: Proceedings of the IEEE/CVF international conference on computer vision. pp. 9650–9660 (2021)
  * [8] Chen, H., Shi, Z.: A spatial-temporal attention-based method and a new dataset for remote sensing image change detection. Remote Sensing 12(10), 1662 (2020)
  * [9] Chen, T., Kornblith, S., Norouzi, M., Hinton, G.: A simple framework for contrastive learning of visual representations. In: International conference on machine learning. pp. 1597–1607. PMLR (2020)
  * [10] Chen, X., He, K.: Exploring simple siamese representation learning. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 15750–15758 (2021)
  * [11] Christie, G., Fendley, N., Wilson, J., Mukherjee, R.: Functional map of the world. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 6172–6180 (2018)
  * [12] Chu, G., Potetz, B., Wang, W., Howard, A., Song, Y., Brucher, F., Leung, T., Adam, H.: Geo-aware networks for fine-grained recognition. In: 2019 IEEE/CVF International Conference on Computer Vision Workshop (ICCVW). pp. 247–254 (2019)
  * [13] Cong, Y., Khanna, S., Meng, C., Liu, P., Rozi, E., He, Y., Burke, M., Lobell, D., Ermon, S.: Satmae: Pre-training transformers for temporal and multi-spectral satellite imagery. Advances in Neural Information Processing Systems 35, 197–211 (2022)
  * [14] Daudt, R.C., Le Saux, B., Boulch, A., Gousseau, Y.: Urban change detection for multispectral earth observation using convolutional neural networks. In: IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium. pp. 2115–2118. Ieee (2018)
  * [15] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. pp. 248–255. Ieee (2009)
  * [16] Everingham, M., Eslami, S.A., Van Gool, L., Williams, C.K., Winn, J., Zisserman, A.: The pascal visual object classes challenge: A retrospective. International journal of computer vision 111, 98–136 (2015)
  * [17] Geospatial, I.N.: Ibm-nasa-geospatial/multi-temporal-crop-classification · datasets at hugging face, https://huggingface.co/datasets/ibm-nasa-geospatial/multi-temporal-crop-classification
  * [18] Grill, J.B., Strub, F., Altché, F., Tallec, C., Richemond, P., Buchatskaya, E., Doersch, C., Avila Pires, B., Guo, Z., Gheshlaghi Azar, M., et al.: Bootstrap your own latent-a new approach to self-supervised learning. Advances in neural information processing systems 33, 21271–21284 (2020)
  * [19] Guo, X., Lao, J., Dang, B., Zhang, Y., Yu, L., Ru, L., Zhong, L., Huang, Z., Wu, K., Hu, D., He, H., Wang, J., Chen, J., Yang, M., Zhang, Y., Li, Y.: Skysense: A multi-modal remote sensing foundation model towards universal interpretation for earth observation imagery (2023)
  * [20] He, K., Chen, X., Xie, S., Li, Y., Dollár, P., Girshick, R.: Masked autoencoders are scalable vision learners. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16000–16009 (2022)
  * [21] He, K., Fan, H., Wu, Y., Xie, S., Girshick, R.: Momentum contrast for unsupervised visual representation learning. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 9729–9738 (2020)
  * [22] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 770–778 (2016)
  * [23] Helber, P., Bischke, B., Dengel, A., Borth, D.: Eurosat: A novel dataset and deep learning benchmark for land use and land cover classification. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12(7), 2217–2226 (2019)
  * [24] Hong, D., Zhang, B., Li, X., Li, Y., Li, C., Yao, J., Yokoya, N., Li, H., Ghamisi, P., Jia, X., Plaza, A., Gamba, P., Benediktsson, J.A., Chanussot, J.: Spectralgpt: Spectral remote sensing foundation model (2024)
  * [25] Huang, X., Schneider, A., Friedl, M.A.: Mapping sub-pixel urban expansion in china using modis and dmsp/ols nighttime lights. Remote Sensing of Environment 175, 92–108 (2016). https://doi.org/10.1016/j.rse.2015.12.042
  * [26] Imhoff, M.L., Lawrence, W.T., Elvidge, C.D., Paul, T., Levine, E., Privalsky, M.V., Brown, V.: Using nighttime dmsp/ols images of city lights to estimate the impact of urban land use on soil resources in the united states. Remote Sensing of Environment 59(1), 105–117 (1997). https://doi.org/10.1016/S0034-4257(96)00110-1
  * [27] Klemmer, K., Rolf, E., Robinson, C., Mackey, L., Rußwurm, M.: Satclip: Global, general-purpose location embeddings with satellite imagery (2023)
  * [28] Li, W., Chen, K., Chen, H., Shi, Z.: Geographical knowledge-driven representation learning for remote sensing images. IEEE Transactions on Geoscience and Remote Sensing 60, 1–16 (2021)
  * [29] Lin, T.Y., Maire, M., Belongie, S., Hays, J., Perona, P., Ramanan, D., Dollár, P., Zitnick, C.L.: Microsoft coco: Common objects in context. In: Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13. pp. 740–755. Springer (2014)
  * [30] Long, Y., Xia, G.S., Li, S., Yang, W., Yang, M.Y., Zhu, X.X., Zhang, L., Li, D.: On creating benchmark dataset for aerial image interpretation: Reviews, guidances and million-aid. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 14, 4205–4230 (2021)
  * [31] Long, Y., Xia, G.S., Zhang, L., Cheng, G., Li, D.: Aerial scene parsing: From tile-level scene classification to pixel-wise semantic labeling (2022)
  * [32] Mall, U., Hariharan, B., Bala, K.: Change-aware sampling and contrastive learning for satellite images. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 5261–5270 (2023)
  * [33] Manas, O., Lacoste, A., Giró-i Nieto, X., Vazquez, D., Rodriguez, P.: Seasonal contrast: Unsupervised pre-training from uncurated remote sensing data. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 9414–9423 (2021)
  * [34] Reed, C.J., Gupta, R., Li, S., Brockman, S., Funk, C., Clipp, B., Keutzer, K., Candido, S., Uyttendaele, M., Darrell, T.: Scale-mae: A scale-aware masked autoencoder for multiscale geospatial representation learning. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 4088–4099 (2023)
  * [35] Rußwurm, M., Klemmer, K., Rolf, E., Zbinden, R., Tuia, D.: Geographic location encoding with spherical harmonics and sinusoidal representation networks (2023)
  * [36] Scheibenreif, L., Hanna, J., Mommert, M., Borth, D.: Self-supervised vision transformers for land-cover segmentation and classification. In: 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW). pp. 1421–1430 (2022)
  * [37] Schmitt, M., Hughes, L.H., Qiu, C., Zhu, X.X.: Sen12ms–a curated dataset of georeferenced multi-spectral sentinel-1/2 imagery for deep learning and data fusion. arXiv preprint arXiv:1906.07789 (2019)
  * [38] Small, C., Pozzi, F., Elvidge, C.D.: Spatial analysis of global urban extent from dmsp-ols night lights. Remote Sensing of Environment 96(3-4), 277–291 (2005). https://doi.org/10.1016/j.rse.2005.02.002
  * [39] Sumbul, G., Charfuelan, M., Demir, B., Markl, V.: Bigearthnet: A large-scale benchmark archive for remote sensing image understanding. In: IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium. pp. 5901–5904. IEEE (2019)
  * [40] Sumbul, G., De Wall, A., Kreuziger, T., Marcelino, F., Costa, H., Benevides, P., Caetano, M., Demir, B., Markl, V.: Bigearthnet-mm: A large-scale, multimodal, multilabel benchmark archive for remote sensing image classification and retrieval [software and data sets]. IEEE Geoscience and Remote Sensing Magazine 9(3), 174–180 (2021)
  * [41] Toker, A., Kondmann, L., Weber, M., Eisenberger, M., Camero, A., Hu, J., Hoderlein, A.P., Şenaras, Ç., Davis, T., Cremers, D., et al.: Dynamicearthnet: Daily multi-spectral satellite dataset for semantic change segmentation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 21158–21167 (2022)
  * [42] Tong, Z., Song, Y., Wang, J., Wang, L.: Videomae: Masked autoencoders are data-efficient learners for self-supervised video pre-training. Advances in neural information processing systems 35, 10078–10093 (2022)
  * [43] Turner, W.: Sensing biodiversity. Science 346(6207), 301–302 (2014)
  * [44] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [45] Wang, R., Gamon, J.A.: Remote sensing of terrestrial plant biodiversity. Remote Sensing of Environment 231, 111218 (2019)
  * [46] Wang, Y., Albrecht, C.M., Braham, N.A.A., Mou, L., Zhu, X.X.: Self-supervised learning in remote sensing: A review. arXiv preprint arXiv:2206.13188 (2022)
  * [47] Wang, Y., Braham, N.A.A., Xiong, Z., Liu, C., Albrecht, C.M., Zhu, X.X.: Ssl4eo-s12: A large-scale multi-modal, multi-temporal dataset for self-supervised learning in earth observation. arXiv preprint arXiv:2211.07044 (2022)
  * [48] Woo, S., Debnath, S., Hu, R., Chen, X., Liu, Z., Kweon, I.S., Xie, S.: Convnext v2: Co-designing and scaling convnets with masked autoencoders. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16133–16142 (2023)
  * [49] Wu, X., Shi, Z., Zou, Z.: A geographic information-driven method and a new large scale dataset for remote sensing cloud/snow detection. ISPRS Journal of Photogrammetry and Remote Sensing 174, 87–104 (2021). https://doi.org/https://doi.org/10.1016/j.isprsjprs.2021.01.023, https://www.sciencedirect.com/science/article/pii/S0924271621000290
  * [50] Yao, F., Lu, W., Yang, H., Xu, L., Liu, C., Hu, L., Yu, H., Liu, N., Deng, C., Tang, D., et al.: Ringmo-sense: Remote sensing foundation model for spatiotemporal prediction via spatiotemporal evolution disentangling. IEEE Transactions on Geoscience and Remote Sensing (2023)
  * [51] Zhang, S., Chen, H., Yang, H., Sun, X., Yu, P.S., Xu, G.: Graph masked autoencoders with transformers. arXiv preprint arXiv:2202.08391 (2022)
  * [52] Zhao, D., Wang, Q., Zhang, J., Bai, C.: Mine diversified contents of multi-spectral cloud images along with geographical information for multi-label classification. IEEE Transactions on Geoscience and Remote Sensing (2023)
  * [53] Zhou, Y., Smith, S.J., Zhao, K., Imhoff, M., Thomson, A., Bond-Lamberty, B., Asrar, G.R., Zhang, X., He, C., Elvidge, C.D.: A global map of urban extent from nightlights. Environmental Research Letters 10(5), 054011 (2015). https://doi.org/10.1088/1748-9326/10/5/054011


# Computer-aided quantization and numerical analysis of superconducting
circuits

Sai Pavan Chitta1, Tianpu Zhao1, Ziwen Huang1,2, Ian Mondragon-Shem1, Jens
Koch1 1 Department of Physics and Astronomy, Northwestern University,
Evanston, Illinois 60208, USA  2 _current address:_ Superconducting Quantum
Materials and Systems Center, Fermi National Accelerator Laboratory (FNAL),
Batavia, IL 60510, USA

###### Abstract

The development of new superconducting circuits and the improvement of
existing ones rely on the accurate modeling of spectral properties which are
key to achieving the needed advances in qubit performance. Systematic circuit
analysis at the lumped-element level, starting from a circuit network and
culminating in a Hamiltonian appropriately describing the quantum properties
of the circuit, is a well-established procedure, yet cumbersome to carry out
manually for larger circuits. We present work utilizing symbolic computer
algebra and numerical diagonalization routines versatile enough to tackle a
variety of circuits. Results from this work are accessible through a newly
released module of the scqubits package.

<EMAIL_ADDRESS>superconducting qubits, circuit quantization

##  1 Introduction

The past two decades have seen dramatic progress in the design and control of
quantum devices based on superconducting circuits [1, 2, 3]. A central part of
this effort is the design of qubits which are robust against decoherence and
can be controlled and measured with high fidelity [4, 5, 6]. The systematic
discovery of new qubits requires exploring the wealth of possible circuits
formed by combining inductors, capacitors, and Josephson junctions as the
basic circuit elements. This has led to the design of a growing number of
qubit types, each with its own advantages [7, 8, 9, 10, 4, 11, 12]. The task
of exploring new circuit designs, however, can be tedious and become
numerically challenging since the Hilbert space dimension needed to describe
the low-energy physics of the system tends to grow exponentially with the
number of elements in the circuit. In order to make progress in the design of
devices of increasing complexity, it is essential to promote tools that
facilitate the efficient analysis of spectra and related key properties of
superconducting circuits.

Several prescriptions have been proposed to approximate the low-energy
Hamiltonian of a superconducting circuit, each depending on the level of
detail needed for an accurate description [13, 14, 15, 16, 17, 18]. For
example, several approaches and software tools model the full geometry and
distributed nature of the circuit [13, 15, 18]. This is computationally
costly, and often works most readily for circuits with only weak
anharmonicity. The most basic, yet quite widely applicable quantization
framework is based on the lumped-element description of a circuit [19, 20, 21,
22, 23, 17]. Despite the fact that this approach is the simplest one, it
nonetheless presents challenges once the number of circuit elements is
increased beyond the bare minimum of two (found, for example, in the transmon
qubit). This has led to the development of software packages that aim to
automate the construction of the Hamiltonian [24, 25, 26, 27].

In this work, we present a streamlined computer-aided approach to quantize
lumped-element systems. By examining the circuit topology, we show how to
classify and systematically detect four different types of degrees of freedom
in a user-specified circuit: free, frozen, periodic and extended degrees of
freedom. The identification of these variables can help reduce the number of
degrees of freedom needed to describe the system and makes manifest the basis
that should be used for an efficient representation of the Hamiltonian. This
functionality has recently been incorporated into the scqubits Python package
[28], leveraging the full array of capabilities already available. We
complement our approach with a hierarchical-diagonalization algorithm. This
further mitigates the circuit-size limitations governed generally by the size
of the Hamiltonian matrix which describes the low-energy spectrum of interest.

For illustration, consider the task of finding the circuit Hamiltonian and
numerical eigenenergies of a KITE circuit [29]. Using a simple input file
`KITE.yaml` that encodes the circuit graph and associated circuit-element
parameters (see Appendix A.1), generating the circuit Hamiltonian and
performing numerical diagonalization is readily achieved by a few lines of
code:

⬇

In: import scqubits as scq

kite = scq.Circuit.from_yaml("KITE.yaml", ext_basis="harmonic")

kite.sym_hamiltonian()

⬇

Out:
$(15.7Q_{1}^{2}+26.4Q_{1}Q_{2}-52.8Q_{1}Q_{3}+26.4Q_{2}^{2}-52.8Q_{2}Q_{3}+52.8Q_{3}^{2})+(0.18\Phi_{1}^{2}+0.72\Phi_{1}\theta_{1}-0.36\Phi_{1}\theta_{2}+0.36\Phi_{1}\theta_{3}+0.18\Phi_{2}^{2}+0.72\Phi_{2}\theta_{1}-0.36\Phi_{2}\theta_{2}+1.9\theta_{1}^{2}-1.44\theta_{1}\theta_{2}+1.18\theta_{1}\theta_{3}+0.36\theta_{2}^{2}-0.36\theta_{2}\theta_{3}+0.295\theta_{3}^{2}-5.9\cos{\left(\theta_{2}\right)}-5.9\cos{\left(\theta_{2}+\theta_{3}\right)})$

Once suitable parameters for Hilbert space truncation are chosen, the
eigenenergies are obtained simply by

⬇

In: kite.eigenvals()

⬇

Out: array([4.17710271, 6.47610652, 7.24813846, 8.47986653, 8.713839,
9.19210006])

The outline of the paper is as follows. In Section 2 we set the stage for the
conceptual discussion of our approach by reviewing the lumped-element
description of superconducting circuits. This includes a simple discussion of
constructing the circuit Lagrangian in terms of node variables. In Section
2.2, we show how to systematically identify the four types of degrees of
freedom in circuits consisting of capacitances, inductances,111Terminology
employed here refers to the ideal circuit elements by “capacitance” and
“inductances”, whereas “capacitor” and “inductor” are understood as real
building blocks that will include parasitic inductance and capacitance. and
Josephson junctions. In Section 4, we describe how the identification of these
variables is implemented in the scqubits package to efficiently represent the
Hamiltonian of a user-provided circuit. Finally, Section 5 concludes with a
summary and outlook of future developments.

##  2  Circuit Lagrangian and classification of variables

###  2.1 Brief review of circuit Lagrangians

The common starting point for the quantization of a superconducting circuit is
the construction of a circuit Lagrangian [20, 30, 31, 32, 21, 33]. In this
paper, we exclusively focus on lumped-element circuits which are composed of
capacitances, inductances, and Josephson tunnel junctions. (For more general
discussions including additional elements such as gyrators, or extending the
discussion beyond the lumped-element regime, consult Refs. [34, 35, 36, 37].)
Although seemingly restrictive, this small set of circuit elements spans a
large family of different circuits, and includes the majority of
superconducting qubits studied in the past, and used today. Each circuit
element forms a branch in the connected circuit network.

Formulation of the Lagrangian depends on the definition of an appropriate set
of generalized coordinates (which we also refer to as variables). In circuits
including Josephson junctions, generalized-flux variables are the most
convenient choice. For each branch $b,$ one defines

$\Phi_{b}(t)=\int_{t_{0}}^{t}d\tau\,V_{b}(\tau),$ (1)

where $t_{0}$ is a reference time222In some references, $t_{0}$ is chosen as
$-\infty$. To avoid questions about convergence of the integral, we choose
$t_{0}$ to be finite and denote a time in the “distant” past when no external
fields were present. and $V_{b}(t)$ is the voltage across branch $b$. The
energy stored in a branch $b$ is given by the

$\displaystyle{\text{{charging}}\atop\text{\raisebox{8.53581pt}{{energy:}}}}\;\;E=\frac{C_{b}\dot{\Phi}^{2}_{b}}{2},\qquad{\text{{inductive}}\atop\text{\raisebox{8.53581pt}{{energy:}}}}\;\;E=\frac{\Phi^{2}_{b}}{2L_{b}},\qquad{\text{{junction}}\atop\text{\raisebox{8.53581pt}{{energy:}}}}\;\;E=-E_{Jb}\cos(2\pi\Phi_{b}/\Phi_{0}),$

for a capacitance $C_{b}$, inductance $L_{b}$, and junction with Josephson
energy $E_{Jb}$, respectively. ($\Phi_{0}=h/2e$ is the flux quantum.)

The collection of all branch variables $\Phi_{b}(t)$ is in general not a set
of independent variables: for each elementary loop $\mathsf{L}_{\ell}$ in the
circuit, Kirchhoff’s loop rule in conjunction with fluxoid quantization [38,
21] imposes a constraint on the corresponding branch fluxes:

$\Phi_{\text{ext}}^{\ell}=\sum_{b\in\mathsf{L}_{\ell}}(\pm)_{b}\Phi_{b}(t).$
(2)

Here, $\Phi_{\text{ext}}^{\ell}$ is the external magnetic flux through the
$\ell$-th loop, and the branch sign $(\pm)_{b}$ is determined according to the
directionality of the branch relative to the chosen loop orientation.

The construction of a spanning tree provides a systematic way to determine a
set of independent branch variables consistent with the constraints Eq. (2). A
spanning tree $\mathsf{T}$ is a subset of connected branches such that every
circuit node is connected to ground by one, and only one, path within the
spanning tree [21] [see the example in Fig. 1(a)].333Spanning trees are not
unique. This reflects the fact that the constraints (2) leave open the choice
of which branch variable in a loop to eliminate. Each branch not part of the
spanning tree is a closure branch. Iterative inclusion of closure branches,
i.e., starting from the spanning tree, adding the first closure branch, then
the second, etc., leads to the closing of one circuit loop at a
time.444Careful analysis of the circuit in Fig. 1 illustrates that there are
cases in which the inclusion of the next closure branch appears to close
multiple loops at a time. However, given the motivating purpose of resolving
the loop constraints, it is always the newly closed loop containing exactly
one closure branch which determines the relevant external loop flux. A more
extensive discussion of this point has been presented in [39] Once the
constraints have been resolved in this way, the circuit Lagrangian can be
expressed in terms of the branch variables $\Phi_{b\in\mathsf{T}}$ belonging
to the spanning tree, with $\Phi_{b\notin\mathsf{T}}$ eliminated by means of
Eq. (2),

$\mathcal{L}=T-V=\sum_{b}\left[\frac{1}{2}C_{b}\dot{\Phi}_{b}^{2}-\frac{1}{2L_{b}}\Phi_{b}^{2}+E_{J,b}\cos(2\pi\Phi_{b}/\Phi_{0})\right]\bigg{|}_{b\notin\mathsf{T}:\text{
replace $\Phi_{b}$ via Eq.\ \eqref{Eq_SimpConstr}}}.$ (3)

Here, notation follows a simple scheme assigning the capacitance $C_{b}$ to
branch $b=(n^{\prime},n)$ if a capacitance indeed connects nodes $n$ and
$n^{\prime}$. Otherwise, the term is omitted from the summation. Analogous
statements hold for inductance and junction contributions to the potential
energy. We note in passing that time-dependent external flux renders the
holonomic constraint (2) among the coordinates time-dependent, and hence also
yields a constraint among the generalized velocities $\dot{\Phi}_{b}$ [40,
37]; for simplicity we limit ourselves to the time-independent case here.

Figure 1: Example of a lumped-element superconducting circuit. (a) One choice
of a spanning tree (marked in red) for the circuit, connecting the ground node
with each other node by a unique path. (Only a subset of labels for branch
variables and external fluxes are shown.) (b) For the same circuit, coloring
here indicates: a capacitively-coupled subcircuit or isolated island (blue), a
junction-coupled subcircuit or superconducting island (red), and an
inductively-coupled subcircuit (green). Labels denote the node variables.

Node variables present a useful alternative to branch variables. The node
variable $\phi_{n}$ at node $n$ corresponds to the generalized flux obtained
by adding up branch fluxes along the spanning tree; specifically, the spanning
tree defines one path $\mathsf{P}_{n}$ connecting the ground node to node $n$.
Then, the node variable $\phi_{n}$ is obtained by summing the branch variables
along that path555Note: while the voltage across branch $b$ is well-defined,
one cannot generally fall back to a description in terms of electric-potential
differences. In the presence of non-zero external magnetic flux, even if
constant in the present, the magnetic field had to be switched on after
$t_{0}$ and be ramped up to its present value. As a result, the integral in
equation (1) necessarily extends over a time interval during which the
electric field was not curl-free, and could not be represented as a gradient
field of the electric potential.:

$\phi_{n}=\sum_{b\in\mathsf{P}_{n}}\Phi_{b},$ (4)

as illustrated in Fig. 1(a). Branch variables within the spanning tree as well
as the generalized flux associated with each closure branch can then be
expressed in terms of node variables via

$\Phi_{b\in\mathsf{T}}=\phi_{n}-\phi_{n^{\prime}},\qquad\Phi_{b\notin\mathsf{T}}=\phi_{n}-\phi_{n^{\prime}}+\Phi_{\text{ext}}^{\ell}.$
(5)

Here $n,n^{\prime}$ are the nodes connected by the branch $b=(n,n^{\prime})$,
and $\Phi_{\text{ext}}^{\ell}$ is the external flux threading the loop $\ell$
formed by the closure branch $b$ and the spanning tree. The circuit Lagrangian
(2) can then be re-expressed in terms of node variables in the general form

$\mathcal{L}(\dot{\boldsymbol{\phi}},\boldsymbol{\phi})=\frac{1}{2}\dot{\boldsymbol{\phi}}^{\intercal}\mathbf{C}\,\dot{\boldsymbol{\phi}}-V(\boldsymbol{\phi}),$
(6)

where $\boldsymbol{\phi}^{\intercal}=(\phi_{1},\ldots,\phi_{N})$, and the
$N\times N$ capacitance matrix $\mathbf{C}$ has non-zero off-diagonal entries
$(\mathbf{C})_{nn^{\prime}}=-C_{b}$ whenever a capacitance $C_{b}$ connects
nodes $n\not=n^{\prime}$; diagonal entries are given by
$(\mathbf{C})_{nn}=-\sum_{m\neq n}(\textbf{C})_{nm}$. The potential energy in
node variables is obtained via Eq. (3) by substitution,

$V(\boldsymbol{\phi})=\sum_{b=(n,n^{\prime})}\left[\frac{1}{2L_{b}}(\phi_{n}^{\prime}-\phi_{n}+\mathfrak{cl}_{b}\,\Phi_{\text{ext}}^{\ell_{b}})^{2}-E_{J,b}\cos\left(\frac{2\pi}{\Phi_{0}}[\phi_{n}^{\prime}-\phi_{n}+\mathfrak{cl}_{b}\,\Phi_{\text{ext}}^{\ell_{b}}]\right)\right],$
(7)

where $\mathfrak{cl}_{b}$ is the indicator for closure branches, i.e.
$\mathfrak{cl}_{b}=1$ for $b\notin\mathsf{T}$ and $\mathfrak{cl}_{b}=0$
otherwise.

In principle, one might now proceed with either the branch- or node-variable
Lagrangian, perform a Legendre transformation to obtain the Hamiltonian, and
subsequently impose canonical quantization rules. However, it is generally
worthwhile (especially where the ultimate aim involves numerical treatment) to
consider a further transformation of variables which will achieve two
important goals: (1) we may reduce the number of relevant variables, for
example by identifying conserved quantities, and (2) we cleanly distinguish
variables associated with different boundary conditions upon quantization. A
suitable variable transformation can thus help make the quantum description of
the circuit both more transparent and numerically efficient.

Table 1:  Overview of the different types of variables. Classification is
stated readily in the Lagrangian framework: a variable is _extended_ if the
potential is confining, and _periodic_ if the potential is periodic along
variable’s axis. A variable (or coordinate) absent from the Lagrangian is
free, while a generalized velocity absent leads to a _frozen_ variable. Each
variable category is associated with specific conditions on how the
corresponding subcircuit (hatched square) connects to the remaining circuit.
The total charge $Q$ on the subcircuit is continuous for extended variables,
quantized in units of $2e$ for periodic variables, and conserved in case of a
free variable.

###  2.2 Classification of variables

Any variable $\theta$, no matter whether branch, node, or collective variable
of neither type, can be classified based on four categories:

  * •

Extended:  a variable $\theta_{e}$ is extended if it describes a direction in
configuration space in which the potential energy is confining, expressed as
$\lim_{\theta_{e}\rightarrow\pm\infty}V(\boldsymbol{\theta})=\infty$.

  * •

Periodic:  a variable $\theta_{p}$ is periodic if it aligns with a direction
in configuration space in which the potential energy is periodic, i.e.,
$V(\boldsymbol{\theta}+\Phi_{0}\mathbf{e}_{p})=V(\boldsymbol{\theta})$. Here,
$\mathbf{e}_{p}$ is the unit vector pointing along the axis associated with
the $\theta_{p}$ coordinate. We exclude the case of a variable in which the
potential energy is constant; this situation is covered by the following
category.

  * •

Free: a variable $\theta_{c}$ is free666Another common, though potentially
misleading, name for this variable type is “cyclic.” if
$\partial_{\theta_{c}}\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})=0$.
This means that $\theta_{c}$ does not appear in the potential energy, and the
corresponding conjugate charge
$Q_{c}=\partial_{\dot{\theta}_{c}}\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})$
is conserved (the subscript $c$ refers to constant charge).

  * •

Frozen: a variable $\theta_{f}$ is frozen if
$\partial_{\dot{\theta}_{f}}\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})=0$.
In this case $\dot{\theta}_{f}$ is missing from the kinetic energy. This
somewhat pathological case only arises when too many parasitic capacitances
are neglected in the lumped element model. In the simplest setting, a frozen
variable leads to an additional constraint given by
$\partial_{\theta_{f}}\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})=0$.

The diagram on the right shows the relation among these four categories. The
_free_ , _periodic_ , and _extended_ categories are mutually exclusive
properties of the potential energy. The _frozen_ category regards the kinetic
energy and as such can occur for each of the former three categories.

This classification holds several advantages when analyzing the Hamiltonian of
a circuit. Recognizing free and frozen variables helps reduce the number of
degrees of freedom needed to describe the system. Furthermore, knowledge of
periodic and extended variables guides a consistent and efficient choice of
basis respecting the appropriate boundary conditions. Neither branch nor node
variables are generally the suitable choice which recognizes the full number
of periodic, free, and frozen degrees of freedom. Failure to find all the
variables of each type makes diagonalization of the quantum Hamiltonian
usually more difficult, and can potentially lead to wrong conclusions about
offset-charge dependence and charge-noise sensitivity of circuit eigenstates.

##  3 Adaptive variable transformation

In the following, we show how to systematically construct a linear variable
transformation which will pinpoint all periodic, free, and frozen degrees of
freedom. Against naive instinct, it is helpful to define this transformation
by expressing the “old” node variables in terms of the new variables, i.e.,
$\phi_{n}=\sum^{N}_{m=1}z_{nm}\theta_{m}$, where $\mathbf{Z}=(z_{nm})$ is an
$N\times N$ invertible transformation matrix.

Figure 2: Variable transformation, and example of a junction-coupled
subcircuit $\mathsf{S}_{p}$ and its complement $\overline{\mathsf{S}}_{p}$.

We determine the $z_{nm}$ by studying the subcircuits of the system, defined
as connected sub-graphs of the circuit. We distinguish three relevant classes
of subcircuits: (1) A subcircuit forms an isolated island if it is connected
to the rest of the circuit exclusively by capacitances; (2) a subcircuit forms
a superconducting island if it is coupled to the remainder of the circuit by
capacitances and at least one junction (but not via any inductances); (3) a
subcircuit is purely inductively coupled if it connects to the rest of the
circuit only through inductances. Examples of each type of subcircuit are
shown in Fig. 1(b).

###  3.1 Periodic and extended variables

For the analysis of periodic variables, consider a subcircuit $\mathsf{S}_{p}$
forming a superconducting island as in the example of Fig. 2. There,
$\mathsf{S}_{p}$ is shaded in gray and couples to its circuit complement
$\overline{\mathsf{S}}_{p}$ exclusively through capacitances and junctions.
Any inductances present in the circuit are entirely inside $\mathsf{S}_{p}$ or
inside $\overline{\mathsf{S}}_{p}$, but not at their interface. Given this
partitioning, we see there is a simple choice of coefficients $z_{n1}$ that
will render the variable $\theta_{1}$ periodic, namely:

$z_{n1}=1\text{ for all }n\in\mathsf{S}_{p},\quad z_{n1}=0\text{ for all
}n^{\prime}\in\overline{\mathsf{S}}_{p}.$ (8)

In other words, a change in the new flux variable $\theta_{1}$ will
simultaneously raise the flux of each node variable inside the superconducting
island by equal amounts, but keep the flux on all nodes outside the
superconducting island unchanged. Since $\theta_{1}$ clearly produces no flux
difference across any inductance, the potential energy must indeed be periodic
with respect to $\theta_{1}$.

The full set of periodic variables is obtained by identifying all
superconducting islands in the circuit. This indicates that the full number of
periodic variables, $N_{p}$, is equal to the number of irreducible
superconducting islands $\mathsf{S}_{p}$ which cannot be decomposed into
smaller superconducting islands. Once all periodic variables are assigned, any
remaining variable $\theta_{e}$ must now appear in the quadratic potential
energy terms generated by the inductances. These terms have the usual
confining property
$\lim_{\theta_{e}\rightarrow\pm\infty}V(\boldsymbol{\theta})=\infty$. As a
result, all these degrees of freedom are extended.

###  3.2 Free variables

Free variables are treated in a manner nearly identical to the periodic case.
Suppose $\theta_{c}=\theta_{1}$ is one of the free variables. In this case,
changes in $\theta_{c}$ cannot produce node-flux differences across any of the
inductances or junctions. Again, partitioning into subcircuits facilitates the
definition of new variables. For free variables, the subcircuit
$\mathsf{S}_{c}$ must now be an _isolated island_ connected to the remainder
of the circuit $\overline{\mathsf{S}}_{c}$ via capacitances exclusively.
Carrying over the definition of $z_{n1}$ to the isolated island case,

$z_{n1}=1\text{ for all }n\in\mathsf{S}_{c},\quad z_{n1}=0\text{ for all
}n\in\overline{\mathsf{S}}_{c},$ (9)

equips us with a free variable: since flux differences generated by
$\theta_{1}$ are across capacitances only, $\theta_{1}$ will be absent from
all potential energy terms associated with inductances and junctions. As a
result, $\partial_{\theta_{c}}V(\boldsymbol{\theta})=0$, so $\theta_{c}$ is
indeed free. The full number of free variables, $N_{c}$, corresponds to the
total number of irreducible isolated islands in the given circuit.

###  3.3 Frozen variables

As the last category, we discuss the case of frozen variables. For a frozen
variable, the relevant subcircuit is of type (3): a subcircuit coupled to the
remaining circuit through inductances only. Once such a subcircuit is
identified, a frozen variable $\theta_{f}$ is readily defined following the
now familiar way: $z$-coefficients for nodes inside and outside the subcircuit
are set to 1 and 0, respectively. As a result, the only flux differences that
occur under a change in $\theta_{f}$ are across inductance terms in the
Lagrangian. Kinetic energy terms due to capacitances remain unchanged such
that $\partial_{\dot{\theta}_{f}}\mathcal{L}=0$, and $\theta_{f}$ is, hence,
indeed a frozen variable. The full number of frozen variables corresponds to
the total number of irreducible, purely inductively coupled subcircuits.

A few remarks about frozen variables are in order. Variables only appear
frozen in situations where too many parasitic capacitances have been dropped
from the description and the capacitance matrix has become singular. A
foolproof, physically accurate but not necessarily efficient way of handling
this situation is to re-admit the omitted parasitic capacitances which
physically are always present. This is not always desirable as it generally
increases the number of degrees of freedom retained in the modeling of the
circuit.

Dropping a parasitic capacitance from the charging-energy contributions is
akin to setting this capacitance to zero. In the mechanically analogous
system, this corresponds to setting the mass of a particle to zero. The $C\to
0$ limit is generally a singular limit [41] that exhibits subtleties
particularly significant when junction capacitances are omitted [34]. We
exclude this case here and focus on the more benign scenario in which one of
the Euler-Lagrange equations turns into a linear holonomic constraint for
$C=0$, see Subsection 3.5.

###  3.4 Full transformation

As an example, consider the circuit in Fig. 1(b), which has three subcircuits.
The corresponding transformation is

$\begin{pmatrix}\phi_{1}\\\ \phi_{2}\\\ \phi_{3}\\\ \phi_{4}\\\ \phi_{5}\\\
\phi_{6}\\\ \phi_{7}\\\ \end{pmatrix}=\begin{pmatrix}0&0&1&\\\\[-6.99997pt]
0&0&0&\vdots\\\ 0&1&0&\\\ 0&1&0&\,\,\cdots\,\,\,z_{nm}\,\,\,\cdots\,\,\\\
0&0&0&\\\\[-6.99997pt] 1&0&0&\vdots\\\
1&0&0&\end{pmatrix}\begin{pmatrix}\theta_{c}\\\ \theta_{\\!f}\\\ \theta_{p}\\\
\theta_{e,1}\\\ \theta_{e,2}\\\ \theta_{e,3}\\\ \theta_{e,4}\end{pmatrix}.$
(10)

This transformation identifies a free, a frozen and a periodic variable. The
remaining variables are extended. The matrix elements corresponding to the
extended variables $\\{\theta_{e,i}\\}_{i=1}^{4}$ are chosen such that the
transformation is invertible. The Lagrangian in terms of the new variables is

$\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})=\tfrac{1}{2}\dot{\boldsymbol{\theta}}^{\intercal}\widetilde{\mathbf{C}}\,\dot{\boldsymbol{\theta}}-V(\boldsymbol{\theta}),$
(11)

where $\widetilde{\mathbf{C}}=\mathbf{Z}^{\intercal}\mathbf{C}\,\mathbf{Z}.$

###  3.5 Elimination of non-dynamical variables

To obtain the circuit Hamiltonian, we first eliminate the frozen and free
variables. Eliminating frozen variables is, in fact, a necessary step, since
the capacitance matrix [Eq. (11)] is singular in the presence of frozen
variables.

To carry out the elimination of frozen variables from the Lagrangian, we use
that
$\partial_{\theta_{\\!f,i}}\mathcal{L}=\frac{d}{dt}\partial_{\dot{\theta}_{\\!f,i}}\mathcal{L}=0,$
where $i=1,\ldots,N_{f}$ labels the inductively-coupled subcircuits. In
general, this set of $N_{f}$ equations are nonlinear due to the cosine terms
in the Lagrangian. As a crucial simplification, we assume that none of the
junction capacitances are omitted in the description. In that case, only
subcircuits with purely linear-inductive coupling to the remainder of the
circuit lead to frozen variables. Thus, the variables $\\{\theta_{\\!f,i}\\}$
can be solved for using the constraints emerging from the Euler-Lagrange
equations. Upon eliminating the $\\{\theta_{\\!f,i}\\}$ from the Lagrangian,
the kinetic energy is expressed in terms of the invertible capacitance matrix
$\widetilde{\mathbf{C}}_{cpe}$ obtained by removing the rows and columns in
$\widetilde{\mathbf{C}}$ that were associated with frozen variables. At this
point, we perform a Legendre transform to obtain the classical Hamiltonian

$H=\boldsymbol{Q}^{\intercal}_{cpe}\cdot\dot{\boldsymbol{\theta}}_{cpe}-\mathcal{L}(\dot{\boldsymbol{\theta}}_{cpe},\boldsymbol{\theta}_{pe}),$

where
$\boldsymbol{Q}_{cpe}=\partial_{\dot{\boldsymbol{\theta}}_{cpe}}\mathcal{L}$
and $\boldsymbol{\theta}_{cpe}$ corresponds to free, periodic and extended
variables.

For every free variable $\theta_{c}$, the conjugate charge $Q_{c}$ is a
conserved quantity. As shown in B, these conserved charges have no impact on
the dynamics of the remaining degrees of freedom, and they may hence be
omitted. The kinetic energy of the Hamiltonian can thus be written in the form
$T=\frac{1}{2}\mathbf{Q}^{\intercal}_{pe}\widetilde{\mathbf{C}}^{-1}_{pe}\mathbf{Q}_{pe},$
where $\widetilde{\mathbf{C}}^{-1}_{pe}$ is the matrix block of
$\widetilde{\mathbf{C}}^{-1}_{cpe}$ that corresponds to the extended and
periodic variables only. Similarly, the vector $\mathbf{Q}_{pe}$ only includes
the non-conserved charges conjugate to the variables
$\boldsymbol{\theta}_{pe}.$

###  3.6 Quantization

Finally, canonical quantization yields the quantum mechanical Hamiltonian.
Quantization is achieved by promoting the coordinate and conjugate momenta to
operators $\hat{\theta}_{n}$ and $\hat{Q}_{n}=-i\hbar\partial_{\theta_{n}}$
satisfying canonical commutation relations. The Hamiltonian finally takes the
form

$\hat{H}=\tfrac{1}{2}(\hat{\mathbf{Q}}_{pe}-\mathbf{Q}_{g})^{\intercal}\widetilde{\mathbf{C}}_{pe}^{-1}(\hat{\mathbf{Q}}_{pe}-\mathbf{Q}_{g})+V(\hat{\boldsymbol{\theta}}_{pe}),$
(12)

where we have introduced a vector $\mathbf{Q}_{g}$ representing off-set
charges. Since periodic and extended degrees of freedom are identified
explicitly, the boundary conditions required for the wavefunction
$\Psi(\boldsymbol{\theta_{pe})}$ can be stated transparently:

$\displaystyle{\text{{Periodic
directions:}}}\;\;\Psi(\boldsymbol{\theta}_{pe}+\mathbf{e}_{p}\Phi_{0})=\Psi(\boldsymbol{\theta}_{pe}),\qquad{\text{{extended
directions:}}}\;\;\lim_{\theta_{e}\rightarrow\pm\infty}\Psi(\boldsymbol{\theta}_{pe})=0.$

The basis choice is further simplified by the distinction between periodic and
extended variables. Periodic degrees of freedom can be represented in the
discrete charge basis, extended variables in the phase basis or an
appropriately chosen oscillator basis. The process of constructing the
Hamiltonian $\hat{H}$ of a given circuit and obtaining its spectrum has been
implemented in the scqubits package [28], as we explain in the following.

##  4 Implementation in scqubits

Figure 3: Flowchart outlining the implementation of computer-aided
quantization and numerical diagonalization in the scqubits package.

###  4.1 Variable transformation

The variable transform discussed in Section 2.2 is implemented
programmatically in order to automate the quantization of custom
superconducting circuits. Several steps described in terms of subcircuits are
equivalently formulated in more compact linear-algebra language applied to the
vector space of node variables.

The algorithm proceeds by dividing the vector space of node variables into
independent subspaces associated with the different variable categories. For
instance, the subspace corresponding to periodic variables is spanned by those
vectors $\mathbf{z}_{p}=(z_{1p},z_{2p},\ldots,z_{Np})$ which are compatible
with the condition of vanishing flux difference across all inductances.
Similar subspaces are constructed for free and frozen degrees of freedom
present in the circuit. The full transformation matrix777Note that the
transformation matrix generated can have any real valued coefficients. The
only condition is that $z_{ip}$ \- $z_{jp}\in\mathbb{Z}$ such that
$i\in\mathsf{S}_{p}$ and $j\notin\mathsf{S}_{p}$, which ensures that the
associated conjugate charge is periodic in a single Cooper pair. In some
instances, having non-integer elements in the transformation matrix can help
to increase the sparsity of the resultant matrix representation of the
Hamiltonian, which can significantly reduce the diagonalization time[17]. is
obtained by combining the vectors spanning the above three subspaces, and
adding additional linearly independent vectors, if needed, to complete the
basis.

Fig. 3 summarizes the steps involved in this algorithm which has been
implemented as a new module extending the scqubits package [28]. The Circuit
class added for this purpose facilitates the generation and manipulation of
symbolic circuit Hamiltonians by leveraging the Python package SymPy.

###  4.2 Numerical diagonalization

To obtain the spectrum of the circuit numerically, its Hamiltonian must be
represented as a matrix with respect to a suitable basis. We choose bases
$\\{|\mu_{i}\rangle\\}$ ($i=1,2\ldots$) for all individual degrees of freedom.
For the periodic degrees of freedom we employ the discrete charge eigenstates
as a convenient (though not necessarily ideal) basis. For extended degrees of
freedom, we either apply simple discretization with multi-point stencils
turning derivatives into finite differences, or employ the harmonic oscillator
basis associated with the potential energy obtained when omitting junction
terms.

In the ordinary diagonalization scheme, the basis of the full Hilbert space is
chosen as the Kronecker product
$|\mu_{1}\rangle\otimes|\mu_{2}\rangle\otimes\ldots$ For each degree of
freedom, individual truncation levels must be monitored to ensure convergence
of the procedure. By contrast, hierarchical diagonalization aims to mitigate
the fact that the dimension of the Hilbert space grows exponentially with the
number of nodes888Details of convergence depend not only on the node number,
but also critically on the specific parameter regime determined by the circuit
parameters.. The Circuit class offers an implementation of the hierarchical
diagonalization procedure following ideas inspired by Kerman [42]. In the
first stage, subsystems are diagonalized individually to obtain low-lying
energy eigenstates $\\{|E_{i}\rangle\\}.$ The basis of the full Hilbert space
used in this case is the Kronecker product
$|E_{1}\rangle\otimes|E_{2}\rangle\otimes\ldots$ Under favorable conditions of
moderate hybridization among subsystems, this procedure lowers the overall
Hilbert space dimension. Once either of the two diagonalization schemes has
been chosen, the full functionality of scqubits for visualizing spectra and
wavefunctions, performing parameter sweeps, etc. is available to explore the
behavior of the custom circuit.

As a concrete illustration of the hierarchical-diagonalization scheme,
consider the fluxonium coupler circuit (see Appendix A.2), composed of five
nodes and four degrees of freedom. Simulation of the full circuit Hamiltonian
generally calls for diagonalization of matrices of size $10^{6}\times 10^{6}$.
In suitable parameter regimes of practical interest, the system can be
meaningfully divided into weakly coupled sub-systems, and hierarchical
diagonalization proves to be a key step that decreases matrix sizes to
$10^{3}\times 10^{3}$ and dramatically reduces the required runtime by two
orders of magnitude.

##  5  Conclusions

We described a computer-aided framework to quantize general lumped-element
superconducting circuit systems. We showed how to systematically identify
periodic, extended, and non-dynamical (free and frozen) degrees of freedom in
a system. Eliminating the non-dynamical degrees of freedom leads to an
efficient representation of the Hamiltonian. This framework is implemented
programmatically in the new module circuit of scqubits Python library, which
can simulate a user-defined circuit described as a graph. Hierarchical
diagonalization was also implemented, which made numerical analysis for larger
circuits accessible.

Future work will further extend the scope of the module by incorporating
coherence-time calculations for custom circuits, handle the quantization and
numerical analysis of circuits in the presence of time-dependent magnetic
flux, and utilize the computation of perturbative corrections for hierarchical
diagonalization discussed in Ref. [42]. Additional opportunities for enhancing
the performance of our code lie in optimizing the choice of variable
transformation which is generally not completely fixed by the protocol
reported here. The remaining freedom may offer routes to increasing numerical
efficiency, especially if variable transformations can reduce explicit
coupling terms [17].

## Note added

During preparation of this manuscript, we were made aware of ongoing work in
Amir Safavi-Naeini’s group with a similar scope. Our two groups did not
discuss details of our projects, but agreed to post preprints on the arXiv
simultaneously.

We thank Andrew Houck, Zlatko Minev, Michel Devoret, Xinyuan You, Daniel Weiss
and Venkat Chandrasekhar for valuable discussions. This material is based upon
work in part supported by the U.S. Department of Energy, Office of Science,
National Quantum Information Science Research Centers, Co-design Center for
Quantum Advantage (C2QA) under contract number DE-SC0012704, and in part by
the Army Research Office under grant W911NF-1910016. For this work, IMS was
funded by the DOE, ZH and SPC were funded by the ARO, and JK was funded by ARO
and DOE. During the writing and editing of the manuscript, ZH was funded by
the U.S. Department of Energy, Office of Science, National Quantum Information
Science Research Centers, Superconducting Quantum Materials and Systems Center
(SQMS) under contract number DE-AC02-07CH11359. We thank the Air Force Office
of Scientific Research for supporting the scqubits package under grant
FA9550-20-1-0271.

##  Appendix A Examples

###  A.1 User-defined transformation for the KITE circuit

The following string encodes the circuit graph for the example KITE circuit in
Section 1. This string can be passed to the method from_yaml directly or
through a text file.

The parameters are (in units of GHz): $E_{C}=2.5$, $E_{L1}=0.23$,
$E_{L2}=0.36$, $E_{J}=5.9$, $E_{CJ}=6.6$. The following snippet sets the
truncation for the Hilbert space, chooses a more efficient transformation
following Ref. [29] and also chooses the closure branches.

⬇

In: trans_mat = np.array([[ -1, -1, -1, 1],

[ -1, -1, 3, 1],

[ -1, 3, -1, 1],

[ 3, -1, -1, 1]])*0.25

closure_branches = [kite.branches[-2], kite.branches[-1]]

kite.configure(transformation_matrix=trans_mat,
closure_branches=closure_branches)

kite.cutoff_ext_1 = 30

kite.cutoff_ext_2 = 30

kite.cutoff_ext_3 = 30

###  A.2 Hierarchical diagonalization for two fluxonia with a tunable coupler

The following string for the input file fluxonium-coupler.yaml encodes the
circuit graph for the fluxonium coupler [43, 44]:

The circuit can be diagonalized using hierarchical diagonalization by the
following lines of code:

⬇

In: import scqubits as scq

fluxonium_coupler = scq.Circuit.from_yaml("fluxonium-coupler.yaml",
ext_basis="harmonic", initiate_sym_calc=False, basis_completion="canonical")

system_hierarchy = [[1], [2], [3], [4]]

closure_branches = [fluxonium_coupler.branches[i] for i in [0, 4, -1]]

fluxonium_coupler.configure(closure_branches=closure_branches,
system_hierarchy=system_hierarchy,

subsystem_trunc_dims = [6, 6, 6, 6])

fluxonium_coupler.$\Phi$1 = 0.5 + 0.01768

fluxonium_coupler.$\Phi$2 = -0.2662

fluxonium_coupler.$\Phi$3 = -0.5 + 0.01768

fluxonium_coupler.cutoff_ext_1 = 110

fluxonium_coupler.cutoff_ext_2 = 110

fluxonium_coupler.cutoff_ext_3 = 110

fluxonium_coupler.cutoff_ext_4 = 110

eigs = fluxonium_coupler.eigenvals()

eigs - eigs[0]

⬇

Out:  array([0. , 0.03559404, 0.05819727, 0.09378676, 4.39927874, 4.43488613])

We have verified that these eigenvalues agree with results obtained
independently by the authors of [43].

##  Appendix B Role of conserved charges in the dynamics of superconducting
circuits

In this appendix, we show that conserved charges do not impact the physics of
the dynamical variables. To simplify the analysis, we will assume that there
is a single variable $\theta_{1}$ that is free. We perform a partial Legendre
transform with respect to the free variable to obtain the Routhian function

$\mathcal{R}(\dot{\boldsymbol{\theta}}_{pe},\boldsymbol{\theta}_{pe};Q_{1})={Q}_{1}\dot{{\theta}}_{1}-\mathcal{L}(\dot{\boldsymbol{\theta}},{\boldsymbol{\theta}}),$
(13)

where
${Q}_{1}=\partial_{\dot{{\theta}}_{pe}}\mathcal{L}(\dot{\boldsymbol{\theta}},\boldsymbol{\theta})$
and $\boldsymbol{\theta}_{pe}$ represents the variables excluding
$\theta_{1}.$ The Routhian can be written in the form

$\mathcal{R}(\dot{\boldsymbol{\theta}}_{pe},\boldsymbol{\theta}_{pe};Q_{1})=-\mathcal{L}_{0}(\dot{\boldsymbol{\theta}}_{pe},\boldsymbol{\theta}_{pe})+\frac{1}{2(\widetilde{\mathbf{C}}_{cpe})_{1,1}}\left(\sum_{i\neq
1}(\widetilde{\mathbf{C}}_{cpe})_{i,1}(\dot{\boldsymbol{\theta}}_{pe})_{i}\right)^{2}+\frac{d}{dt}G(Q_{1},t),$
(14)

where
$\mathcal{L}_{0}(\dot{\boldsymbol{\theta}}_{pe},\boldsymbol{\theta}_{pe})$ is
the Lagrangian of the system without the free degree of freedom, and

$G(Q_{1},t)=-\frac{Q_{1}}{(\widetilde{\mathbf{C}}_{cpe})_{1,1}}\left(\sum_{i\neq
1}(\widetilde{\mathbf{C}}_{cpe})_{i,1}(\boldsymbol{\theta}_{pe})_{i}\right)+\frac{Q_{1}^{2}}{2(\widetilde{\mathbf{C}}_{cpe})_{1,1}}t.$
(15)

The dynamics of the variables $\boldsymbol{\theta}_{pe}$ is determined by the
Euler-Lagrange equations
$\frac{d}{dt}\partial_{\dot{\boldsymbol{\theta}}_{pe}}\mathcal{R}=\partial_{\boldsymbol{\theta}_{pe}}\mathcal{R}.$
Since the conserved charge $Q_{1}$ enters through $G(Q_{1},t)$ as a total time
derivative in the Routhian, the value of $Q_{1}$ has no impact on the
dynamics. As a result, after obtaining the classical Hamiltonian, the charge
$Q_{1}$ can be set to zero. The generalization to more than one free variable
can be carried out straightforwardly.

## References

## References

  * [1] Gambetta J M, Chow J M and Steffen M 2017 npj Quantum Inf. 3 2
  * [2] Wendin G 2017 Rep. Prog. Phys. 80 106001
  * [3] Blais A, Grimsmo A L, Girvin S M and Wallraff A 2021 Rev. Mod. Phys. 93(2) 025005
  * [4] Brooks P, Kitaev A and Preskill J 2013 Phys. Rev. A 87(5) 052306
  * [5] Gyenis A, Di Paolo A, Koch J, Blais A, Houck A A and Schuster D I 2021 Phys. Rev. X Quantum 2(3) 030101
  * [6] Zhang H, Chakram S, Roy T, Earnest N, Lu Y, Huang Z, Weiss D K, Koch J and Schuster D I 2021 Phys. Rev. X 11(1) 011010
  * [7] Nakamura Y, Pashkin Y A and Tsai J S 1999 Nature 398 786–788
  * [8] Orlando T P, Mooij J E, Tian L, van der Wal C H, Levitov L S, Lloyd S and Mazo J J 1999 Phys. Rev. B 60(22) 15398–15413
  * [9] Koch J, Yu T M, Gambetta J, Houck A A, Schuster D I, Majer J, Blais A, Devoret M H, Girvin S M and Schoelkopf R J 2007 Phys. Rev. A 76 042319
  * [10] Manucharyan V E, Koch J, Glazman L and Devoret M 2009 Science 326 113–116
  * [11] Smith W C, Kou A, Xiao X, Vool U and Devoret M H 2020 npj Quantum Inf. 6 8
  * [12] Gyenis A, Mundada P S, Di Paolo A, Hazard T M, You X, Schuster D I, Koch J, Blais A and Houck A A 2021 Phys. Rev. X Quantum 2(1) 010339
  * [13] Nigg S E, Paik H, Vlastakis B, Kirchmair G, Shankar S, Frunzio L, Devoret M H, Schoelkopf R J and Girvin S M 2012 Phys. Rev. Lett. 108(24) 240502
  * [14] Bourassa J, Beaudoin F, Gambetta J M and Blais A 2012 Phys. Rev. A 86(1) 013814
  * [15] Solgun F, Abraham D W and DiVincenzo D P 2014 Phys. Rev. B 90(13) 134504
  * [16] Smith W C, Kou A, Vool U, Pop I M, Frunzio L, Schoelkopf R J and Devoret M H 2016 Phys. Rev. B 94(14) 144507
  * [17] Ding D, Ku H S, Shi Y and Zhao H H 2021 Phys. Rev. B 103(17) 174501
  * [18] Minev Z K, Leghtas Z, Mundhada S O, Christakis L, Pop I M and Devoret M H 2021 npj Quantum Inf. 7 131
  * [19] Yurke B and Denker J S 1984 Phys. Rev. A 29(3) 1419–1437
  * [20] Devoret M H 1995 Quantum fluctuations in electrical circuits Quantum Fluctuations, Les Houches, Session LXIII ed Reynaud S, Giacobino E and Zinn-Justin J (Elsevier Sci.) chap 10
  * [21] Vool U and Devoret M 2017 Int. J. Circuit Theory Appl. 45 897–934
  * [22] Burkard G, Koch R H and DiVincenzo D P 2004 Phys. Rev. B 69 064503
  * [23] Burkard G 2005 Phys. Rev. B 71(14) 144511
  * [24] Klots A R 2013 Superquantpackage https://github.com/andreyklots/SuperQuantPackage
  * [25] Parra Rodríguez A 2014 Software Tool to Analyse Superconducting Qubits Ph.D. thesis Queensland University of Technology
  * [26] Gely M F and Steele G A 2020 New Journal of Physics 22 013025
  * [27] Aumann P, Menke T, Oliver W D and Lechner W 2021 arXiv:2106.05342
  * [28] Groszkowski P and Koch J 2021 Quantum 5 583
  * [29] Smith W C, Villiers M, Marquet A, Palomo J, Delbecq M R, Kontos T, Campagne-Ibarcq P, Douçot B and Leghtas Z 2022 Phys. Rev. X 12(2) 021002
  * [30] Burkard G, Koch R H and DiVincenzo D P 2004 Phys. Rev. B 69(6) 064503
  * [31] Yurke B and Denker J S 1984 Phys. Rev. A 29(3) 1419–1437
  * [32] Girvin S M 2011 Circuit QED: superconducting qubits coupled to microwave photons (Oxford University Press)
  * [33] Rasmussen S, Christensen K, Pedersen S, Kristensen L, Bækkegaard T, Loft N and Zinner N 2021 Phys. Rev. X Quantum 2 040204
  * [34] Rymarz M, Bosco S, Ciani A and DiVincenzo D P 2021 Phys. Rev. X 11
  * [35] Egusquiza I L and Parra-Rodriguez A 2022 arXiv:2203.06167
  * [36] Parra-Rodriguez A 2021 arXiv:2104.09410
  * [37] Riwar R P and DiVincenzo D P 2022 npj Quantum Inf. 8 36
  * [38] Tinkham M 2004 Introduction to superconductivity (Courier Corporation)
  * [39] Thomas R E, Rosa A J and Toussaint G J 2016 The analysis and design of linear circuits (John Wiley & Sons)
  * [40] You X, Sauls J A and Koch J 2019 Phys. Rev. B 99 174512
  * [41] Bender C M and Orszag S A 1999 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory (Springer)
  * [42] Kerman A J 2020 arXiv:2010.14929 [quant-ph]
  * [43] Weiss D K, Zhang H, Ding C, Schuster D I and Koch J High-fidelity entangling gates for fluxonium qubits via flux modulation of a tunable coupler bulletin of the American Physical Society, APS March Meeting 2022, Chicago. Number 5 in Session T41 (2022) URL https://meetings.aps.org/Meeting/MAR22/Session/T41.5
  * [44] Zhang H, Ding C, Weiss D K, Schuster D I and Koch J Fast control of low-frequency fluxonium qubits bulletin of the American Physical Society, APS March Meeting 2022, Chicago. Number 4 in Session T41 (2022) URL https://meetings.aps.org/Meeting/MAR22/Session/T41.4