# Wavelets Are All You Need for Autoregressive Image Generation

Wael Mattar, Idan Levy, Nir Sharon and Shai Dekel

(28th June 2024)

Abstract. In this paper, we take a new approach to autoregressive image
generation that is based on two main ingredients. The first is wavelet image
coding, which allows to tokenize the visual details of an image from coarse to
fine details by ordering the information starting with the most significant
bits of the most significant wavelet coefficients. The second is a variant of
a language transformer whose architecture is re-designed and optimized for
token sequences in this ‘wavelet language’. The transformer learns the
significant statistical correlations within a token sequence, which are the
manifestations of well-known correlations between the wavelet subbands at
various resolutions. We show experimental results with conditioning on the
generation process.

##  1 Introduction

The generation of high-resolution visual information is certainly one of the
most remarkable achievements of modern-age artificial intelligence. One of the
prominent methods is diffusion-based models (1, 2, 3, 4, 5). In essence,
diffusion models attempt to learn inversions of ill-posed operators, such as
additive Gaussian noise, blurring, etc., so an image may be generated from
random noisy or blurry seeds. One then enforces various conditions on images
created as the time steps of the inversion process so that the final generated
image may correspond to a given text prompt.

Another line of research is designing autoregressive models, in an attempt to
leverage on the powerful Large Language Models (LLMs) (6, 7, 8). These models
(9, 10) provide different methods to convert the image pixel representation to
a series of visual tokens and then apply generative language techniques.

In this paper, we refine this line of research and provide a mathematically
robust approach to the autoregressive image generation process. To this end,
we reach out to a classic technique in image processing, specifically, wavelet
image coding (11, 12, 13). Wavelets (14, 15, 16) are one of the main tools of
modern approximation theory for nonlinear, adaptive approximation. The various
wavelet transforms provide the means to transform an image into a
representation that captures the essence of the visual information in a sparse
way. Typically, the significant wavelet coefficients are a small fraction of
the coefficients and represent important edge and texture information, while
the insignificant coefficients with small absolute values are associated with
smooth regions of the image. The goal of wavelet image compression is then to
efficiently store the information of only the significant coefficients. In
fact, the underlying method of the popular JPEG image compression algorithm
(17), invented in the 80s, contains many elements of wavelet coding, where a
local Discrete Cosine Transform, a precursor of wavelets, is used. However, in
this paper, we leverage the progressive wavelet compression technique, a more
advanced form of image compression. It creates a bit-stream where every bit
corresponds to the next most important piece of visual information. Since we
are generating images rather than decoding them from a compressed file, there
is no need to create actual binary bit-streams, and using a ‘wavelet language’
of a limited number of tokens is sufficient.

Thus, our new approach to autoregressive image generation is based on two main
ingredients. The first is progressive wavelet image coding, which allows to
tokenize the visual information of an image from coarse to fine details. This
is achieved using only 7 tokens, by ordering the information starting with the
most significant bits of the most significant wavelet coefficients. The second
ingredient is a variant of an NLP decoder-only transformer (6, 7, 8) whose
architecture was re-designed and optimized for token sequences in this
‘wavelet language’. The transformer learns the significant statistical
correlations within a token sequence, which are the manifestations of well-
known correlations between the wavelet subbands at various resolutions (18,
19, 20). During inference, this allows the generation of visually meaningful
images from an initial random seed generated from sampled from the
distribution of the scaling function coefficients at the lowest resolution.

Using the wavelet autoregressive approach, where the ‘wavelet language’
contains only 7 tokens, provides many attractive features. The length of the
token sequences during training or inference can be flexible, where longer
sequences imply more detailed or higher-resolution images. Guiding the
generative process using a class affiliation or text prompting is easily
achieved by concatenating the corresponding vector representations to the
tokens’ vector representation of dimension 7. Stochastic control using simple
transformer inference techniques, allows to create from one textual prompt a
diversity of corresponding images. Furthermore, since each token is associated
with the local support in the image domain of the corresponding wavelet, one
can switch the guidance during the generative process to allow different
prompting for different regions.

Our paper is organized as follows. We begin with a review of related work in
Section 2. In Section 3, we review wavelet image coding and explain how one
may extract from the classical theory the ability to tokenize the visual
information of images. In Section 4, we focus on components of language
transforms that we redesign to serve our special wavelet language. We further
provide several methods that can be used to direct the generation process
under certain conditions: class label and/or textual prompt. In Section 5, we
provide experimental results. Finally, in Section 6, we discuss possible
future applications of our method, such as multi-modality generation and
compositions of blobs.

##  2 Related work

In this section, we first review the current state of the art in image
generation. We then review some methods that apply wavelets as a frequency
decomposition backbone for various aspects of style transfer, acceleration,
and optimization of existing image generation methods, etc.

Currently, many commercial solutions apply diffusion-based models (1, 2, 3, 4,
5, 21). In essence, diffusion models learn inversions of ill-posed operators,
such as additive Gaussian noise, blurring, etc., so images may be generated
from random noisy or blurry seeds. One then enforces various conditions on
images created through the time steps of the inversion process so that the
final generated image may correspond to a given text prompt.

The methods of VQGAN (9) and DALL-E (10) along with (22, 23) utilize a visual
tokenizer to discretize images into grids of 2D tokens, which are then
flattened to a 1D sequence for autoregressive learning, mirroring the process
of sequential language modeling. For example, in (10) a discrete variational
autoencoder is trained to compress each $256\times 256$ RGB image into a
$32\times 32$ grid of image tokens, where each such token can assume 8192
possible values. This creates a relatively short context sequence of $1024$
tokens, but with a vocabulary of 8192 word tokens. The TiTok method, recently
introduced in (24), shows how to compose Vision Transformers with the Vector-
Quantization method to arrive at an autoregressive method that may use only 32
tokens. In comparison, our method uses only 7 tokens for any image resolution
and any level of fine detail generation.

In contrast, to the raster-scan “next-token prediction”, the method of (25)
provides autoregressive learning on images as “next-scale prediction” or
“next-resolution prediction”

Some methods, such as (26, 27), use wavelets as means for frequency
decomposition representations for image inpainting, style transfer, and
generative adversarial network methods. Some works propose to use wavelets as
part of diffusion methods (28, 29) to speed up the diffusion approach by
applying the denoising process in the wavelet regime.

##  3 Elements of Wavelet Image Coding

In this section, we review some elements of wavelet image coding (11, 12, 13)
that we use for our generative method. Essentially, we are interested in the
process that takes an image in its raw pixel form as input and generates a
sequence of tokens that capture its visual details. The structure of the
sequence from coarse to fine details is achieved by ordering the information
starting with the most significant bits of the most significant wavelet
coefficients. En par with wavelet coding, we also have a goal to create token
sequences that are as short as possible. This creates shorter contexts for the
transformer decoder and improves its performance. However, in our generative
application, we are not concerned with efficient encoding of the token
sequence to a compressed stream of bits. This is a subtopic of information
theory and is typically implemented using arithmetic encoders.

###  3.1 Wavelet Transforms

A univariate wavelet system (14, 16) is a family of real functions
$\left\\{{\psi_{j,k}\thinspace:\thinspace\thinspace(j,k)\in{\mathbb{Z}}^{2}}\right\\}$
in $L_{2}({\mathbb{R}})$, built by dilating and translating a unique mother
wavelet function $\psi$

$\psi_{j,k}(x):=2^{-j/2}\psi(2^{-j}x-k),$

where the mother wavelet typically has compact support (or fast decay) and has
$r$ vanishing moments

$\int_{\mathbb{R}}{x^{k}\psi(x)dx=0},\qquad k=0,1\ldots,r-1.$ (1)

Wavelet systems can be constructed to serve a basis of $L_{2}(\mathbb{R})$. To
facilitate applications, one then also constructs a dual $\tilde{{\psi}}$ of
$\psi$, where
$\langle\psi_{j,k},\tilde{\psi}_{j^{\prime},k^{\prime}}\rangle=\delta_{j,j^{\prime}}\delta_{k,k^{\prime}}$,
so that for each $f\in L_{2}(\mathbb{R})$,

$f=\sum\limits_{j,k}{\langle f,\tilde{{\psi}}_{j,k}\rangle\psi_{j,k}}.$

For special choices of $\psi$ , the set $\left\\{{\psi_{j,k}}\right\\}$ forms
an orthonormal basis for $L_{2}({\mathbb{R}})$ and then,
$\psi=\tilde{{\psi}}$.

Usually, one starts the construction of a wavelet system from a Multi-
Resolution Analysis (MRA) generated by a scaling function $\varphi\in
L_{2}\left({\mathbb{R}}\right)$ that satisfies a two-scale equation

$\varphi=\sum\limits_{k}{a_{k}\varphi\left({2\cdot-k}\right)}.$

One then sets

$V_{j}=\overline{span}\left\\{{\varphi_{j,k}:=2^{{-j}\mathord{\left/{\vphantom{{-j}2}}\right.\kern-1.2pt}2}\varphi\left({2^{-j}\cdot-k}\right):k\in{\mathbb{Z}}}\right\\},\quad
j\in\mathbb{Z},$

which implies (under certain mild conditions)

$....V_{2}\subset V_{1}\subset V_{0}\subset V_{-1}\subset V_{-2}...,\quad\cap
V_{j}=\left\\{0\right\\},\quad\cup_{j}V_{j}=L_{2}({\mathbb{R}}).$

Again, to facilitate applications, one may also construct a dual
$\tilde{{\varphi}}$ of $\varphi$, where
$\langle\varphi_{0,k},\tilde{\varphi}_{0,k^{\prime}}\rangle=\delta_{k,k^{\prime}}$,
so that for each $f\in V_{j}$,

$f=\sum\limits_{k}{\langle f,\tilde{{\varphi}}_{j,k}\rangle\varphi_{j,k}}.$

Equipped with the MRA, one then proceeds to construct the wavelet $\psi$ such
that $W_{j}:=\overline{span}\left\\{{\psi_{j,k}:k\in{\mathbb{Z}}}\right\\}$,
with $V_{j+1}+W_{j+1}=V_{j}$. A classic example for an orthonormal MRA and
wavelet system where $\varphi=\tilde{\varphi}$ and $\psi=\tilde{\psi}$, are
the Haar scaling function and Haar wavelet

$\varphi\left(x\right):=\left\\{{{\begin{array}[]{*{20}c}{1,\hfill}&{x\in\left[{0,1}\right],\hfill}\\\
{0,\hfill}&{\mbox{else}.\hfill}\\\
\end{array}}}\right.\quad\psi\left(x\right):=\left\\{{{\begin{array}[]{*{20}c}1\hfill&{x\in\left[{0,\frac{1}{2}}\right),\hfill}\\\
{-1\hfill}&{x\in\left[{\frac{1}{2},1}\right],\hfill}\\\
0\hfill&\mbox{else.}\hfill\\\ \end{array}}}\right.$

The bivariate Haar system (see below) is a good choice when working with
piecewise constant images, such as the MNIST handwritten digits (30). For some
of our experiments, we use a famous wavelet system from the
Cohen–Daubechies–Feauveau (CDF) family of wavelets (14), which is sometimes
termed bior4.4 ($r=4$ in (1)) or [9,7] in the signal processing community (the
supports of the scaling functions and wavelets, as well as the lengths of the
associated filters, are 9 and 7). The generating functions of the bior4.4 are
depicted in Figure 1.

Figure 1: The CDF [9,7] wavelet system (figure reproduced from (31)).

The wavelet model can be easily generalized to any dimension, via a tensor
product of the wavelet and the scaling functions. Assume that the univariate
dual scaling functions $\varphi,\tilde{\varphi}$ and dual wavelets
$\psi,\tilde{\psi}$, are given. Then, a wavelet bivariate basis is constructed
using three types of basic wavelets

$\psi^{1}(x_{1},x_{2}):=\varphi(x_{1})\psi(x_{2}),\quad\psi^{2}(x_{1},x_{2}):=\psi(x_{1})\varphi(x_{2}),\quad\psi^{3}(x_{1},x_{2}):=\psi(x_{1})\psi(x_{2}),$
$\tilde{\psi}^{1}(x_{1},x_{2}):=\tilde{\varphi}(x_{1})\tilde{\psi}(x_{2}),\quad\tilde{\psi}^{2}(x_{1},x_{2}):=\tilde{\psi}(x_{1})\tilde{\varphi}(x_{2}),\quad\tilde{\psi}^{3}(x_{1},x_{2}):=\tilde{\psi}(x_{1})\tilde{\psi}(x_{2}).$

The bivariate wavelet transform of $f\in L_{2}({\mathbb{R}}^{2})$, in terms of
the bivariate wavelet tensor basis

$\psi^{e}_{j,k}:=2^{-j}\psi^{e}(2^{-j}\cdot-k),\quad\tilde{\psi}^{e}_{j,k}:=2^{-j}\tilde{\psi}^{e}(2^{-j}\cdot-k),\quad
e=1,2,3,j\in\mathbb{Z},k\in\mathbb{Z}^{2},$

is then

$f=\sum_{e=1,2,3,j\in\mathbb{Z},k\in\mathbb{Z}^{2}}{\langle
f,\tilde{\psi}^{e}_{j,k}\rangle\psi_{j,k}^{e}}.$

The bivariate wavelet decomposition can thus be interpreted as a signal
decomposition in a set of three spatially oriented frequency subbands:
$LH(e=1)$ detects horizontal edges; $HL$ ($e=2)$ detects vertical edges and
$HH$ ($e=3)$ detects diagonal edges.

Under the assumption that $\psi$ and $\tilde{\psi}$ are compactly supported
(or have fast decay), a wavelet coefficient $\langle
f,\tilde{\psi}^{e}_{j,k}\rangle$ at a scale $j$ represents the information
about the function in the spatial region in the neighborhood of
$2^{j}k,\thinspace k\in{\mathbb{Z}}^{2}$. At the next finer scale $j-1$, the
information about this region is represented by four wavelet coefficients,
which are described as the children of $\langle
f,\tilde{\psi}^{e}_{j,k}\rangle$. This leads to a natural tree structure
organized in a quad tree structure of each of the three subband types as shown
in Figure 2. As $j$ decreases, the child coefficients add finer and finer
details into the spatial regions occupied by their ancestors.

Figure 2: Wavelet coefficient tree structure across the subbands (MRA
decomposition).

In image processing, one uses the Discrete Wavelet Transform (DWT). It works
by initially assuming that the image pixels
$\\{f_{k}=f_{k_{1},k_{2}}\\}_{k1,k2=1}^{M}$ are good approximants of the
projections on the shifts of the dual scaling function with the underlying
function $f$

$f_{k}\approx\langle f,\tilde{\varphi}_{0,k}\rangle.$

With these coefficients as input, one uses the DWT to compute coefficients
down to some predefined low-resolution $m$. For simplicity, we may assume that
$M=2^{m}$ and that we use the DWT to compute

$\\{\langle f,\tilde{\varphi}_{m,k}\rangle\\},\quad\\{\langle
f,\tilde{\psi}^{e}_{j,k}\rangle\\},\quad 1\leq j\leq m,\quad e=1,2,3.$ (2)

Wavelet representations are considered very efficient for image compression
(11, 12, 13). The edge information typically constitutes a small portion of a
typical image, while the dual wavelet coefficients have a large absolute value
only if edges intersect the support of the corresponding dual wavelets.
Consequently, the image can be approximated well using a few significant
wavelet coefficients. A clear statistical structure also follows: large/small
values of wavelet coefficients tend to propagate through the scales of the
quadtrees depicted in Figure 2. As an example, a sparse wavelet representation
of a $512\times 512$ fishing boat image and a compressed version of it are
shown in Figure 3, where the compression algorithm JPEG2000 is based on the
sparse representation. The Figure clearly depicts that the significant wavelet
coefficients (coefficients with relatively large absolute values) are located
on strong edges of the image.

(a) Fishing boat image.

(b) 15267 significant coefficients.

(c) Compressed image 1:17.

Figure 3: Image compression based on sparse wavelet approximation.

###  3.2 Embedded Wavelet Tokenization

The sparse wavelet representation (2) of an image provides the perfect
infrastructure for the generation of embedded coding representations (11, 12,
13). Embedded coding is similar in spirit to binary finite precision
representations of real numbers. For each digit added to the right, more
precision is added. Yet, the “encoding” can cease at any time and provide the
“best” representation of the real number achievable within the framework of
the binary digit representation. Similarly, the embedded coder can cease at
any time and provide the “best” representation of an image achievable within
its framework. Embedded coding streams can be generated from wavelet
representations by ordering the information on the wavelet representation
starting with the most significant bits of the most significant coefficients.
that is, the coefficients with the largest absolute value.

First, for simplicity of notation, using (2), denote for $I=(i_{1},i_{2})$,
$1\leq i_{1},i_{2}\leq 2$

$\alpha_{I}=\langle f,\tilde{\varphi}_{m,I}\rangle.$

We also map the coefficients $\\{\langle f,\tilde{\psi}^{e}_{j,k}\rangle\\}$
for $1\leq j\leq m$, $e=1,2,3$, based on their location $I=(i_{1},i_{2})$,
$3\leq i_{1},i_{2}\leq M$, in the coefficient matrix

$\alpha_{I}\leftarrow\langle f,\tilde{\psi}^{e}_{j,k}\rangle.$

We note in passing that one may assume that the low resolution scaling
function coefficients from (2) are known during training and are randomly
sampled from some distribution during image generation and therefore need not
be part of the tokenization.

####  3.2.1 Encoding a wavelet representation into a token sequence

We now show how to process the numeric representation of the coefficients,
from most significant to least significant and ‘encode’ it into a relatively
compact series of tokens. The representation using the series of tokens should
be invertible. That is, one can ‘decode’ the sequence back to the wavelet
representation.

To this end, assuming the image pixels are normalized to the range $[0,1]$,
one can show that for an image of dyadic dimension $[M,M]=[2^{m},2^{m}]$

$\max_{I}|\alpha_{I}|\leq 2^{m-1}.$ (3)

Assuming for simplicity that all images of a given dataset have the same
dyadic dimensions $[M,M]$, then this bound holds for all of their wavelet
representations. Our first option is to initialize a threshold $T=2^{m-2}$ and
begin scanning the wavelet coefficients of the image, in a predetermined order
(see below) for significance, with the goal of reporting only those
coefficients for which the following holds

$T\leq|\alpha_{I}|<2T.$

Our second option, is to compute separately for each image in the dataset

$\tilde{m}:=\lceil{\log_{2}\max_{I}|\alpha_{I}|}\rceil,$ (4)

and then initialize for the specific image $T=2^{\tilde{m}-1}$. In this
scenario, we store and use the parameter $\tilde{m}$ for each image in the
training set along with its sequence of tokens.

We also maintain a matrix of approximated wavelet coefficients
$\\{\tilde{\alpha}_{I}\\}$ which we initialize with zeros. Once we complete
the processing at a given bit plane, we update $T\leftarrow T/2$ and repeat
the process. At each bit-plane we report the significant coefficients that
were just uncovered in this bit-plane using a token ‘NowSignificantNeg’ if the
coefficient is negative or a token ‘NowSignificantPos’ if it is positive. At
the time of uncovering, we modify the approximation of the coefficient
$\tilde{\alpha}_{I}$ to have the absolute value $3T/2$, with the reported
sign. Next, we add a token to represent the coefficient’s next significant
bit, ‘NextAccuracy0’ if the coefficient satisfies
$|\alpha_{I}|\leq|\tilde{\alpha}_{I}|$ or the token ‘NextAccuracy1’ if
$|\alpha_{I}|>|\tilde{\alpha}_{I}|$. The approximation $\tilde{\alpha}_{I}$ is
updated accordingly by subtracting or adding $T/4$ (depending on the sign of
the coefficient and the accuracy bit type).

Let us demonstrate with an example. Assume T=16 and $\alpha_{I}=-17.45$.
Therefore, the coefficient is first uncovered in the current bit plane. When
we arrive at the index $I$, we report a ‘NowSignificantNeg’ token for this
coefficient, providing it with a temporary approximation
$\tilde{\alpha}_{I}=-24$, which lies in the middle of the segment
$[-T,-2T]=[-16,-32]$. Next, since in fact $|\alpha_{I}|\leq 24$, we report a
token ‘NextAccuracy0’ to represent the coefficient’s next significant bit,
providing an updated approximation $\tilde{\alpha}_{I}=-20$, which lies in the
middle of the segment $[-T,-3T/2]=[-16,-24]$, leading to a better
approximation of the ground truth value.

In case a coefficient has been uncovered in any of the previous bit-planes and
is already known to be significant, we only add one of the tokens
‘NextAccuracy0’ if $|\alpha_{I}|\leq|\tilde{\alpha}_{I}|$ or ‘NextAccuracy1’
if $|\alpha_{I}|>|\tilde{\alpha}_{I}|$. We then update the approximation
$\tilde{\alpha}_{I}$ by subtracting or adding $T/4$ (depending on the sign of
the coefficient and the accuracy bit type).

Assuming the bit-plane scanning order of the coefficients is fixed, one then
only needs to add the token ‘Insignificant’ to provide a valid invertible
tokenization process. One simply scans the coefficients in the fixed order and
uses their true known value to test and apply one of three possibilities:

  1. (i)

$|\alpha_{I}|\geq 2T$: The coefficient has already been reported as
significant in a previous bit-plane. Therefore one reports the token
‘NextAccuracy0’ or ‘NextAccuracy1’ depending on the test
$|\alpha_{I}|<|\tilde{\alpha}_{I}|$.

  2. (ii)

$T\leq|\alpha_{I}|<2T$: First report the token ‘NowSignificantNeg’ or
‘NowSignificantPos’ depending on the sign and then report the token
‘NextAccuracy0’ or ‘NextAccuracy1’.

  3. (iii)

$|\alpha_{I}|<T$: report ‘Insignificant’.

The process described above, although completely sufficient for invertible
tokenization, potentially creates long sequences. Specifically, it does not
take into consideration the local correlations among ‘neighboring’
insignificant wavelet coefficients. Due to the sparsity property of the
wavelet transform, during the scanning process, many of the ‘Insignificant’
coefficients form local groups. Moreover, there are correlations between local
groups of insignificant coefficients of the same subband type across
resolutions in the manner of the quad-tree structure of Figure 2. Image
compression algorithms such as the EZW (11) or SPIHT (12), are based on
statistical zero tree models that try to capture these correlations across
resolutions. As we shall later see, for image generation, we actually rely on
the powerful capabilities of the transformer models to learn correlation
patterns of the ‘wavelet language’ of the given dataset. However, we do ‘ease
the burden’ off the transformers significantly by utilizing the structure of
the groups of insignificant coefficients to reduce the size of the token
sequences, thereby creating shorter contexts.

To this end, we add two additional tokens for groups of insignificant
coefficients: ‘Group4x4’ and ‘Group2x2’ and modify the scanning process to
visit the coefficients based on groups of $4\times 4$. The first token is used
in locations where the scan is at an index $(4l_{1},4l_{2})$, for some
integers $l_{1},l_{2}$. If at the current bit plane, all the 16 coefficients
with indices $I=(i_{1},i_{2})$, $4l_{1}\leq i_{1}\leq 4(l_{1}+1),4l_{2}\leq
i_{2}\leq 4(l_{2}+1)$, are still insignificant, we issue the token ‘Group4x4’
and the tokenization process continues to the next group of $4\times 4$
coefficients. However, if any of the coefficients of the $4\times 4$ group
becomes significant in the current bit-plane, the group breaks down to 4
groups of $2\times 2$. If a group of $2\times 2$ is still composed of
insignificant coefficients at the current bit-plane, we add a token
‘Group2x2’. If a group of $2\times 2$ breaks down, then each coefficient from
the group is reported individually as being ‘Insignificant’ or one of
‘NowSignificantNeg’, ‘NowSignificantPos’. The scanning process keeps track of
which groups broke up, so that only necessary and informative tokens are
generated. We summarize the seven tokens and their roles below

  1. (i)

‘Group4x4’ – At the index $(4l_{1},4l_{2})$, the group of 16 coefficients
$\\{\alpha_{I}\\}$, $4l_{1}\leq i_{1}\leq 4l_{1}+4$, $4l_{2}\leq i_{2}\leq
4l_{2}+4$, are still insignificant, $|\alpha_{I}|<T$.

  2. (ii)

‘Group2x2’ – At the index $(2l_{1},2l_{2})$, the group of 4 coefficients
$\\{\alpha_{I}\\}$, $2l_{1}\leq i_{1}\leq 2l_{1}+2$, $2l_{2}\leq i_{1}\leq
2l_{2}+2$, are still insignificant, $|\alpha_{I}|<T$.

  3. (iii)

’NowSignificantNeg’, ’NowSignificantPos’ – At the current location $I$, the
coefficient satisfies $T\leq|\alpha_{I}|<2T$. If the coefficient was part of a
group of insignificant coefficients at the previous bit-plane, the group is
now automatically dissolved.

  4. (iv)

‘Insignificant’ – At the current location $I$, the coefficient is still
insignificant and satisfies $|\alpha_{I}|<T$. If the coefficient was part of a
group of insignificant coefficients at the previous bit-plane, the group is
now automatically dissolved.

  5. v)

‘NextAccuracy0’, ‘NextAccuracy1’ – At the current location $I$, the
coefficient has already been reported to be significant since it satisfies
$|\alpha_{I}|\geq T$. Here, we improve the accuracy of its approximation using
one of these tokens, depending on the test
$|\alpha_{I}|<|\tilde{\alpha}_{I}|$.

The bit-plane scan is carried out in two nested loops; the outer loop proceeds
from low resolution to high resolution, each time traversing the three types
of wavelet subbands. The inner loop traverses the $4\times 4$ blocks. Figure 4
illustrates the outer and inner scanning patterns.

(a) Subband scanning order.

(b) Scanning order of $4\times 4$ blocks.

Figure 4: A sketch illustrating the outer and inner scanning orders.

Figure 5 exemplifies the tokenization algorithm of an image from the MNIST
dataset (30). The image was padded with zeros to be of dimensions $M\times
M=32\times 32$, with $m=5$. The bottom row of the figure shows the tokens and
their locations on the wavelet image for the first three bit planes. To make
the process clearer, we explicitly write the resulted sequence of tokens for
the first bit plane shown in Figure 5(d).

$\displaystyle\\{$ ‘Insignificant’, ‘Insignificant’, ‘NowSignificantPos’,
‘Insignificant’, ‘Insignificant’, ‘Insignificant’, ‘NowSignificantNeg’,
‘Insignificant’, ‘Group2x2’, ‘Group2x2’, ‘Insignificant’, ‘Insignificant’,
‘Insignificant’, ‘NowSignificantNeg’, ‘Group2x2’,‘Group2x2’, ‘Group2x2’,
$\displaystyle\text{`Group4x4'},\dots,\text{`Group4x4'}\\}$

The token sequences of the second and third bit-planes follow the same
scanning pattern. Eventually, the three sequences are concatenated in the
natural order to form the final sequence which describes the three bit planes
wavelet image appearing in Figure 5(c).

There is a very important hyper-parameter which is the choice of the smallest
threshold at the final bit-plane. Through this hyper-parameter, the wavelet
representation provides us with a very robust and stable trade-off of fine
detail generation and length of token sequences. Choosing a final threshold
provides very consistent control over visual quality relating to: “Visually
Lossless”, “High”, “Medium”, “Low”, etc. This is en par with the quality
settings in digital cameras, which in turn lead to a selection of the
corresponding quantization tables of the JPEG algorithm generating the
compressed images.

(a) $32\times 32$ padded MNIST image.

(b) Wavelet Transform.

(c) Three bit planes of coefficients.

(d) First bit plane.

(e) Second bit plane.

(f) Third bit plane.

Figure 5: Depiction of the tokenization process. On the top left and middle, a
$32\times 32$ padded MNIST image and its wavelet transform. On the top right,
the wavelet approximation generated by the first three bit-planes. The bottom
row illustrates the tokens and their locations on the $32\times 32$ grid,
where, ‘NowSignificantNeg’ and ‘NowSignificantPos’ tokens are annotated with
orange “$-$” and blue “$+$” and signs respectively. The tokens ‘NextAcurracy0’
and ‘NextAccuracy1’ are marked with green down and red up triangles. The
purple dots represent ‘Insignificant’ coefficients and the brown and pink
squares represent the ‘Group2x2’ and ‘Group4x4’ zero block tokens.

####  3.2.2 Decoding the token sequence into an approximate wavelet
representation

The tokenization process described in the previous subsection can be easily
inverted back to an approximate wavelet representation. Moreover, any initial
sub-sequence can be inverted to provide a possibly coarser approximation. We
initialize a matrix of size $M\times M$ of the approximated wavelet
coefficients $\\{\tilde{\alpha}_{I}\\}$ with zeros and begin the scanning
process with the first bit-plane. Based on (3) or (4), we know how to
initialize the first bit-plane with the initial threshold $T=2^{m-2}$ or
$T=2^{\tilde{m}-1}$. We then process the token sequence and update the
approximated coefficients using the corresponding ‘significant’ and ‘bit
accuracy’ tokens. If for any given reason, the sequence of tokens terminates,
we have the best possible approximated coefficients $\\{\tilde{\alpha}_{I}\\}$
from which we can obtain an approximated image by applying the inverse DWT.
Our decoding process relies on the assumption that the token sequence is
valid. For example, a ‘Group4x4’ token cannot appear while the decoder scan
position is at a location of indices not divisible by 4. It is obvious how to
achieve this in the context of image coding. However, during an image
generation process, this needs to be enforced using the conditional next-token
inference described in Subsection 4.6.

##  4 The Generative Wavelet Transformer

Assume that for a certain dataset of images, we have established the
translation of the visual information of each image to a sequence of tokens
encapsulating the visual information from coarse to fine details as explained
in Subsection 3.2. We assume that within the sequences, distinct patterns and
relations exist between the tokens. For example, the wavelet coefficients
$\\{\langle f,\tilde{\psi}^{e}_{j,k}\rangle\\}$ of wavelets
$\\{\psi^{e}_{j,k}\\}$ whose support intersects with a certain portion of an
edge of the image, will be significant and aligned across scales in a tree-
like structure as per Figures 2 and 3(b). At the same time, coefficients of
wavelets whose support intersects with a smooth area of the image will be
insignificant and they appear in local groups. This leads to the intuition
that the powerful transformers created over the last few years (7) are able to
learn the patterns of the ‘wavelet language’ and to generate them from some
random seeds during inference.

Interestingly, as we experienced through empirical experimentation, one can
actually take off-the-shelf pre-trained transformers such as T5 (32) or
DistilBert (33) and fine-tune them to our purpose of learning the wavelet
language with no additional modifications and with some reasonable generative
results. That is, even though the pre-trained models were trained on a set of
more than 20,000 tokens and on text datasets, one is able to fine-tune them on
wavelet-based token sequences of images, containing only 7 tokens. However,
some classical components of transformers are not aligned with the ‘wavelet
language’ and should be replaced. At the same time, there are some very useful
components and techniques that come with language models that can be also
leveraged successfully in the context of image generation.

In this section, we describe how we modified the architecture of the
DistilGPT2 model (33) (we used the code from HuggingFace (34) as a starting
point) to optimize it to align with the wavelet-based image generation method.
This obviously requires training the modified model from scratch.

###  4.1 Token vector representation

Typically, in the standard scenarios of spoken languages, transformers apply a
‘pre-processing’ learnable transform to tokens to convert them to vector
representations. The idea is that similar words should be converted to vectors
with some proximity, which intuitively serve as better input for the
transformer’s neural network. However, as explained in Subsection 3.2.1, our
wavelet dictionary includes only seven tokens that have very distinctive and
different roles. Therefore, the simple transformation of the tokens to the
one-shot encoding of the standard basis of dimension 7 is probably a better,
if not optimal choice. Thus, the initial vector representation of a token is:
$\text{`Group4x4'}\rightarrow(1,0,0,0,0,0,0)$,
$\text{`Group2x2'}\rightarrow(0,1,0,0,0,0,0)$, etc. This means that in the
transformers we modified, we removed the ‘token $\rightarrow$ vector’
learnable transformation.

###  4.2 Initial bit-plane threshold

Recall that we have two options: to use a uniform initial bit-plane threshold
for all images in the dataset derived from (3), or to use an adaptive initial
threshold for each image of the training set using (4). In the latter case, we
need to inform the transformer, per image, which initial threshold the token
sequence is associated with. We do this as follows: assume a given dataset has
$l$ possible values for $\tilde{m}$ in (4) (e.g., $l=4$ for the MNIST dataset,
see Figure 6). Then, we concatenate a one-shot encoding of dimension $l$ of
the initial threshold parameter of the given image to each vector
representation of each token.

For image generation, one may sample randomly from the distribution of $l$
possible initial thresholds. In the case that the image generation is
conditioned on a certain class (see Subsection 4.4), one can sample from the
conditional distribution of the possible thresholds of the specific class.

Figure 6: Distribution of $\log_{2}$ of the initial thresholds for the 70,000
MNIST images with the Haar wavelet transform.

###  4.3 Positional encoding

In classic transformer architectures (7), one adds the positional encoding
$v_{p}(t)$ of the position $t$ to the token’s vector representation
$v_{e}(x(t))$. Learnt positional embedding applied a learned transform
$t\rightarrow v_{p}(t)$. Some transformers use hard-coded mapping of the
position. Assuming the vector embedding dimension is $d_{e}$ and the maximum
length of a sequence is $l_{\max}$, then

$v_{p}(t)(2i-1)=\sin(t/l_{\max}^{2i/d_{e}}),\quad
v_{p}(t)(2i)=\cos(t/l_{\max}^{2i/d_{e}}),\quad 1\leq i\leq d_{e}/2.$

In our scenario of the wavelet language, the position of a token in a sequence
is $(bp,I)$, where $bp$ is the enumeration of the bit-plane and
$I=(i_{1},i_{2})$ is the index of the current coefficient $\alpha_{I}$ in the
scan order. Therefore, we concatenate to the vector representation of a token
from Subsection 4.1, a vector component of dimension 3 with the location of
the token $(bp,i_{1},i_{2})$.

###  4.4 Generative guidance

It is obviously critical for any image generation method to allow guidance of
the generative process by placing a condition on the class type of the
generated image or a text prompt that describes it. Some image generation
models apply a joint embedding space for text and images for this purpose. One
such method is to used a pretrained model such as CLIP (35) that maps text and
images to a joint embedding space. The CLIP contains an image encoder $f$ and
a caption encoder $g$, that during training over pairs of images with captions
$\\{(x,c)\\}$, optimizes a contrastive cross-entropy loss that encourages high
dot-products $\langle f(x),g(c)\rangle$ in the joint embedding space. Thus,
any image generation method, can use the vector embedding of the given text
prompt $c$ to guide the generative process by conditioning the image embedding
$f(x)$ to be highly correlated with the embedding of the textual prompt.

In our case, since we converted the problem of image generation to a ‘wavelet-
language’ generation, we can apply ‘text’-type prompting methods. Having
access to a joint embedding text-image space allows us to train using the
vector representation of the image training set. Then, at image generation, we
use the vector representation of the given text prompt to guide the generative
process. There are very simple ways of using these vector representations. We
choose to concatenate them to the vector representation of each token and its
position (as explained above). For example, as shown in Section 5, for the
image datasets MNIST or FashionMNIST with 10 classes, it is easy to
concatenate a vector of length 10 representing the class of the image. In the
case where we wish to guide the generative process using a textual prompt, we
may concatenate the CLIP vector embedding (35) of the textual prompt to each
token vector representation. As we discuss in Subsection 6.2, we hope this
approach to guiding the generative process can be generalized to composition
of blobs (36), where a given guiding vector of a blob is used only at
positions of the scan where the support of the corresponding wavelet
intersects the blob.

###  4.5 Initialization of the generative process

Since the guidance of the generative process (Subsection 4.4) is applied
through the concatenation of vector representations to each token vector
representation, in some cases, the initialization becomes a minor issue. For
example, when training on MNIST and generating digits, one can get away with a
simple random choice from the subset: ‘Insignificant’, ‘NowSignificantNeg’ or
‘NowSignifiantPos’ for the first token and from there the transformer will
generate a valid token sequence which is converted to an adequate image of a
digit from the pre-selected class.

A more robust method is as follows. Suppose we wish to generate a handwritten
digit from a certain digit class. Let $\\{f_{s}\\}_{s\in S}$ be the subset of
MNIST images from that specific digit class and let

$\\{\langle f_{s},\tilde{\varphi}_{m,k}\rangle\\},\quad s\in S,\quad
k=(k_{1},k_{2}),\quad 1\leq k_{1},k_{2}\leq 2,$ (5)

be the subset of low-resolution coefficients of these images defined by (2).
Let $N(v,\Sigma)$ be the fourth-dimensional normal distribution, approximated
by the subset (5). We then sample from $N(v,\Sigma)$, a random group of four
low-resolution coefficients. Now, the token representation of these
coefficients can serve as a basis for a robust initialization of the
generative process of the required digit. In the case where the guidance is
provided by a vector representation of some text-prompt, one can create the
normal distribution using a subset of $K$-nearest neighbors in the image
vector representation space.

Once some random seeding allows us to initialize the token sequence, we may
introduce as much diversity as required using the methods of Subsection 4.7 so
that even using the same seed may generate various images corresponding to the
given guidance.

###  4.6 Conditional next token inference

In Greedy generative mode, the next selected token $x(t)$, $1\leq x(t)\leq 7$,
at location $t$, is the token for which the transformer assigns the highest
probability from $(p_{1}(t),...,p_{7}(t))$. As described in Subsection 4.7
below, there are various alternative methods to control the output of the
transformer. However, since each generative token inference step is a
statistical event, it may occur that the next predicted token is not valid at
the current position of the wavelet bit-plane scan. To overcome this, we
ensure any selected token satisfies the conditions below relating to the
context and the current position in the scan. For example, when using the
Greedy method, one simply picks from the subset of tokens satisfying the
conditions below, the one with the highest probability.

  1. (i)

‘Group4x4’ - The scan is at an index $(4l_{1},4l_{2})$ and the group has not
yet dissolved.

  2. (ii)

‘Group2x2’ - The scan is at an index $(2l_{1},2l_{2})$ and the group has not
yet dissolved.

  3. (iii)

’NowSignificantNeg’, ’NowSignificantPos’ - At the current location $I$, the
coefficient $\alpha_{I}$ is still insignificant, possibly as part of a group
of insignificant coefficients.

  4. (iv)

‘Insignificant’ - At the current location $I$, the coefficient $\alpha_{I}$ is
still insignificant, possibly as part of a group of insignificant
coefficients.

  5. v)

‘NextAccuracy0’, ‘NextAccuracy1’ - At the current location $I$, the
coefficient $\alpha_{I}$ has already been reported to be significant.

###  4.7 Controlling the degree of generative diversity during inference

Since we are applying a language transformer model we may use various simple
stochastic mechanisms to control the generative process during inference and
allow a diversity of possible images to be generated from a single prompt.
Some of the available stochastic methods are: Beam search with multinomial
sampling, Top-$k$ and Top-$p$. In our experiments, we tested the latter two:

  1. (i)

Top-$k$ sampling - The Top-$k$ inference method (37) filters the $k$ most
likely next words first and then samples from the probability mass that is
redistributed among only those $k$ next words. GPT2 adopted this sampling
scheme, which was one of the reasons for its success in story generation. In
Figure 8 below, we see a diversity of sandals generated by guiding the model
with the vector representation of the corresponding FashionMNIST ‘sandal’
class and using the Top-$2$ method. We see that using $k=2$ is sufficient to
move the generative process from a deterministic process to a sufficiently
diverse stochastic process, yet with output that fits the class description.

  2. (ii)

Top-$p$ sampling- In Top-$p$ sampling or nucleus sampling, the selection pool
for the next token is determined by the cumulative probability of the most
probable tokens. Setting a threshold $p$, the model includes just enough of
the most probable tokens so that their combined probability reaches or exceeds
this threshold. Again, the distribution mass is redistributed among these
tokens and then the next token is sampled using this distribution. In Figure 7
we see different examples of the digits ‘3’ and ‘8’ generated using the
Top-$0.6$ method.

##  5 Experimental results

We conducted experiments on the MNIST and FashionMNIST datasets. Here are some
details:

  * •

The images in both datasets were padded with zeros to $M\times M=32\times 32$
and normalized to have values within $[0,1]$.

  * •

We used the Haar wavelet basis for the MNIST images and the bior4.4 wavelet
basis for the FashionMNIST.

  * •

The images were tokenized with a final threshold of $T=2^{-3}$ for MNIST and
$T=2^{-4}$ for FashionMNIST.

  * •

The maximal token sequence lengths were 1742 for MNIST and 3098 for
FashionMNIST.

  * •

We trained two separate distillgpt2 models from scratch on the two datasets.
As for the training configurations, both training sessions had batch size 4,
learning rate 0.0004, and weight decay 0.01.

  * •

Models were trained on an NVIDIA A100 GPU with 80GB; MNIST occupied around
22GB while FashionMNIST occupied 61GB. Both models were trained for a few
days.

Results with different controlling methods appear below in Figures 7 and 8.

Figure 7: Digits generated with Top-$p=0.6$ along with a depiction of the
generated wavelet coefficients.

Figure 8: Sandals generated with Top-$k=2$ along with a depiction of the
generated wavelet coefficients.

More generated images for different classes of MNIST and FashionMNIST appear
in the following figures.

Figure 9: More MNIST results.

Figure 10: More FashionMNIST results.

##  6 Discussion and future work

In this paper, we introduced a novel method for image generation that is based
on elements of wavelet image coding and NLP transformers. Unfortunately, our
research group does not have access to sufficient computational resources at
the moment, so this work serves as a first modest proof of concept. Indeed,
the wavelet representation is a powerful tool in image processing that can
serve as a basis for many image generation functionalities. Here, we list some
directions that we will consider for future work.

###  6.1 Generation of color images at high resolution and with fine details

In our experiments, we only generated small grayscale images. We provide here
some details on how the method can be generalized:

  1. (i)

Color images - For color images (or even spectral images), we may adopt a
well-known paradigm from image compression. For improved performance, one may
transform input images in the $RGB$ color space to the $YCbCr$ color space.
The $Y$ component is the luminance component, essentially the image’s
grayscale part. The other two components, $Cb$ and $Cr$, capture the color
information of the image. Typically, the luminance component carries most of
the visual information, and thus also, its encoding is usually the significant
part of an encoded image. In image coding, one usually encodes separately each
of the three channels. Our method can then be generalized to color images by
applying the DWT and the tokenization process separately to each color
channel.

  2. (ii)

Image resolution - Support of any image resolution simply translates to more
wavelet coefficients and longer token sequences. Obviously, this requires
larger transformer models that can support longer contexts and more
computational resources. The method should also generalize well, just as
wavelet image coding is being applied to Gigabyte images.

  3. (iii)

Generating fine details - Using our wavelet model, finer details are captured
at higher bit-planes. The choice of the final threshold of the final bit-plane
provides excellent and very consistent control over the amount of detail one
wishes to generate (see Subsection 3.2.1). This quantization technique is at
the heart of the JPEG algorithms and translates to very specific modes in
digital cameras that can be set to: ”Visually Lossless”, “High”, “Medium”,
etc. This exact form of control also applies to wavelets but, unfortunately,
is not the default mode of operation in JPEG2000. Obviously, to generate finer
details, one needs to train the transformer on longer token sequences, again
requiring more computational resources.

###  6.2 Support for generation of compositions of blobs

In many cases, one wishes to apply fine-grained control of compositional text-
to-image generation, where certain locations in the image, marked perhaps with
bounding boxes or ellipses, receive different textual descriptions (36). One
possible method to accomplish this using the wavelet generative approach is to
apply the transformer in evaluation mode and apply the vector representation
of the blob’s textual prompt as described in Subsection 4.4 whenever the bit-
plane scan is at indices of wavelet coefficients whose support intersects the
blob.

###  6.3 Multi-modal generation

The ability to represent an image’s visual information as a sequence of tokens
presents an attractive possibility of merging the wavelet-based tokens with
other language tokens to create a unified multi-modal transformer.

## Funding

N. Sharon is partially supported by the NSF-BSF award 2019752. W. Mattar is
partially supported by The Nehemia Levtzion Scholarship for Outstanding
Doctoral Students from the Periphery (2023). N. Sharon and W. Mattar are
partially supported by the DFG award 514588180.

## References

  * Dhariwal and Nichol (2021) Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in neural information processing systems_ , 34:8780–8794, 2021.
  * Saharia et al. (2022) Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily L Denton, Kamyar Ghasemipour, Raphael Gontijo Lopes, Burcu Karagol Ayan, Tim Salimans, et al.  Photorealistic text-to-image diffusion models with deep language understanding.  _Advances in neural information processing systems_ , 35:36479–36494, 2022.
  * Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 6840–6851. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper_files/paper/2020/file/4c5bcfec8584af0d967f1ab10179ca4b-Paper.pdf.
  * Rombach et al. (2022) Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 10684–10695, 2022.
  * Ramesh et al. (2022) Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.  Hierarchical text-conditional image generation with clip latents.  _arXiv preprint arXiv:2204.06125_ , 1(2):3, 2022.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Phuong and Hutter (2022) Mary Phuong and Marcus Hutter.  Formal algorithms for transformers, 2022.  URL https://arxiv.org/abs/2207.09238.
  * Tay et al. (2022) Yi Tay, Mostafa Dehghani, Dara Bahri, and Donald Metzler.  Efficient transformers: A survey.  _ACM Computing Surveys_ , 55(6):1–28, 2022.
  * Esser et al. (2021) Patrick Esser, Robin Rombach, and Bjorn Ommer.  Taming transformers for high-resolution image synthesis.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 12873–12883, 2021.
  * Ramesh et al. (2021) Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen, and Ilya Sutskever.  Zero-shot text-to-image generation.  In _International conference on machine learning_ , pages 8821–8831. Pmlr, 2021.
  * Shapiro (1993) Jerome Shapiro.  Embedded image coding using zerotrees of wavelet coefficients.  _IEEE Transactions in signal processing_ , 41(12):3445–3462, 1993.
  * Said and Pearlman (1996) Amir Said and William Pearlman.  A new, fast, and efficient image codec based on set partitioning in hierarchical trees.  _IEEE Transactions on Circuits and Systems for Video Technology_ , 6(3):243–250, 1996.
  * Taubman and Marcellin (2002) David Taubman and Michael Marcellin.  _JPEG2000: Image Compression Fundamentals, Standards and Practice, 2nd edition_.  Springer, 2002.
  * Daubechies (1992) Ingrid Daubechies.  _Ten Lectures on Wavelets_.  SIAM, 1992.
  * DeVore (1998) Ronald DeVore.  Nonlinear approximation.  _Acta Numerica_ , 7:51–150, 1998.
  * Mallat (2009) Stephan Mallat.  _A Wavelet tour of signal processing, the sparse way_.  Academic Press, 2009.
  * Wallace (1992) G.K. Wallace.  The jpeg still picture compression standard.  _IEEE Transactions on Consumer Electronics_ , 38(1):xviii–xxxiv, 1992.  doi: 10.1109/30.125072.
  * Mihcak et al. (1999) M Kivanc Mihcak, Igor Kozintsev, and Kannan Ramchandran.  Spatially adaptive statistical modeling of wavelet image coefficients and its application to denoising.  In _1999 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings. ICASSP99 (Cat. No. 99CH36258)_ , volume 6, pages 3253–3256. IEEE, 1999.
  * Kivanc Mihcak et al. (1999) M. Kivanc Mihcak, I. Kozintsev, K. Ramchandran, and P. Moulin.  Low-complexity image denoising based on statistical modeling of wavelet coefficients.  _IEEE Signal Processing Letters_ , 6(12):300–303, 1999.  doi: 10.1109/97.803428.
  * Buccigrossi and Simoncelli (1999) Robert W Buccigrossi and Eero P Simoncelli.  Image compression via joint statistical characterization in the wavelet domain.  _IEEE transactions on Image processing_ , 8(12):1688–1701, 1999.
  * Nichol et al. (2022) Alexander Quinn Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob Mcgrew, Ilya Sutskever, and Mark Chen.  GLIDE: Towards photorealistic image generation and editing with text-guided diffusion models.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 16784–16804. PMLR, July 2022.  URL https://proceedings.mlr.press/v162/nichol22a.html.
  * Wang et al. (2023) Xinlong Wang, Wen Wang, Yue Cao, Chunhua Shen, and Tiejun Huang.  Images speak in images: A generalist painter for in-context visual learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 6830–6839, 2023.
  * Lee et al. (2022) Doyup Lee, Chiheon Kim, Saehoon Kim, Minsu Cho, and Wook-Shin Han.  Autoregressive image generation using residual quantization.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 11523–11532, 2022.
  * Yu et al. (2024) Qihang Yu, Mark Weber, Xueqing Deng, Xiaohui Shen, Daniel Cremers, and Laing-Chieh Chen.  An image is worth 32 tokens for reconstruction and generation.  _arXiv preprint arXiv:2406.07550_ , 2024.
  * Tian et al. (2024) Keyu Tian, Yi Jiang, Zehuan Yuan, Bingyue Peng, and Liwei Wang.  Visual autoregressive modeling: Scalable image generation via next-scale prediction.  _arXiv preprint arXiv:2404.02905_ , 2024.
  * Zhu et al. (2023) Qing Zhu, Xiumei Li, Junmei Sun, and Huang Bai.  Wdig: a wavelet domain image generation framework based on frequency domain optimization.  _EURASIP Journal on Advances in Signal Processing_ , 2023(1):66, 2023.
  * Yu et al. (2021) Yingchen Yu, Fangneng Zhan, Shijian Lu, Jianxiong Pan, Feiying Ma, Xuansong Xie, and Chunyan Miao.  Wavefill: A wavelet-based generation network for image inpainting.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pages 14114–14123, 2021.
  * Phung et al. (2023) Hao Phung, Quan Dao, and Anh Tran.  Wavelet diffusion models are fast and scalable image generators.  In _2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 10199–10208, 2023.  doi: 10.1109/CVPR52729.2023.00983.
  * Guth et al. (2022) Florentin Guth, Simon Coste, Valentin De Bortoli, and Stephane Mallat.  Wavelet score-based generative modeling.  _Advances in Neural Information Processing Systems_ , 35:478–491, 2022.
  * Deng (2012) Li Deng.  The mnist database of handwritten digit images for machine learning research [best of the web].  _IEEE signal processing magazine_ , 29(6):141–142, 2012.
  * Cerná and Finěk (2011) Dana Cerná and Václav Finěk.  Discrete CDF 9 / 7 wavelet transform for finite-length signals, 2011.  URL https://api.semanticscholar.org/CorpusID:208013335.
  * Chung et al. (2024) Hyung Won Chung, Le Hou, Shayne Longpre, Barret Zoph, Yi Tay, William Fedus, Yunxuan Li, Xuezhi Wang, Mostafa Dehghani, Siddhartha Brahma, Albert Webson, Shixiang Shane Gu, Zhuyun Dai, Mirac Suzgun, Xinyun Chen, Aakanksha Chowdhery, Alex Castro-Ros, Marie Pellat, Kevin Robinson, Dasha Valter, Sharan Narang, Gaurav Mishra, Adams Yu, Vincent Zhao, Yanping Huang, Andrew Dai, Hongkun Yu, Slav Petrov, Ed H. Chi, Jeff Dean, Jacob Devlin, Adam Roberts, Denny Zhou, Quoc V. Le, and Jason Wei.  Scaling instruction-finetuned language models.  _Journal of Machine Learning Research_ , 25(70):1–53, 2024.  URL http://jmlr.org/papers/v25/23-0870.html.
  * Sanh (2019) V Sanh.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  In _Proceedings of Thirty-third Conference on Neural Information Processing Systems (NIPS2019)_ , 2019.
  * (34) Huggingface Repository.  HuggingFace DistilGPT2.  https://huggingface.co/distilbert/distilgpt2, 2019.  Accessed: 2024-06-27.
  * Radford et al. (2021) Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al.  Learning transferable visual models from natural language supervision.  In _International conference on machine learning_ , pages 8748–8763. PMLR, 2021.
  * Nie et al. (2024) Weili Nie, Sifei Liu, Morteza Mardani, Chao Liu, Benjamin Eckart, and Arash Vahdat.  Compositional text-to-image generation with dense blob representations, 2024.
  * Fan et al. (2018) Angela Fan, Mike Lewis, and Yann Dauphin.  Hierarchical neural story generation.  In Iryna Gurevych and Yusuke Miyao, editors, _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 889–898, Melbourne, Australia, July 2018. Association for Computational Linguistics.  doi: 10.18653/v1/P18-1082.  URL https://aclanthology.org/P18-1082.


11institutetext: Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300
RA Leiden, The Netherlands 22institutetext: European Southern Observatory,
Karl-Schwarzschild-Str. 2, 85748, Garching, Germany 33institutetext: SRON
Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The
Netherlands44institutetext: INAF – Istituto di Radioastronomia, via P. Gobetti
101, 40129 Bologna, Italy55institutetext: Department of Astrophysical
Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA

# Complex diffuse emission in the $z=0.52$ cluster PLCK G004.5-19.5

J. G. Albert 11 C. Sifón 1155 A. Stroe ESO Fellow2211 F. Mernier 3311 H. T.
Intema 11 H. J. A. Röttgering 11 G. Brunetti 44

We present radio observations of the galaxy cluster PLCK G004.5-19.5
($z=0.52$) using the Giant Metrewave Radio Telescope at 150 MHz, 325 MHz, and
610 MHz. We find an unusual arrangement of diffuse radio emission in the
center and periphery of the cluster, as well as several radio galaxies with
head-tail emission. A patch of peripheral emission resembles a radio relic,
and central emission resembles a radio halo. Reanalysis of archival _XMM-
Newton_ X-ray data shows that PLCK G004.5-19.5 is disturbed, which has a known
correlation with the existence of radio relics and halos. Given that the
number of known radio halos and radio relics at $z>0.5$ is very limited, PLCK
G004.5-19.5 is an important addition to understanding merger-related particle
acceleration at higher redshifts.

###### Key Words.:

galaxies: clusters: individual (PLCK G004.5-19.5) – galaxies: clusters:
intracluster medium – large-scale structure of universe – radiation
mechanisms: non-thermal – X-rays: galaxies: clusters

##  1 Introduction

Galaxy clusters are the most massive gravitationally bound structures in the
Universe with masses and volumes of the order of
$10^{14-15}\leavevmode\nobreak\ \mathrm{M}_{\sun}$ and
$100\leavevmode\nobreak\ \text{Mpc}^{3}$. Hierarchical mergers between
clusters at the intersections of the cosmic web, with relative velocities near
$1000\leavevmode\nobreak\ \mathrm{km}\,\mathrm{s}^{-1}$, can release
gravitational energy of the order of $10^{63}\leavevmode\nobreak\
\mathrm{erg}$. A fraction of this energy is dissipated into shocks and
turbulence, which in turn accelerate cosmic ray electrons (CRe) and hadrons in
the intracluster medium (ICM). The shock acceleration might occur via
diffusive shock acceleration (DSA; for a comprehensive review see Drury 1983;
Malkov & O’C Drury 2001; Kang & Ryu 2013). Turbulence can reaccelerate
relativistic particles via second-order Fermi mechanisms (for an updated
review see Brunetti & Jones 2014). The relativistic CRe give rise to
synchrotron radiation in the presence of magnetic fields, which are typically
of the order of $\mu\mathrm{G}$ (Roettiger et al. 1999).

Emission tracing the shock fronts is referred to as radio relics, which are
typically diffuse, elongated, and polarised sources located near the cluster
periphery. Central unpolarized, diffuse radio sources are known as radio halos
and likely result from turbulent acceleration. Clusters hosting radio relics
and halos have a disturbed morphology (Cassano et al. 2010b), but not all
merging clusters exhibit halos, suggesting that halos originate from a
hierarchy of complex mechanisms (Donnert et al. 2013; Brunetti 2016).

The current catalogues of radio relics and halos are far from complete, mainly
because of the limited sensitivities of existing radio surveys. Nuza et al.
(2012) performed an analysis of MareNostrum, a high-resolution cosmological
simulation, and modelled the abundance of radio relics across redshift. By
assuming a DSA efficiency, they estimate $\gtrsim 100$(800) to be found at
$0.5<z<1$ by the upcoming Tier-1 survey of LOFAR (van Haarlem et al. 2013) at
60 MHz (120 MHz). Currently, there are only four relics observed beyond
$z=0.5$ (Bonafede et al. 2012, MACS J1149.5+2223 at $z=0.54$; Bonafede et al.
2009; van Weeren et al. 2009, MACS J0717.5+3745 at $z=0.55$; Lindner et al.
2014; Botteon et al. 2016, El Gordo at $z=0.87$ and Riseley et al. 2017, MACS
J0025.4-1222 at $z=0.586$).

The number of observed halos beyond $z=0.5$ has been limited by the
sensitivity of low-frequency radio observations. It is estimated that there
are about 500–1000 radio halos for $z<0.6$ that should be found at 150 MHz by
LOFAR Tier-1 Cassano et al. (2006). This is because inverse Compton (IC)
losses (scaling with $(1+z)^{4}$) become more efficient than DSA at higher
redshifts (Cassano et al. 2006). Thus IC losses are able to suppress high-
energy CRe, giving rise to an increasing fraction of radio halos with very
steep spectra, which are only bright at low frequencies (e.g. Brunetti et al.
2008).

One as yet unanswered question is how efficient DSA and turbulence are at
electron acceleration. In particular, it is unclear how the two mechanisms
evolve over cosmic time, and what their impact is on the underlying magnetic
field. Furthermore, we eventually wish to understand how shocks and turbulence
affect host galaxies following a merger.

We present galaxy cluster PLCK G004.5-19.5, which lies at a redshift
($z=0.516$) with few known radio relics and halos; this makes it a tantalizing
target to search for diffuse emission. It was discovered by the Planck
satellite through the Sunyaev-Zel’dovich (SZ) effect and confirmed with X-ray
observations (Planck Collaboration et al. 2011, 2014). This cluster is hot,
$10.2\pm 0.5$ keV, and very massive, $M^{SZ}_{500}=(10.4\pm 0.7)\times
10^{14}\leavevmode\nobreak\ \mathrm{M}_{\sun}$, and it hosts strong lensing
arcs (Sifón et al. 2014). Initial follow-up studies in low-resolution low-
frequency archival data found that PLCK G004.5-19.5 hosted strong radio
emission (Sifón et al. 2014). However, since the radio sources were
unresolved, the nature of the radio emission could not be solved.

In this paper we offer new high-resolution 150 MHz, 325 MHz, and 610 MHz
observations of the cluster PLCK G004.5-19.5 using the Giant Metrewave Radio
Telescope (GMRT), and a reanalysis of XMM-Newton data for a morphology study
of the X-ray emission. In Section 2 we explain our reduction method and how we
corrected XMM-Newton X-ray data of PLCK G004.5-19.5 for background and
vignetting. In Section 3 we analyse the X-ray and radio emission.

In this work, we assume a flat $\Lambda$CDM cosmology with $h=0.70$, and
$\Omega_{\rm m}=0.30$, which gives an angular scale at $z=0.516$ of
$6.2\leavevmode\nobreak\ $kpc/″.

##  2 Observations and data reduction

###  2.1 GMRT observations of PLCK G004.5-19.5

We carried out a set of observations on PLCK G004.5-19.5 using the GMRT at 150
MHz, 325 MHz (PI: A. Stroe, Project: 27_051), and 610 MHz (PI: C. Sifón,
Project: 25_036) with total times on source of 445 min., 390 min., and 798
min. We performed 15 min. observations on flux calibrators 3C 48 and 3C 286 at
the start and end of the observation session, and 10 min. interleaved
observations of the nearby and bright phase calibrator J1924-292 that has a
flat spectrum amplitude (Healey et al. 2007).

We removed the initial radio frequency interference (RFI) of the obvious
corrupt time and frequency data by flagging by hand. We then used the source
peeling and atmospheric modelling pipeline (SPAM; for a full description see
Intema et al. 2009; Intema 2014). SPAM performs antenna delay, bandpass,
phase, and amplitude calibration, and then multiple rounds of self-calibration
and further flagging of bad data. It then iteratively subtracts all sources
except for bright calibrators, and calibrates direction dependently. It then
spatially fits a model to the direction-dependent calibration products and
derives interpolated ionospheric corrections.

For our choice of initial sky model for self-calibration we used a bootstrap
method. We first used the 150 MHz model from the TIFR GMRT Sky Survey
Alternative Data Release (TGSSADR; Intema et al. (2016)) to start self-
calibration for our 150 MHz observations, and then used our deeper 150 MHz
model as the initial model for our 325 MHz self-calibration, likewise using
the 325 MHz image for the 610 MHz initial sky model. The TGSSADR provides an
excellent starting model for self-calibration in this lower frequency range
and is available to the public111http://tgssadr.strw.leidenuniv.nl/.

We performed multi-scale (multiscale = (0,4,8,32) pixels) multi-frequency
deconvolution with CASA on the SPAM-calibrated products with a slightly
uniform Briggs weighting (robust$=-0.3$, Briggs 1995) and four pixels per
beam. We also corrected for a non-coplanar array using $w$-projection
(Cornwell et al. 2008). The multi-scale clean avoids the issues related to
representing resolved diffuse emission with a sum of point-like sources, and
this choice of weighting balances beam uniformity in side lobes with beam
resolution. The highest minimum baseline of the three frequencies is
$120\leavevmode\nobreak\ \lambda$, corresponding to a highest sensitivity
scale of $\sim 20\leavevmode\nobreak\ \arcmin$.

The background root-mean-squared noises for the individual radio maps (Figs. 1
and 2), $\sigma_{\rm rms}$, were calculated as the standard deviation of the
residual after Gaussian source fitting and subtraction using Python Blob
Detection and Source Measurement (PyBDSM) v1.8. Our reduction yielded
$\sigma_{\rm rms}=1.40$ mJy beam-1 at 150 MHz, $\sigma_{\rm rms}=120$ $\mu$Jy
beam-1 at 325 MHz, and $\sigma_{\rm rms}=90$ $\mu$Jy beam-1 at 610 MHz. The
resolutions of the maps are $43.3\leavevmode\nobreak\ \arcsec\times
18.9\leavevmode\nobreak\ \arcsec$, $17.5\leavevmode\nobreak\ \arcsec\times
9.5\leavevmode\nobreak\ \arcsec$, and $7.2\leavevmode\nobreak\ \arcsec\times
4.9\leavevmode\nobreak\ \arcsec$, respectively. We note that the side lobes
from a nearby bright source ($S_{1.4\leavevmode\nobreak\ {\rm GHz}}=230$ mJy
from the NRAO VLA Sky Survey) cause a spatially varying background noise in
our target field.

Figure 1: 150 MHz contours (red) and 610 MHz contours (blue) overlaid on the
corrected _XMM-Newton_ EPIC image, with radio beams in the bottom left corner.
Contours are at $(5,10,20,40,80,160)\sigma_{\rm rms}$ levels. Thinner dashed
contours mark $-5\sigma_{\rm rms}$ (they are only noticeable for 150 MHz).
Black crosses show GMOS spectroscopic cluster members within the redshift
range $z=0.516\pm 0.015$. $R_{500}$ is shown by the black circle. The 150 MHz
beam is $43.3\leavevmode\nobreak\ \arcsec\times 18.9\leavevmode\nobreak\
\arcsec$, and 610 MHz is $7.2\leavevmode\nobreak\ \arcsec\times
4.9\leavevmode\nobreak\ \arcsec$.

Figure 2: Same legend as in Figure 1 except that the red contours are 325 MHz.
The 325 MHz beam is $17.5\leavevmode\nobreak\ \arcsec\times
9.5\leavevmode\nobreak\ \arcsec$, and 610 MHz is $7.2\leavevmode\nobreak\
\arcsec\times 4.9\leavevmode\nobreak\ \arcsec$.

###  2.2 XMM-Newton reduction

PLCK G004.5-19.5 was observed by XMM-Newton on 23 March, 2010 (ObsID:
0656201001), for a total duration of 14.2 ks. We reduced the EPIC MOS 1, MOS
2, and pn data using the XMM Science Analysis Software (SAS) v14. We ran the
standard pipeline tasks emproc and epproc for MOS and pn data, respectively.
For each of the three detectors, we filtered the data for solar flare events,
following the same method as described in Mernier et al. (2015). In the 10–12
keV band of MOS and the 12–14 keV band of pn, we fitted count-rate histograms
binned in 100 s intervals with a Poissonian curve. We excluded all the time
intervals in which the count-rate lay above $\mu+2\sqrt{\mu}$, where $\mu$ is
the mean of the Poissonian distribution. We repeated the procedure for the
0.3–10 keV band in the three instruments with histograms binned in 10 s
intervals, since De Luca & Molendi (2004) reported that soft flare events can
also affect soft X-ray bands. As recommended by the calibration reports, we
kept only the highest quality events (flag$=0$) for the three instruments.
Only single events were are allowed in pn (pattern$=0$), while we allowed
single, double, triple, and quadruple events in MOS (pattern$\leq 12$).

We extracted the MOS 1, MOS 2, and pn images of PLCK G004.5-19.5 within four
distinct energy bands (0.3–2 keV, 2–4.5 keV, 4.5–7.5 keV, and 7.5–12 keV), and
extracted similar “background” images, taken from stacked filter-wheel closed
observations222XMM-Newton SOC website (http://xmm.esac.esa.int).. We scaled
the average count-rate of the “background” images to the average count-rate of
the EPIC images within 10–12 keV, where negligible cluster emission is
expected. The EPIC images in each band, and for each instrument, were then
subtracted from this instrumental background and corrected for vignetting
effects, after dividing them by their respective exposure maps generated using
eexpmap.

The background- and vignetting-corrected images were then combined into one
full EPIC image (Figs 1–Fig. 3).

Figure 3: XMM-Newton EPIC (MOS+pn) background- and exposure-map corrected
image of PLCK G004.5-19.5 (0.3–12 keV). The point-spread function of the EPIC
cameras is $\approx 6\leavevmode\nobreak\ \arcsec$. Overlaid are the
extraction regions (wedges) and point-source masks (circles) used for the
radial profile measurement.

###  2.3 Spectroscopic data

We complement our X-ray and radio analysis with spectroscopic measurements of
galaxy redshifts derived from two sets of observations. In addition to the six
redshifts published by Sifón et al. (2014), we retrieved spectroscopic
observations performed with the Focal reducer and low-dispersion spectrograph
(FORS2) on the Very Large Telescope (VLT; Program ID: 090.A-0925(A), PI: H.
Böhringer) as part of a spectroscopic follow-up of Planck $z>0.5$ cluster
candidates. We reduced these data with the ESO Reflex pipeline, which performs
standard reduction steps such as bias subtraction, flat-field correction and
wavelength calibration (Freudling et al. 2013). We measured galaxy redshifts
by cross-correlating the resulting spectra with template galaxy spectra from
the Sloan Digital Sky Survey using the rvsao software (Kurtz & Mink 1998).
Combining there two data sets, we find 16 galaxies in the redshift range
0.50–0.53 (corresponding to a range of 3,000 $\mathrm{km}\,\mathrm{s}^{-1}$)
that are isolated in redshift space and have a velocity dispersion of $\sim
1200\leavevmode\nobreak\ \mathrm{km}\,\mathrm{s}^{-1}$. The cluster redshift
from the combined 16 galaxies is $z_{\mathrm{cl}}=0.519$.

##  3 Results and discussion

###  3.1 X-ray structure

Our processed EPIC image of PLCK G004.5-19.5 (Fig. 3) suggests that the
cluster is disturbed. Assuming $z=0.516$, we find an offset $\sim
100\leavevmode\nobreak\ \mathrm{kpc}$ between the X-ray peak and the X-ray
centroid, estimated from a circular aperture of $R_{\rm
ap}=500\leavevmode\nobreak\ \mathrm{kpc}$ on the X-ray peak.

The degree of cluster disturbance can be provided in a more quantitative way,
by estimating at least two morphological parameters (described in detail in
Cassano et al. 2010b).

  1. 1.

The centroid shift, $w$ (Poole et al. 2006; Maughan et al. 2008), is estimated
from fitting 2D $\beta$-profiles in a series of $N$ circular apertures
$\xi_{i}R_{\text{ap}}$ (with $0.05\leq\xi_{i}\leq 1$) centred on the X-ray
peak, and can be expressed as

$w=\frac{1}{R_{\text{ap}}}\times\bigg{[}\frac{1}{N-1}\sum_{i=1}^{N}\,(\Delta_{i}-\langle\Delta\rangle)^{2}\bigg{]}^{1/2},$
(1)

where $\Delta_{i}-\langle\Delta\rangle$ is the difference between the X-ray
peak and the centroid estimated in the $i$th aperture. Although Cassano et al.
(2010b) proceeded from $0.05R_{\rm ap}$ to $R_{\rm ap}$ in 5% steps, here the
larger PSF of EPIC does not allow us to explore the image with such precision.
Therefore, here we chose larger aperture radii of $0.2R_{\rm ap}$, $0.4R_{\rm
ap}$, $0.6R_{\rm ap}$, $0.8R_{\rm ap}$, and $R_{\rm ap}$.

  2. 2.

The concentration parameter, $c$ (Santos et al. 2008), is defined as the ratio
of the peak over the ambient surface brightness $S$:

$c=\frac{S(r<100\leavevmode\nobreak\
\mathrm{kpc})}{S(r<500\leavevmode\nobreak\ \mathrm{kpc})}.$ (2)

Based on our EPIC image, we find $w=0.058\pm 0.007$ and $c=0.079\pm 0.004$.
Referring to Fig. 1a of Cassano et al. (2010b), it appears that PLCK
G004.5-19.5 is situated within the lower-right quadrant where disturbed galaxy
clusters are found that host radio halos. This suggests that PLCK G004.5-19.5
is disturbed, although we note that PLCK G004.5-19.5 lies further in redshift
than all of the clusters used in the derivation of the relation.

We note a small nodule of previously unidentified X-ray emission to the south-
east of the main X-ray profile in Fig. 6. A clustering of red elliptical
galaxies, as well as a spectroscopic member, suggests that this nodule is at
the cluster redshift. Assuming the nodule is purely thermal emission, with
$kT$ ranging from (3–10) keV, we find unabsorbed luminosities in the (0.1–2.4)
keV energy band, within a circular region of $31\arcsec$ (192 kpc) and
excluding a point-like source and main cluster X-ray emission333The south-east
nodule falls between two CCD chips for the EPIC pn instrument, therefore we
were unable to calculate its luminosity.,

$\displaystyle L_{X}^{\rm MOS1}=$ $\displaystyle(6.7\pm 1.4)\times
10^{43}\leavevmode\nobreak\ \mathrm{erg}\,s^{-1}$ $\displaystyle L_{X}^{\rm
MOS2}=$ $\displaystyle(5.8\pm 1.3)\times 10^{43}\leavevmode\nobreak\
\mathrm{erg}\,s^{-1},$

These luminosities are consistent with typical values for small galaxy groups
and sub-clusters (Reiprich & Böhringer 2002). Taking the same aperture,
$31\arcsec$, we calculate a main-cluster luminosity of $(9.2\pm 1.4)\times
10^{44}\leavevmode\nobreak\ \mathrm{erg}\,s^{-1}$. This is 14–16 times
brighter than the luminosity of south-eastern nodule.

In Fig. 4 we show the radial X-ray (0.3–12) keV surface brightness profile of
PLCK G004.5-19.5 in the south-east and north-west, centred on the X-ray
brightness peak (see the mask in Fig. 3). We see a significant uniform over-
brightness in the south-east relative to north-west. This supports that the
ICM is disturbed. We also note an excess located 160 ″to the south-east that
is associated with the X-ray nodule.

Figure 4: Radial surface brightness profile of XMM-Newton EPIC background- and
exposure-map corrected image of PLCK G004.5-19.5 (0.3–12) keV in the NW and SE
directions, centred on the X-ray brightness peak.

###  3.2 Radio emission

We are able to detect several diffuse and compact radio sources, which are
labelled in bold in Figures 1 and 2. We performed Gaussian source extraction
using PyBDSM v1.8.7 with a detection threshold of $5\sigma_{\rm rms}$. For
point sources lacking a detection at one frequency, we measured the flux
within an aperture of one beam. For one resolved source, we measured the flux
within a common aperture, and we discuss this procedure below. The flux
measurements are given in Table 5, with nondetection measurements labelled
accordingly.

We take the flux variance to be the sum of the measurement variance and a
systematic variance equal to 15%444Suggestion of H. T. Intema. of the flux,
that is, $\sigma_{S}^{2}=\sigma_{\rm rms}^{2}A/A_{\rm beam}+(0.15S)^{2}$,
where $S$, $A$, and $A_{\rm beam}$ are the measured flux, source area, and
beam area.

We modelled the flux with two parameters, the spectral index $\alpha$ and
$S_{\rm 1.4\leavevmode\nobreak\ GHz}$, according to

$\displaystyle S(\nu)=S_{\rm 1.4\leavevmode\nobreak\ GHz}\left(\frac{\nu}{\rm
1.4\leavevmode\nobreak\ GHz}\right)^{\alpha}.$ (3)

We used the Metropolis-Hastings algorithm to sample the posterior distribution
of the model. This is a common straightforward way to sample from the
posterior distribution of a model given only knowledge of the a priori
distribution and a method of computing the observables. We report the
resulting posterior model parameters and uncertainties in Table 5 with
Gaussian $\alpha$ and log-normal $S_{\rm 1.4\leavevmode\nobreak\ GHz}$. The
resulting model for each source with $1\sigma$ bands is shown in Fig. 5.

Table 1: Integrated spectral fluxes, spectral index, and rest-frame $P_{\rm
1.4\leavevmode\nobreak\ GHz}$ with $\Lambda$CDM and $z=0.516$.

Source | $S_{150\leavevmode\nobreak\ \rm{MHz}}$ | $S_{325\leavevmode\nobreak\ \rm{MHz}}$ | $S_{610\leavevmode\nobreak\ \rm{MHz}}$ | $\alpha$ | $S_{1.4\leavevmode\nobreak\ \rm{GHz}}$ | $P_{1.4\leavevmode\nobreak\ \rm{GHz}}$ | Nature
---|---|---|---|---|---|---|---
| (mJy) | (mJy) | (mJy) |  | (mJy) | ($10^{24}$W/Hz) |
C1+2 | $65\pm 10$ | $8\pm 1$ |  $4.2\pm 0.7$bb$b$Both resolved at 610 MHz; $S_{610\leavevmode\nobreak\ \rm{MHz}}^{\rm C1}=2.4\pm 0.2$ mJy and $S_{610\leavevmode\nobreak\ \rm{MHz}}^{\rm C2}=1.9\pm 0.2$ mJy. | $-1.9\pm 0.4$ | $0.6^{+0.4}_{-0.2}$ | $0.6^{+0.3}_{-0.2}$ | lobes?
NW1 | $160\pm 30$ | $62\pm 9$ | $46\pm 7$ | $-0.9\pm 0.2$ | $19^{+7}_{-5}$ | $13^{+3}_{-3}$ | relic?
NW2 | $150\pm 20$ | $75\pm 10$ | $47\pm 7$ | $-0.8\pm 0.2$ | $23^{+7}_{-5}$ | $15^{+3}_{-3}$ | head-tail
H |  $10\pm 2$cc$c$Resolved source, measurement within $17\leavevmode\nobreak\ \arcsec$ circular aperture. |  $4.0\pm 0.7$cc$c$Resolved source, measurement within $17\leavevmode\nobreak\ \arcsec$ circular aperture. |  $1.2\pm 0.5$cc$c$Resolved source, measurement within $17\leavevmode\nobreak\ \arcsec$ circular aperture. | $-1.2\pm 0.4$ | $0.7^{+0.7}_{-0.3}$ | $0.5^{+0.4}_{-0.2}$ | halo?
E | $55\pm 10$ | $22\pm 3$ | $8\pm 1$ | $-1.4\pm 0.2$ | $2.4^{+0.7}_{-0.6}$ | $2.0^{+0.4}_{-0.4}$ | head-tail
X1 | $40\pm 10$ | $4.6\pm 0.7$ |  $0.5\pm 0.1$aa$a$Point-source, peak flux below detection limit $5\sigma_{\rm rms}$, and aperture flux measurement. | $-3.1\pm 0.3$ | $0.04^{+0.02}_{-0.01}$ | $0.07^{+0.02}_{-0.02}$ | AGN?
X2 |  $10\pm 2$aa$a$Point-source, peak flux below detection limit $5\sigma_{\rm rms}$, and aperture flux measurement. | $5.5\pm 0.9$ | $4.4\pm 0.7$ | $-0.5\pm 0.2$ | $2.7^{+1.0}_{-0.7}$ | $1.6^{+0.4}_{-0.3}$ | point?
S | $30\pm 10$ | $15\pm 2$ | $9\pm 1$ | $-0.8\pm 0.2$ | $4^{+2}_{-1}$ | $3.0^{+0.8}_{-0.6}$ | head-tail
555

Figure 5: The measured fluxes, $S(\nu)$, for the sources with are in blue. The
mean posterior model fluxes are with orange dashed lines with $1\sigma$ shaded
intervals.

We labelled three sources S, X1, and X2 near $R_{500}$. Source S has an
extended radio morphology that is indicative of a tailed radio galaxy, where
the tail direction implies its projected relative velocity with respect to the
ICM (see for e.g. Miley et al. 1974). The tail of source S implies that it is
moving away from the projected cluster centre. It has a spectral index of
$-0.8\pm 0.2$, which is in line with a radio lobe.

Sources X1 and X2 are situated above the X-ray nodule of the SE sub-cluster,
and they have quite differing spectral indices of $-3.1\pm 0.3$ and $-0.5\pm
0.2$, respectively. Figure 6 shows the optical region around the sources. Red
elliptical galaxies cluster strongly around source X1. This suggests that X1
is associated with one of the giant red elliptical galaxies that make up the
sub-cluster, and its steep spectral index suggests that it might be a merger-
reactivated AGN.

Figure 6: X-ray (white contours) sub-nodule region with _Gemini_ gri-colour
image, 320 MHz contours (red), and GMOS spectroscopic sources as cyan circles.

To the east of the main cluster we have labelled an extended source, E, which
has a tailed morphology similar to source S, and also a spectroscopically
confirmed cluster member source near the radio surface brightness peak. The
spectral index of source E is $-1.4\pm 0.2$, which suggests that there is some
significant electron ageing. This suggests that there might be some merger-
induced reacceleration of aged electrons from a previous era of activity, or
alternatively, some flux contamination from the central cluster region. The
orientation of the tail implies that source E is moving towards the projected
centre of the main cluster.

Two compact sources, labelled C1 and C2, are located between source E and the
cluster centre. At 610 MHz, both patches are resolved and have similar fluxes
of $2.4\pm 0.2$ mJy and $1.9\pm 0.2$ mJy, respectively. However, they are
unresolved at 150 MHz and 325 MHz. We therefore calculated the integrated
spectral index of the sum of the two sources (in Table 5 as C1+2), and found
it to be very steep, $-1.9\pm 0.4$. There are spectroscopic cluster members
near the two brightness peaks of C1 and C2, although it is unclear whether
they are counterparts, as there is are also galaxies between C1 and C2 that
are not confirmed as cluster members. The steep spectral index supports that
C1 and C2 are lobes of a galaxy.

To the north-west of the main cluster are two extended sources labelled NW1
and NW2. Figure 7 shows a zoom of this region with 610 MHz contours over the
_Gemini_ gri-colour image.

Figure 7: _Gemini_ gri-colour image with 610 MHz contours at
$(5,10,20,40,80,160)\sigma_{\rm rms}$, and GMOS spectroscopic sources as cyan
circles. We also show a strong lensing arc south of NW1. The brightest cluster
galaxy is the southernmost cyan circle. X-ray contours are white.

Source NW2 has an orientation perpendicular to the X-ray profile, with the
largest linear scale (LLS) of 40″(250 kpc). Its southernmost 610 MHz contour
is 62″(380 kpc) away from the brightest cluster galaxy (BCG;
$19\mathrm{h}17\mathrm{m}05.09\mathrm{s}$, $-33\mathrm{d}31\arcmin
20.9\arcsec$). There is a slight brightening at the southern tip and a nearby
spectroscopically confirmed cluster member (a giant red elliptical galaxy,
shown by the cyan circle in Fig. 7). The integrated spectral index of NW2 is
$-0.8\pm 0.2$.

Source NW1 is oriented along the tangent of the X-ray profile, with LLS of
44″(270 kpc). Its nearest 610 MHz contour is 28″(170 kpc) from the BCG. The
surface brightness along its length is uniform, with little variation compared
to the brightness of source NW2. The integrated spectral index is $-0.9\pm
0.2$.

In the central region, at 150 MHz and 325 MHz, there is a diffuse bridge of
radio emission between the north-western and eastern parts of the cluster,
which we have labelled H because its central location and extended morphology
suggest that it is a radio halo. This source is noticeably extended in a
region containing other compact sources (C1 and C2), therefore we did not use
Gaussian source extraction to measure fluxes for H. After fitting all other
sources with Gaussians and subtracting them from the field we took the flux
within a $17\arcsec$ circular aperture.

We chose this aperture for two reasons. The first is that we required the
aperture to fit completely within the region between C1+2 and NW1 to avoid
flux contamination. In the 150 MHz map, the angular separation between C1+2
and NW1 is approximately 23″because of the lower resolution. At 325 MHz the
separation is approximately 45″. The second reason is that we wish to have at
least one beam of collecting area at all frequencies. The aperture area is
$900\leavevmode\nobreak\ \mathrm{arcsec}^{2}$, which is equivalent to 1 beam
at 150 MHz, 5 beams at 325 MHz, and 23 beams at 610 MHz. These criteria fixed
the circular aperture size to 17″. The spectral index of this extended source
is measured to be $-1.2\pm 0.4$.

###  3.3 Nature of the diffuse emission

In this disturbed galaxy cluster, the extended morphology and tangential
orientation of source NW1 is in line with a radio relic tracing a shock front,
and its spectral index is consistent with recently accelerated electrons. The
X-ray data exposure time, 14.2 ks, is too short to detect a counterpart
temperature jump, and there are no polarisation data. The projected distance
of NW1 from the BCG is relatively small for a radio relic (Vazza et al. 2012),
since shock structure is typically found farther out in the cluster periphery.
A possible explanation is that we are viewing this system from a partially
head-on perspective. In addition to the creation of relics and halos, shocks
in the ICM can revive fossil radio plasma from a previous episode of activity
(Enßlin & Gopal-Krishna 2001) which can have filamentary structure. On the
other hand, we might expect the spectra of NW1 to be steeper if this were the
case (see e.g. sources in van Weeren et al. 2011).

Source NW2 has a morphology similar to sources E and S, with a tail direction
implying it is an infalling cluster member. The spectroscopic cluster member
could be related to the slight brightening in the southern tip of NW2. We
remark that the red elliptical galaxy at the northern tip of NW2 is likely
also a cluster member that might likewise contribute to the extended emission.

The diffuse source H is consistent with a radio halo, in line with the
evidence that these sources are found in dynamically disturbed systems.
Current data suggest that the halo has a spectrum of $-1.2\pm 0.4$, although
higher sensitivity observations are required to derive better constraints.
Measuring halo spectra at this redshift is very important since at $z>0.5$ IC
losses are expected to dominate. Simple modelling predicts that only one-third
of very massive $\sim 1-2\times 10^{15}M_{\odot}$ merging systems at redshift
$\sim 0.5$ (such as PLCK G004.5-19.5) are predicted to have spectra flatter
than $\alpha=-1.5$, with the majority of cases generating halos with much
steeper spectra (Cassano et al. 2010a).

Radio halos can also, if less frequently, have filamentary structure (e.g.
Govoni et al. 2005). Since NW1 has a rather small projected distance on the
sky from the BCG, and since we have no polarisation data and not enough data
for an X-ray shock, we cannot rule out that NW1 is part of the radio halo that
has survived IC losses. However, we rate this as unlikely since we expect IC
losses, which scale as $(1+z)^{4}$, to heavily diminish the flux at high
frequencies, which leads to the USS halo nature. The brightness of NW1 at 610
MHz indeed seems to argue against this. Furthermore, Sifón et al. (2014) found
very bright emission in the region of NW1 in archival 1.4 GHz NVSS. Given that
PLCK G004.5-195 is at a mean redshift of 0.516, we would not expect remnant
halo flux to be bright at higher frequencies.

##  4 Conclusion

We have performed high-resolution 610 MHz, 325 MHz, and 150 MHz observations
of the cluster PLCK G004.5-19.5 with the GMRT. We found a complex arrangement
of central and peripheral diffuse emission, as well as several head-tail
galaxies.

We reprocessed archival _XMM-Newton_ data and found morphological parameters
that imply that PLCK G004.5-19.5 is a disturbed system (Cassano et al. 2010b).
Furthermore, a previously unidentified south-eastern sub-cluster is observed
with X-ray unabsorbed luminosity $(6.7\pm 1.4)\times
10^{43}\leavevmode\nobreak\ \mathrm{erg}\,s^{-1}$, approximately a factor
14–16 lower than the luminosity of them main cluster in the same aperture.

Most remarkably, we found a partially face-on radio relic (NW1) and centrally
located radio halo (H) with spectral index ($-1.2\pm 0.4$). Our interpretation
of this finding favours the dynamically young state suggested by the X-ray
morphology parameters, as well as the morphology of the radio emission.

PLCK G004.5-19.5 is the fifth cluster beyond $z=0.5$ known to host a radio
relic, and it is one of the few known instances of a ($z>0.5$) radio halo.

###### Acknowledgements.

J. A., A. S., and H. R. acknowledge financial support from an NWO top subsidy
(614.001.006). C. S. acknowledges support from the European Research Council
under FP7 grant number 279396. We thank the staff of the GMRT, who have made
these observations possible. The GMRT is run by the National Centre for Radio
Astrophysics of the Tata Institute of Fundamental Research. This work is based
on observations obtained with XMM-Newton, an ESA science mission with
instruments and contributions directly funded by ESA member states and the USA
(NASA). The SRON Netherlands Institute for Space Research is supported
financially by NWO, the Netherlands Organisation for Space Research. G. B.
acknowledges partial support from grant PRIN-INAF2014. J. A. thanks R. Cassano
for private discussion on use of morphological parameters at higher redshifts.

## References

  * Bonafede et al. (2012) Bonafede, A., Brüggen, M., van Weeren, R., et al. 2012, MNRAS, 426, 40
  * Bonafede et al. (2009) Bonafede, A., Feretti, L., Giovannini, G., et al. 2009, A&A, 503, 707
  * Botteon et al. (2016) Botteon, A., Gastaldello, F., Brunetti, G., & Kale, R. 2016, MNRAS[arXiv:1607.04641]
  * Briggs (1995) Briggs, D. 1995, PhD thesis, New Mexico Institute of Mining Technology
  * Brunetti (2016) Brunetti, G. 2016, Plasma Physics and Controlled Fusion, 58, 014011
  * Brunetti et al. (2008) Brunetti, G., Giacintucci, S., Cassano, R., et al. 2008, Nature, 455, 944
  * Brunetti & Jones (2014) Brunetti, G. & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007
  * Cassano et al. (2010a) Cassano, R., Brunetti, G., Röttgering, H. J. A., & Brüggen, M. 2010a, A&A, 509, A68
  * Cassano et al. (2006) Cassano, R., Brunetti, G., & Setti, G. 2006, MNRAS, 369, 1577
  * Cassano et al. (2010b) Cassano, R., Ettori, S., Giacintucci, S., et al. 2010b, ApJ, 721, L82
  * Cornwell et al. (2008) Cornwell, T. J., Golap, K., & Bhatnagar, S. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 647
  * De Luca & Molendi (2004) De Luca, A. & Molendi, S. 2004, A&A, 419, 837
  * Donnert et al. (2013) Donnert, J., Dolag, K., Brunetti, G., & Cassano, R. 2013, MNRAS, 429, 3564
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973
  * Enßlin & Gopal-Krishna (2001) Enßlin, T. A. & Gopal-Krishna. 2001, A&A, 366, 26
  * Freudling et al. (2013) Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96
  * Govoni et al. (2005) Govoni, F., Murgia, M., Feretti, L., et al. 2005, A&A, 430, L5
  * Healey et al. (2007) Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61
  * Intema (2014) Intema, H. T. 2014, SPAM: Source Peeling and Atmospheric Modeling, Astrophysics Source Code Library
  * Intema et al. (2016) Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2016, ArXiv e-prints [arXiv:1603.04368]
  * Intema et al. (2009) Intema, H. T., van der Tol, S., Cotton, W. D., et al. 2009, A&A, 501, 1185
  * Kang & Ryu (2013) Kang, H. & Ryu, D. 2013, ApJ, 764, 95
  * Kurtz & Mink (1998) Kurtz, M. J. & Mink, D. J. 1998, PASP, 110, 934
  * Lindner et al. (2014) Lindner, R. R., Baker, A. J., Hughes, J. P., et al. 2014, ApJ, 786, 49
  * Malkov & O’C Drury (2001) Malkov, M. A. & O’C Drury, L. 2001, Reports on Progress in Physics, 64, 429
  * Maughan et al. (2008) Maughan, B. J., Jones, C., Forman, W., & Van Speybroeck, L. 2008, ApJS, 174, 117
  * Mernier et al. (2015) Mernier, F., de Plaa, J., Lovisari, L., et al. 2015, A&A, 575, A37
  * Miley et al. (1974) Miley, G. K., van der Laan, H., & Wellington, K. J. 1974, in IAU Symposium, Vol. 58, The Formation and Dynamics of Galaxies, ed. J. R. Shakeshaft, 109
  * Nuza et al. (2012) Nuza, S. E., Hoeft, M., van Weeren, R. J., Gottlöber, S., & Yepes, G. 2012, MNRAS, 420, 2006
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A29
  * Planck Collaboration et al. (2011) Planck Collaboration, Aghanim, N., Arnaud, M., et al. 2011, A&A, 536, A9
  * Poole et al. (2006) Poole, G. B., Fardal, M. A., Babul, A., et al. 2006, MNRAS, 373, 881
  * Reiprich & Böhringer (2002) Reiprich, T. H. & Böhringer, H. 2002, ApJ, 567, 716
  * Riseley et al. (2017) Riseley, C. J., Scaife, A. M. M., Wise, M. W., & Clarke, A. O. 2017, A&A, 597, A96
  * Roettiger et al. (1999) Roettiger, K., Stone, J. M., & Burns, J. O. 1999, ApJ, 518, 594
  * Santos et al. (2008) Santos, J. S., Rosati, P., Tozzi, P., et al. 2008, A&A, 483, 35
  * Sifón et al. (2014) Sifón, C., Menanteau, F., Hughes, J. P., Carrasco, M., & Barrientos, L. F. 2014, A&A, 562, A43
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * van Weeren et al. (2011) van Weeren, R. J., Röttgering, H. J. A., & Brüggen, M. 2011, A&A, 527, A114
  * van Weeren et al. (2009) van Weeren, R. J., Röttgering, H. J. A., Brüggen, M., & Cohen, A. 2009, A&A, 505, 991
  * Vazza et al. (2012) Vazza, F., Brüggen, M., Gheller, C., & Brunetti, G. 2012, MNRAS, 421, 3375