revtex4-2Repair the float figurec

††thanks: Equal contributions.††thanks: Equal contributions.

# Pareto optimal fronts of kinetic proofreading

Davide Chiuchiu  Biological Complexity Unit, Okinawa Institute of Science and
Technology Graduate University, Onna, Okinawa 904-0495, Japan  Shrabani Mondal
Biological Complexity Unit, Okinawa Institute of Science and Technology
Graduate University, Onna, Okinawa 904-0495, Japan  Department of Chemistry,
Physical Chemistry Section, Jadavpur University, Kolkata, 700032, India
Simone Pigolotti<EMAIL_ADDRESS>Biological Complexity Unit, Okinawa
Institute of Science and Technology Graduate University, Onna, Okinawa
904-0495, Japan

###### Abstract

Biological processes such as DNA replication, RNA transcription, and protein
translation operate with remarkable speed and accuracy in selecting the right
substrate from pools of chemically identical molecules. This result is
obtained by non-equilibrium reactions that dissipate chemical energy. It is
widely recognized that there must be a trade-off between speed, error, and
dissipation characterizing these systems. In this paper, we quantify the
trade-off between speed, error, and dissipation using tools from mathematical
optimization theory. We characterize the Pareto optimal front for a
generalization of Hopfield’s kinetic proofreading model, which is a
paradigmatic example of biological error correction. We find that models with
more proofreading steps are characterized by better trade-offs. Furthermore,
we numerically study scaling relations between speed, accuracy, and
dissipation on the Pareto front.

Subject Areas: Biological Physics, Statistical Physics

##  I Introduction

In living cells, information encoded in the DNA is constantly transcribed into
RNA, which in turn is translated into proteins. Moreover, this information
must be reliably copied into new DNA before a cell division occurs. These
processes are fundamental in biology; their speed and accuracy have a profound
impact on the organism fitness [1, 2, 3]. Evolution must therefore have shaped
these processes to achieve high performance.

It has long been recognized that the accuracy of these reactions, let alone
their speed, can not be achieved close to thermodynamic equilibrium [4]. This
means that information-replicating enzymes must necessarily catalyze non-
equilibrium chemical reactions. A paradigmatic example is the kinetic
proofreading scheme independently proposed by Hopfield [5] and Ninio [6].
Kinetic proofreading can lead to an accuracy significantly higher than the
equilibrium one, at the cost of dissipating chemical energy. This idea
profoundly impacted our understanding of error correction in biology.

The concept of kinetic proofreading sparked an interest in experimentally
characterizing the reaction networks of replicating enzymes and measuring
their kinetic rates [7, 8, 9]. These studies have revealed that these networks
usually include several intermediate steps and are more complex than the
kinetic proofreading model. Unfortunately, these intermediate states are
usually very short-lived [7] and, therefore, hard to observe experimentally.

Alternative approaches attempt to characterize the performance of replicating
enzymes without relying on knowledge of the underlying kinetic details. Often,
these approaches invoke some form of optimality principle [10]. For example,
speed and accuracy of in vitro translation are simultaneously affected by
altering the concentration of magnesium ions in the assay [11, 12]. This
observation has led to the idea that a trade-off exists between accuracy and
speed. Besides speed and accuracy, energy dissipation is another important
property that biological systems can optimize [13]. Several studies have
focused on trade-offs among speed, accuracy, and dissipation in biological
error correction [14, 15, 16, 17, 18, 19]. However, these trade-offs might
depend on the choice of the parameter being tuned [20, 21]. More fundamental
bounds on speed, error, and dissipation of replicating enzymes are set by the
second law of thermodynamics [22, 23, 24]. But despite their theoretical
interest, these bounds are usually far from the operating regimes of
replicating enzymes. This suggests that the formulation of optimality
principles and trade-offs requires concrete implementations of the reaction
networks. Several studies investigated tradeoffs between pairs of observables
[25, 26]. In particular, a recent study [27] has theoretically derived a
trade-off between error and dissipation. However, the modeling assumptions
made in [27] are such that the process speed can be varied independently of
error and dissipation and is therefore not subject to a tradeoff.

In this paper, we numerically study trade-offs between speed, error, and
dissipation using the concept of a Pareto front. We consider, as a
paradigmatic example, a multi-step generalization of the original Hopfield
model. We find that the model performance significantly improves at increasing
the number of intermediate proofreading steps. We also characterize scaling
relations between speed, error, and dissipation on the Pareto front.

##  II Pareto optimal front

We illustrate the idea of a Pareto front with an example in which we
simultaneously optimize the speed $v$, the error rate $\eta$, and the
dissipation per incorporated monomer $\Delta\sigma$, as a function of the
kinetic rates. Kinetic constraints determine feasible combinations of speed,
error, and dissipation. An example of such constraints is the fact that the
reaction network is driven out of thermodynamic equilibrium by a finite energy
budget. The feasible combinations $(v,\eta,\Delta\sigma)$ form a set, that we
call $F$. As an illustration, the set of feasible combinations for a
generalization of the Hopfield model is shown in Figure 1(a). The Pareto front
$P$ is the subset of $F$ made up of the optimal configurations, i.e., those in
which one observable improves only at the expense of other observables [28].
Here, by “improvement of observable” we mean an increase of the speed and a
decrease of the dissipation and the error. Formally, a triplet
$(v,\eta,\Delta\sigma)$ belongs to the Pareto front if there is no other
triplet $(v^{\prime},\eta^{\prime},\Delta\sigma^{\prime})\in F$ such that
$v^{\prime}\geq v$, $\eta^{\prime}\leq\eta$, and
$\Delta\sigma^{\prime}\leq\Delta\sigma$. Feasible combinations of observables
values that do not lie on the Pareto optimal front constitute sub-optimal
solutions where a suitable parameter change can improve some observables
without penalizing others.

For simplicity of illustration, we show feasible configurations and Pareto
fronts for pairs of observables, see Figures 1(b)–(d). The definitions of
feasible configurations and Pareto fronts for pairs of observables are
obtained by marginalization from the three-dimensional case. For example, a
pair $(v,\eta)$ is a feasible speed-error pair if $(v,\eta,\Delta\sigma)\in F$
for some $\Delta\sigma$. The corresponding Pareto front is the set of feasible
configurations $(v,\eta)$ such that there are no other pairs
$(v^{\prime},\eta^{\prime})$ such that $v^{\prime}\geq v$ and
$\eta^{\prime}\leq\eta$. Pareto fronts for pairs of observables are
represented as red curves in Figures 1(b)–(d).

The Pareto front constitutes only a part of the boundary of the set of
feasible configurations. In the example of Figure 1(c), the boundary on the
right side of the figure is not part of the Pareto front since dissipation is
not optimal along this line.

Figure 1: Feasible values of speed, error, and dissipation for the generalized
Hopfield network with four proofreading steps. (a) Feasible configuration in
the speed-error-dissipation phase space. Two-dimensional projections of the
speed-error-dissipation phase space for (b) speed-error, (c) speed-
dissipation, and (d) dissipation-error. The Pareto optimal front is
represented by a red curve. Model details and parameter values are discussed
in Sections III and IV.

##  III Generalized Hopfield model

### Definition of the model

Figure 2: Generalized Hopfield model with $N$ intermediate steps. a) Energy
landscape of the generalized Hopfield model. The depth $\Delta E_{x}$ of the
first valley represents the binding energy of monomer $x$. The following
states are characterized by positive energy increments $\Delta
E_{{\mathrm{conf}}}^{i}$. Forward reactions are driven out of equilibrium by
energies $\Delta E_{{\mathrm{hyd}}}^{i}$ originating from hydrolysis. We
assume that neither $\Delta E_{{\mathrm{conf}}}^{i}$ nor $\Delta
E_{{\mathrm{hyd}}}^{i}$ depend on the monomer type. The energy of the final
state is not shown because we assume that the last reaction is fully
irreversible. b) Chemical reaction network of the generalized Hopfield model.
The enzyme goes through $N$ intermediate steps before incorporating monomer
$x$. The reactions leading to the intermediate states $m_{2}\ldots m_{N}$ are
assisted by an external nonequilibrium driving. In addition, $N-1$
proofreading pathways permit to reject the candidate monomers from the
intermediate states $m_{2}\ldots m_{N}$. Rate constants and energy values that
are kept fixed, rather than varied in the optimization process, are shown in
red.

We study a generalized version of the Hopfield model [5] including multiple
proofreading steps. A replicative enzyme $E$ can bind to free monomers of type
$x\in\\{r,w\\}$, where $r$ represents a right monomer and $w$ a wrong one. The
enzyme may then reject the monomer based on its binding free energy $\Delta
E_{x}$. In the incorporation pathway, the enzyme sequentially adopts $N$
configurations $m_{i}^{x}$, $i=1\dots N$ before the final incorporation of
monomer $x$, where $i=1$ denotes the initial binding, see Figure 2(a).

We assume that the free energy changes $\Delta E^{i}_{{\mathrm{conf}}}$ when
going from configuration $i$ to configuration $i+1$ do not depend on the
monomer type. Here and in the following, all energies are expressed in units
of the thermal energy $k_{B}T$, where $k_{B}$ is the Boltzmann constant and
$T$ is the temperature of the environment. From the configurations $i=2\dots
N$, monomer $x$ can be discarded via a proofreading reaction that brings the
replicative enzyme back to its initial free state. Proofreading rates depend
on the energy of intermediate states and the monomer type. Similar to the
original Hopfield model, the enzyme consumes chemical energy from hydrolysis
to assist the configuration changes. The enzyme absorbs this energy from a
single hydrolysis event. However, we assume that this energy can be
partitioned among the multiple steps leading to each of the intermediate
states. Thus, the energy $\Delta E_{{\mathrm{hyd}}}^{i}$ consumed in going
from the $(i-1)$th to the the $i$th configuration satisfies

$\sum_{i=1}^{N-1}\Delta E_{{\mathrm{hyd}}}^{i}=\Delta
E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}$ (1)

where $\Delta E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}$ is the (fixed) total energy
budget that the enzyme gains from one hydrolysis event. We also assume that
all reactions except the final incorporation are reversible. The reaction
rates are determined by the energy landscape in Figure 2(a).

Specifically, the concentrations $[m_{i}^{x}]$ of the intermediate states
evolve as

$\displaystyle\begin{aligned} \frac{d[m_{1}^{x}]}{dt}=&k^{0}+e^{\Delta
E_{{\mathrm{conf}}}^{1}}k^{1}[m_{2}^{x}]\\\ &\,-\left(e^{\Delta
E_{hyd}^{1}}\,k^{1}+e^{-\Delta E_{x}}\,k^{0}\right)[m_{1}^{x}]\\\
\frac{d[m_{i}^{x}]}{dt}=&\,k^{i,p}+e^{\Delta
E_{{\mathrm{hyd}}}^{i-1}}\,k^{i-1}[m_{i-1}^{x}]+e^{\Delta
E_{{\mathrm{conf}}}^{i}}k^{i}[m_{i+1}^{x}]\\\ &\,-\left(e^{\Delta
E_{{\mathrm{hyd}}}^{i}}\,k^{i}+\,e^{[-\Delta E_{x}+\sum_{j=1}^{i-1}\Delta
E_{{\mathrm{conf}}}^{j}]}k^{i,p}\right)[m_{i}^{x}]\\\ &+\left(e^{\Delta
E_{{\mathrm{conf}}}^{i-1}}k^{i-1}\right)[m_{i}^{x}]\\\ &\mbox{for }i=2\ldots
N-1\\\ \frac{d[m_{N}^{x}]}{dt}=&\,k^{N-1,p}+e^{\Delta
E_{{\mathrm{hyd}}}^{N-1}}\,k^{N-1}[m_{N-1}^{x}]\\\ &\,-\left(k^{N}+e^{[-\Delta
E_{x}+\sum_{j=1}^{N-1}\Delta
E_{{\mathrm{conf}}}^{j}]}k^{N-1,p}\right)[m_{N}^{x}]\\\ &+\left(e^{\Delta
E_{{\mathrm{conf}}}^{N-1}}k^{N-1}\right)[m_{N}^{x}]\,,\end{aligned}$ (2)

see Figure 2(b). The binding energy $\Delta E_{x}$ depends on the monomer
type, leading to different rejection rates for different monomer types. For
simplicity, we assume that the concentrations of enzyme $[E]$ and monomers
$[x]$ are maintained constant, and we implicitly incorporate them in the
reaction rates involving free enzymes or monomers.

We compute the steady values of the concentrations $[m_{i}^{x}]$ from Eq. (2).
We then express the error $\eta$, the polymerization speed $v$, and the
dissipation per incorporated monomer $\Delta\sigma$ as

$\displaystyle\eta=$
$\displaystyle\,\frac{[m_{N}^{w}]}{[m_{N}^{w}]+[m_{N}^{r}]}\,,$ (3)
$\displaystyle v=$
$\displaystyle\,k^{N}\left([m_{N}^{r}]+[m_{N}^{w}]\right)\,,$
$\displaystyle\Delta\sigma=$
$\displaystyle\,\frac{\sigma^{r}+\sigma^{w}}{v}\,,$

where

$\displaystyle\sigma^{x}=$ $\displaystyle k^{0}\left(1-e^{-\Delta
E_{x}}\left[m_{1}^{x}\right]\right)\log\left(\frac{1}{e^{-\Delta
E_{x}}\left[m_{1}^{x}\right]}\right)$ (4)
$\displaystyle+\sum_{i=1}^{N-1}k^{i}\left(e^{\Delta
E_{{\mathrm{hyd}}}^{i}}\left[m_{i}^{x}\right]-e^{\Delta
E_{{\mathrm{conf}}}^{i}}\left[m_{i-1}^{x}\right]\right)$
$\displaystyle\times\log\left[\frac{e^{\Delta
E_{{\mathrm{hyd}}}^{i}}\left[m_{i}^{x}\right]}{e^{\Delta
E_{{\mathrm{conf}}}^{i}}\left[m_{i-1}^{x}\right]}\right]$
$\displaystyle+\sum_{i=1}^{N-1}k^{i,p}\left(1-e^{[-\Delta
E_{x}+\sum_{j=1}^{i}\Delta
E_{{\mathrm{conf}}}^{j}]}\left[m_{i}^{x}\right]\right)$
$\displaystyle\times\log\left[\frac{1}{e^{[-\Delta E_{x}+\sum_{j=1}^{i}\Delta
E_{{\mathrm{conf}}}^{j}]}\left[m_{i}^{x}\right]}\right]$

is the average entropy production rate depending on the monomer type. This
definition of the average entropy production rate does not include the
contribution of the irreversible final incorporation step, as discussed in
Ref. [20].

We take as free parameters the rate constants $k_{1},k_{2},\dots,k_{N}$. We
assume that the initial binding rate $k_{0}$ is determined by diffusion and
monomer concentration and thus can not be optimized. We therefore set
$k_{0}=1$. Further free parameters are the hydrolysis energies at each step
$\Delta E_{{\mathrm{hyd}}}^{i}$, $i=1\dots N-1$, and the intrinsic
proofreading rates $k^{2,p}\dots k^{N,p}$. The total number of free parameters
is thus $3N-2$, which reduces to $3N-3$ by considering that the hydrolysis
energies are constrained by Eq. (1).

In the optimization, we also impose that the direction of the average flux
associated with each reaction is consistent with a process of polymer
synthesis:

$\displaystyle 1-e^{-\Delta E_{x}}[m_{1}^{x}]\geq$ $\displaystyle 0$ (5)
$\displaystyle e^{\Delta E_{{\mathrm{hyd}}}^{i}}k^{i}[m_{i}^{x}]-e^{\Delta
E_{{\mathrm{conf}}}^{i}}k^{i}[m_{i+1}^{x}]\geq$ $\displaystyle 0$
$\displaystyle\mbox{for }i=1\ldots N-1$ $\displaystyle k^{i,p}-e^{[-\Delta
E_{x}+\sum_{j=1}^{i}\Delta E_{{\mathrm{conf}}}^{j}]}k^{i,p}[m_{i+1}^{x}]\leq$
$\displaystyle 0$ $\displaystyle\mbox{for }i=1\ldots N-1\,.$

Additionally, we require that the enzyme has a fixed energy budget for each
incorporated monomer via Eq. (1), and we impose that

$\frac{k^{i+1,p}}{k^{i,p}}<1\quad\mbox{for }i=1\ldots N-1$ (6)

so that the states $m_{i}^{x}$ have increasing energy as shown in Figure 2.

Figure 3: Additional proofreading steps improve Pareto optimal fronts. a)
Three-dimensional Pareto front between $\eta$, $v$ and $\Delta\sigma$ for
different numbers of steps $N$. Each blue point corresponds to the speed and
error of particular sets of rate constants. All the blue points are bounded by
the trade-off curves shown in (b). Holes in the plot indicate that the
optimization algorithm did not find a local minimum consistent with the
definition of a Pareto front. Panels (b), (c), and (d) present trade-off
curves between pairs of observables in $\\{\eta,v,\Delta\sigma\\}$. The color
dashed lines in plots with $\eta$ correspond to the minimum error and the
black dashed lines correspond to the equilibrium error. The black line in
panel (d) represents the scaling law expressed by Eq. (12). Parameters are
$\Delta E_{r}=6.5$ and $\Delta E_{w}=5.21$, which are the averages DNA binding
energies for cognate and non-cognate pairings from ab-initio calculations
[29], and $\Delta E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=100$, which is
sufficiently large value to recover results from the ideal Hopfield model.
Results for $\Delta E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=25$ and $\Delta
E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=12.5$ are qualitatively similar, see
Appendix A (Figure 7).

### Optimization algorithm

We find the optimal solutions by a computational scheme based on multi-
functional optimization [28]. In particular, we systematically explore the
parameter space using genetic algorithms [30] to identify feasible values of
the observables that are compatible with the constraints and find the Pareto
optimal front between $\eta$, $v$ and $\Delta\sigma$ as a function of the free
parameters. We employ this method for different numbers of proofreading steps.

In practice, the coefficients of the linear system satisfied by the
concentrations $[m_{i}^{x}]$ can be of very different magnitude, which might
cause numerical instabilities. To counter this problem, we rescale Eqs. (2),
(3), and (5) by performing the change of variables

$\displaystyle K_{\rightarrow}^{1}=$ $\displaystyle\frac{k^{1}}{k^{0}}$ (7)
$\displaystyle K_{\rightarrow}^{i}=$
$\displaystyle\frac{k^{i}}{e^{\sum_{j=1}^{i-1}\Delta
E_{{\mathrm{conf}}}^{j}}k^{i-1,p}}\quad\mbox{for }i=2\ldots N$ $\displaystyle
K_{\leftarrow}^{i}=$ $\displaystyle\frac{e^{\Delta
E_{{\mathrm{conf}}}^{i}}k^{i}}{e^{\sum_{j=1}^{i}\Delta
E_{{\mathrm{conf}}}^{j}}k^{i,p}}\quad\mbox{for }i=1\ldots N-1$
$\displaystyle[M_{1}^{x}]=$ $\displaystyle[m_{1}^{x}]$
$\displaystyle[M_{i}^{x}]=$
$\displaystyle\frac{[m_{i}^{x}]}{e^{-\sum_{j=1}^{i-1}\Delta
E_{{\mathrm{conf}}}^{j}}}\quad\mbox{for }i=2\ldots N.$

With this prescription, the forward and the reverse rate constants are
rescaled by the corresponding reverse proofreading rate constants, but without
including the contributions from the binding energy and the hydrolysis energy.
The concentrations are scaled by their corresponding Boltzmann factor, i.e.,
with their equilibrium value. We remark that the rescaling defined in Eq. (7)
is invertible, and the rescaled equations have therefore the same number of
free parameters as the original equations.
More details about the rescaled form of the chemical equations are presented
in the Appendix A.

##  IV Results

### Pareto fronts

We compute the Pareto optimal fronts for the generalized Hopfield model with
different numbers of intermediate steps, see Figure 3(a). Theory [5, 6]
predicts that the error ranges from the minimum error

$\eta_{H}=\frac{1}{1+e^{N(\Delta E_{r}-\Delta E_{w})}}$ (8)

up to the equilibrium error

$\eta_{eq}=\frac{1}{1+e^{\Delta E_{r}-\Delta E_{w}}}.$ (9)

In all cases, the Pareto front between pairs of variables is a monotonic curve
by definition, see Figure 3(b), (c), and (d). Speed-error and dissipation-
error trade-off curves are defined in the error range
$\eta\in[\eta_{eq},\eta_{H}]$ as expected. In particular, at the equilibrium
error the speed attains its maximum, and the dissipation tends to zero, see
Figures 3(b), (c). The Pareto fronts for different $N$ appear to tend to a
common limit $v\approx 1$ for $\eta\rightarrow\eta_{eq}$. In this limit, the
reaction speed is limited by the rate of monomer binding $k^{0}$, while all
other reaction steps are much faster. We note that the estimated dissipation
per step steeply drops to zero for $\eta\rightarrow\eta_{eq}$, see Figure
3(c). We expect the exact dependence of $\Delta\sigma$ on $\eta$ to be smooth
around $\eta_{eq}$, as in exactly solvable models (see, e.g., [13,
sartori2013kinetic]. We therefore expect the abrupt jump in Figure 3(c) to be
due to our finite numerical resolution and our choice of representing results
in log scale.

The Pareto fronts between $\eta$ and $v$ and between $\eta$ and $\Delta\sigma$
substantially improve with the number of proofreading steps, see Figure 3(b)
and 3(c), respectively. Instead, the optimal front between $v$ and
$\Delta\sigma$ barely depends on the number of proofreading steps, see Figure
3(d). To understand this observation, we start from the full three-dimensional
Pareto fronts represented in Figure 7(a). If we now minimize $\Delta\sigma$ at
fixed $v$ along these fronts, we find that the minimum is attained for large
values of $\eta$ (i.e., close to $\eta_{eq}$). In this limit, proofreading is
effectively inactivated, as its action would slow down the process and cause
additional dissipation. The monomer incorporation pathway is then reduced to a
linear one, for which varying the number of steps at fixed chemical driving
has little effect on the speed-dissipation Pareto front.

Figure 4: Scaling relation between error speed, and error dissipation in two-
dimensional Pareto fronts for different $N$. Plotting as a function of
$(\eta-\eta_{H})$, we obtain similar curves for the Pareto fronts between (a)
$\eta$ and $v$ and (c) $\eta$ and $\Delta\sigma$. Panels (b) and (d) show data
collapses according to Eqs.(10) and (11), respectively. In both cases, we find
a data collapses for small error rates. In all curves, we fix $\Delta
E_{{\mathrm{hyd}}}=100$.

### Scaling Laws

We study the behavior of the speed and dissipation on the Pareto front for
errors close to the minimum error $\eta_{H}$ given by Eq. (8). Plotting
velocity and dissipation as a function of $(\eta-\eta_{H})$ leads to a good
data collapse, see Figure 4(b, d). However, a scrutiny of the data collapse in
(a) and (c) shows that the slope is different for different $N$. In fact, we
find that, for small errors, the behaviors of the velocity and the dissipation
are well described by scaling relations

$v\sim\left(\frac{\eta-\eta_{H}}{\eta_{H}}\right)^{N}$ (10)

and

$\Delta\sigma\sim\left(\frac{\eta-\eta_{H}}{\eta_{H}}\right)^{1-N},$ (11)

see Figure 4(b) and 4(d), respectively.

Figure 5: Higher hydrolysis energy $\Delta E_{{\mathrm{hyd}}}$ yields better
Pareto fronts. (a), (b) Trade-offs between pairs of observables in
$\\{\eta,v,\Delta\sigma\\}$ for different $\Delta E_{{\mathrm{hyd}}}$. In all
cases, the Pareto fronts improve at higher hydrolysis energy up to a certain
limit when the $\Delta E_{{\mathrm{hyd}}}$ (e.g $\Delta
E_{{\mathrm{hyd}}}=12.5$) is not sufficient for the enzyme to achieve the high
energy state. All curves in this figure are for $N=5$.

Our data suggest a scaling law between $v$ and $\Delta\sigma$ of the form

$v\sim(\Delta\sigma)^{\gamma}$ (12)

with $\gamma=1.85\pm 0.05$, at least for a range of small values of
$\Delta\sigma$, see Figure 3.

Figure 6: Comparison of error-speed trade-off curves for symmetric and general
Hopfield models.

Both speed-error and dissipation-error Pareto fronts improve with $\Delta
E_{{\mathrm{hyd}}}$ and seems to saturate for large $\Delta
E_{{\mathrm{hyd}}}$, see Figure 5. In particular, an increase from $\Delta
E_{{\mathrm{hyd}}}=12.5$ to $\Delta E_{{\mathrm{hyd}}}=25$ leads to
significant improvements. More extensive numerical results on the Pareto
fronts for $\Delta E_{{\mathrm{hyd}}}=12.5$ are presented in Appendix B. On
the other hand, a further increase from $\Delta E_{{\mathrm{hyd}}}=25$ to
$\Delta E_{{\mathrm{hyd}}}=100$ does not lead to appreciable changes. This
result shows that $\Delta E_{{\mathrm{hyd}}}=25$ (a value on the order of the
free energy of ATP hydrolysis) can be already considered as very large, at
least for the values of $N$ that we considered.

### Symmetric choice of parameters.

The rescaling of the parameters introduced in Eq. (7) serves to have rates of
comparable order of magnitude. In this section, we explore what happens if we
assume that the rescaled optimal parameters are independent of the reaction
step. In particular, we assume that $K^{i}_{\to}=K_{\to}$ and
$K^{i}_{\leftarrow}=K_{\leftarrow}$. We also assume that the system consume
equal amount of hydrolysis energy at each step, $\Delta
E^{i}_{{\mathrm{hyd}}}=\Delta E_{{\mathrm{hyd}}}$ for $i=1,\dots N$. As a
result, we are left with only two free parameters, $K_{\to}$ and
$K_{\leftarrow}$, irrespective of $N$. We find that, the Pareto fronts
obtained under this approximation are similar to the general one, see Figure
6. In particular, for $N=2$, the Pareto fronts between speed and error in the
symmetric and general model are indistinguishable in our numerical
simulations. The difference between the Pareto fronts in the two models
increases with the number of steps.

##  V Conclusions

In this work, we numerically studied the Pareto optimal fronts between speed,
error rate, and dissipation in biological error correction. Our algorithm,
based on multi-object optimization theory, permits to reconstruct the Pareto
front with high accuracy.

Our results demonstrate aspects of error corrections that would be difficult
to reveal with other theoretical methods. For example, we find that the
error/speed and error/dissipation trade-offs strongly depend on the number of
proofreading steps. In contrast, the trade-off between speed and dissipation
is rather insensitive to the number of steps. In certain regimes, we found
that speed, error rate, and dissipation along the Pareto front are
characterized by non-trivial scaling laws. Finally, we have found that a
rescaling of the chemical equations, that we introduce to obtain more stable
numerical results, suggests a symmetry assumption on the optimal rates that
leads to a good approximation of the Pareto front. This result is a
potentially useful hint for future theoretical approaches.

We focused our study on a generalization of the Hopfield model, where the
discrimination in the main incorporation pathway is only present in backward
rates. In many real error correction networks, such as those implemented by
ribosomes, this discrimination strategy is complemented by forward
discrimination, i.e., discrimination based on energy barrier differences [7,
15, 25]. The presence of forward discrimination leads to more complex
scenarios, where the choice of parameters to be optimized and global
constraints might play a delicate role. Understanding these cases using the
numerical approach outlined in our work is an interesting venue of
investigation for future studies.

##  Appendix A Computational details of error-speed-dissipation trade-off

The rescaled version of Eq. (2) at steady state are expressed by

$\displaystyle\begin{aligned} 0=&1+\left(-e^{-\Delta E_{x}}-e^{\Delta
E_{{\mathrm{hyd}}}^{1}}K_{\rightarrow}^{1}\right)M_{1}^{x}+K_{\rightarrow}^{1}M_{2}^{x}\\\
0=&1+e^{\Delta E_{{\mathrm{hyd}}}^{i-1}}K_{\leftarrow}^{i-1}[M_{i-1}^{x}]\\\
&+\left(-e^{-\Delta E_{x}}-K_{\leftarrow}^{i-1}-e^{\Delta
E_{{\mathrm{hyd}}}^{i}}K_{\rightarrow}^{i}\right)[M_{i}^{x}]\\\
&+K_{\rightarrow}^{i}[M_{i+1}^{x}]\\\ &\mbox{for }i=2\ldots N-1\\\
0=&1+e^{\Delta E_{{\mathrm{hyd}}}^{N-1}}K_{\leftarrow}^{N-1}[M_{N-1}^{x}]\\\
&+\left(-e^{-\Delta
E_{x}}-K_{\leftarrow}^{N-1}-K_{\rightarrow}^{N}\right)[M_{N}^{x}]\\\
\end{aligned}$ (13)

Eqs. (13) are linear in the rescaled variables $[M_{i}^{x}]$. To compute the
steady concentrations, we express Eq. (13) in matrix form

$\mathbf{P}\cdot\mathbf{M}=\mathbf{Q}$ (14)

where $\mathbf{M}=\\{[M_{1}^{x}]...[M_{N}^{x}]\\}$, $\mathbf{Q}$ is a constant
vector with all entries equal to $-1$, and $\mathbf{P}$ is the matrix of
coefficients of the linear system. For example, for $N=2$, the matrix
$\mathbf{P}$ reads

$\mathbf{P}=\begin{pmatrix}-e^{-\Delta E_{x}}-e^{\Delta
E_{{\mathrm{hyd}}}^{1}}K_{\rightarrow}^{1}&K_{\rightarrow}^{1}\\\ e^{\Delta
E_{{\mathrm{hyd}}}^{1}}K_{\leftarrow}^{1}&-e^{-\Delta
E_{x}}-K_{\leftarrow}^{1}-K_{\rightarrow}^{2}.\end{pmatrix}$ (15)

Next, we invert the rescaled relations to express the original parameter as a
function of the rescaled parameters $K_{\rightarrow}^{i}$ and
$K_{\leftarrow}^{i}$. From Eq. (7) we obtain the inverse relations

$k_{\rightarrow}^{i,p}=k_{\rightarrow}^{0}\prod_{j=1}^{i}\frac{K_{\rightarrow}^{j}}{K_{\leftarrow}^{j}}$
(16)

$k^{1}=k^{0}K_{\rightarrow}^{1}$ (17)

$\displaystyle k_{\rightarrow}^{i}=$ $\displaystyle\,e^{\sum_{j=1}^{i-1}\Delta
E_{{\mathrm{conf}}}^{j}}K_{\leftarrow}^{i}k_{\rightarrow}^{0}\prod_{j=1}^{i}\frac{K_{\rightarrow}^{j}}{K_{\leftarrow}^{j}}$
(18) $\displaystyle\mbox{for }i=2\ldots N-1$ $\displaystyle k^{N}=$
$\displaystyle\,k^{0}e^{\sum_{j=1}^{i-1}\Delta
E_{{\mathrm{conf}}}^{j}}\frac{\prod_{j=1}^{N}K_{\rightarrow}^{j}}{\prod_{j=1}^{N-1}K_{\leftarrow}^{j}}$

We can now write the observables $\eta$, $v$, and $\Delta\sigma$ in terms of
the rescaled rates and concentrations:

$\displaystyle\eta=$
$\displaystyle\,\frac{[M_{N}^{w}]}{[M_{N}^{w}]+[M_{N}^{r}]}$ (19)
$\displaystyle v=$
$\displaystyle\,k^{0}\frac{\prod_{j=1}^{N}K_{\rightarrow}^{j}}{\prod_{j=1}^{N-1}K_{\leftarrow}^{j}}\left([M_{N}^{r}]+[M_{N}^{w}]\right)$
$\displaystyle\Delta\sigma=$ $\displaystyle\,\frac{\sigma^{r}+\sigma^{w}}{v}$

where

$\displaystyle\sigma_{n}^{x}=$ $\displaystyle
k^{0}\left(\sum_{i=1}^{N-1}e^{\Delta
E_{{\mathrm{hyd}}}^{i}}K_{\rightarrow}^{i}\log\left[\frac{e^{-\Delta
E_{{\mathrm{hyd}}}^{i}}M_{i}^{x}}{M_{i+1}^{x}}\right]\right)$ (20)
$\displaystyle\times\left(\left(M_{i}^{x}-e^{\Delta
E_{{\mathrm{hyd}}}^{i}}M_{i+1}^{x}\right)\prod_{j=1}^{i-1}\frac{K_{\rightarrow}^{j}}{K_{\leftarrow}^{j}}\right)$
$\displaystyle+k^{0}\left(\sum_{i=0}^{N-1}-\log[\frac{e^{-\Delta
E_{x}}}{M_{i+1}^{x}}](-1+e^{-\Delta E_{x}}M_{i+1}^{x}\right)$
$\displaystyle\times\prod_{j=1}^{i}\frac{K_{\rightarrow}^{j}}{K_{\leftarrow}^{j}}$

Finally, we feed Eq. (14) and Eq. (19) to the MATLAB function gamultiobj to
obtain three-dimensional Pareto fronts for $\eta$, $v$, and $\Delta\sigma$.
The function gamultiobj uses genetic algorithms [30] and numerical evaluations
of Eq. (19) to find the combination of free parameters corresponding to the
dominant configurations of $\eta$, $v$, and $\Delta\sigma$.

The rescaled problem can still generate badly scaled matrices in some rare
instances, which we further mitigate by implementing arbitrary precision
computations with the Advanpix multi-precision package for MATLAB [31]. In
this way, we efficiently obtain numerically stable results.

##  Appendix B Results for different values of the hydrolysis energy.

We computed the Pareto front by setting the hydrolysis energy $\Delta
E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=12.5$ and $\Delta
E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=25$ (see Figure 7) Figure 7 is qualitative
similar to the case for $\Delta E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=100$ (see
Figure 3) and agree with our observation that additional proofreading steps
improve Pareto fronts. However, for $\Delta
E_{{\mathrm{hyd}}}^{{\mathrm{tot}}}=12.5$ Figure 7(a-d), minimum errors
obtained from the simulation are substantially larger than the minimum errors
defined in Eq. 8, in particular for $N>3$.

Figure 7: Pareto optimal fronts (a-d) for $\Delta E_{{\mathrm{hyd}}}=12.5$ and
(e-h) for $\Delta E_{{\mathrm{hyd}}}=25$ (a),(e) 3-Dimensional Pareto optimal
front between $\eta$, $v$ and $\Delta\sigma$ for different number of steps
$N$. (b), (c), (d), (f), (g), and (h) marginalized trade-off curves between
pairs observables in $\\{\eta,v,\Delta\sigma\\}$. Dashed lines in plots with
$\eta$ correspond to the ideal error of the Hopfield model with $N$
intermediate steps.

###### Acknowledgements.

We thank Deepak Bhat for useful discussions. SP was supported by JSPS KAKENHI
Grant Number JP18K03473 and by the Okawa Foundation (Grant Number 21-01).

## References

  * Zaher and Green [2009] H. S. Zaher and R. Green, Fidelity at the molecular level: Lessons from protein synthesis, Cell 136, 746 (2009).
  * Johnson [1993] K. A. Johnson, Conformational coupling in dna polymerase fidelity, Annual Review of Biochemistry 62, 685 (1993).
  * Savir and Tlusty [2013a] Y. Savir and T. Tlusty, The ribosome as an optimal decoder: A lesson in molecular recognition, Cell 153, 471 (2013a).
  * Pauling [1957] L. Pauling, The probability of errors in the process of synthesis of protein molecules (1957).
  * Hopfield [1974] J. J. Hopfield, Kinetic proofreading: A new mechanism for reducing errors in biosynthetic processes requiring high specificity, Proceedings of the National Academy of Sciences 71, 4135 (1974).
  * Ninio [1975] J. Ninio, Kinetic amplification of enzyme discrimination, Biochimie 57, 587 (1975).
  * Rodnina  _et al._ [2017] M. V. Rodnina, N. Fischer, C. Maracci, and H. Stark, Ribosome dynamics during decoding, Philosophical Transactions of the Royal Society B: Biological Sciences 372, 20160182 (2017).
  * Dieckman  _et al._ [2010] L. M. Dieckman, R. E. Johnson, S. Prakash, and M. T. Washington, Pre-steady state kinetic studies of the fidelity of nucleotide incorporation by yeast dna polymerase delta, Biochemistry 49, 7344 (2010).
  * Goodman  _et al._ [1993] M. F. Goodman, S. Creighton, L. B. Bloom, J. Petruska, and D. T. A. Kunkel, Biochemical basis of dna replication fidelity, Critical Reviews in Biochemistry and Molecular Biology 28, 83 (1993).
  * Savir and Tlusty [2013b] Y. Savir and T. Tlusty, The ribosome as an optimal decoder: A lesson in molecular recognition, Cell 153, 471 (2013b).
  * Johansson  _et al._ [2012] M. Johansson, J. Zhang, and M. Ehrenberg, Genetic code translation displays a linear trade-off between efficiency and accuracy of trna selection, Proceedings of the National Academy of Sciences 109, 131 (2012).
  * Vorstenbosch  _et al._ [1996] E. Vorstenbosch, T. Pape, M. Rodnina, B. Kraal, and W. Wintermeyer, The g222d mutation in elongation factor tu inhibits the codon-induced conformational changes leading to gtpase activation on the ribosome., The EMBO journal 15, 6766 (1996).
  * Bennett [1979] C. H. Bennett, Dissipation-error tradeoff in proofreading., Bio Systems 11 2-3, 85 (1979).
  * Savageau and Freter [1979] M. A. Savageau and R. R. Freter, Energy cost of proofreading to increase fidelity of transfer ribonucleic acid aminoacylation, Biochemistry 18, 3486 (1979).
  * Ehrenberg and Blomberg [1980] M. Ehrenberg and C. Blomberg, Thermodynamic constraints on kinetic proofreading in biosynthetic pathways, Biophysical Journal 31, 333 (1980).
  * Murugan  _et al._ [2012] A. Murugan, D. A. Huse, and S. Leibler, Speed, dissipation, and error in kinetic proofreading, Proceedings of the National Academy of Sciences 109, 12034 (2012).
  * Wong  _et al._ [2018] F. Wong, A. Amir, and J. Gunawardena, Energy-speed-accuracy relation in complex networks for biological discrimination, Phys. Rev. E 98, 012420 (2018).
  * Yu  _et al._ [2020] Q. Yu, J. D. Mallory, A. B. Kolomeisky, J. Ling, and O. A. Igoshin, Trade-offs between speed, accuracy, and dissipation in trnaile aminoacylation, The Journal of Physical Chemistry Letters 11, 4001 (2020).
  * Hartich  _et al._ [2015] D. Hartich, A. C. Barato, and U. Seifert, Nonequilibrium sensing and its analogy to kinetic proofreading, New Journal of Physics 17, 055026 (2015).
  * Rao and Peliti [2015] R. Rao and L. Peliti, Thermodynamics of accuracy in kinetic proofreading: dissipation and efficiency trade-offs, Journal of Statistical Mechanics: Theory and Experiment 2015, P06001 (2015).
  * Pigolotti and Sartori [2016] S. Pigolotti and P. Sartori, Protocols for copying and proofreading in template-assisted polymerization, Journal of Statistical Physics 162, 1167 (2016).
  * Sartori and Pigolotti [2015] P. Sartori and S. Pigolotti, Thermodynamics of error correction, Phys. Rev. X 5, 041039 (2015).
  * Seifert [2018] U. Seifert, Stochastic thermodynamics: From principles to the cost of precision, Physica A: Statistical Mechanics and its Applications 504, 176 (2018).
  * Chiuchiú  _et al._ [2019] D. Chiuchiú, Y. Tu, and S. Pigolotti, Error-speed correlations in biopolymer synthesis, Phys. Rev. Lett. 123, 038101 (2019).
  * Banerjee  _et al._ [2017] K. Banerjee, A. B. Kolomeisky, and O. A. Igoshin, Elucidating interplay of speed and accuracy in biological error correction, Proceedings of the National Academy of Sciences 114, 5183 (2017).
  * Mallory  _et al._ [2019] J. D. Mallory, A. B. Kolomeisky, and O. A. Igoshin, Trade-offs between error, speed, noise, and energy dissipation in biological processes with proofreading, The Journal of Physical Chemistry B 123, 4718 (2019).
  * Yu  _et al._ [2022] Q. Yu, A. B. Kolomeisky, and O. A. Igoshin, The energy cost and optimal design of networks for biological discrimination, Journal of The Royal Society Interface 19, 20210883 (2022).
  * Miettinen [2012] K. Miettinen,  _Nonlinear Multiobjective Optimization_ , International Series in Operations Research & Management Science (Springer US, 2012).
  * Šponer  _et al._ [1996] J. Šponer, J. Leszczynski, and P. Hobza, Structures and energies of hydrogen-bonded dna base pairs. a nonempirical study with inclusion of electron correlation, The Journal of Physical Chemistry 100, 1965 (1996).
  * Deb [2001] K. Deb,  _Multi-Objective Optimization using Evolutionary Algorithms_ , Wiley Interscience Series in Systems and Optimization (Wiley, 2001).
  * for MATLAB 4.7.0. [13589] M. C. T. for MATLAB 4.7.0.13589, (Advanpix LLC., Yokohama, Japan).


# The Einstein ring GAL-CLUS-022058s: a Lensed Ultrabright Submillimeter
Galaxy at z=1.4796

A. Díaz-Sánchez  Departamento Física Aplicada, Universidad Politécnica de
Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain
H. Dannerbauer Instituto de Astrofísica de Canarias, Vía Láctea, La Laguna
38200, Spain
Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco
Sánchez, La Laguna, 38200, Spain  N. Sulzenauer Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
S. Iglesias-Groth Instituto de Astrofísica de Canarias, Vía Láctea, La Laguna
38200, Spain
Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco
Sánchez, La Laguna, 38200, Spain  R. Rebolo Instituto de Astrofísica de
Canarias, Vía Láctea, La Laguna 38200, Spain
Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco
Sánchez, La Laguna, 38200, Spain  Consejo Superior de Investigaciones
Cientícas, Spain

(Received 2021 February 19; Revised 2021 June 25; Accepted 2021 June 27)

###### Abstract

We report an ultra-bright lensed submillimeter galaxy at $z_{spec}=1.4796$,
identified as a result of a full-sky cross-correlation of the AllWISE and
Planck compact source catalogs aimed to search for bright submillimeter
galaxies at $z\sim 1.5-2.8$. APEX/LABOCA observations of the candidate galaxy
reveal a source with flux, S${}_{870\mu m}=54\pm 8$ mJy. The position of the
APEX source coincides with the position of the AllWISE mid-IR source, and with
the Einstein ring GAL-CLUS-022058s, observed with the HST. Archival VLT/FORS
observations reveal the redshift of this Einstein ring, $z_{spec}=1.4796$, and
detection of the CO(5-4) line at $z_{spec}=1.4802$ with APEX/nFLASH230
confirms the redshift of the submillimeter emission. The lensed source appears
to be gravitationally magnified by a massive foreground galaxy cluster lens at
$z=0.36$. We use Lenstool to model the gravitational lensing, which is near to
a “fold arc” configuration for an elliptical mass distribution of the central
halo, where four images of the lensed galaxy are seen; the mean magnification
is $\mu_{\rm L}=18\pm 4$. We have determined an intrinsic rest-frame infrared
luminosity of $L_{\rm IR}\approx 10^{12}L_{\sun}$ and a likely star formation
rate of $\sim 70-170$ $M_{\sun}\ yr^{-1}$. The molecular gas mass is $M_{\rm
mol}\sim 2.6\times 10^{10}M_{\sun}$ and the gas fraction is $f=0.34\pm 0.07$.
We also obtain a stellar mass log$(M_{\ast}/M_{\sun})=10.7\pm 0.1$ and a
specific star formation rate log$(sSFR/Gyr^{-1})=0.15\pm 0.03$. This galaxy
lies on the so-called main sequence of star-forming galaxies at this redshift.

galaxies: active — galaxies: clusters: general — galaxies: evolution —
galaxies: starburst — gravitational lensing: strong — submillimeter: galaxies

2021 September 23

##  1 Introduction

Submillimeter-selected galaxies (see Casey et al., 2014, for a review) allow
the study of the key phases of the formation and evolution of the most massive
galaxies in the distant universe. The peak of the submillimeter galaxy (SMG)
redshift distribution is at $z\approx 1.7-2.8$ (Chapman et al., 2005)
coinciding with the highest cosmic star-formation rate density. Furthermore,
the contribution to the star-formation history of the most luminous, dusty
star-forming galaxies steeply increases during earlier time (Le Floc’h et al.,
2005). Forming stars in compact starbursts with extremely high rates (see e.g.
Greve et al., 2005), exceeding thousand solar masses per year, could explain
the rapid build-up of present-day giant ellipticals (see Ivison et al., 2013).
Surveys with the Herschel space observatory, such as H-ATLAS (Negrello et al.,
2010), HerMES (Mitchell-Wynne et al., 2012), and with the South Pole Telescope
(SPT; Vieira et al., 2013; Weiß et al., 2013) have provided an increasing
number of the most luminous of these dusty star-forming galaxies. But due to
their compact sizes, resolved observations of their star-forming Interstelar
medium (ISM) remained limited, even with contemporary mm/submm facilities (see
e.g. Dye et al. (2015)). In addition, high dust attenuation at optical and
near-infrared (NIR) bands impede the identification of SMGs (Dannerbauer et
al., 2002) and further analysis of gas lines in SMGs at these wavelengths.

Strong gravitational lensing, especially with galaxy clusters as deflectors,
enables efficient detection of SMGs (Negrello et al., 2010). The Herschel and
the Planck space mission have provided detections of some of the brightest
lensed SMGs at the whole sky. In Sun et al. (2021) they carried out ALMA
observations of lensed SMGs selected by Herschel, where 29 sources with $z\sim
1.0-3.2$ in 20 fields of massive galaxy clusters have been detected, with a
median total IR luminosity $\mu_{L}L_{\rm IR}=10^{12.92\pm 0.07}L_{\odot}$.
Follow up observations of 24 strongly lensed star-forming galaxies between
$z\sim 1.0-3.5$ from the Planck satellite have been presented in Harrington et
al. (2021) with $\mu_{L}L_{\rm IR}\sim 10^{13-14.6}L_{\odot}$. The Cosmic
Eyebrow (Díaz-Sánchez et al., 2017) is one of the brightest SMGs in the sky,
it is a lensed dusty star-forming high redshift galaxy at the peak epoch of
the star formation in the Universe, at $z=2.04$. Observations of this galaxy
with IRAM/NOEMA at the redshifted CO(3-2) line led to the brightest detection
of this line in a SMG confirming its submillimeter nature and extreme
molecular gas properties (Dannerbauer et al., 2019). The source magnifications
permit resolved investigations of giant molecular clouds at a scale of
$\sim$100 pc. Together with high-precision lensing analysis, observations of
the molecular media of lensed galaxies allow to test models on the process of
star-formation in the early universe (Cava et al., 2018; Dessauges-Zavadsky et
al., 2019). Evidence of significant star-forming clumps down to a spatial
scale of 80 pc has been found in the Cosmic Snake at z=1.04 (Dessauges-
Zavadsky et al., 2019) but it is lacking in the Cosmic Eyelash at z= 2.3
(Ivison et al., 2020). Although non-lensed SMGs with high emission in the
submillimeter are obscured in the NIR and optical bands, lensing magnification
also helps to study the rest-frame UV and optical lines (see e.g. Acharyya et
al., 2019).

Based on data from the Planck and WISE space missions, we report here the
discovery of WISE J022057.56-383311.4, a magnified SMG at $z=1.48$. VLT
archival data have allowed us to obtain the redshift, and follow up
observations with APEX have confirmed the optical redshift and the strong
submillimiter fluxes observed with Planck coming from the molecular gas. The
paper is structured as follows: we first summarize our method for searching
lensed SMGs in section 2. In section 3, we present the observations and data
reduction for APEX observations, archival HST observations, archival VLT
spectroscopy, archival GEMINI spectroscopy, archival radio sources, and
archival Chandra observations. We analyze the data in section 4, we first
study the galaxy cluster properties, then we model the gravitational lensing,
after that the spectral energy density (SED) is presented and some properties
of the SMG are derived, and finally the CO lines are studied. In section 5, we
discuss the physical properties of the SMG, and in section 6 some conclusions
are presented. We adopt a flat $\Lambda$CDM cosmology from Planck
Collaboration et al. (2016) with $H_{0}=68$ km s-1 Mpc-1, $\Omega_{m}=0.31$
and $\Omega_{\Lambda}=1-\Omega_{m}$.

##  2 The search

As described in Díaz-Sánchez et al. (2017), we carried out a cross-matching
between the AllWISE111http://wise2.ipac.caltech.edu/docs/release/allwise/ and
Planck222http://pla.esac.esa.int/pla/#home full-sky compact source catalogs in
order to find bright lensed SMGs in the full-sky. First we built a full-sky
selection of galaxies verifying the color criteria $W1-W2>0.8$, $W2-W3<2.4$,
$W3-W4>3.5$, taken from the Cosmic Eyelash referenced SED (Swinbank et al.,
2010), and detection at $S/N>5$ in the four bands of the AllWISE catalog, we
limited the sample to galactic latitudes $|b|\geq 20^{o}$. Only SMGs at the
peak epoch of the star formation in the Universe ($z=1.5-2.8$) are expected to
fulfill these criteria (Iglesias-Groth et al., 2017). Then, we requested the
AllWISE selected objects to have a Planck source detected within 5 arcmin with
submillimeter flux ratios consistent with those expected for SMGs at
$z=1.5-3$. We removed the sources which clearly showed contamination by
galactic dust emission. This led to the identification of two sources WISE
J132934.18+224327.3, nicknamed Cosmic Eyebrow, and WISE J022057.56-383311.4.
This second source is the most likely counterpart of the submillimeter Planck
source PCCS2 857 G249.95-68.09 and a very promising candidate as an ultra-
bright SMG.

A subsequent search in the HST333http://archive.stsci.edu/ archive provided
images of a lensed galaxy at the same region where WISE J022057.56-383311.4 is
located. This source is within 30 arcsec from the centroid of the Planck
source and shows a spectacular strongly gravitational lensed galaxy with arc
structures larger than 25 arcsec, the contours in the WISE image are elongated
along the main arc structure. A massive galaxy cluster at $z=0.36$ is lensing
this galaxy. This galaxy cluster and the lensed SMG have been observed
previously with HST, VLT/FORS2, GEMINI-S and Chandra.

##  3 Observations and data reduction

Figure 1: (a) APEX/LABOCA signal map at 870 $\mu$m (345 GHz), we show a field
of $5\times 5$ arcmin2, north is up; east is left. Contour levels are
$3\sigma$, $3.5\sigma$, $4\sigma$, $4.5\sigma$, $4.8\sigma$, and $4.9\sigma$.
The beam size is also shown bottom-left in white. (b) APEX/nFLASH230 spectrum
centered at the redshifted CO(5-4) emission line frequencies of the lensed
source, obtaining a spectroscopic redshift $z=1.4802$. On the vertical axis,
two representations of flux density are shown. Equally, the horizontal axis
provides information on the channel frequency and radial velocity in the rest-
frame of the source. Additionally, careful noise reduction is performed via
the rms selection function and the noise content per binned channel is shown
as a gray fringe along the velocity axis. S/N in units of confidence level
$\sigma$ are also plotted across the radial velocity axis. Systemics
uncertainty from the antenna noise, at a constant level among a single
observing session, is shown in a box as $3\sigma$ error bars at the right
border, referring to the 90 % quantile of the line’s surface density value.
The data were binned in frequency to 60 km s-1 per channel. (c) APEX/SEPIA180
spectrum centered at the redshifted CO(4-3) emission line frequencies from the
lensed source, at redshift $z=1.4802$. The data were further binned in
frequency, to 180 km s-1 per channel, to try to increase the contrast

(same notation as for the central panel).

###  3.1 APEX Observations

We have observed this galaxy with APEX in order to confirm the origin of the
strong submillimeter fluxes and the spectroscopic redshift of this emission.
For that we have carried out photometric and spectroscopic observations with
APEX. We now describe these observations and the data reduction.

####  3.1.1 APEX photometry

Our observations were carried out between 18 and 19 April 2018, program ID
0101.B-0560(A) (PI Iglesias-Groth) using the bolometer LABOCA on APEX (Siringo
et al., 2009). We observed in mostly excellent weather conditions (PWV 0.5 mm
and 0.9 mm). The data were reduced using the CRUSH
444https://www.sigmyne.com/crush/ package, which is a free data reduction
software for specific astronomical imaging arrays. It is especially designed
for use with ground-based or air-borne (sub)millimeter wave cameras.
Individual maps were co-added (noise-weighted) and the final map was beam
smoothed, resulting in a spatial resolution with full width at half maximum
(FWHM) of 19.5 arcsec. The total source observing time of the used data for
this analysis is 8400 seconds and the average rms across the field is $8$ mJy
beam-1. A submillimeter source is detected with coordinates consistent with
those of WISE J022057.56-383311.4 which is partially responsible of the
observed Planck submillimeter fluxes. We measure a APEX/LABOCA flux of
$S_{870\mu m}=54\pm 8$ mJy. The differences with the expected Planck
measurements (see Fig. 6) are due to the flux of the Sunyaev-Zeldovich ($SZ$)
effect of a near galaxy cluster PSZ2 G249.80-68.11, at 4.5 arcmin north from
our source at $z=0.228$ (Cruddace et al., 2002). In fact, the position of the
source coincides with the arc structures of the strong lensed galaxy seen in
the HST images as can be seen in Figure 2(a). In Figure 1(a) we see the flux
detection at S/N$\sim 6.6$.

Figure 2: (a) $1\times 1$ arcmin2 HST/WFC3 RGB-image of GAL-CLUS-022058s with
the filters, blue F555W, green F775W and red F160W. North is up; east is left.
Yellow contours are from the APEX/LABOCA, contour levels are for $3\sigma$,
$3.5\sigma$, $4\sigma$, $4.5\sigma$, $4.8\sigma$, and $4.9\sigma$. (b) RGB-
image of GAL-CLUS-022058, which shows 3 slits from the VLT/FORS MOS (marked in
magenta with their measured redshifts). The long slit (marked in blue) is from
GEMINI GMOS observations. The obtained redshifts from the long slit are shown
beside the galaxies in green.

####  3.1.2 APEX spectroscopy

In order to obtain the redshift of the submillimeter emission we observed the
submillimeter source with APEX program ID E-0105.20N6 (PI Díaz-Sánchez). We
first observed with APEX/SEPIA180 (Belitsky et al., 2018) the CO(4-3) line at
186 GHz on 19 October 2020 with an exposure time of 47.2 minutes and PWV 2.1
mm. The CO(5-4) line was observed with APEX/nFLASH230555http://www.apex-
telescope.org/ns/nflash/ at 232 GHz the days 23 October, 7 and 8 December 2020
with a total exposure time of 57.7 minutes and PWV 2.1-2.6 mm. The total time
spent in this program was 5 hours including overheads.

Data reduction for the APEX spectra was performed with the software
GILDAS666https://www.iram.fr/IRAMFR/GILDAS with software package CLASS (Pety,
2005). After collecting all daily observing scans, possible line signals are
masked and polynomial baseline correction of order $n=5$ within an observed-
frame frequency window of $\Delta\nu\approx 1.5$ GHz centered on the signals
was performed. Close to the CO(5-4) emission line, a spurious instrumental
feature at $\nu_{obs}=232.1$ GHz is present in some of the scans. Thus
contaminated channels were blanked. In order to suppress excess noise from
inhomogeneous weather conditions across the daily observations, we ran a
simple CLASS script rmse_selection_function that was specifically developed
for blind spectroscopic line search by our group and is publicly available on
GitHub777https://github.com/NiSZR/rmse_selection_function. It works by
optimally selecting individual scans according to the ranked baseline rms
noise contribution in order to maximize the signal-to-noise content per
coadded spectrum. Selected scans are averaged and the frequency channels are
re-binned to source rest-frame 60 km s-1 to compromise between sampling of the
line profile and S/N-level for the full line detection. The flux conversion
factor for APEX/SEPIA180 (nFLASH230) is 38.4$\pm$2.8 (35.0$\pm$3.0) Jy/K. The
total uncertainty of the flux calibration is assumed to be less than 10 %.
Final flux calculations and visualization (see Figure 1(b,c)) were performed
in Python with modules NumPy (Virtanen et al., 2019) and astropy (Astropy
Collaboration et al., 2018).

We have detected the CO(5-4) line emission with $\mathrm{S/N}=6$ obtaining a
redshift of $z=1.4802$ with velocity-integrated line intensity of
$I_{CO(5-4)}=21\pm 4$ Jy km s-1 (Figure 1(b)). The FWHM line width for CO(5-4)
is $235\pm 45$ km s-1. We also have a tentative detection of CO(4-3) line with
upper limit velocity-integrated line intensity of $I_{CO(4-3)}\approx 44$ Jy
km s-1 at $3\sigma$ (Figure 1(c)). At this redshift, the flux for line CO(4-3)
is shifted towards the edge of the 2 mm atmospheric window, causing strong,
asymmetric noise contamination across the bandwidth.

###  3.2 Archival HST Optical and NIR imaging

We have identified the arc structure near to WISE J022057.56-383311.4 in
archival HST images taken in 2015, proposal ID 13756 (PI Saurabh Jha), aimed
at obtaining HST/WFC3 optical and infrared imaging of a large and nearly-
complete Einstein ring. The optical images were taken using HST/WFC3/UVIS with
filters F555W and F775W, with two exposures in each filter and a total
exposure time of 1227 seconds and 1270 seconds respectively in 2015-05-15. The
NIR image was taken with HST/WFC3/IR in filter F160W, with four exposures and
a total exposure time of 2812 seconds in 2015-05-16. We have taken the raw
data and processed them with the WFC3 pipeline (Gennaro, 2018) and IRAF
package (Figure 2). There is not yet any scientific publication related this
archival data.

###  3.3 Archival VLT Spectroscopy

Figure 3: Top: The spectra for the two MOS slits on the arc structure, taken
with two grisms (GRIS_600I+25 & GRIS_600z+23), they are ordered as in the
bottom panel. The doublet of OII $\lambda$3727, $\lambda$3729 Å is detected in
the two slits and for the two grisms, yielding a spectroscopic redshift
$z=1.479$. Bottom: The spectra are shown as a function of the rest-frame
wavelength, the sum of all spectra is also presented. Blue lines are for the
double of OII.

VLT/FORS2 observations were taken in 2008, ESO archive program ID
081.A-0693(A) (PI Graham Smith). For this cluster they designed a multi-slit
mask containing 25 slits of 1.0 arcsec width to include as many of the
candidate multiple systems as possible, with two tilted slits in the long arc
structures (Figure 2(b)). Then they carried out 3600 seconds of exposure time
with grism GRIS_600I+25, and 2700 seconds of exposure using grism name
GRIS_600z+23.

We have processed and calibrated these archival data using IRAF package and
our own specific code written in C++. Fig. 3 shows the final results for the
two slits and two grism in the spectral region where the doublet of OII
$\lambda$3727, $\lambda$3729 Å is detected. We have detected this doublet in
the two slits on the arc structure and for the two grisms. From this doublet
we obtain a spectroscopic redshift $z=1.479$ for the arc structures, which is
consistent with the spectroscopic redshift obtained from APEX observations.
For Slit 1 and 2 the doublet is detected in the whole region with significant
flux, it is impossible distinguish if it comes from the bluest or the reddest
clumps on the arc structures due to their separation (less than 1 arcsec).
Considering the sum of the flux from both slits and grisms, we have $S/N>10$
for the doublet of OII with a total exposure time of 1.75 hours, no more lines
are detected in the observed range $\lambda_{obs}$=6700-10600 Å (or
$\lambda_{rest}$=2700-4270 Å). We have measured the FWHM line width for the
OII emission line (corrected from the instrumental profile) and we obtained
$120\pm 20$ km s-1 for slit 1 and $145\pm 20$ km s-1 for slit 2.

We have also measured redshifts for 19 more galaxies on mask slits. We show
the redshift of these galaxies in Table 1. For galaxies with emission lines we
obtained the redshift from lines H$\alpha$, NII $\lambda$6585Å, and SII
$\lambda$6718, $\lambda$6733 Å, and for the galaxies without emission lines we
used mainly the absorption lines Na $\lambda$5896 Å, and Mg $\lambda$5177Å;
the H$\beta$ line was used when available. For a low redshift galaxy at
$z=0.00009$, it was obtained from the triplet of CaII $\lambda$8498,
$\lambda$8542, $\lambda$8662 Å, absorption lines. A redshift of a galaxy of
the cluster is found in the literature (Guzzo et al., 2009).

Figure 4: (a) $1\times 1$ arcmin2 X-ray image smoothed with a 2 arcsec FWHM
Gaussian. North is up; east is left. The first line of the contours is at
5$\sigma$ increasing lineally to 42$\sigma$, 14 levels with 2.6$\sigma$
separation among them. The radius of the circle is 25 arcsec and it is used to
made aperture photometry. (b) $1\times 1$ arcmin2 HST/WFC3 RGB-image of GAL-
CLUS-022058s, the X-ray contours are superimposed. North is up; east is left.
Table 1: Spectroscopic redshift of GAL-CLUS-022058s

RA | Dec |  $z_{spec}$aafootnotemark:  | Classbbfootnotemark:  | Telescope | Typeccfootnotemark:
---|---|---|---|---|---
02:20:57.27 | -38:33:12.2 | 1.47935 | SMG | VLT | E
02:20:58.36 | -38:33:10.7 | 1.47998 | SMG | VLT | E
02:20:59.50 | -38:32:52.0 | 0.20030 | G | VLT | A
02:20:59.44 | -38:32:27.4 | 0.36418 | CG | VLT | A
02:20:58.51 | -38:31:40.2 | 0.12377 | G | VLT | E
02:21:02.34 | -38:32:44.0 | 0.00009 | G | VLT | A
02:20:59.96 | -38:31:17.7 | 0.23106 | G | VLT | E
02:21:00.70 | -38:31:04.6 | 0.41227 | G | VLT | A
02:20:59.87 | -38:30:13.1 | 0.22670 | G | VLT | E
02:21:07.26 | -38:32:55.3 | 0.17080 | G | VLT | E
02:21:02.93 | -38:30:33.8 | 0.16500 | G | VLT | A
02:21:10.23 | -38:33:00.3 | 0.35784 | CG | VLT | E
02:21:06.18 | -38:30:56.6 | 0.44211 | G | VLT | A
02:21:08.99 | -38:31:31.5 | 0.36089 | CG | VLT | A
02:21:09.05 | -38:30:59.5 | 0.23022 | G | VLT | A
02:21:11.46 | -38:31:22.0 | 0.17824 | G | VLT | E
02:21:12.54 | -38:31:27.8 | 0.23761 | G | VLT | E
02:20:49.00 | -38:31:04.3 | 0.35851 | CG | VLT | A
02:20:52.37 | -38:33:09.6 | 0.36473 | CG | VLT | A
02:20:54.46 | -38:35:23.8 | 0.36066 | CG | VLT | A
02:20:49.60 | -38:34:33.7 | 0.35551 | CG | VLT | A
02:20:56.88 | -38:33:06.0 | 0.35732 | CG | GEMINI-S | A
02:20:57.73 | -38:33:03.1 | 0.35967 | BCG | GEMINI-S | E
02:20:58.71 | -38:33:00.7 | 0.36574 | CG | GEMINI-S | A
02:20:59.36 | -38:32:58.4 | 4.30949 | G | GEMINI-S | E
02:20:53.79 | -38:27:38.1 | 0.36057 | CG | Literature | E

Note. — aSpectroscopic redshift. bSMG=Submillimeter Galaxy, CG= Cluster
Galaxy, BCG=Brightest Cluster Galaxy, G=Galaxy. cE=Emission, A=Absorption.

###  3.4 Archival Gemini Spectroscopy

Spectroscopic follow-up of the cluster of galaxies, using GMOS long slit with
GEMINI-South in 2008, have been found in the GEMINI archive, archival program
ID GS-2008A-Q-5 (PI Graham Smith). The observation were taken with GMOS-S,
long slit of 1 arcsec width, grating B600+G5323, and 2400 seconds of exposure
time. The long slit used in the observations can be seen in Figure 2(b). We
have processed these data and we have obtained the redshift of three galaxies
of the cluster including the brightest cluster galaxy (BCG). The redshifts are
shown in Table 1. Moreover, a galaxy with the redshift $z=4.309$ has also been
detected.

For the BCG we obtain the redshift from the emission lines [OII]
$\lambda$3726, $\lambda$3729 Å, [OIII] $\lambda$5008 Å, and H$\beta$, and the
absorption lines K, H, and G. For the other two bright galaxies in the long
slit we determine the redshift from the absorption lines K, H, and G. For the
galaxy with the redshift $z=4.309$ the $Ly\alpha$ emission line is detected
with $S/N\sim 15$ in 20 minutes of exposure time for two different set-ups of
the instrument, and it is also detected in the HST/F160W band with AB-
magnitude of 25.20, which is magnified a factor $\mu_{L}\sim 5$ (calculated
with our lensing model in section 4.2).

###  3.5 Archival X-Ray Observations

Archival Chandra observations were taken in 2008, archival program ID 9411 (PI
Graham Smith), they observed a massive cluster of galaxies RXCJ0220.9-3829 at
redshift $z=0.2287$ near of our galaxy cluster, at $\sim 4.5$ arcmin north. We
can see in these observations X-ray emission coming from the galaxy cluster
Fig. 4. This X-ray emission have been mistaken with a point like source in
Wang et al. (2016) with a flux of $F=5.79\times 10^{-14}$ mW m-2 but in the
catalog given in Evans & Evans (2017) there is a warning saying this is an
extended source but they give the flux with an aperture of 2.4 arcsec which is
$F=4.01\times 10^{-14}$ mW m-2. We have taken the processed image in the
energy band 0.5-7.0 keV, from the Chandra X-ray
Center888https://cda.harvard.edu/chaser/. In Fig. 4(a) we show the smoothed
image with a 2 arcsec FWHM Gaussian. We have made a 25 arcsec aperture radius
photometry (Díaz-Sánchez et al., 2007), circle in green color in Fig. 4(a),
and obtained a flux of $F=(3.6\pm 0.2)\times 10^{-13}$ mW m-2. We show in Fig.
4(b) the X-ray contours superimposed to the RGB-image.

##  4 Analysis

###  4.1 Galaxy Cluster

The lensed galaxy is gravitationally magnified by a massive foreground galaxy
cluster at redshift $z=0.36$. We have obtained the spectroscopic redshift of
11 galaxies of the cluster. We have estimated the velocity dispersion of the
cluster along the line-of-sight from the 11 galaxy members of the cluster with
available spectroscopic redshifts (Table 1), which is $\sigma_{1\rm D}=684$ km
s-1. In order to estimate the size of the cluster we calculate the radius
$R_{200}$ , which approximates the virial radius, from $\sigma_{1\rm D}$
(Díaz-Sánchez et al., 2007) and we obtain $R_{200}=1.39$ Mpc. We use the
relation given in reference Saro et al. (2013) between the mass $M_{\rm cl}$
and $\sigma_{1\rm D}$ to calculate the cluster mass, which gives $M_{\rm
cl}=(3.3\pm 0.5)\times 10^{14}M_{\sun}$. Assuming that the X-ray emission
comes from the cluster of galaxies its X-ray luminosity is $L_{X}=(3.5\pm
0.2)\times 10^{43}$ erg s-1, which is compatible with its $M_{200}$ (Eckmiller
et al., 2011; Capasso et al., 2020). We have measured the FWHM line width for
two emission lines of the BCG, for the emission line [OII] $\lambda\lambda
3726,3728$ we obtain $373.2\pm 20$ km s-1, for the emission line [OIII]
$\lambda 5008$ we obtain $379.5\pm 20$ km s-1 and for H$\beta$ $396.1\pm 30$
km s-1. These are narrow emission lines and we discard Type I AGN for the BCG.

Figure 5:  (a) $28\times 22$ arcsec2 RGB-image of GAL-CLUS-022058s with the
HST/WFC3 filters. North is up; east is left. We mark the four images of the
lensed galaxy we have used to calculate the lensing model (identifying clumps
in these images). This image is taken from the image release as picture of the
week by ESA/Hubble, where the authors named this object GAL-CLUS-022058s, also
nicknamed “Molten Ring” a . (b) $30\times 30$ arsec2 RGB-image of GAL-
CLUS-022058s. The clumps used to model the gravitational lens in the four
images of the galaxy are marked in yellow and there are four families ($a$,
$b$, $c$, and $d$). Critical lines obtained from the lensing model for
$z=1.48$ are drawn in blue.

a https://www.spacetelescope.org/images/potw2050a/

###  4.2 Lens modeling

We have used the public software
Lenstool999https://projets.lam.fr/projects/lenstool/wiki (Kneib et al., 1993;
Jullo et al., 2007) to derive the lensing model. This software assumes a
parametric model for the distribution of dark matter and performs a mass
reconstruction of the foreground cluster. This model is constrained using the
location of the multiple images identified in the cluster. We used a simple
model with a single cluster-scale mass component, as well as individual
galaxy-scale mass components centered on each cluster member. For each
component, we used a dual pseudo-isothermal elliptical mass distribution
(dPIE, also known as a truncated PIEMD (Limousin et al., 2005)). Using
SExtractor (Bertin & Arnouts, 1996) and visual inspection of the HST/WFPC3
images, we have chosen multiple clumps with similar colors on the arc
structure as family members of multiply lensed background components of the
lensed galaxy at $z=1.479$. In selecting these image families we are supported
by a preliminary lens model, which is based on the primary arcs (Sharon et
al., 2020). The initial model is successful in predicting the approximate
locations of possible counter-images, in our case it has also predicted an
elliptical mass distribution for the central halo. Using the theoretical
results for this halo mass distribution we have identified four main images of
the elliptical lensed galaxy. The lensing configuration is near to the so
called “fold arc” (Kneib & Natarajan, 2011) (see Fig. 5(a)). The western image
is also distorted by the halo of a galaxy of the cluster.

Four families of four clumps were selected assuming redshift $z=1.479$ for all
of them (see Fig. 5(b) and Table 2). Family $a$ is in the center of the
elliptical galaxy and it is the reddest clump in the galaxy. Family $b$ is
near to the central part and an arm of the galaxy and it has green-yellow
color. And finally families $c$ and $d$ are on an arm of the galaxy and their
colors are the bluest ones. For the eastern image of the galaxy we have
measured the spectroscopic redshift with slit 2 (see section 3.3), as the flux
coming from the doublet of OII is in the whole region of the slit where the
clumps are present, we have the spectroscopic redshift for the families of
this image ($a.2$, $b.2$, $c.2$, and $d.2$). For the low central image of the
lensed galaxy we have measured the spectroscopic redshift with slit 1 (see
section 3.3) and we have only flux coming from the doublet of OII for the
central region of the galaxy ($a.3$) and for a blue arm ($c.3$).

We checked the model with the reconstruction of the arclets, and found good
agreement between the positions for image clumps and the reconstruction for
each component of the families (the rms between the original and the
resconstructed positions is 0.41 arcsec). The data are best fitted with a
cluster-scale potential of ellipticity $e=0.172$, position angle PA =102.08o,
normalization $\sigma_{\rm PIEMD}=796$ km s-1, core radius $r_{\rm
core}=35.82$ kpc and cut radio $r_{\rm cut}=1500$ kpc. The Einstein radius is
$\theta_{e}\approx 10$ arcsec. The enclosed mass within $R_{200}=1.39$ Mpc is
$M_{200}=3.8\times 10^{14}M_{\sun}$, in concordance with the $M_{200}$
obtained in section 4.1. We can see the critical lines for $z=1.48$ in Fig.
5(b), the magnification for each clump is given in Table 2, obtained from the
output of Lenstool software. The mean magnification in the center of the
galaxy (family $a$) is $\mu_{\rm L}=18\pm 4$ calculated from the four images
of family $a$.

Table 2: List of Lensing Constrains for the Lensing Model.

ID aafootnotemark:  | RA | Dec |  $\mu_{\rm L}$bbfootnotemark:  |  $z_{spec}$ccfootnotemark:
---|---|---|---|---
a.1 | 35.2404838 | -38.5494009 | $16\pm 4$ |
a.2 | 35.2432099 | -38.5529587 | $15\pm 3$ | 1.47998
a.3 | 35.2387213 | -38.5534089 | $25\pm 6$ | 1.47935
a.4 | 35.2365832 | -38.5512788 | $15\pm 3$ |
b.1 | 35.2404524 | -38.5491923 | $16\pm 4$ |
b.2 | 35.2436567 | -38.5525034 | $12\pm 2$ | 1.47998
b.3 | 35.2399095 | -38.5535043 | $17\pm 4$ |
b.4 | 35.2364258 | -38.5510711 | $12\pm 2$ |
c.1 | 35.2400105 | -38.5490930 | $17\pm 4$ |
c.2 | 35.2439491 | -38.5522620 | $10\pm 2$ | 1.47998
c.3 | 35.2393579 | -38.5532990 | $17\pm 4$ | 1.47935
c.4 | 35.2366492 | -38.5508668 | $14\pm 3$ |
d.1 | 35.2397692 | -38.5491489 | $18\pm 4$ |
d.2 | 35.2441465 | -38.5522513 | $9\pm 1$ | 1.47998
d.3 | 35.2384417 | -38.5529615 | $24\pm 6$ |
d.4 | 35.2368640 | -38.5509108 | $18\pm 3$ |

Note. — aThe letter refers to the name of the group of images coming from the
same region in the source galaxy. The number refers to the number of the image
seen in the image plane. bMagnification. cSpectroscopic redshift.

###  4.3 SED

Figure 6: Multi-wavelength spectral energy distribution of GAL-CLUS-022058s.
We have used the MAGPHYS code to fit the SED (da Cunha et al., 2008). The
difference between the Planck measurements and the SED fit are due to the flux
of the $SZ$ effect coming from a nearby galaxy cluster, PSZ2 G249.80-68.11, at
4.5 arcmin north from our source, and with $z=0.228$ (Cruddace et al., 2002).
The fluxes from the compact source (CS) plus the ones from the SZ source are
marked with red circles.

We plot the SED (Fig. 6) and list the photometry in Table 3. We have also
found data for our SMG in the GALEX NUV band (Bianchi et al., 2011). We have
calculated the upper limit for 60 and 100 $\mu$m from IRAS Sky Survey Atlas
101010http://irsa.ipac.caltech.edu/applications/IRAS/ISSA/. We have performed
the photometry for WISE J022057.56-383311.4 using the WISE channels 3 and 4
from the WISE images. As there is no emission from any galaxy of the cluster
for these channels, we have taken the whole measured flux in the region of the
arc structures. For channel 4 there is only 0.4 magnitudes of difference with
respect to the catalog magnitude and for channel 3 the difference is 1.2
magnitudes. For the other 2 channels we have taken the magnitudes given in the
catalog for this source on the arc structure, in these channels another source
on the BCG is detected and probably the flux for these bands is underestimated
due to resolution and deblending process. For the HST images we have performed
photometric measurements and taken the detected flux with S/N $>3$ on the arc
structures.

This SED is consistent with a source at redshift $z=1.479$ and brighter than
other Planck/Herschel SMGs found so far (Frye et al., 2019). There is a
cluster of galaxies with $SZ$ emission, PSZ2 G249.80-68.11, at 4.5 arcmin
north from our source with $z=0.228$ (Cruddace et al., 2002). The flux
detected by Planck is also coming from this source. We use the MAGPHYS code
that allows us to fit simultaneously the ultraviolet-to-FIR SED, so we can
constrain physical parameters (da Cunha et al., 2008). We assume that the
submillimeter flux is coming from a region near to the center of the lensed
galaxy and we take the mean amplification for family $a$ as the submillimeter
amplification, which is $\mu_{\rm L}=18$. The fit is shown in Figure 6. From
this fit we have the following physical parameters corrected for lensing
amplification and image multiplicity, stellar mass
log$(M_{\ast}/M_{\sun})=10.7\pm 0.1$, specific star formation rate (sSFR)
log$(sSFR/Gyr^{-1})=0.15\pm 0.03$, average V-band dust attenuation
$A_{V}=2.30\pm 0.01$, luminosity-averaged dust temperature $T_{\rm
dust}/K=36\pm 2$, mass-weighted age log$(age_{M}\ yr^{-1})=8.84\pm 0.02$, dust
luminosity log$(L_{\rm dust}/L_{\sun})=12.0\pm 0.1$, dust mass log$(M_{\rm
dust}/M_{\sun})=8.12\pm 0.09$, and star formation rate (SFR)
log$(SFR/M_{\sun}\ yr^{-1})=1.83\pm 0.02$. We also calculate the SFR from the
dust luminosity and the expression given in Kennicutt (1998) $SFR\sim 172$
$M_{\sun}\ yr^{-1}$, which assumes a Salpeter IMF, assuming a Chabrier IMF
would possibly be a factor of $\sim$ 1.8 lower (Casey et al., 2014).

The sSFR agrees with the so-called main sequence of star-forming galaxies with
a stellar mass of $10^{11}$ M☉ at this redshift (Sun et al., 2021), and if we
compare the sSFR with the dust-to-stellar mass ratio (also referred to as the
specific dust mass), this correlation also matches with previous high redshift
SMGs (Sun et al., 2021). The SFR dependence with the stellar mass is
consistent with the main sequence behavior (Caliendo et al., 2021; Whitaker et
al., 2014; Danielson et al., 2017).

Table 3: Ultraviolet to Submillimeter photometry of GAL-CLUS-022058s

Wavelength($\mu$m) | Flux (mJy)a,ba,bfootnotemark:  | Observatory/Instrument
---|---|---
0.23157 | 0.009 $\pm 0.003$ | $GALEX/NUV$
0.5308 | 0.02 $\pm 0.01$ | $HST/WFC3$
0.7647 | 0.03 $\pm 0.02$ | $HST/WFC3$
1.5369 | 0.42 $\pm 0.02$ | $HST/WFC3$
3.4 | 0.76 $\pm 0.02$ | $WISE$
4.6 | 0.89 $\pm 0.02$ | $WISE$
12 | 2.45 $\pm 0.08$ | $WISE$
22 | 9.9 $\pm 0.7$ | $WISE$
60 |  $<$ 200c | $IRAS$
100 |  $<$ 300c | $IRAS$
350 | 1400 $\pm 600$d | $Planck$
550 | 800 $\pm 300$d | $Planck$
850 |  $<$ 200c | $Planck$
870 | 54 $\pm 8$ | APEX/LABOCA
1400 |  $<$ 50c | $Planck$

Note. — aUncorrected for lensing amplification. bTotal arc flux. cUpper flux
limit. d $SZ$ \+ Compact Source fluxes.

###  4.4 CO lines

Figure 1(b,c) shows the CO line detections at the position of our SMG. We have
identified two lines CO(5-4) and CO(4-3). The binned $\Delta V=60$ km s-1
CO(5-4) line is detected with flux density $S_{CO(5-4)}=40\pm 6$ mJy,
velocity-integrated line intensity $I_{CO(5-4)}=21\pm 4$ Jy km s-1, FWHM line
width $235\pm 45$ km s-1 and redshift $z=1.4802$. The binned $\Delta V=180$ km
s-1 CO(4-3) tentative detection line would have $S_{CO(4-3)}\approx 93$ mJy
and $I_{CO(4-3)}\approx 44$ Jy km s-1, at the redshift $z=1.4802$. This last
line is located close to a strong atmospheric H2O line at 183 GHz. Accounting
for the low $S/N\sim 3$ of this line measurement, the inferred line intensity
should be considered as an upper-limit.

For the CO(5-4) line we derive a luminosity $\mu_{\rm
L}L^{\prime}_{CO(5-4)}=1.0\pm 0.2\times 10^{11}$ K km s-1 pc2, which is
uncorrected from lensing amplification and image multiplicity. If we compare
this luminosity with the also uncorrected IR luminosity, log$(\mu_{\rm
L}L_{\rm IR}/L_{\sun})=13.85$, it agrees with the expectation for high
redshift SMGs (Carilli & Walter, 2013). The luminosity of this line is high
and it can be compared with the ones listed in Harrington et al. (2021), 5
from the 24 SMGs listed in this reference have a lower CO(5-4) luminosity,
$\mu_{\rm L}L^{\prime}_{CO(5-4)}$.

Bothwell et al. (2013) calculated average flux ratios for a sample of 40 non-
lensed luminous SMGs at comparable redshift. Demanding similar subthermal
conditions as in this reference and small differential lensing effects, the
velocity-integrated line intensity for CO(1-0) line can be derived,

$I_{CO(1-0)}\approx\frac{I_{CO(5-4)}}{25\langle r_{2,1}\rangle\langle
r_{3,2}\rangle\langle r_{4,3}\rangle\langle r_{5,4}\rangle}=14.7{\rm\ Jy\ km\
s^{-1}}.$ (1)

We use the CO luminosity to molecular gas mass conversion factor $\alpha_{\rm
CO}=0.8$ $M_{\sun}$ pc-2 (K km s-1)-1, a commonly used value for merger-
induced star formation (Solomon & Vanden Bout, 2005). Including a factor of
1.36 to account for helium, and correcting for lensing amplification and image
multiplicity, we derive Mmol $\approx 2.6\times 10^{10}M_{\sun}$. This value
agrees with the main sequence relation of the star-forming galaxies for SFR-
Mmol (Sargent et al., 2014). For galaxies on the main sequence, the gas mass
conversion factor $\alpha_{\rm CO}$ is larger than the one used here, up to
$\sim 5$ larger (Magnelli et al., 2012). For the 24 Planck lensed SMGs listed
in Harrington et al. (2021) with similar gas mass, they found great
dispersion, $\alpha_{\rm CO}=0.8-4$. On the other hand, the conversion factor
given in equation 1 also has a great uncertainty. Nevertheless, the $M_{\rm
mol}$ value obtained here agrees with the one calculated from the dust
continuum in the next section, but it could be between Mmol
$\approx(2.6-13)\times 10^{10}M_{\sun}$.

##  5 Discussion

Cross-matching between the AllWISE and Planck full-sky compact source catalogs
we have identified a bright high redshift lensed SMG. It is lensed by a galaxy
cluster at $z=0.36$ with a mean amplification factor in the center of the
galaxy of $\mu_{\rm L}\approx 18$, four images of the lensed galaxy have been
identified. The redshifts derived from the submillimeter spectra and the
optical one agree with each other very well, indicating that the submillimeter
flux is coming from the arc structure at $z\approx 1.48$. The uncorrected IR
luminosity is log$(\mu_{\rm L}L_{\rm IR}/L_{\sun})=13.85$ and it is among the
brightest high redshift lensed galaxies found with Planck (Harrington et al.,
2021). The intrinsic luminosity is $L_{\rm IR}\approx 10^{12}$ $L_{\sun}$,
which gives a SFR of $\sim 70-170$ $M_{\sun}\ yr^{-1}$. We assume that the
submillimeter flux is coming from a region near to the center of the lensed
galaxy and we take the mean amplification for family $a$ as the submillimeter
amplification. The amplification depends on the position of the flux in the
image plane, see Table 2, if it was near to the critical lines then the
amplification would be much larger. We correct for lensing magnification and
image multiplicity, and with $T_{\rm dust}=36$ K and the observed 850 $\mu m$
dust continuum, we use Equation (16) in Scoville et al. (2016) to derive
M${}_{\rm mol}=1.7\pm 0.9\times 10^{10}M_{\sun}$. This value agrees with the
main sequence relation of the star-forming galaxies for SFR-Mmol (Sargent et
al., 2014) and with the value obtained from the CO(5-4) line calculated in
section 4.4. We have derived physical properties from the SED using MAGPHYS.
Comparing the stellar mass with the star formation rate, we find that this
galaxy lies on the main sequence of the star-forming galaxies at $z\sim 1.5$
(Caliendo et al., 2021; Whitaker et al., 2014; Danielson et al., 2017). The
sSFR also agrees with the main sequence of star-forming galaxies at this
redshift (Sun et al., 2021; Dudzevičiūtė et al., 2020). The comparison between
the sSFR and the dust-to-stellar mass ratio also matches with previous high
redshift SMGs (Sun et al., 2021). From the luminosity for the CO(5-4) line,
$\mu_{\rm L}L^{\prime}_{CO(5-4)}$, and the uncorrected IR luminosity $\mu_{\rm
L}L_{\rm IR}$, this galaxy compares well with the main sequence of star-
forming galaxies (Carilli & Walter, 2013). The gas fraction is $f=M_{\rm
mol}/(M_{\rm mol}+M_{\ast})=0.34\pm 0.07$, which agrees with the expectations
at this redshift for SMGs (Dudzevičiūtė et al., 2020; Bothwell et al., 2013).

We have measured different widths for the [OII] and C(5-4) lines, one is
detected in the optical and the other one in the submillimeter regime, so it
is possible that these emissions would come from different regions in the
galaxy, i.e. the optical line would come from the spiral arms and the
submillimeter line from the center of the galaxy near to an AGN. The
submillimeter emission would come from the interaction of dense ISM around a
central supermassive black hole, so it is wider than the optical one.

##  6 Conclusion

We have found an ultra-bright SMG at $z_{spec}=1.4796$ cross-matching the
AllWISE and Planck catalogs. Our APEX/LABOCA observations give the position of
the submillimeter flux which coincides with the position of a multiple lensed
galaxy, the Einstein ring GAL-CLUS-022058s, observed with HST. We have
obtained the spectroscopic redshift from archival VLT/FORS observations and
confirmed the redshift of the submillimeter emission with our APEX/nFLASH230
observations. The lensed source appears to be gravitationally magnified by a
massive foreground galaxy cluster lens at $z=0.36$, we have identified four
images of the lensed galaxy. We have used Lenstool to model the gravitational
lensing which is near to a “fold arc” configuration for an elliptical mass
distribution of the central halo; the mean magnification is $\mu_{\rm L}=18\pm
4$. We have determined an intrinsic rest-frame infrared luminosity of $L_{\rm
IR}\approx 10^{12}L_{\sun}$ and a likely SFR of $\sim 70-170$ $M_{\sun}\
yr^{-1}$. The molecular gas mass is $M_{\rm mol}\sim 2.6\times
10^{10}M_{\sun}$ and the gas fraction is $f=0.34\pm 0.07$. We also obtain a
stellar mass log$(M_{\ast}/M_{\sun})=10.7\pm 0.1$ and a specific SFR
log$(sSFR/Gyr^{-1})=0.15\pm 0.03$. This galaxy lies on the so-called main
sequence of the star-forming galaxies for several parameters. Finally we
emphasize that a large amount of data used here are from archive. Our new
observations have been carried out with APEX telescope.

This publication is based on data acquired with the Atacama Pathfinder
Experiment (APEX) under programmes ID 105.20N6 and ID 0101.B-0560(A). APEX is
a collaboration between the Max-Planck-Institut fur Radioastronomie, the
European Southern Observatory, and the Onsala Space Observatory. We thank
Carlos De Breuck, and the APEX team for their support. Based on observations
collected at the European Southern Observatory under ESO programme ID
081.A-0693(A). Based on observations made with the NASA/ESA Hubble Space
Telescope, and obtained from the Hubble Legacy Archive, which is a
collaboration between the Space Telescope Science Institute (STScI/NASA), the
Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian
Astronomy Data Centre (CADC/NRC/CSA). This publication makes use of data
products from the Wide-field Infrared Survey Explorer. Based on observations
obtained at the international Gemini Observatory program ID GS-2008A-Q-5. The
scientific results reported in this article are based in part on data obtained
from the Chandra Data Archive. This work has been partially funded by projects
“Participation in the NISP instrument and preparation for the scientific
exploitation of Euclid”, PID2019-110614GB-C22/AEI/10.13039/501100011033, &
PID2019-110614GB-C21 financed by the “Agencia Estatal de Investigación” (AEI-
MCINN). H.D. acknowledges financial support from the Spanish Ministry of
Science, Innovation and Universities (MICIU) under the 2014 Ramón y Cajal
program RYC-2014-15686 and under the AYA2017-84061-P, co-financed by FEDER
(European Regional, Development Funds), and in addition, from the Agencia
Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN)
under grant (La evolución de los cíumulos de galaxias desde el amanecer hasta
el mediodía cósmico) with reference
(PID2019-105776GB-I00/DOI:10.13039/501100011033). N.S. acknowledges a research
stay at the Instituto de Astrofísica de Canarias, funded by “Short-term grant
abroad” (KWA) programme of the University of Vienna.

## References

  * Acharyya et al. (2019) Acharyya, A., Kewley, L. J., Rigby, J. R., et al. 2019, MNRAS, 488, 5862, doi: 10.1093/mnras/stz1987
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Belitsky et al. (2018) Belitsky, V., Lapkin, I., Fredrixon, M., et al. 2018, A&A, 612, A23, doi: 10.1051/0004-6361/201731458
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bianchi et al. (2011) Bianchi, L., Herald, J., Efremova, B., et al. 2011, Ap&SS, 335, 161, doi: 10.1007/s10509-010-0581-x
  * Bothwell et al. (2013) Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047, doi: 10.1093/mnras/sts562
  * Caliendo et al. (2021) Caliendo, J. N., Whitaker, K. E., Akhshik, M., et al. 2021, ApJ, 910, L7, doi: 10.3847/2041-8213/abe132
  * Capasso et al. (2020) Capasso, R., Mohr, J. J., Saro, A., et al. 2020, MNRAS, 494, 2736, doi: 10.1093/mnras/staa888
  * Carilli & Walter (2013) Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105, doi: 10.1146/annurev-astro-082812-140953
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * Cava et al. (2018) Cava, A., Schaerer, D., Richard, J., et al. 2018, Nature Astronomy, 2, 76, doi: 10.1038/s41550-017-0295-x
  * Chapman et al. (2005) Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772, doi: 10.1086/428082
  * Cruddace et al. (2002) Cruddace, R., Voges, W., Böhringer, H., et al. 2002, ApJS, 140, 239, doi: 10.1086/324519
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595, doi: 10.1111/j.1365-2966.2008.13535.x
  * Danielson et al. (2017) Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2017, ApJ, 840, 78, doi: 10.3847/1538-4357/aa6caf
  * Dannerbauer et al. (2019) Dannerbauer, H., Harrington, K., Díaz-Sánchez, A., et al. 2019, AJ, 158, 34, doi: 10.3847/1538-3881/aaf50b
  * Dannerbauer et al. (2002) Dannerbauer, H., Lehnert, M. D., Lutz, D., et al. 2002, ApJ, 573, 473, doi: 10.1086/340762
  * Dessauges-Zavadsky et al. (2019) Dessauges-Zavadsky, M., Richard, J., Combes, F., et al. 2019, Nature Astronomy, 3, 1115, doi: 10.1038/s41550-019-0874-0
  * Díaz-Sánchez et al. (2017) Díaz-Sánchez, A., Iglesias-Groth, S., Rebolo, R., & Dannerbauer, H. 2017, ApJ, 843, L22, doi: 10.3847/2041-8213/aa79ef
  * Díaz-Sánchez et al. (2007) Díaz-Sánchez, A., Villo-Pérez, I., Pérez-Garrido, A., & Rebolo, R. 2007, MNRAS, 377, 516, doi: 10.1111/j.1365-2966.2007.11643.x
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828, doi: 10.1093/mnras/staa769
  * Dye et al. (2015) Dye, S., Furlanetto, C., Swinbank, A. M., et al. 2015, MNRAS, 452, 2258, doi: 10.1093/mnras/stv1442
  * Eckmiller et al. (2011) Eckmiller, H. J., Hudson, D. S., & Reiprich, T. H. 2011, A&A, 535, A105, doi: 10.1051/0004-6361/201116734
  * Evans & Evans (2017) Evans, J. D., & Evans, I. N. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 512, Astronomical Data Analysis Software and Systems XXV, ed. N. P. F. Lorente, K. Shortridge, & R. Wayth, 309
  * Frye et al. (2019) Frye, B. L., Pascale, M., Qin, Y., et al. 2019, ApJ, 871, 51, doi: 10.3847/1538-4357/aaeff7
  * Gennaro (2018) Gennaro, M. 2018, Baltimore: STScI, 4.0
  * Greve et al. (2005) Greve, T. R., Bertoldi, F., Smail, I., et al. 2005, MNRAS, 359, 1165, doi: 10.1111/j.1365-2966.2005.08979.x
  * Guzzo et al. (2009) Guzzo, L., Schuecker, P., Böhringer, H., et al. 2009, A&A, 499, 357, doi: 10.1051/0004-6361/200810838
  * Harrington et al. (2021) Harrington, K. C., Weiss, A., Yun, M. S., et al. 2021, ApJ, 908, 95, doi: 10.3847/1538-4357/abcc01
  * Iglesias-Groth et al. (2017) Iglesias-Groth, S., Díaz-Sánchez, A., Rebolo, R., & Dannerbauer, H. 2017, MNRAS, 467, 330, doi: 10.1093/mnras/stx041
  * Ivison et al. (2020) Ivison, R. J., Richard, J., Biggs, A. D., et al. 2020, MNRAS, 495, L1, doi: 10.1093/mnrasl/slaa046
  * Ivison et al. (2013) Ivison, R. J., Swinbank, A. M., Smail, I., et al. 2013, ApJ, 772, 137, doi: 10.1088/0004-637X/772/2/137
  * Jullo et al. (2007) Jullo, E., Kneib, J. P., Limousin, M., et al. 2007, New Journal of Physics, 9, 447, doi: 10.1088/1367-2630/9/12/447
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kneib et al. (1993) Kneib, J. P., Mellier, Y., Fort, B., & Mathez, G. 1993, A&A, 273, 367
  * Kneib & Natarajan (2011) Kneib, J.-P., & Natarajan, P. 2011, A&A Rev., 19, 47, doi: 10.1007/s00159-011-0047-3
  * Le Floc’h et al. (2005) Le Floc’h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169, doi: 10.1086/432789
  * Limousin et al. (2005) Limousin, M., Kneib, J.-P., & Natarajan, P. 2005, MNRAS, 356, 309, doi: 10.1111/j.1365-2966.2004.08449.x
  * Magnelli et al. (2012) Magnelli, B., Saintonge, A., Lutz, D., et al. 2012, A&A, 548, A22, doi: 10.1051/0004-6361/201220074
  * Mitchell-Wynne et al. (2012) Mitchell-Wynne, K., Cooray, A., Gong, Y., et al. 2012, ApJ, 753, 23, doi: 10.1088/0004-637X/753/1/23
  * Negrello et al. (2010) Negrello, M., Hopwood, R., De Zotti, G., et al. 2010, Science, 330, 800, doi: 10.1126/science.1193420
  * Pety (2005) Pety, J. 2005, in SF2A-2005: Semaine de l’Astrophysique Francaise, ed. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani, 721
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Sargent et al. (2014) Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19, doi: 10.1088/0004-637X/793/1/19
  * Saro et al. (2013) Saro, A., Mohr, J. J., Bazin, G., & Dolag, K. 2013, ApJ, 772, 47, doi: 10.1088/0004-637X/772/1/47
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83, doi: 10.3847/0004-637X/820/2/83
  * Sharon et al. (2020) Sharon, K., Bayliss, M. B., Dahle, H., et al. 2020, ApJS, 247, 12, doi: 10.3847/1538-4365/ab5f13
  * Siringo et al. (2009) Siringo, G., Kreysa, E., Kovács, A., et al. 2009, A&A, 497, 945, doi: 10.1051/0004-6361/200811454
  * Solomon & Vanden Bout (2005) Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677, doi: 10.1146/annurev.astro.43.051804.102221
  * Sun et al. (2021) Sun, F., Egami, E., Rawle, T. D., et al. 2021, ApJ, 908, 192, doi: 10.3847/1538-4357/abd6e4
  * Swinbank et al. (2010) Swinbank, A. M., Smail, I., Longmore, S., et al. 2010, Nature, 464, 733, doi: 10.1038/nature08880
  * Vieira et al. (2013) Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Nature, 495, 344, doi: 10.1038/nature12001
  * Virtanen et al. (2019) Virtanen, P., Gommers, R., Burovski, E., et al. 2019, scipy/scipy: SciPy 1.3.2, v1.3.2, Zenodo, doi: 10.5281/zenodo.3533894
  * Wang et al. (2016) Wang, S., Liu, J., Qiu, Y., et al. 2016, ApJS, 224, 40, doi: 10.3847/0067-0049/224/2/40
  * Weiß et al. (2013) Weiß, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88, doi: 10.1088/0004-637X/767/1/88
  * Whitaker et al. (2014) Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104, doi: 10.1088/0004-637X/795/2/104