# Compressing multidimensional weather and climate data into neural networks

Langwen Huang
Department of Computer Science
ETH Zurich
<EMAIL_ADDRESS>
&Torsten Hoefler
Department of Computer Science
ETH Zurich
<EMAIL_ADDRESS>
###### Abstract

Weather and climate simulations produce petabytes of high-resolution data that
are later analyzed by researchers in order to understand climate change or
severe weather. We propose a new method of compressing this multidimensional
weather and climate data: a coordinate-based neural network is trained to
overfit the data, and the resulting parameters are taken as a compact
representation of the original grid-based data. While compression ratios range
from 300× to more than 3,000×, our method outperforms the state-of-the-art
compressor SZ3 in terms of weighted RMSE, MAE. It can faithfully preserve
important large scale atmosphere structures and does not introduce significant
artifacts. When using the resulting neural network as a 790× compressed
dataloader to train the WeatherBench forecasting model, its RMSE increases by
less than 2%. The three orders of magnitude compression democratizes access to
high-resolution climate data and enables numerous new research directions.

##  1 Introduction

Numerical weather and climate simulations can produce hundreds of terabytes to
several petabytes of data (Kay et al., 2015; Hersbach et al., 2020) and such
data are growing even bigger as higher resolution simulations are needed to
tackle climate change and associated extreme weather (Schulthess et al., 2019;
Schär et al., 2019). In fact, kilometer-scale climate data are expected to be
one of, if not the largest, scientific datasets worldwide in the near future.
Therefore it is valuable to compress those data such that ever-growing
supercomputers can perform more detailed simulations while end users can have
faster access to them.

Data produced by numerical weather and climate simulations contains
geophysical variables such as geopotential, temperature, and wind speed. They
are usually stored as multidimensional arrays where each element represents
one variable evaluated at one point in a multidimensional grid spanning space
and time. Most compression methods (Yeh et al., 2005; Lindstrom, 2014; Lalgudi
et al., 2008; Liang et al., 2022; Ballester-Ripoll et al., 2018) use an auto-
encoder approach that compresses blocks of data into compact representations
that can later be decompressed back into the original format. This approach
prohibits the flow of information between blocks. Thus larger blocks are
required to achieve higher compression ratios, yet larger block sizes also
lead to higher latency and lower bandwidth when only a subset of data is
needed. Moreover, even with the largest possible block size, those methods
cannot use all the information due to computation or memory limitations and
cannot further improve compression ratio or accuracy.

Figure 1: Diagram of the neural network structure: Coordinates are transformed
into 3D coordinates and fed to the Fourier feature layer (green). Then they
flow into a series of fully connected blocks (light blue) where each block
consists of two feed forward layers (dark blue) with batch norms (orange) and
a skip connection. Solid lines indicate flows of data in both compression
(training) and decompression (inference) processes, and dash lines correspond
to the compression-only process.

We present a new lossy compression method for weather and climate data by
taking an alternative view: we compress the data through training a neural
network to surrogate the geophysical variable as a continuous function mapping
from space and time coordinates to scalar numbers (Figure 1). The input
horizontal coordinates are transformed into three-dimensional Cartesian
coordinates on the unit sphere where the distance between two points is
monotonically related to their geodesic distance, such that the periodic
boundary conditions over the sphere are enforced strictly. The resulting
Cartesian coordinates and other coordinates are then transformed into Fourier
features before flowing into fully-connected feed forward layers so that the
neural network can capture high-frequency signals (Tancik et al., 2020). After
the neural network is trained, its weights are quantized from float32 into
float16 to get an extra 2× gain in the compression ratio without losing much
accuracy. This method allows compressing a large data block and squeezes
redundant information globally while still having low latency and high
bandwidth when retrieving decompressed values. It performs well in the high
compression ratio regime: it can achieve 300× to more than 3,000× in
compression ratio with better quality than the state-of-the-art SZ3 compressor
(Liang et al., 2022) in terms of weighted root mean squared error (RMSE) and
mean absolute error (MAE). The compression results preserve large-scale
atmosphere structures that are important for weather and climate predictions
and do not introduce artifacts, unlike SZ3. The neural network produced by our
method can be used as an IO-efficient dataloader for training machine learning
models such as convolutional neural networks (CNNs) and only lead to a limited
increase in test errors compared with training with original data. Although
equipped with Fourier features, error analysis shows that our method poorly
captures rarely occurring sharp signals such as hurricanes: it tends to smooth
out extreme values inside the center of a hurricane. It also suffers from slow
training speed, especially in the low compression ratio regime.

Besides high compression ratios, our method provides extra desirable features
for data analysis. Users can access any subset of the data with a cost only
proportional to the size of the subset because accessing values with different
coordinates is independent and thus trivially parallelizable on modern
computation devices such as GPUs or multi-core CPUs. In addition, the
functional nature of the neural network provides ”free interpolation” when
accessing coordinates that do not match existing grid points. Both features
are impossible or impractical to implement on traditional compressors with a
matching compression ratio.

###  1.1 Related work

#### Compression methods for multidimensional data

Existing lossy compression methods compress multidimensional data by
transforming it into a space with sparse representations after truncating, and
then quantize and optionally entropy-encode the resulting sparse data. As one
of the most successful methods in the area, SZ3 (Liang et al., 2022) finds
sparse representations in the space of locally spanning splines. It provides
an error-bounded compression method and can achieve a 400× compression ratio
in our test dataset. TTHRESH (Ballester-Ripoll et al., 2018) compresses data
by decomposing it into lower dimensional tensors. It performs well on
isotropic data such as medical images and turbulence data with compression
ratios of around 300×, but it is not the case for heavily stratified weather
and climate data where it either fails to compress or results in poor
compression ratios of around 3×. ZFP (Lindstrom, 2014) is another popular
compression method that provides a fixed-bitrate mode by truncating blocks of
orthogonally transformed data. It provides only low compression ratios that
usually do not exceed 10×. While there are lossless methods for
multidimensional data (Yeh et al., 2005), they cannot achieve compression
ratios more than 2× because scientific data stored in multidimensional
floating point arrays rarely has repeated bit patterns (Zhao et al., 2020).
SimFS (Girolamo et al., 2019) is a special lossless compression method that
compresses the simulation data by storing simulation checkpoints every N time
steps and, whenever needed, recreating the discarded result by re-running the
simulation starting from the nearest saved time step. It can reach arbitrary
high compression ratios by simply discarding the respective amount of
intermediate results, but reconstruction (decompression) may be expensive
because it always produces all the prognostic variables of a simulation. Table
1 summarizes the performance of the methods above when compressing weather and
climate data.

#### Neural representation

Using neural networks to represent multidimensional objects is gaining
popularity after the “spectral bias” of the vanilla MLP is resolved by Fourier
features (Tancik et al., 2020) or SIREN (Sitzmann et al., 2020). Previous
works of applying neural representation for compression focus on compressing
3D scenes (Mildenhall et al., 2021), images (Dupont et al., 2021), videos
(Chen et al., 2021), 3D volumetric data (Lu et al., 2021), and isotropic fluid
data (Pan et al., 2022). Instead of optimizing the compression ratio, it is
also possible to optimize training time and learned details by adding a
trainable embedding layer (Müller et al., 2022; Barron et al., 2021). Our work
differs because the target data is defined on a manifold instead of the
Euclidean space. Also, weather and climate data are typically larger,
smoother, and much more anisotropic. COIN++ (Dupont et al., 2022) and MSCN
(Schwarz & Teh, 2022) applied meta-learning methods to compress time slices of
climate data into latent vectors of several bits. The two methods can achieve
an asympotic compression ratio of 3,000× assuming infinite number of time
slices, yet the actual compression ratios are typically below 50× considering
the size of neural network weights.

Table 1: Comparison of methods for multidimensional compressing weather and climate data with the same RMSE bound. The decompression speed is measured in both continuous and random access patterns. SimFS’ decompression speed depends heavily on the compression ratio. Method | Comp. ratio | Comp. speed | Decomp. speed
---|---|---|---
|  |  | Cont. | Rand.
SimFS (Girolamo et al., 2019) | Arbitrary |  |  |
ZFP (Lindstrom, 2014) |  $<$ 10× |  |  |
TTHRESH (Ballester-Ripoll et al., 2018) |  $<$ 300× |  |  |
SZ3 (Liang et al., 2022) |  $<$ 400× |  |  |
NN (Ours) | 300× - 3,000× |  |  |

##  2 Method

A geophysical variable in weather and climate dataset can be seen as a scalar
function of four inputs $f(t,p,\psi,\phi)$ with time $t$, pressure $p$,
latitude $\psi$ and longitude $\phi$. Note that pressure is used in
meteorology as the coordinate for vertical levels while the height of an air
parcel with constant pressure becomes a variable. We propose using a neural
network with weights $\bm{\theta}$ as a compact representation of the function
$\hat{f}(t,p,\psi,\phi;\bm{\theta})$. As Figure 1 shows, the neural network is
trained with $N$ random samples of the coordinates
$\\{(t_{i},p_{i},\psi_{i},\phi_{i})|i\in[1,N]\\}$ as inputs, and evaluations
of the function $\\{f(t_{i},p_{i},\psi_{i},\phi_{i})|i\in[1,N]\\}$ at those
coordinates as targets, and use the mean squared error (MSE) as the loss
function. In practice, given a weather and climate data stored in a four-
dimensional array, one can implement the function $f(t,p,\psi,\phi)$ by
returning the interpolated value at the input coordinate. Since the aim of the
method is to “overfit” the neural network to the data, there is no need for
generalization to unseen inputs. Therefore we use all the data available in
the training process. After training, the weights $\bm{\theta}$ are quantized
from float32 into float16 and then stored as the compressed representation of
the original data. In the decompression process, the stored weights are cast
back to float32 and loaded into the neural network along with coordinates of
the data to produce decompressed values at those input coordinates. We remark
that more complex quantization (Gholami et al., 2021) and sparsification
(Hoefler et al., 2021) schemes could be used for further compressing of neural
network parameters that we explain in the following.

###  2.1 Neural network structure

Table 2: Hyper-parameters Symbol | Explanation
---|---
LR | Learning rate
$c_{t}$ | Scaling factor for time
$c_{p}$ | Scaling factor for pressure
$m$ | Number of Fourier features
$\sigma$ | The standard deviation of Fourier features
$d$ | Number of FCBlocks
$w_{i}$ | Width of the i-th FCBlock

We use the fully connected MLP with skip connections as the backbone of the
neural network (Figure 1). It contains a sequence of fully connected blocks
(FCBlocks), and each FCBlock contains two fully connected feed forward layers
with batch norms and a skip connection in the end. We use the Gaussian Error
Linear Unit (GELU) (Hendrycks & Gimpel, 2016) as the activation function in
feed forward layers to make the neural network a smooth function. The widths
of feed forward layers $m$ and the number of FCBlocks $d$ are configurable
hyperparameters (Table 2) that influence the compression ratio.

Firstly, the input coordinates are normalized and then transformed into
Fourier features. During the normalization process, time and pressure
coordinates are normalized by two constants:
$\hat{t}=t/c_{t},\hat{p}=p/c_{p}$. Latitude and longitude coordinates are
transformed into 3D Cartesian XYZ coordinates on the unit sphere:
$x=\cos(\psi)\cos(\phi),y=\cos(\psi)\sin(\phi),z=\sin(\psi)$ such that two
points close to each other on the sphere are also close to each other on the
XYZ coordinates. This transformation can enforce the periodic property along
longitude strictly. Then, the five coordinates
$\bm{v}=[\hat{t},\hat{p},x,y,z]$ are transformed into Fourier features
$\bm{\gamma}(\bm{v})=[\cos(\bm{B}v^{T})^{T},\sin(\bm{B}v^{T})^{T}]$ where
$\cos$ and $\sin$ are applied element-wise. $\bm{B}$ is a blocked diagonal
random matrix with $3m$ rows and 5 columns whose non-zero values are
independently and randomly sampled from a Gaussian distribution with 0 mean
and standard deviation $\sigma$. The $\bm{B}$ matrix is blocked in a way that
there are only interactions among $x,y$, and $z$:

$\bm{B}^{3m\times 5}=\begin{bmatrix}\bm{B}_{t}^{m\times 1}&&\\\
&\bm{B}_{p}^{m\times 1}&\\\ &&\bm{B}_{xyz}^{m\times 3}\end{bmatrix}.$

The values of the matrix $\bm{B}$ are fixed for a neural network and setting
it to trainable does not lead to a noticeable difference according to Tancik
et al. (2020).

As weather and climate data is anisotropic, some geophysical variables may
vary several magnitudes along the pressure and latitude dimensions. In that
case, the neural network only captures variations along that dimension while
ignoring variations along all other dimensions. Therefore, we scale the neural
network output to appropriate dynamic ranges at that coordinate. The scaling
function can be computed by interpolating mean and range along pressure and
latitude dimensions.

###  2.2 Sampling method

As we train the neural network using a stochastic optimization method, it is
important to feed the inputs randomly as required by the stochastic
optimization method to have a good convergence rate (Bottou, 2004). To achieve
that, we randomly sample coordinates for latitude, longitude, and time within
the data definition region. The target values on those coordinates are then
interpolated such that coordinates are not restricted to grid point
coordinates of the data. This makes our method independent from the input data
grid and helps the neural network learn the underlying smooth function from
the data instead of remembering values at each grid point and outputting
meaningless values anywhere else. We use a low discrepancy quasi-random
generator (Owen, 1998) to reduce the occurrence of closely clustered
coordinates as well as reduce regions that are not sampled in each batch which
also helps the neural network learn uniformly over the sampling region. In
addition, to avoid biased sampling over polar regions, the latitude $\psi$ and
longitude $\phi$ are sampled according to the metrics on the sphere (Figure
1): $\psi=\frac{1}{2}\pi-\cos^{-1}(1-2y),\phi=2\pi x$ where $(x,y)$ is drawn
from the low discrepancy quasi-random generator within range
$[0,1]\times[0,1]$. For pressure coordinates, we simply sample the exact grid
points because they are usually very sparse with only tens of grid points, and
interpolating on them introduces too much error.

##  3 Experimental Evaluation

In this section, we test our method by varying model structure and data
resolution, apply it to the training dataset of a CNN, and compare the
resulting test errors. The data used for experiments are extracted from the
ERA5 dataset (Hersbach et al., 2020) and interpolated to regular longitude-
latitude grids. When computing the mean of errors such as RMSE and MAE over
the longitude-latitude grid, we assign weights to each grid point that are
proportional to the cosine of latitudes to balance the denser polar regions
and the thinner equator region (see Appendix A.1). Throughout the experiments,
we fix the number of FCBlocks $d$ to 12, the number of Fourier features $m$ to
128, and the standard deviation of Fourier features $\sigma$ to 1.6. We use
the Adam optimizer (Kingma & Ba, 2014) with a learning rate of 3e-4. The
neural networks are trained for 20 epochs where each has 1,046,437,920
randomly sampled coordinates (matching the total number of grid points in
dataset 1) that are split into batches of 389,880 (matching the size of 3
horizontal slices in dataset 1) coordinates. All the experiments are performed
on a single compute node with 3 NVIDIA RTX 3090 GPUs.

###  3.1 Compressing data with different resolutions

Table 3: Comparison of compressing data with different spatial and temporal resolutions. Dataset | 1 | 2 | 3 | 4
---|---|---|---|---
Variable | Geopotential
Time period | 2016 | 2016 | 1979-2018 | 1979-2018
Spatial res. | 0.50∘ | 0.25∘ | 5.625∘ | 2.8125∘
Temporal res. | 24h | 24h | 1h | 1h
Pressure levels | 11 | 11 | 1 | 1
Data shape | (366,11,361,720) | (366,11,721,1440) | (350640,1,32,64) | (350640,1,64,128)
Data size | 3.9GB | 15.6GB | 2.7GB | 10.7GB

We test the performance of our method by compressing geopotential data from
the ERA5 dataset with various sizes and resolutions (Table 3). Geopotential is
selected because it is the most important variable for weather and climate
predictions (Rasp et al., 2020). We also compress the temperature data with
the same formats with resulting error plots in Appendix A.4. With datasets 1
and 2, we test the performance for compressing 4D high spatial resolution data
following the settings in Grönquist et al. (2021). Both datasets have 11
pressure levels which are meteorologically important (Ashkboos et al., 2022).
The year of 2016 is selected arbitrarily. We expect a similar result when
selecting a different year as shown in the error plots of compressing datasets
with four different years (1998, 2004, 2010, 2016) in Appendix A.3. On the
other hand, we test the performance for compressing high temporal resolution
data on the 500hPa pressure level with datasets 3 and 4. Datasets 3 and 4 also
can be used for training machine learning models for weather prediction under
the WeatherBench (Rasp et al., 2020) framework. Note that 0.25∘ is the highest
spatial resolution and 1 hour is the highest temporal resolution available in
the ERA5 dataset.

Table 4: Compression ratios of our method under different datasets and widths. |  | Dataset
---|---|---
Width | Model size | 1 | 2 | 3 | 4
512 | 13.86MB | 287.92 | 1,150.07 | 197.58 | 790.32
256 | 3.81MB | 1,046.91 | 4,181.83 | 718.43 | 2,873.73
128 | 1.15MB | 3,482.97 | 13,912.56 | 2,390.16 | 9,560.64
64 | 0.40MB | 9,913.17 | 39,597.75 | 6,802.84 | 27,211.37
32 | 0.18MB | 22,372.67 | 89,366.74 | 15,353.09 | 61,412.36

We use neural networks with different widths from 32 to 512 to create a
spectrum of compression ratios (Table 4). Because larger models take more than
20 hours to train on our devices which are not practical , we only test our
method on the high compression ratio regime ($>$200×). As a comparison, we
compress the same data using SZ3 (Liang et al., 2022) with max absolute error
tolerances from 10 to 2,000 and resulting compression ratios ranging from 10×
to 700×. Additionally, We tested ZFP (Lindstrom, 2014) and TTHRESH (Ballester-
Ripoll et al., 2018) but they cannot reach compression ratios larger than 10×
and thus are not included in the following analysis. For each compression
result, we compute weighted RMSE, weighted MAE, max absolute error, and
0.99999-quantile of absolute error. The error is measured in the unit of
$\text{m}^{2}/\text{s}^{2}$ (the same as geopotential). We also analyze error
distribution and perform a case study for outputs on Oct 5th 2016 when
hurricane Matthew emerged over the Caribbean Sea.

#### Error analysis

Figure 2: Weighted RMSE (upper left), weighted MAE (upper right),
0.99999-quantile of absolute error (lower left), and max absolute error (lower
right) between the compressed and original data.

As shown in Figure 2, our method can achieve at least 10 times more
compression ratios with the same weighted RMSE and MAE compared to SZ3. The
gap between the two methods is even larger in higher error regions where our
method can achieve more than 10,000× while SZ3 is capped at 600×. Comparing
the results between datasets 1 and 2, 3 and 4 where the only spatial
resolution differs, our method gets a free 4× gain in compression ratio
without hurting any of the error measurements. This indicates our method can
capture the real information underneath the continuous grid data. However, in
terms of max absolute error, it is worse than SZ3 even with similar
compression ratios due to different error distributions, as shown in Figure 3
(left). The majority of errors in our method concentrate near 0, but there
exists a small fraction that is relatively large. That small fraction is
around 0.001% according to the 0.99999-quantile of absolute error in Figure 2
where the two methods have similar errors at similar compression ratios.
Meanwhile, SZ3 distributes errors more evenly within the error tolerance. In
terms of the geographic distribution of errors, we calculate the mean and
standard deviation of errors over time and pressure levels from compressed
dataset 2 using our method (Figure 3 right). The mean values of errors of
approximately -10 are distributed evenly over the globe with exceptions over
mountainous areas. The standard deviation of errors shows meridional
difference where it grows from the equator to polar regions.

Figure 3: Left: Histogram of compression errors when compressing dataset 2
using our method with a width of 512 and SZ3 with an absolute error tolerance
of 1000. Right: mean and standard deviation of compression errors of our
method.

#### Case study

Geopotential at 500hPa is one of the most important variables for weather
prediction. Thus we select from Dataset 2 the 500hPa slice on Oct 5th 2016
when hurricane Matthew emerged over the Caribbean Sea and compare the
compression results of the two methods with the original data. We use the
biggest model with a compression ratio of 1,150× for our approach. For SZ3, we
select the one with an absolute error tolerance of 1,000 and a compression
ratio of 358× to match the 0.99999-quantile of absolute error.

In the global plot (Figure 4), our method clearly preserves the original data
better despite having a much higher compression ratio: locations and shapes of
minimum and maximum areas are nearly identical to the original ones, and the
data is perfectly continuous over the $0^{\circ}$ longitude line. The only
noticeable difference is that some high-frequency details are smoothed out. On
the SZ3 side, the average value over the equator region deviates from the
original data, the structure of minimum and maximum areas is broken due to
numerous artifacts, and there is discontinuity crossing the $0^{\circ}$
longitude line. Moreover, we can easily identify from data compressed by our
method two subtropical highs over the North Atlantic and the North Pacific
with slightly smoothed shapes compared to the original. The two subtropical
highs are important for predictions of precipitations and tracks of tropical
cyclones (hurricanes/typhoons). But from the SZ3 compressed data, we can only
identify one over the Pacific while the other breaks into small pieces, which
is meteorologically incorrect.

Figure 4: Global plots of geopotential at 500hPa on Oct 5th 2016 from original
dataset 2 (left), compressed data using our method with a width of 512
(center), and SZ3 compressed data with an absolute error tolerance of 1,000
(right). See Appendix A.2 for plots with varying compression ratios.

In a closer look at the hurricane region (Figure 5), our method failed to
capture the extreme values in the center of the hurricane although a smoothed
version of the surrounding structure is preserved. SZ3 better preserves
extreme value in the hurricane center but largely deviates from the original
in other regions. As the error histogram in Figure 3 shows, the errors of our
method in the center of the hurricane fall into the ”tail” of the distribution
that rarely occurs. Considering the hurricane center is just a few pixels
wide, our method can provide a lot more details with fewer data compared with
widely used lower resolution data that average out those extreme value
regions.

,

Figure 5: Local plots of geopotential at 500hPa on Oct 5th 2016 with hurricane
Matthew in the center. The plots use data from original dataset 2 (left),
compressed data using our method with a width of 512 (center), and SZ3
compressed data with an absolute error tolerance of 1,000 (right). See
Appendix A.2 for plots with varying compression ratios.

###  3.2 Compressing CNN training data

Training predictive machine learning models on large scientific datasets such
as weather and climate data is often I/O bound (Dryden et al., 2021). Our
method of compressing the dataset can potentially accelerate the training
process by replacing the traditional dataloader with a trained neural network
that generates training data on demand. Yet, it is unclear how point-wise
error estimates such as RMSE or MAE of the training data relate to changes in
test error of the trained neural network. Here, we apply our method to the
input data of a CNN to show its effectiveness in this scenario. We choose the
WeatherBench (Rasp et al., 2020) baseline model as the CNN. We use
geopotential at 500hPa in datasets 3 and 4 (Table 3) and temperature at 850hPa
with the same range and shape. We select the biggest available neural network
with a width of 512 to compress the geopotential and temperature data,
respectively. For SZ3, we set the absolute error tolerances to 900
$\text{m}^{2}/\text{s}^{2}$ for geopotential and 5.0 ${}^{\circ}\text{C}$ for
temperature to match the 0.99999-quantile of absolute error with our method.
The outcomes are measured from the weighted RMSE error of CNN calculated using
the test set.

Table 5: Weighted RMSE test error of CNN trained with original dataset and compressed dataset. C.R. represents compression ratio. | C.R. | Weighted RMSE error (test set)
---|---|---
Dataset 3 |  | Geopotential at 500 hPa ($\text{m}^{2}/\text{s}^{2}$) | Temperature at 850 hPa (${}^{\circ}\text{C}$)
Original |  | 632.9 | 2.906
NN compressed | 198× | 637.3 (+0.7%) | 2.944 (+1.3%)
SZ3 compressed | 71× | 650.6 (+2.8%) | 2.985 (+2.7%)
Dataset 4 |  |  |
Original |  | 688.8 | 2.834
NN compressed | 790× | 697.3 (+1.2%) | 2.888 (+1.9%)
SZ3 compressed | 106× | 702.9 (+2.0%) | 2.887 (+1.9%)

#### Results

The results in Table 5 shows that our method leads to small increases in
weighted RMSE test errors of not more than 1.9% while achieving high
compression ratios of 198× and 790× for dataset 3 and 4. On the other hand,
SZ3 results in larger test errors and lower compression ratios of 71× and 106×
probably due to the additional artifacts observed in Figure 4 and 5. The high
compression ratio and limited increase in test errors of our method make it
possible to feed machine learning models with much more data which will
eventually surpass the model trained with original data in performance.

###  3.3 Ablation study

We perform an ablation study by selectively disabling some features of the
model and comparing the resulting weighted RMSE loss between compressed data
and original data. We select output scaling, XYZ transform, Fourier feature,
skip connection, and batch norm as the targets because they are the special
parts added to the fully connected MLP to compress data effectively. We also
experiment changing the activation function from GELU to ReLU. We choose
dataset 1 in Table 3 as the inputs and a neural network with a width of 512
which results in a compression ratio of 288×.

Table 6:  Weighted RMSE between neural network compressed data and original data. In each experiment, parts of the model are disabled or replaced.Scaling means output scaling, XYZ means XYZ transform. | Scaling | GELU | XYZ | Fourier feature | Skip connection | Batch norm | WRMSE
---|---|---|---|---|---|---|---
1 |  |  |  |  |  |  | 818
2 | ✓ |  |  |  |  |  | 420
3 | ✓ |  | ✓ |  |  |  | 395
4 | ✓ |  | ✓ | ✓ |  |  | 347
5 | ✓ |  |  | ✓ |  |  | 346
6 | ✓ | ✓ |  |  |  |  | 303
7 | ✓ |  | ✓ |  | ✓ |  | 292
8 | ✓ | ✓ | ✓ |  | ✓ |  | 254
9 | ✓ | ✓ | ✓ |  | ✓ | ✓ | 158
10 | ✓ | ✓ | ✓ | ✓ | ✓ |  | 78
11 | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | 76

#### Results

As Table 6 shows, all of the features combined resulted in the lowest error.
Output scaling, Fourier feature, skip connections and batch norm can indeed
reduce the compression error compared with experiments without them. The GELU
activation function can help reduce the error compared with ReLU when combined
with output scaling and skip connection. The XYZ transform part does not have
much impact on the final result, however, it is important to enforce the
physically correct hard constraint on the boundary condition of output data.

##  4 Limitations and future works

Despite having high compression ratios and relatively low error, our method
suffers from long compression/training time. It takes approximately 8 hours to
train the network on dataset 2 (Table 3) with 4 NVIDIA RTX 3090 GPUs, whereas
SZ3 takes just several minutes to compress with 32 CPU cores. The slow
training process also limits us from using bigger models to get lower errors.
Our method also struggles to capture small-scale extreme values as discussed
in Section 3.1. While error analysis shows that it rarely occurs, such events
can relate to very important weather processes like hurricanes.

We are going to investigate potential ways to integrate techniques like
trainable embeddings (Müller et al., 2022; Barron et al., 2021) and meta-
learning (Finn et al., 2017; Nichol & Schulman, 2018) to accelerate training
and reduce extreme errors. In addition, we will expand our method to compress
broader weather and climate data such as satellite images and weather station
measurements into neural networks hoping to act as a complement to the
traditional data assimilation in the numerical weather prediction pipeline
(Bonavita et al., 2016).

##  5 Conclusion

We present a new method of compressing multidimensional weather and climate
data into neural network weights. We test our method on subsets of the ERA5
dataset with different resolutions, and compare the results with the best
available compressor SZ3 (Liang et al., 2022). Our method can achieve less
weighted RMSE and MAE and similar 0.99999-quantile of absolute errors compared
with SZ3 at compression ratios from 300× to 3,000×. While its max absolute
error is worse than SZ3, the high errors are shown to be rare occurring (less
than 0.001%). Through the case study on the hurricane Matthew event, our
method can reproduce high-quality data that preserve the large-scale structure
at a compression ratio of 1,150×, but it struggled to capture the extreme
values in the hurricane center. It is also effective in replacing neural
networking training data. Considering the compression ratio can be projected
to 3,000× or even more on high-resolution data, our method has the potential
to compress the whole ECMWF archive with hundreds of petabytes into tens of
terabytes of data, fitting on a single hard drive, that can drive machine
learning and other scientific applications.

#### Reproducibility Statement

We provide the source code implementing our method as well as instructions to
reproduce the experiments in Section 3. The data we used is also publicly
available through the download script in the source code. The source code is
available in https://github.com/spcl/NNCompression .

#### Ethics Statement

Our work can help reduce power consumption due to reduced data movement
whenever loading the compressed data. It can make more weather and climate
data easily available, which helps people better understand and mitigate
climate change and severe weather.

#### Acknowledgements

This project received EuroHPC-JU funding under grant MAELSTROM, No. 955513. We
thank the Swiss National Supercomputing Center (CSCS) for providing computing
resources.

## References

  * Ashkboos et al. (2022) Saleh Ashkboos, Langwen Huang, Nikoli Dryden, Tal Ben-Nun, Peter Dueben, Lukas Gianinazzi, Luca Kummer, and Torsten Hoefler.  Ens-10: A dataset for post-processing ensemble weather forecast.  _arXiv preprint arXiv:2206.14786_ , 2022.
  * Ballester-Ripoll et al. (2018) Rafael Ballester-Ripoll, Peter Lindstrom, and Renato Pajarola.  Tthresh: Tensor compression for multidimensional visual data.  6 2018.  doi: 10.1109/TVCG.2019.2904063.  URL http://arxiv.org/abs/1806.05952http://dx.doi.org/10.1109/TVCG.2019.2904063.
  * Barron et al. (2021) Jonathan T Barron, Ben Mildenhall, Matthew Tancik, Peter Hedman, Ricardo Martin-Brualla, and Pratul P Srinivasan.  Mip-nerf: A multiscale representation for anti-aliasing neural radiance fields.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 5855–5864, 2021.
  * Bonavita et al. (2016) Massimo Bonavita, Elias Hólm, Lars Isaksen, and Mike Fisher.  The evolution of the ecmwf hybrid data assimilation system.  _Quarterly Journal of the Royal Meteorological Society_ , 142(694):287–303, 2016.
  * Bottou (2004) Léon Bottou.  Stochastic learning.  In Olivier Bousquet and Ulrike von Luxburg (eds.), _Advanced Lectures on Machine Learning_ , Lecture Notes in Artificial Intelligence, LNAI 3176, pp. 146–168. Springer Verlag, Berlin, 2004.  URL http://leon.bottou.org/papers/bottou-mlss-2004.
  * Chen et al. (2021) Hao Chen, Bo He, Hanyu Wang, Yixuan Ren, Ser Nam Lim, and Abhinav Shrivastava.  Nerv: Neural representations for videos.  _Advances in Neural Information Processing Systems_ , 34:21557–21568, 2021.
  * Dryden et al. (2021) Nikoli Dryden, Roman Böhringer, Tal Ben-Nun, and Torsten Hoefler.  Clairvoyant Prefetching for Distributed Machine Learning I/O.  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC21)_. ACM, 11 2021\.
  * Dupont et al. (2021) Emilien Dupont, Adam Goliński, Milad Alizadeh, Yee Whye Teh, and Arnaud Doucet.  Coin: Compression with implicit neural representations.  _arXiv preprint arXiv:2103.03123_ , 2021.
  * Dupont et al. (2022) Emilien Dupont, Hrushikesh Loya, Milad Alizadeh, Adam Golinski, Y Whye Teh, and Arnaud Doucet.  Coin++: Neural compression across modalities.  _Transactions on Machine Learning Research_ , 2022(11), 2022\.
  * Finn et al. (2017) Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _International conference on machine learning_ , pp. 1126–1135. PMLR, 2017.
  * Gholami et al. (2021) Amir Gholami, Sehoon Kim, Zhen Dong, Zhewei Yao, Michael W. Mahoney, and Kurt Keutzer.  A survey of quantization methods for efficient neural network inference, 2021.  URL https://arxiv.org/abs/2103.13630.
  * Girolamo et al. (2019) Salvatore Di Girolamo, Pirmin Schmid, Thomas Schulthess, and Torsten Hoefler.  Simfs: A simulation data virtualizing file system interface.  1 2019.  doi: 10.48550/arxiv.1902.03154.  URL https://arxiv.org/abs/1902.03154.
  * Grönquist et al. (2021) Peter Grönquist, Chengyuan Yao, Tal Ben-Nun, Nikoli Dryden, Peter Dueben, Shigang Li, and Torsten Hoefler.  Deep learning for post-processing ensemble weather forecasts.  _Philosophical Transactions of the Royal Society A_ , 379(2194):20200092, 2021.
  * Hendrycks & Gimpel (2016) Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  _arXiv preprint arXiv:1606.08415_ , 2016.
  * Hersbach et al. (2020) Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al.  The era5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.
  * Hoefler et al. (2021) Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste.  Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks, 2021.  URL https://arxiv.org/abs/2102.00554.
  * Kay et al. (2015) Jennifer E Kay, Clara Deser, A Phillips, A Mai, Cecile Hannay, Gary Strand, Julie Michelle Arblaster, SC Bates, Gokhan Danabasoglu, James Edwards, et al.  The community earth system model (cesm) large ensemble project: A community resource for studying climate change in the presence of internal climate variability.  _Bulletin of the American Meteorological Society_ , 96(8):1333–1349, 2015.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Lalgudi et al. (2008) H.G. Lalgudi, Ali Bilgin, Michael Marcellin, and Mariappan Nadar.  Compression of multidimensional images using jpeg2000.  _Signal Processing Letters, IEEE_ , 15:393 – 396, 02 2008\.  doi: 10.1109/LSP.2008.922285.
  * Liang et al. (2022) Xin Liang, Kai Zhao, Sheng Di, Sihuan Li, Robert Underwood, Ali M Gok, Jiannan Tian, Junjing Deng, Jon C Calhoun, Dingwen Tao, et al.  Sz3: A modular framework for composing prediction-based error-bounded lossy compressors.  _IEEE Transactions on Big Data_ , 2022.
  * Lindstrom (2014) Peter Lindstrom.  Fixed-rate compressed floating-point arrays.  _IEEE Transactions on Visualization and Computer Graphics_ , 20, 08 2014.  doi: 10.1109/TVCG.2014.2346458.
  * Lu et al. (2021) Yuzhe Lu, Kairong Jiang, Joshua A. Levine, and Matthew Berger.  Compressive neural representations of volumetric scalar fields.  4 2021.  doi: 10.48550/arxiv.2104.04523.  URL http://arxiv.org/abs/2104.04523.
  * Mildenhall et al. (2021) Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Communications of the ACM_ , 65(1):99–106, 2021\.
  * Müller et al. (2022) Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller.  Instant neural graphics primitives with a multiresolution hash encoding.  1 2022.  doi: 10.48550/arxiv.2201.05989.  URL https://arxiv.org/abs/2201.05989.
  * Nichol & Schulman (2018) Alex Nichol and John Schulman.  Reptile: a scalable metalearning algorithm.  _arXiv preprint arXiv:1803.02999_ , 2(3):4, 2018\.
  * Owen (1998) Art B Owen.  Scrambling sobol’and niederreiter–xing points.  _Journal of complexity_ , 14(4):466–489, 1998\.
  * Pan et al. (2022) Shaowu Pan, Steven L. Brunton, and J. Nathan Kutz.  Neural implicit flow: a mesh-agnostic dimensionality reduction paradigm of spatio-temporal data.  4 2022.  doi: 10.48550/arxiv.2204.03216.  URL https://arxiv.org/abs/2204.03216.
  * Rasp et al. (2020) Stephan Rasp, Peter D Dueben, Sebastian Scher, Jonathan A Weyn, Soukayna Mouatadid, and Nils Thuerey.  Weatherbench: a benchmark data set for data-driven weather forecasting.  _Journal of Advances in Modeling Earth Systems_ , 12(11):e2020MS002203, 2020.
  * Schulthess et al. (2019) Thomas Schulthess, P. Bauer, Oliver Fuhrer, Torsten Hoefler, C. Schaer, and N. Wedi.  Reflecting on the goal and baseline for exascale computing: a roadmap based on weather and climate simulations.  _Computing in Science and Engineering (CiSE)_ , 21(1), Jan. 2019.  ISSN 1521-9615.
  * Schwarz & Teh (2022) Jonathan Richard Schwarz and Yee Whye Teh.  Meta-learning sparse compression networks.  _arXiv preprint arXiv:2205.08957_ , 2022.
  * Schär et al. (2019) Christoph Schär, Oliver Fuhrer, Andrea Arteaga, Nikolina Ban, Christophe Charpilloz, Salvatore Di Girolamo, Laureline Hentgen, Torsten Hoefler, Xavier Lapillonne, David Leutwyler, Katherine Osterried, Davide Panosetti, Stefan Rüdisühli, Linda Schlemmer, Thomas Schulthess, Michael Sprenger, Stefano Ubbiali, and Heini Wernli.  Kilometer-scale climate models: Prospects and challenges.  _Bulletin of the American Meteorological Society_ , 100(12), Dec. 2019.  Early Online Release.
  * Sitzmann et al. (2020) Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  6 2020.  doi: 10.48550/arxiv.2006.09661.  URL http://arxiv.org/abs/2006.09661.
  * Tancik et al. (2020) Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  6 2020.  doi: 10.48550/arxiv.2006.10739.  URL https://arxiv.org/abs/2006.10739.
  * Yeh et al. (2005) Pen-Shu Yeh, Philippe Armbruster, Aaron Kiely, Bart Masschelein, Gilles Moury, Christoph Schaefer, and Carole Thiebaut.  The new ccsds image compression recommendation.  In _2005 IEEE Aerospace Conference_ , pp. 4138–4145. IEEE, 2005\.
  * Zhao et al. (2020) Kai Zhao, Sheng Di, Xin Lian, Sihuan Li, Dingwen Tao, Julie Bessac, Zizhong Chen, and Franck Cappello.  Sdrbench: Scientific data reduction benchmark for lossy compressors.  In _2020 IEEE International Conference on Big Data (Big Data)_ , pp. 2716–2724. IEEE, 2020.

##  Appendix A Appendix

###  A.1 Weights for RMSE and MAE on the regular longitude-latitude grid

Both RMSE and MAE in this paper need to compute the average over the
horizontal dimensions. The pointwise squared and absolute error can be seen as
a function $f(\psi,\phi)$ mapping from latitude $\psi$, longitude $\phi$ to a
scalar as the error that location. Ideally, the mean of the pointwise error
over the sphere $\Omega$ is defined as:

$\int_{\Omega}{f(\psi,\phi)dS}/\int_{\Omega}{1dS}$ (1)

. Both integrations can be decomposed into sums of integrations over
individual cells $\Omega_{ij}$ in the regular longitude-latitude grid:

$\sum_{i,j}{\left(\int_{\Omega_{ij}}{f(\psi,\phi)dS}\right)}/\sum_{i,j}{S_{ij}}$
(2)

where $i$ denotes the cell index along the latitude dimension, $j$ denotes the
index along the longitude dimension, $S_{ij}$ means the area of cell
$\Omega_{ij}$. By approximating integrals using the product of cell center
value and cell areas, we have:

$\sum_{i,j}{\left(f(\psi_{i},\phi_{j})S_{ij}\right)}/\sum_{i,j}{S_{ij}}$ (3)

. If the grid size is large enough (which is the case for datasets in Table
3), the cell area can be further simplified:

$S_{ij}=r\cos(\psi_{i})\Delta\psi\Delta\phi$ (4)

where $r$ is the radius of the sphere, $\Delta\psi$ and $\Delta\phi$ are the
height and length of each cell. Canceling out $r$, $\Delta\psi$ and
$\Delta\phi$ as they are constant values in the regular longitude-latitude
grid, the formula (1) becomes:

$\sum_{i,j}{f(\psi_{i},\phi_{j})\cos(\psi_{i})}/\sum_{i,j}{\cos(\psi_{i})}$
(5)

. It can be seen as the sum of pointwise errors $f(\psi_{i},\phi_{j})$
weighted by normalized cosine of latitudes
$\cos(\psi_{i})/\sum_{i}{\cos(\psi_{i})}$.

###  A.2 Plots of geopotential on Oct 5th 2016

In this section, we add extra plots to complement the case study in Section
3.1.

Figure 6: Global plots of geopotential at 500hPa on Oct 5th 2016 from
compressed dataset 2 using our method with a range of widths. CR means
compression ratio. Figure 7: Local plots of geopotential at 500hPa on Oct 5th
2016 with hurricane Matthew in the center from compressed dataset 2 using our
method with a range of widths. CR means compression ratio. Figure 8: Global
plots of geopotential at 500hPa on Oct 5th 2016 from SZ3 compressed dataset 2
with a range of absolute error tolerances. CR means compression ratio. Figure
9: Local plots of geopotential at 500hPa on Oct 5th 2016 with hurricane
Matthew in the center from SZ3 compressed dataset 2 with a range of absolute
error tolerances. CR means compression ratio.

###  A.3 Error plots for compression experiments with different years

Figure 10: Weighted RMSE (upper left), weighted MAE (upper right),
0.99999-quantile of absolute error (lower left), and max absolute error (lower
right) between the compressed and original geopotential data. Similar datasets
with four different years (1998, 2004, 2010, 2016) are compressed.

###  A.4 Error plots for compression experiments with temperature data

Figure 11: Weighted RMSE (upper left), weighted MAE (upper right),
0.99999-quantile of absolute error (lower left), and max absolute error (lower
right) between the compressed and original temperature data.


[a]Ciaran Hughes

# Theory Overview of Heavy Exotic Spectroscopy

###### Abstract

This proceeding broadly overviews the current landscape of heavy exotic
spectroscopy. Such work includes the composition of certain $X$, $Y$, and $Z$
states, and proceeds to discuss tetraquarks made exclusively of four quarks.
The slides for this talk can be found at the following link.

##  1 What Defines a State to be Exotic?

In order to review exotic states, it is first necessary to know what defines a
state as exotic. By definition, being exotic is being non-conventional. So,
what do we mean by a conventional state111A state is rigorously defined to be
a pole singularity of the S-matrix.? A conventional state is defined to be any
state which we understand “well enough”222This definition also applies to
nuclear physics, where there also exists exotic nuclei.. Currently, mesons and
baryons are the only states deemed to be conventional, since they are
phenomenologically understood well enough in a potential model [1]. Everything
else, such as four-quark states, are classified as exotic. This definition of
conventional is ever changing, and if we understood a new class of hadrons
“well enough”, we would include this new class under the conventional
umbrella.

For illustrative purposes, the current (as of 2020) status of exotic states in
the charmonium sector is shown in Fig. 1 and 2 of Richard Lebed’s talk [2].
Taking all exotic states from all sectors, there have been 44 experimentally
observed exotic candidates, and 15 experimentally established333See the
non-$q\bar{q}$ mesons 2020 review from the PDG. Also, PDG defines an
established signal as being seen by two independent experiments with greater
than $5\sigma$ precision. exotic candidates. Colloquially, these exotic
candidates are called $XYZ$ states, but the PDG has defined a naming
convention444See Naming Scheme for Hadrons 2020 review from the PDG..

###  1.1 What Can The Exotic States Be?

Gell-Mann knew in 1964 that color confinement allowed a plethora of states
like $\bar{q}gq$, $\bar{q}\bar{q}qq$, $\bar{q}qqqq$, etc, where $q$ represents
a quark, and $g$ an excited gluon degree of freedom. For current exotic
states, there are four exotic configurations largely being considered:

  1. 1.

Hybrids. These $\bar{q}gq$ configurations exist when there is an excited gluon
in an octet representation which combines with a $\bar{q}q$ in an octet
representation. Phenomenological models include the constituent gluon model
and the flux tube model.

  2. 2.

Molecules. A state composed of two or more constituent hadrons, typically
bound by a Yakawa like force.

  3. 3.

Hadro-Quarkonium. A constituent hadron core with a quark and an anti-quark
cloud orbiting the core.

  4. 4.

Compact Tetraquarks. A four quark state composed of a tightly arranged $qq$
diquark and a separate $\bar{q}\bar{q}$ anti-diquark.

The goal of exotic spectroscopy is to correctly match which exotic
configurations map onto a specific experimental exotic signal. Certain models
start from the proposition that exotic signals arise from a single particular
exotic configuration, e.g, are solely a compact tetraquark. However, a single
configuration frequently explains only certain aspects of the data, but not
all. This has given rise to debate about which single configuration is most
apt, for example to describe the $X(3872)$. Such debate disappears if we start
from the view point that the most general solution consists of an admixture of
multiple configurations. It is then not so surprising that exotic candidates
exhibit complicated phenomena which can be explained by the mixture of
configurations.

##  2 The $\chi_{c1}(3872)$ aka $X(3872)$

The $X(3872)$ was the first exotic candidate discovered, in 2003 by BELLE, and
consequently is the most well known and studied. Its quantum numbers have been
established to be $J^{PC}=1^{++}$, $Q=0$, and $I^{G}=0^{+}$. Its decay modes
include $\mathcal{B}(\pi^{+}\pi^{-}J/\psi)>3.2\%$. Since charmonium states
like the $J/\psi$ have suppressed annihilation effects, the $X(3872)$ must
have at least $\bar{c}c$ valence components. Other modes include
$\mathcal{B}(D^{0}\bar{D^{*}})>30\%$, and
$\mathcal{B}(\rho^{0}J/\psi)\sim\mathcal{B}(\omega J/\psi)$. The isospin
violations in this latter decay enter during the decay process, and are caused
[5] by the 8 MeV mass difference between the $D^{0}\bar{D^{*}}$ and
$D^{+}\bar{D^{*-}}$ components of the $X(3872)$, i.e., isospin violations do
not arise in the state itself, but during the decay process.

One amazing piece of information about the $X(3872)$ is that its binding
energy lies exactly on the $D^{0}\bar{D^{*}}$ threshold, with a binding energy
of $0\pm 180$ keV. It also has a very narrow width of $\Gamma<1.2$ MeV. Being
so close to the two-particle S-wave threshold indicates cusp effects can
occur, however a pure cusp explanation of the data has been ruled out [6].
Both a virtual state and bound state are compatible with experimental data
[6]. In a potential model which only includes the $\bar{c}c$ component, the
nearest state, the $\chi_{c1}(2P)$, is over 100 MeV too high [7]. However,
when the $\bar{c}c$ is coupled to the two meson $D\bar{D}^{*}$ component, then
the energy of the $\chi_{c1}(2P)$ drops and is largely in line with the
experimental $X(3872)$ [7]. This is the reason that the $X(3872)$ is
considered to be an exotic candidate. Lattice QCD studies find a bound state
pole with a binding energy of $\mathcal{O}(10)$ MeV at unphysical pion masses,
and indicate that only $\bar{c}c$ and $D\bar{D}$ are important contributions
to the wavefunction. They disfavor a diquark interpretation. Still, lattice
QCD calculations with physical pion mass could find a different binding energy
and pole structure.

The most likely model of the $X(3872)$ is an appreciable mixture of a
$\bar{c}c$ and a $D\bar{D}$. There has been an interesting proposal that
highlights that we can experimentally measure the binding energy an order of
magnitude more precisely by using a triangle singularity with low energy
$D^{*0}\bar{D}^{0*}$s [8].

##  3 The $\psi(4230)$ aka $Y(4230)$ aka $Y(4260)$ State

In 2020, the PDG replaced the previously known $Y(4260)$ (found by BELLE in
2013) with the $Y(4230)$ (found by BESIII in 2017), due to the latter being
more precise. The quantum numbers of $Y$ states are defined to be
$J^{PC}=1^{--}$, $Q=0$, and $I^{G}=0^{-}$. The $Y(4230)$ decays to
$\pi^{+}\pi^{-}J/\psi$. Because charmonium states are annihilation suppressed,
the $Y(4230)$ must have a minimal $\bar{c}c$ valence component. However, the
most notable statement about the $Y(4230)$ is its lack of open-charm decays
into $D\bar{D}$ states - unlike excited charmonium states555Charmonium states
above threshold decay into open-charm thresholds $2-3$ orders of magnitude
more frequently than to hidden-charm thresholds.. The lack of open-charm
decays implies that the $Y(4230)$ has to contain more than just $\bar{c}c$,
making it an exotic candidate.

Any model for the $Y(4230)$ needs to explain the lack of open-charm decays.
One model suggests that the $Y(4230)$ is a $D\bar{D}_{1}(2420)$ molecular
state. The needed $65$ MeV binding energy is within reach of a potential model
[9], making it a viable option. As molecular states decay through their
constituents, this explains the $Y(4230)=D\bar{D}_{1}(2420)\to D\pi^{+}D^{*-}$
decay, as well as the $Y(4230)\to\gamma X(3872)$ decay (caused by a triangle
singularity). However, by studying the $\pi\pi$ and $K\bar{K}$ distributions
in the final states of $Y(4230)$ decays, [10] finds that the $Y(4230)$ needs
to have a sizable but not dominant $SU(3)$-flavor octet component. This means
that the $Y(4230)$ needs a $SU(3)$-flavor singlet component.

The $Y(4230)$ is also consistent with a hybrid scenario. Here, lattice QCD
energies of hybrids are around 180 MeV too high, but if systematic errors were
included then the masses are roughly consistent [1]. Additionally, with
potentials extracted from the lattice, pNRQCD finds a hybrid energy consistent
with the $Y(4230)$ [1]. Due to hybrid dynamics, decays of hybrids into S-wave
open charm are forbidden, making the S-P wave $D\bar{D}_{1}$ channel the
dominant decay process. Moreover, heavy quark spin symmetry is broken more in
hybrids than in quarkonia, occurring at
$\mathcal{O}(\Lambda_{\text{QCD}}/m_{Q})$. This could help explain the
observation of the $Y(4230)$ decaying into both the heavy quark spin triplet
$\pi^{+}\pi^{-}J/\psi$ and the spin singlet $\pi^{+}\pi^{-}h_{c}$ at
comparable rates.

Given the current experimental data, the most likely model for the $Y(4230)$
is a mixture of a $D\bar{D}_{1}$ molecule and a hybrid, causing a bound state
pole.

##  4 The $Z_{c}(3900)$ and $Z_{b}(10610)$ States

The $Z_{b}(10610)$ has decay modes $\mathcal{B}((B\bar{B}^{*})^{+})=86\%$,
$\mathcal{B}(\Upsilon(nS)\pi^{+})\sim 3\%$. For the $Z_{c}(3900)$ we have
$\mathcal{B}((D\bar{D}^{*})^{\pm})/\mathcal{B}(J/\psi(nS)\pi^{\pm})=6.2$.
Again, due to quarkonium annihilation effects being suppressed, and the $Z$
states being charged, the $Z$ states must have a minimal valence contribution
consisting of four quarks $\bar{Q}Qq_{1}q_{2}$. $Z$ states are defined by the
quantum numbers $Q=\pm 1,0$, $I^{G}=1^{+}$, and $J^{PC}=1^{+-}$.

For the $Z_{b}(10610)$, its large branching fraction into the S-wave
$B\bar{B}^{*}$ threshold is explained by its small binding energy of $3(3)$
MeV relative to that threshold. In fact, recent lattice QCD work [11]
extracted the potential between the $B$ and $\bar{B}^{*}$ and found sizable
attraction in the potential for small $r$. Certain parameterisations of the
potential allow for a virtual state - which would be identified as the
$Z_{b}(10610)$. It should not be surprising then that the $Z_{b}(10610)$ has
been identified as a virtual state when examining the experimental data by
[12]. As a shallow virtual state, it would be molecular [4], and potentially
be a mixture of multiple molecular components including $B\bar{B}^{*}$,
$\pi\Upsilon$, etc.

Similarly, the $Z_{c}(3900)$ is 13 MeV lower than the S-wave $D\bar{D}^{*}$
threshold. A lattice QCD [1] calculation includes diquark and two meson type
operators but does not find evidence of a bound state or narrow resonance.
However, [13] indicates that this lattice work is consistent with a virtual
state or broad resonance, which is also consistent with experimental data. To
distinguish between these cases, smaller bin sizes and better energy
resolution are needed in experimental data.

For the spatial distribution, the $Z_{c}(3900)$ is most likely a multi-
molecular system, where the mixing between the $\pi
J/\psi(\rho\eta_{c})-D\bar{D}^{*}$ is as important as the diagonal parts of
the potential [14].

##  5 Four Quark States Containing Two or More Heavy Quarks

As we have seen, the most understood exotic states are likely a mixture of
very different components, e.g., the $X(3872)$ is likely a mixture of
$\bar{c}c$ and $D\bar{D}^{*}$. This mixing can, and does, cause complicated
experimental phenomena, and this in turn makes it difficult to exclude regions
in the theoretical space of models. As such, we should first attempt to find
an exotic state that is solely composed of a single component and understand
this state fully. Afterwards, we can use this understanding when describing
the multi-component exotic candidates.

To make full use of our intuition, we should focus on bound states rather than
virtual or resonance states. A bound state is likely when at least two of the
quarks are heavy, as then there is a possibility of the formation of diquark
or anti-diquark pairs due to the non-relativistic behavior. The $QQ$ diquark
can be in a $\bar{3}/6$ representation with an attractive/repulsive potential.
Phenomenologically, the attractive $QQ$ diquark may be deep enough in the
heavy-quark potential to cause binding. For a heavy-light diquark $Qq$, there
is no known mechanism derived from QCD that can cause this structure to exist,
and so I will not discuss it further.

###  5.1 Bound $\bar{b}\bar{b}bb$ States

The above reasoning prompted my coauthors and I to search for a bound
$\bar{b}\bar{b}bb$ state using lattice QCD. We searched in three channels
($J^{PC}=0^{++}$, $1^{+-}$, and $2^{++}$) and used a full basis of two-meson
and diquark anti-diquark interpolating operators. We did not find evidence for
any bound state below the lowest S-wave thresholds, as shown in Fig. 9 of
[15]. When a signal is not found, a bound needs to be set on the likelihood
that you missed that signal. This is frequently not done in lattice QCD
calculations, but should be. In Fig. 11 [15] we show the probability that we
would have missed a bound state at a specific energy within our statistical
precision. As such, we have ruled out the possibility of a $\bar{b}\bar{b}bb$
compact tetraquark bound state to $5\sigma$.

###  5.2 Resonance $\bar{c}\bar{c}cc$ States

In 2020, LHCb found evidence [16] of a state which decays into 2$J/\psi$. Its
mass was around 700 MeV above the 2$J/\psi$ threshold, and was called the
$X(6900)$. Consequently, it would be composed of $\bar{c}\bar{c}cc$. Being
above threshold, this state is a resonance. A compact tetraquark with (anti-)
diquark constituents is the only likely model for the state due to the non-
relativistic behavior of charm quarks in this energy region. The $X(6900)$ is
very exciting as it is the first clear evidence for a state that can be
explained by a diquark model which can be connected to QCD.

There are a few noticeable features in Fig. 2 of the LHCb data. First, there
are three bumps, and second, there is a dip around 6.8 GeV. All these features
need to be explained by a model. As we have discussed, two body S-wave
thresholds are important for exotic states. Some important S-wave thresholds
are the $2\chi_{c0}$, $\chi_{c1}\chi{c0}$, and the $\Xi_{cc}\bar{\Xi}_{cc}$
open flavour threshold.

As these compact tetraquark states should be narrow, the broad width of the
experimental bumps hint that multiple states are contributing to the signal.
LHCb have not yet performed an amplitude analysis which could distinguish
nearby $J^{PC}$ states. What most models seem to agree on is that the
$X(6900)$ is some combination of multiple 2S $0^{++}$ states. Then the 1S
$0^{++}$ state could either be the first bump around $6.5$ GeV [17], or
alternatively could be below the $2J/\psi$ threshold (where LHCb does not have
data). If the 1S $0^{++}$ is below the $2J/\psi$ threshold, then the first
bump would be some combination of 1P $0^{-+}$ states [18]. The dip is
explained by destructive interference from the $2\chi_{c0}$ threshold becoming
accessible. Further, such states are likely in the bottom sector.

###  5.3 Bound $bb\bar{u}\bar{d}$ State

A $J^{P}=1^{+}$, $I=0$ exotic state composed of $bb\bar{u}\bar{d}$ quarks has
been predicted to be bound within a multitude of different theoretical
frameworks. Fig. 8 of [19] illustrates this consensus between lattice NRQCD
calculations, potentials extracted from lattice QCD, and phenomenological
model calculations. For a review talk of this state, and similar
$Q_{1}Q_{2}\bar{q_{3}}\bar{q_{4}}$ exotic states, see Anthony Francis’ talk.

In fact, there is a straightforward intuitive understanding why this state is
bound. First, start with the two $b$-quarks. Assume they are in the infinitely
heavy quark mass limit. Then they arrange themselves into an attractive
$\bar{3}$ diquark. The $\bar{u}$ and $\bar{d}$ then form a light quark cloud
that screens the $bb$ interaction. This is sufficient to form deep binding
around $100$ MeV. Adding in finite $b$-quark mass corrections does not change
this intuitive picture. Given the robustness of this prediction, a
confirmation is needed from experiment, potentially via the
$bb\bar{u}\bar{d}\to\Xi^{0}_{bc}\bar{p}$ or $\to B^{-}D^{+}\pi^{-}$ processes.
If found, this would be the first bound state that is composed exclusively of
four quarks, and would be a very useful leverage point in our understanding of
exotics.

###  5.4 The Bound $D_{s0}(2317)$ State

The $D_{s0}(2317)$ has many similar features to the $X(3872)$. First, prior to
its discovery, it was expected to be the $j=\frac{1}{2}^{+}$ state composed of
$c\bar{s}$ quarks, above threshold, and very broad through its decay to $DK$
[7]. However, experimentally, the $D_{s0}(2317)$ state is below the
$D^{0}K^{+}$ threshold and very narrow, in contrast to quark model
expectations. This is what makes it an exotic candidate. Secondly, when quark
potential model calculations include the two-meson $DK$ couplings to the
$c\bar{s}$ component, the eigenstate mass dramatically shifts downwards and in
line with experimental results. Multiple lattice QCD calculations [20] have
been performed, taking into account the various systematical errors, and find
a bound state pole nearby the experimental mass. These lattice QCD results
find that both the $c\bar{s}$ and two-meson $DK$ components are important, but
the diquark operators are not. For these reasons, the $D_{s0}(2317)$ is a
prototype to understanding the $X(3872)$. Yet, the $D_{s0}(2317)$ is easier to
study theoretically. We can use our understanding of the $D_{s0}(2317)$ as
input into understanding the $X(3872)$.

It is useful to ask [21] what perturbations can be performed to the
$D_{s0}(2317)$ system in order to understand exotics more quantitatively, and
the $X(3872)$ specifically? With lattice QCD, we can explore how exotics
change (both their mass, pole type, and residue) as we vary the quark mass - a
task that experiment cannot perform. In this way, we can supplement
experimental data with lattice QCD. Specifically, [21] advocates examining the
$D_{s0}\to D~{}K$ process as a function of strange quark mass $m_{s}$, and
using leading order heavy-quark effective field theory to describe heavy-light
systems, and leading order chiral perturbation theory for light-light
pseudoscalar mesons. If we take $m_{s}\to m_{s}-\epsilon$, then $M_{D}$
reduces by $-\epsilon$, while $M_{K}^{2}$ changes by $-B\epsilon$, with $B$
the leading order $\chi$-PT constant. Applying this logic to $D_{s0}\to
D~{}K$, we see that we can cause $D_{s0}$ to sit right at threshold, and even
to lie above threshold. Studying how the $D_{s0}$ changes would quantitatively
illuminate the role of the S-wave threshold effects on exotic states. It would
also be interesting to see the bound state pole move in the complex plane and
become a virtual/resonance state as the threshold effects change as a function
of $m_{s}$. We could even move the binding energy to be exactly that of the
$X(3872)$, and understand this arch-typical state more quantitatively.

Such a calculation is computationally feasible also. All the expensive pieces
can be reused, and only the cheap unphysical strange quark propagators would
need to be recomputed. Additionally, the systematic errors would be minimal.
Consequently, such a program is imminently needed and has few minimal
downsides.

##  6 Conclusion

Taken together, we have explored the latest exotic state developments. This
has included the experimentally established $X(3872)$, $Y(4230)$,
$Z_{c}(3900)$, $Z_{b}(10610)$, the recently discovered all charm tetraquark
$X(6900)$, and the $D_{s0}(2317)$. We also discussed the non-existence of a
bound all-bottom tetraquark, how lattice QCD can be used to study slight
perturbation of the $D_{s0}(2317)$ in order to quantitatively understand
S-wave threshold effects, and how imminent results are needed from experiments
to discover the bound $bb\bar{u}\bar{d}$ state.

We did not get to discuss all of the interesting exotic states currently in
the literature. However, given the physics we have learned, we can generalise
our knowledge to a multitude of other states. For example, the $Z_{b}(10650)$
is likely the $B^{*}\bar{B}^{*}$ spin-partner of the $Z_{b}(10610)$.
Similarly, the $Z_{c}(4020)$ is likely the $D^{*}\bar{D}^{*}$ spin-partner of
the $Z_{c}(3900)$. The heavy flavour equivalent of the $X(3872)$, the $X_{b}$,
has likely not been seen as the $\chi_{b1}(2P)$ is below the open-bottom
threshold. This is in contrast to the charm case, which highlights the
importance of the $\chi_{c1}(2P)$ for the $X(3872)$. The pentaquark $P_{c}$ is
likely a molecule of $\bar{c}c$ and a nucleon. LHCb have released preliminary
results about a exotic flavour $cs\bar{u}\bar{d}$ state, where the lattice QCD
HadSpec collaboration finds suggestions of this state in the $I=0$,
$J^{P}=0^{+}$ channel.

There are other established states that all models fit accurately. More
experimental data is needed to nullify some of the models describing these
states. The exotic candidates in this class are the $Y(4360)$ (which could be
the $D_{1}D^{*}$ partner of the $Y(4230)$), the $Y(4660)$ (which could be the
$D_{s}D_{s1}$strange partner of the $Y(4230)$), and the $Z(4430)$.

The future of exotic spectroscopy is bright. With the BELLE and BES upgrades,
the new PANDA experiment, and the continuous measurements from the LHC,
further data continues to be collected and analysed. Such work will continue
to find more exotic candidates, and further shed light into the exotic ways
Nature operates at the smallest of scales.

## References

  * [1] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C.E. Thomas et al., _The $XYZ$ states: experimental and theoretical status and perspectives_, _Phys. Rept._ 873 (2020) 1 [1907.07583].
  * [2] R. Lebed, “The Spectroscopy of Exotics in the Dynamical Diquark Model.” 2020.
  * [3] F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao and B.-S. Zou, _Hadronic molecules_ , _Rev. Mod. Phys._ 90 (2018) 015004 [1705.00141].
  * [4] I. Matuschek, V. Baru, F.-K. Guo and C. Hanhart, _On the nature of near-threshold bound and virtual states_ , 2007.05329.
  * [5] Z.-Y. Zhou and Z. Xiao, _Comprehending Isospin breaking effects of $X(3872)$ in a Friedrichs-model-like scheme_, _Phys. Rev. D_ 97 (2018) 034011 [1711.01930].
  * [6] LHCb collaboration, _Study of the lineshape of the $\chi_{c1}(3872)$ state_, _Phys. Rev. D_ 102 (2020) 092005 [2005.13419].
  * [7] J. Ferretti, G. Galatà and E. Santopinto, _Interpretation of the X(3872) as a charmonium state plus an extra component due to the coupling to the meson-meson continuum_ , _Phys. Rev. C_ 88 (2013) 015207 [1302.6857].
  * [8] F.-K. Guo, _Novel Method for Precisely Measuring the $X(3872)$ Mass_, _Phys. Rev. Lett._ 122 (2019) 202002 [1902.11221].
  * [9] X.-K. Dong, Y.-H. Lin and B.-S. Zou, _Prediction of an exotic state around 4240 MeV with $J^{PC}=1^{-+}$ as C-parity partner of Y(4260) in molecular picture_, _Phys. Rev. D_ 101 (2020) 076003 [1910.14455].
  * [10] Y.-H. Chen, L.-Y. Dai, F.-K. Guo and B. Kubis, _Nature of the $Y(4260)$: A light-quark perspective_, _Phys. Rev. D_ 99 (2019) 074016 [1902.10957].
  * [11] S. Prelovsek, H. Bahtiyar and J. Petkovic, _$Z_{b}$ tetraquark channel from lattice QCD and Born-Oppenheimer approximation_, _Phys. Lett. B_ 805 (2020) 135467 [1912.02656].
  * [12] Q. Wang, V. Baru, A. Filin, C. Hanhart, A. Nefediev and J.-L. Wynen, _Line shapes of the $Z_{b}(10610)$ and $Z_{b}(10650)$ in the elastic and inelastic channels revisited_, _Phys. Rev. D_ 98 (2018) 074023 [1805.07453].
  * [13] M. Albaladejo, P. Fernandez-Soler and J. Nieves, _$Z_{c}(3900)$ : Confronting theory and lattice simulations_, _Eur. Phys. J. C_ 76 (2016) 573 [1606.03008].
  * [14] P.G. Ortega, J. Segovia, D.R. Entem and F. Fernández, _The $Z_{c}$ structures in a coupled-channels model_, _Eur. Phys. J. C_ 79 (2019) 78 [1808.00914].
  * [15] C. Hughes, E. Eichten and C.T.H. Davies, _Searching for beauty-fully bound tetraquarks using lattice nonrelativistic QCD_ , _Phys. Rev. D_ 97 (2018) 054505 [1710.03236].
  * [16] LHCb collaboration, _Observation of structure in the $J/\psi$-pair mass spectrum_, _Sci. Bull._ 65 (2020) 1983 [2006.16957].
  * [17] M. Karliner and J.L. Rosner, _Interpretation of structure in the di- $J/\psi$ spectrum_, _Phys. Rev. D_ 102 (2020) 114039 [2009.04429].
  * [18] J.F. Giron and R.F. Lebed, _Simple spectrum of $c\bar{c}c\bar{c}$ states in the dynamical diquark model_, _Phys. Rev. D_ 102 (2020) 074003 [2008.01631].
  * [19] L. Leskovec, S. Meinel, M. Pflaumer and M. Wagner, _Lattice QCD investigation of a doubly-bottom $\bar{b}\bar{b}ud$ tetraquark with quantum numbers $I(J^{P})=0(1^{+})$_, _Phys. Rev. D_ 100 (2019) 014503 [1904.04197].
  * [20] G.K. Cheung, C.E. Thomas, D.J. Wilson, G. Moir, M. Peardon and S.M. Ryan, _$DK$ $I=0,$ $D\bar{K}$ $I=0,1$ scattering and the $D_{s0}^{\ast}(2317)$ from lattice QCD_, 2008.06432.
  * [21] E. Eichten and C. Hughes, _Exploring S-Wave Threshold Effects in QCD: A Heavy-Light Approach_ , _Phys. Lett. B_ 802 (2020) 135250 [1911.02024].