# On the first eigenvalue of the generalized laplacian

Julián Fernández Bonder and Ariel M. Salort  Departamento de Matemática, FCEN
– Universidad de Buenos Aires, and Instituto de Cálculo – CONICET $0+\infty$
building, Ciudad Universitaria (1428), Buenos Aires, Argentina.
<EMAIL_ADDRESS><EMAIL_ADDRESS>page:http://mate.dm.uba.ar/$\sim$jfbonder, http://mate.dm.uba.ar/$\sim$asalort

###### Abstract.

In this work we investigate the energy of minimizers of Rayleigh-type
quotients of the form

$\frac{\int_{\Omega}A(|\nabla u|)\,dx}{\int_{\Omega}A(|u|)\,dx}.$

These minimizers are eigenfunctions of the generalized laplacian defined as
$\Delta_{a}u=\mathop{\text{div}}\left(a(|\nabla u|)\frac{\nabla u}{|\nabla
u|}\right)$ where $a(t)=A^{\prime}(t)$ and the Rayleigh quotient is comparable
to the associated eigenvalue. On the function $A$ we only assume that it is a
Young function but no $\Delta_{2}$ condition is imposed.

Since the problem is not homogeneous, the energy of minimizers is known to
strongly depend on the normalization parameter
$\alpha=\int_{\Omega}A(|u|)\,dx$. In this work we precisely analyze this
dependence and show differentiability of the energy with respect to $\alpha$
and, moreover, the limits as $\alpha\to 0$ and $\alpha\to\infty$ of the
Rayleigh quotient.

The nonlocal version of this problem is also analyzed.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Elementary estimates on the energy function
  4. 4 Continuity and differentiability of the energy function
  5. 5 Sharp asymptotic behavior when $\alpha\to 0$ and $\alpha\to\infty$
  6. 6 The behavior of $E(\alpha)/\alpha$ without the $\Delta_{2}-$condition
  7. 7 Nonlocal case

##  1\. Introduction

Let $A\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ be convex function superlinear at
0 and at $\infty$. Given an open and bounded set $\Omega\subset\mathbb{R}^{n}$
and a fixed normalization parameter $\alpha>0$, we consider the energy
function $E(\alpha)$ defined as

$E(\alpha)=\inf\left\\{\int_{\Omega}A(|\nabla u|)\,dx\colon u\in
C_{c}^{\infty}(\Omega),\int_{\Omega}A(|u|)\,dx=\alpha\right\\}.$ (1.1)

Without any further growth assumption on $A$, in [8, 9] it is proved that the
minimization problem $E(\alpha)$ is solvable, i.e., there exists a function
$u_{\alpha}$ for which the infimum in $E(\alpha)$ is achieved. Furthermore,
through a generalized version of the Lagrange’s Multiplier Theorem, in [9]
(see also [22]) it is proved that $u_{\alpha}$ satisfies in the weak sense the
equation

$\displaystyle\begin{cases}-\Delta_{a}u=\lambda(\alpha)\frac{a(|u|)}{|u|}u&\quad\text{
in }\Omega\\\ u=0&\quad\text{ in }\partial\Omega,\end{cases}$ (1.2)

where the eigenvalue $\lambda(\alpha)$ is the Lagrange multiplier,
$a(t)=A^{\prime}(t)$ and
$\Delta_{a}u:=\mathop{\text{div}}\left(\frac{a(|\nabla u|)}{|\nabla u|}\nabla
u\right)$ is the so-called _generalized Laplacian_ or the _$a-$ Laplacian
operator_.

Existence of solutions to (1.1) and eigenpairs to (1.2) assuming diverse
growth behaviors on $A$ were widely studied, see for instance [8, 9, 17, 18,
20, 22]. These problems, typically are non homogeneous and in general, the
eigenvalue $\lambda(\alpha)$ does not coincide with the ratio
$E(\alpha)/\alpha$. One known exception arises when the problem is
homogeneous, that is, when $A(t)=t^{p}$ for some $p>1$ and hence (1.2) turns
out to be the well-known eigenvalue problem for the Dirichlet $p-$Laplace
operator. In this case the first eigenvalue of such a problem admits the
variational characterization

$\lambda(1)=\lambda(\alpha)=\frac{E(\alpha)}{\alpha}=\inf_{0\neq u\in
C_{c}^{\infty}(\Omega)}\frac{\int_{\Omega}|\nabla
u|^{p}dx}{\int_{\Omega}|u|^{p}\,dx}$

and it is independent of the normalization parameter $\alpha>0$. For further
information, see for instance [12, 13, 14].

Although, in general, eigenvalues and energy may not coincide, when the Young
function $A$ satisfies the $\Delta_{2}-$condition (see Section 2 for the
definition), it can be readily proved that the following comparison estimate

$c\lambda(\alpha)\leq\frac{E(\alpha)}{\alpha}\leq C\lambda(\alpha),$ (1.3)

holds for some constants $c,C>0$ depending only on $A$. Moreover, in Section 3
the following elementary inequality is proved to hold:

$\min\left\\{\alpha^{p-1},\alpha^{\frac{1}{p}-1}\right\\}E(1)\leq\frac{E(\alpha)}{\alpha}\leq\max\left\\{\alpha^{p-1},\alpha^{\frac{1}{p}-1}\right\\}E(1),$
(1.4)

where $p>1$ is a suitable constant depending only on $A$. However, while this
estimate is straightforward, it lacks the depth necessary to provide
meaningful insights into the behavior of $E(\alpha)$ as $\alpha\downarrow 0$
and $\alpha\uparrow\infty$. This lack of information serves as the primary
motivation for this paper.

Our first results concern the regularity of the curve $E(\alpha)$. In Theorem
4.5, without any additional assumptions regarding the behavior of $A$, we
prove that the energy $E(\alpha)$ is a Lipschitz continuous function of the
parameter $\alpha>0$. Moreover, assuming a technical but fairly natural
condition on $A$ we are able to obtain its differentiability. More precisely,
assuming that the map $t\mapsto ta(t)$ is convex, in Theorem 4.6 we prove that
$E(\alpha)$ is differentiable for any $\alpha>0$ and the following formulas
hold:

$\displaystyle\frac{dE}{d\alpha}(\alpha)=\lambda(\alpha),\qquad\frac{d}{d\alpha}\left(\frac{E(\alpha)}{\alpha}\right)=\frac{1}{\alpha}\left(\lambda(\alpha)-\frac{E(\alpha)}{\alpha}\right),$

In these expressions, the eigenvalue $\lambda(\alpha)$ and the energy quotient
$E(\alpha)/\alpha$ coincide if and only if both quantities are constants
independent of $\alpha$. This is precisely the scenario where $A$ is
homogeneous.

In Section 5 we analyze the behavior of the ratio $E(\alpha)/\alpha$ as
$\alpha\downarrow 0$ and $\alpha\uparrow\infty$. Our main result establishes
that the behavior of that quotient can be described in terms of the the
Matuszewska-Orlicz functions associated to the Young function $A$, which are
defined as follows

$M_{0}(t)=\limsup_{\tau\to 0}\frac{A(\tau t)}{A(\tau)}\qquad\text{and}\qquad
M_{\infty}(t)=\limsup_{\tau\to\infty}\frac{A(\tau t)}{A(\tau)},$

(see Section 2 for further information and properties of these function). We
highlight that the behavior of the energy quotient only takes into account the
local behavior of the Young function $A$ near $0$ and near $\infty$. More
precisely, in Theorem 5.1 we prove that if $A$ satisfies the $\Delta_{2}$
condition near $i=0,\infty$, then

$\limsup_{\alpha\to
i}\frac{E(\alpha)}{\alpha}=\inf_{\tau>0}\frac{E_{i}(\tau)}{\tau}$

where

$E_{i}(\tau)=\inf\left\\{\int_{\Omega}M_{i}(|\nabla u|)\,dx\colon u\in
C_{c}^{\infty}(\Omega),\,\int_{\Omega}M_{i}(|u|)\,dx=\tau\right\\}.$

This result can be refined. We prove that under a fairly general additional
assumption on $A$, there exist numbers $p_{i}\in[1,\infty)$ for which
$M_{i}(t)=t^{p_{i}}$, and these result in

$\lim_{\alpha\to i}\frac{E(\alpha)}{\alpha}=\lambda_{p_{i}},$

where $\lambda_{p}$ denotes the first eigenvalue of the $p-$Laplacian with
Dirichlet boundary conditions in $\Omega$.

Although in the previous result the assumption of the $\Delta_{2}$ condition
on $A$ seems restrictive, in fact this is somehow the interesting case of
study since when this condition is removed, in broad terms, the Matuszewska-
Orlicz functions associated to the Young function $A$ becomes trivial:

$M_{i}(t)=\begin{cases}0&\text{if }t<1\\\ 1&\text{if }t=1\\\ \infty&\text{if
}t>1\\\ \end{cases}$

for $i=0,\infty$. In this case, under suitable assumptions on the domain, the
limit of $E(\alpha)/\alpha$ can be characterized. More precisely, if
$\Omega\subset\mathbb{R}^{n}$ is such that the inner radius is larger than 1,
in Theorem 6.1 we obtain that

$\lim_{\alpha\to i}\frac{E(\alpha)}{\alpha}=0.$

Applying the same technique, we can extend our approach to handle the nonlocal
version of (1.1) within the framework of fractional Orlicz-Sobolev spaces of
order $s\in(0,1)$ introduced in [6]. Specifically, for a fixed $\alpha>0$ we
analyze the minimization problem defined by

$E^{s}(\alpha)=\inf\left\\{\iint_{\mathbb{R}^{2n}}A\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\,d\nu_{n}\colon
u\in C^{\infty}_{c}(\Omega),\int_{\Omega}A(|u|)\,dx=\alpha\right\\}$ (1.5)

where we denote $d\nu_{n}:=|x-y|^{-n}dxdy$. As in the local case, without any
growth assumption on the Young function $A$, in [21] it is proved the
existence of a function $u^{s}_{\alpha}$ for which the infimum in (1.5) is
attained. Moreover, this function satisfies in the weak sense the equation

$\displaystyle\begin{cases}(-\Delta_{a})^{s}u=\lambda^{s}(\alpha)\frac{a(|u|)}{|u|}u&\quad\text{
in }\Omega\\\ u=0&\quad\text{ in }\partial\Omega,\end{cases}$ (1.6)

where the fractional $a-$Laplacian of order $s\in(0,1)$ is defined as

$(-\Delta_{a})^{s}u:=\text{p.v.}\,\displaystyle\int_{\mathbb{R}^{n}}a\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{dy}{|x-y|^{n+s}}$

and $\lambda^{s}(\alpha)$ is the Lagrange’s multiplier. Problems (1.5) and
(1.6) (and further variants) were widely studied, see for instance [3, 4, 5,
6, 7, 19, 21].

### Organization of the paper

After this introduction, in Section 2 we review the results on Young functions
and Orlicz and Orlicz-Sobolev spaces needed in the rest of the article. Most
results in this section are well known and in the cases where we were not able
to find a precise reference, a full proof is provided. Then, we include a
section (Section 3) where, assuming the $\Delta_{2}-$condition on $A$, we
deduce some elementary estimates on the energy $E(\alpha)$ and on the
eigenvalue $\lambda(\alpha)$. In particular, we obtain the estimates (1.3) and
(1.4).

Section 4 is devoted to the study of the regularity properties of the energy
function $E(\alpha)$. The main result being Theorem 4.6.

The behavior of the energy quotient $E(\alpha)/\alpha$ as $\alpha\to 0$ and
$\alpha\to\infty$ is analyzed in Sections 5 and 6. In Section 5 we analyze the
more interesting case where $A$ verifies the $\Delta_{2}-$condition either at
0 or at $\infty$ and in Section 6 the case in which $A$ does not verifies the
$\Delta_{2}-$condition.

Finally, in Section 7 we state our results for the nonlocal case.

##  2\. Preliminaries

In this section we collect all of the preliminaries needed in this work.
Almost all of the results are well known and we try to provide references in
all the cases. There are a couple of new, or not so well known results and in
those cases full proofs are provided.

The standard reference for this section is the book [11].

###  2.1. Young functions

A function $A\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ is called a _Young
function_ if it can be written as

$A(t)=\int_{0}^{t}a(\tau)\,d\tau,$

where $a\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ is nondecreasing, right
continuous, $a(t)>0$ if $t>0$, $a(0)=0$ and $a(t)\to\infty$ as $t\to\infty$.

From this definition, it can readily be checked that $A$ is convex, $A(0)=0$
and it is superlinear at $0$ and at $\infty$, i.e.

$\lim_{t\to 0}\frac{A(t)}{t}=0,\quad\lim_{t\to\infty}\frac{A(t)}{t}=\infty.$

Observe that the convexity of the Young function immediately gives that

$A(\tau t)\leq\tau A(t)\quad\text{ for }0<\tau<1,\qquad A(\tau t)\geq\tau
A(t)\quad\text{ for }\tau>1.$ (2.1)

Given a Young function $A$, an important associated function is the so-called
_complementary function_ , denoted by $\bar{A}$, that is defined as

$\bar{A}(t):=\sup\\{\tau t-A(\tau)\colon\tau\geq 0\\}.$

The complementary function $\bar{A}$ is also a Young function.

Observe that $\bar{A}$ is the optimal Young function in the Young-type
inequality

$\tau t\leq\bar{A}(\tau)+A(t)$

for all $\tau,t\geq 0$. It is also easy to see that $\bar{\bar{A}}=A$.

We now define the class

$\Delta_{2}^{\infty}:=\left\\{A\colon\mathbb{R}_{+}\to\mathbb{R}_{+}\colon\text{Young
function such that }\limsup_{t\to\infty}\frac{A(2t)}{A(t)}<\infty\right\\}$

and

$\Delta_{2}^{0}:=\left\\{A\colon\mathbb{R}_{+}\to\mathbb{R}_{+}\colon\text{Young
function such that }\limsup_{t\to 0}\frac{A(2t)}{A(t)}<\infty\right\\}.$

It is easy to see that $A\in\Delta_{2}^{\infty}$ if and only if there exists
constants $C_{\infty}\geq 2$ and $T_{\infty}>0$ such that

$A(2t)\leq C_{\infty}A(t)\quad\text{ for }t\geq T_{\infty},$

and $A\in\Delta_{2}^{0}$ if and only if there exists $C_{0}\geq 2$ and
$T_{0}>0$ such that

$A(2t)\leq C_{0}A(t)\quad\text{ for }t\leq T_{0}.$

Finally, we define the class

$\Delta_{2}=\Delta_{2}^{0}\cap\Delta_{2}^{\infty}$

and $A\in\Delta_{2}$ if and only if

$\sup_{t>0}\frac{A(2t)}{A(t)}<\infty.$

In [11, Theorem 4.1] it is shown that $A\in\Delta_{2}^{\infty}$ (or
$A\in\Delta_{2}^{0}$, or $A\in\Delta_{2}$) if and only if there exists $p>1$
such that

$\frac{ta(t)}{A(t)}\leq p\quad\text{for all }t\geq T_{\infty}\text{ (or $t\leq
T_{0}$, or $t>0$)}.$ (2.2)

In order to analyze the behavior of a Young function $A$ at 0 and at $\infty$,
we recall the Matuszewska-Orlicz functions that are defined as

$M_{0}(t)=\limsup_{\tau\to 0}\frac{A(\tau t)}{A(\tau)}\qquad\text{and}\qquad
M_{\infty}(t)=\limsup_{\tau\to\infty}\frac{A(\tau t)}{A(\tau)},$ (2.3)

These functions were introduced by Matuszewska and Orlicz in [16] and the
reader can consult with the book [15] for the results presented here. For
$i=0,\infty$, we have that $M_{i}$ is convex and sub-multiplicative, that is
$M_{i}\colon\mathbb{R}_{+}\to[0,\infty]$ verify

$M_{i}(0)=0,\quad M_{i}(1)=1,\quad M_{i}(\tau t)\leq M_{i}(\tau)M_{i}(t).$

When $A\in\Delta_{2}^{i}$, it follows that $M_{i}(t)<\infty$ for every $t>0$.

In the case where $A\not\in\Delta_{2}^{i}$ for some $i=0,\infty$, the
Matuszewska function $M_{i}$ turns out to be trivial. We were not able to find
this result in the literature, so we include a proof of this fact.

######  Proposition 2.1.

Let $A$ be a Young function. Then, if $A\not\in\Delta_{2}^{i}$, $i=0,\infty$,
we have that

$M_{i}(t)=\begin{cases}1&\text{if }t=1\\\ \infty&\text{if }t>1.\end{cases}$

Moreover, if $A$ satisfies the slightly stronger assumption

$\liminf_{t\to i}\frac{A(2t)}{A(t)}=\infty,$ (2.4)

then $M_{i}(t)=0$ if $0<t<1$.

###### Proof.

We prove the result for $i=\infty$, the other case being analogous.

For $t>1$ we argue as in [11, Theorem 4.1], to obtain

$A(st)\geq\int_{s}^{st}a(\tau)\,d\tau\geq(t-1)sa(s).$

Then,

$\limsup_{s\to\infty}\frac{sa(s)}{A(s)}\leq\frac{1}{t-1}\limsup_{s\to\infty}\frac{A(st)}{A(s)}=\frac{1}{t-1}M_{\infty}(t).$
(2.5)

Therefore $A\not\in\Delta_{2}^{\infty}$ implies that $M_{\infty}(t)=\infty$
for $t>1$.

Observe that $M_{\infty}(1)=1$ by definition.

Now, if (2.4) holds, then the limit in (2.5) exists. And so, for $0<t<1$, we
have

$\frac{A(st)}{A(s)}=\frac{1}{\frac{A(st\frac{1}{t})}{A(st)}}.$

Therefore the result follows by taking the limit in the former expression. ∎

As we mentioned, if $A\in\Delta_{2}^{i}$, then $M_{i}(t)<\infty$ for every
$t>0$. In the case where the limit in the definition of $M_{i}$ (2.3) exists,
these Matuszewska function are power functions. This is the content of the
next result.

######  Proposition 2.2.

Assume $A\in\Delta_{2}^{i}$ for some $i=0,\infty$. Assume moreover that the
limit in (2.3) exists. Then, there exists $p_{i}\in[1,\infty)$ such that

$M_{i}(t)=t^{p_{i}}.$

###### Proof.

First observe that, by (2.1), it follows that

$M_{i}(t)\leq t\text{ for }0<t<1\quad\text{and}\quad M_{i}(t)\geq t\text{ for
}t>1.$ (2.6)

Now, the proof follows just observing that if the limit in (2.3) exists, then
$M_{i}(t)$ is multiplicative, i.e. $M_{i}(st)=M_{i}(s)M_{i}(t)$ and
$M_{i}(1)=1$.

Then, if we define $v_{i}\colon\mathbb{R}\to\mathbb{R}$ by
$v_{i}(t)=\ln(M_{i}(e^{t}))$, the function $v_{i}$ is additive, i.e.
$v_{i}(s+t)=v_{i}(s)+v_{i}(t)$ for any $s,t\in\mathbb{R}$. It is a well known
fact that measurable additive functions are linear, therefore there exists
$p_{i}\in\mathbb{R}$ such that $v_{i}(t)=p_{i}t$ from where we can deduce that
$M_{i}(t)=t^{p_{i}}$. This fact, together with (2.6), implies that $p_{i}\geq
1$. ∎

####  2.1.1. Examples

######  Example 2.3.

Let $p>1$, and assume that

$A(t)\sim t^{p}\quad\text{when }t\ll 1.$

Then $A\in\Delta_{2}^{0}$ and $M_{0}(t)=t^{p}$.

If we assume that

$A(t)\sim t^{p}\quad\text{when }t\gg 1,$

then $A\in\Delta_{2}^{\infty}$ and $M_{\infty}(t)=t^{p}$.

As a special case we have the $(p,q)-$Laplacian: if $1<p<q<\infty$

$A(t)=\frac{t^{p}}{p}+\frac{t^{q}}{q},\quad\text{then
}A\in\Delta_{2}\quad\text{and }M_{0}(t)=t^{p},\quad M_{\infty}(t)=t^{q}.$

######  Example 2.4.

Let $\alpha\geq 0$ and $p\geq 1$. Then if

$A(t)\sim t^{p}\ln^{\alpha}t\quad\text{when }t\ll 1,$

then $A\in\Delta_{2}^{0}$ and $M_{0}(t)=t^{p}$. If

$A(t)\sim t^{p}\ln^{\alpha}t\quad\text{when }t\gg 1,$

then $A\in\Delta_{2}^{\infty}$ and $M_{\infty}(t)=t^{p}$.

As a special case, if

$A(t)=\frac{t^{p}}{p}\ln^{\alpha}(1+t^{r}),$

for some $r>0$, then $A\in\Delta_{2}$ and $M_{0}(t)=t^{p+r\alpha}$,
$M_{\infty}(t)=t^{p}$.

######  Example 2.5.

For $\alpha,\beta\geq 0$ and $p\geq 1$, assume that

$A(t)\sim t^{p}\ln^{\alpha}(t)\ln^{\beta}(\ln t)\quad\text{when }t\gg 1.$

Then $A\in\Delta_{2}^{\infty}$ and $M_{\infty}(t)=t^{p}$.

######  Example 2.6.

For $n\in\mathbb{N}$, define $A(t)=e^{t}-\sum_{k=0}^{n-1}\frac{t^{k}}{k!}$.
Then $A\in\Delta_{2}^{0}$, but $A\not\in\Delta_{2}^{\infty}$. Moreover,

$M_{0}(t)=t^{n}\quad\text{and}\quad M_{\infty}(t)=\begin{cases}0&\text{if
}0<t<1\\\ 1&\text{if }t=1\\\ \infty&\text{if }t>1.\end{cases}$

######  Example 2.7.

For $\alpha>0$, assume that

$A(t)\sim e^{-t^{-\alpha}}\quad\text{when }t\ll 1.$

Then $A\not\in\Delta_{2}^{0}$ and $M_{0}(t)=\begin{cases}0&\text{if }0<t<1\\\
1&\text{if }t=1\\\ \infty&\text{if }t>1.\end{cases}$

######  Example 2.8.

Assume that

$A(t)\sim e^{e^{t}}\quad\text{when }t\gg 1.$

Then, $A\not\in\Delta_{2}^{\infty}$ and $M_{\infty}(t)=\begin{cases}0&\text{if
}0<t<1\\\ 1&\text{if }t=1\\\ \infty&\text{if }t>1.\end{cases}$

Observe that all of the examples in this subsection satisfy the hypotheses of
Proposition 2.1 or Proposition 2.2.

###  2.2. Orlicz and Orlicz-Sobolev spaces

As we mentioned in the introduction on this section, the main reference for
Orlicz spaces is the book [11]. As for Orlicz-Sobolev spaces, the reader can
consult, for instance, with [10].

Fractional order Orlicz-Sobolev spaces, as we will use them here, were
introduced in [6] and then further analyze by several authors. The results
used in this paper are found in [1, 2] and in [5].

####  2.2.1. Orlicz spaces

Given $\Omega\subset\mathbb{R}^{n}$ a bounded domain and a Young function $A$,
the Orlicz class is defined as

$\mathcal{L}^{A}(\Omega):=\left\\{u\in
L^{1}_{\text{loc}}(\Omega)\colon\int_{\Omega}A(|u|)\,dx<\infty\right\\}.$

Then the Orlicz space $L^{A}(\Omega)$ is defined as the linear hull of
$\mathcal{L}^{A}(\Omega)$ and is characterized as

$L^{A}(\Omega)=\left\\{u\in L^{1}_{\text{loc}}(\Omega)\colon\text{ there
exists }k>0\text{ such that
}\int_{\Omega}A\left(\frac{|u|}{k}\right)\,dx<\infty\right\\}.$

In general the Orlicz class is strictly smaller than the Orlicz space, and
$\mathcal{L}^{A}(\Omega)=L^{A}(\Omega)$ if and only if
$A\in\Delta_{2}^{\infty}$.

The space $L^{A}(\Omega)$ is a Banach space when it is endowed, for instance,
with the Luxemburg norm, i.e.

$\|u\|_{L^{A}(\Omega)}=\|u\|_{A}:=\inf\left\\{k>0\colon\int_{\Omega}A\left(\frac{|u|}{k}\right)\,dx\leq
1\right\\}.$

This space $L^{A}(\Omega)$ turns out to be separable if and only if
$A\in\Delta_{2}^{\infty}$.

An important subspace of $L^{A}(\Omega)$ is $E^{A}(\Omega)$ that it is defined
as the closure of the functions in $L^{A}(\Omega)$ that are bounded. This
space is characterized as

$E^{A}(\Omega)=\left\\{u\in
L^{1}_{\text{loc}}(\Omega)\colon\int_{\Omega}A\left(\frac{|u|}{k}\right)\,dx<\infty\text{
for all }k>0\right\\}.$

This subspace $E^{A}(\Omega)$ is separable, and we have the inclusions

$E^{A}(\Omega)\subset\mathcal{L}^{A}(\Omega)\subset L^{A}(\Omega)$

with equalities if and only if $A\in\Delta_{2}^{\infty}$.

The following duality relation holds

$(E^{A}(\Omega))^{*}=L^{\bar{A}}(\Omega),$

where the equality is understood via the standard duality pairing. Observe
that this automatically implies that $L^{A}(\Omega)$ is reflexive if and only
if $A,\bar{A}\in\Delta_{2}^{\infty}$.

In the sequel, we will need a well known lemma (see for instance [11]), but we
include the details for completeness.

######  Lemma 2.9.

Let $A$ be a Young function and $u\in L^{A}(\Omega)$. Then

$\|u\|_{A}\leq\max\left\\{1;\int_{\Omega}A(|u|)\,dx\right\\}.$

###### Proof.

We can assume that $u\in\mathcal{L}^{A}(\Omega)$, otherwise there is nothing
to prove. Also, we can assume that $\|u\|_{A}>1$. Then,

$\int_{\Omega}A(|u|)\,dx=\int_{\Omega}A\left(\frac{|u|}{\|u\|_{A}}\|u\|_{A}\right)\,dx\geq\|u\|_{A},$

where we have used (2.1) and the definition of the Luxemburg norm. ∎

####  2.2.2. Orlicz-Sobolev spaces

Once the Orlicz spaces are defined, one can easily define the Orlicz-Sobolev
spaces in the usual way, namely

$W^{1}L^{A}(\Omega)=\\{u\in L^{A}(\Omega)\colon\partial_{x_{i}}u\in
L^{A}(\Omega),i=1,\dots,n\\}$

and

$W^{1}E^{A}(\Omega)=\\{u\in E^{A}(\Omega)\colon\partial_{x_{i}}u\in
E^{A}(\Omega),i=1,\dots,n\\}.$

These two spaces coincide only if $A\in\Delta_{2}^{\infty}$ and in this case
we will denote the Orlicz-Sobolev space as

$W^{1,A}(\Omega)=W^{1}L^{A}(\Omega)=W^{1}E^{A}(\Omega).$

These spaces are endowed with the natural norm

$\|u\|_{W^{1}L^{A}(\Omega)}=\|u\|_{1,A}=\|u\|_{A}+\|\nabla u\|_{A}.$

With this norm, $W^{1}L^{A}(\Omega)$ is a Banach space and
$W^{1}E^{A}(\Omega)$ is a closed subspace.

In order to work with Dirichlet boundary conditions it is necessary to define
the Orlicz-Sobolev functions that vanish at the boundary. So, we define
$W^{1}_{0}L^{A}(\Omega)$ as the closure of $C^{\infty}_{c}(\Omega)$ with
respect to the topology $\sigma(W^{1}L^{A}(\Omega),W^{1}E^{\bar{A}}(\Omega))$
and $W^{1}_{0}E^{A}(\Omega)$ as the closure of $C^{\infty}_{c}(\Omega)$ in
norm topology.

Again, when $A\in\Delta_{2}^{\infty}$ these spaces coincide and will be
denoted as

$W^{1,A}_{0}(\Omega)=W^{1}_{0}L^{A}(\Omega)=W^{1}_{0}E^{A}(\Omega).$

These spaces are reflexive if and only if both
$A,\bar{A}\in\Delta_{2}^{\infty}$. In any other case, what we have is that
there exists a Banach space $X$ such that $X^{*}=W^{1}_{0}L^{A}(\Omega)$. This
Banach space $X$ can be easily characterized in terms of distributions, but it
will not be used in this article.

In [8] the following is proved.

######  Proposition 2.10.

Let $A$ be a Young function. Then the embedding $W^{1}_{0}L^{A}(\Omega)\subset
L^{A}(\Omega)$ is compact.

####  2.2.3. Fractional Orlicz-Sobolev spaces

The construction of the fractional order Orlicz-Sobolev spaces is almost the
same as the Orliz-Sobolev spaces. We will make only a sketch and the reader
can consult with [5].

Given a fractional parameter $s\in(0,1)$, we define the Hölder quotient of a
function $u\in L^{A}(\Omega)$ as

$D^{s}u(x,y)=\frac{u(x)-u(y)}{|x-y|^{s}}.$

Then, the fractional Orlicz-Sobolev space of order $s$ is defined as

$W^{s}L^{A}(\mathbb{R}^{n}):=\\{u\in L^{A}(\mathbb{R}^{n})\colon D^{s}u\in
L^{A}(\mathbb{R}^{2n},d\nu_{n})\\},$

where $d\nu_{n}=|x-y|^{-n}dxdy$ and

$W^{s}E^{A}(\mathbb{R}^{n}):=\\{u\in E^{A}(\mathbb{R}^{n})\colon D^{s}u\in
E^{A}(\mathbb{R}^{2n},d\nu_{n})\\}.$

Since we are now working in the whole space, in order for these spaces to
coincide we require that $A\in\Delta_{2}$. In this case, we denote

$W^{s,A}(\mathbb{R}^{n})=W^{s}L^{A}(\mathbb{R}^{n})=W^{s}E^{A}(\mathbb{R}^{n}).$

As in the previous case, we have that $W^{s}L^{A}(\mathbb{R}^{n})$ is
reflexive if and only if $A,\bar{A}\in\Delta_{2}$.

In these spaces the norm considered is

$\|u\|_{W^{s}L^{A}(\mathbb{R}^{n})}=\|u\|_{s,A}=\|u\|_{A}+\|D^{s}u\|_{A,d\nu_{n}}.$

Again, with this norm, $W^{s}L^{A}(\mathbb{R}^{n})$ is a Banach space and
$W^{s}E^{A}(\mathbb{R}^{n})$ is a closed subspace.

The space $W^{s}_{0}L^{A}(\Omega)$ is then defined as the closure of
$C^{\infty}_{c}(\Omega)$ with respect to the topology
$\sigma(W^{s}L^{A}(\mathbb{R}^{n}),W^{s}E^{\bar{A}}(\mathbb{R}^{n}))$ and
$W^{s}_{0}E^{A}(\Omega)$ as the closure of $C^{\infty}_{c}(\Omega)$ in norm
topology.

Also, in [5] it is shown that for any $A$, there exists a Banach space $X$
such that $X^{*}=W^{s}_{0}L^{A}(\Omega)$ and that the extension of Proposition
2.10 to the fractional case holds:

######  Proposition 2.11.

Let $A$ be a Young function. Then the embedding $W^{s}_{0}L^{A}(\Omega)\subset
L^{A}(\Omega)$ is compact.

###  2.3. A review on the existence of minimizers for the energy function

Given a Young function $A$ without any further assumption on its behavior, in
[8, Theorem 3.1] (see also [9, Theorem 4.1]) it is proved that for any
election of the normalization $\alpha>0$, the minimization problem $E(\alpha)$
is solvable:

######  Proposition 2.12.

Let $E(\alpha)$ be the energy function defined in (1.1). Then for every
$\alpha>0$, there exists a function $u_{\alpha}\in W^{1}_{0}L^{A}(\Omega)$
such that $E(\alpha)=\int_{\Omega}A(|\nabla u_{\alpha}|)\,dx$ and
$\int_{\Omega}A(|u_{\alpha}|)\,dx=\alpha$.

As mentioned in the introduction, if $A,\bar{A}\in\Delta_{2}$ it is
straightforward to see that by means of the Lagrange’s multiplier rule there
exists $\lambda(\alpha)$ such that

$\int_{\Omega}a(|\nabla u_{\alpha}|)\frac{\nabla u_{\alpha}\cdot\nabla
v}{|\nabla
u_{\alpha}|}\,dx=\lambda(\alpha)\int_{\Omega}a(|u_{\alpha}|)\frac{u_{\alpha}v}{|u_{\alpha}|}\,dx\qquad\text{
for all }v\in C_{c}^{\infty}(\Omega).$ (2.7)

That is, $u_{\alpha}$ solves the eigenvalue problem (1.2) in the weak sense.

By using a generalized version of the Lagrange’s multiplier rule, in [9,
Theorem 4.2] the previous claim is proved to hold without any growth
assumption on $A$:

######  Proposition 2.13.

Let $\alpha>0$ and let $u_{\alpha}$ be a minimizer of $E(\alpha)$. Then there
exists $\lambda(\alpha)>0$ for which $u_{\alpha}$ solves (2.7).

As for the fractional problem, this was tackled in [21] (see also [5] where
higher order eigenvalues were also studied) and the following result was
otained.

######  Proposition 2.14.

For any $\alpha>0$, there exists a minimizer $u_{\alpha}^{s}\in
W^{s}_{0}L^{A}(\Omega)$ of the energy function $E^{s}(\alpha)$ defined in
(1.5). Moreover, there exists $\lambda^{s}(\alpha)$ such that $u_{\alpha}^{s}$
is a weak solution of (1.6), that is

$\iint_{\mathbb{R}^{2n}}a(|D^{s}u_{\alpha}^{s}|)\frac{D^{s}u_{\alpha}^{s}D^{s}v}{|D^{s}u_{\alpha}^{s}|}\,d\nu_{n}=\lambda^{s}(\alpha)\int_{\Omega}a(|u_{\alpha}^{s}|)\frac{u_{\alpha}^{s}v}{|u_{\alpha}^{s}|}\,dx,$

for all $v\in C^{\infty}_{c}(\Omega)$.

##  3\. Elementary estimates on the energy function

In this section we will assume that the Young function $A$ verifies the
$\Delta_{2}-$condition either at 0 or at $\infty$ and deduce some immediate
bounds on the energy $E(\alpha)$ and on the eigenvalue $\lambda(\alpha)$.

To this end, let us consider $u\in C^{\infty}_{c}(\Omega)$ to be fixed such
that

$\int_{\Omega}A(|u|)\,dx=1.$

Then, since $A$ is continuous and nondecreasing, if we define the function

$\phi(r)=\phi_{A}^{u}(r)=\int_{\Omega}A(r|u|)\,dx,$ (3.1)

it follows that $\phi\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$ is continuous,
nondecreasing, $\phi(0)=0$, $\phi(1)=1$ and $\phi(\infty)=\infty$. Hence, for
any $\alpha>0$, there exists $r=r_{\alpha}>0$ such that
$\phi(r_{\alpha})=\alpha$.

Recall that if $\alpha\leq 1$, then $r_{\alpha}\leq 1$ and that if $\alpha\geq
1$ then $r_{\alpha}\geq 1$.

We want to estimate $r_{\alpha}$ in terms of $\alpha$. To this end, we recall
an easy lemma regarding Young functions

######  Lemma 3.1.

Let $A$ be a Young function such that $A\in\Delta_{2}^{0}$ (or
$A\in\Delta_{2}^{\infty}$) and let $p>1$ be defined in (2.2).

Then it holds that $\tau^{p}A(t)\leq A(\tau t)\leq\tau A(t)$ for $0<\tau<1$
and $t<T_{0}$ (or $\tau A(t)\leq A(\tau t)\leq\tau^{p}A(t)$ for $\tau>1$ and
$t>T_{\infty}$).

Hence, if $A\in\Delta_{2}$, then

$\min\\{\tau,\tau^{p}\\}A(t)\leq A(\tau
t)\leq\max\\{\tau,\tau^{p}\\}A(t),\quad\text{for $\tau>0$}.$

The proof of this lemma is elementary and is omitted.

With the help of Lemma 3.1 we can obtain a simple estimate on $r_{\alpha}$,
namely

######  Lemma 3.2.

Let $A\in\Delta_{2}$ and $u\in C^{\infty}_{c}(\Omega)$. Let
$\phi=\phi_{A}^{u}$ be defined in (3.1), $\alpha>0$ and let $r_{\alpha}>0$ be
such $\phi(r_{\alpha})=\alpha$. Then

$\min\\{\alpha,\alpha^{\frac{1}{p}}\\}\leq
r_{\alpha}\leq\max\\{\alpha,\alpha^{\frac{1}{p}}\\},$

where $p$ is the constant defined in (2.2).

###### Proof.

It is immediate from the definition of $r_{\alpha}$ and Lemma 3.1. ∎

Next, we combine Lemmas 3.1 and 3.2 to obtain a bound for the energy
$E(\alpha)$.

######  Lemma 3.3.

Given $\alpha>0$, the energy of level $\alpha$ is defined by (1.1). Then, if
$A\in\Delta_{2}$, we have that

$\min\\{\alpha^{p},\alpha^{\frac{1}{p}}\\}E(1)\leq
E(\alpha)\leq\max\\{\alpha^{p},\alpha^{\frac{1}{p}}\\}E(1)$

###### Proof.

Immediate from the previous lemmas. ∎

Next, we want to compare the energy $E(\alpha)$ to the eigenvalue
$\lambda(\alpha)$ of problem (1.2).

######  Lemma 3.4.

Let $A$ be a Young function such that $A\in\Delta_{2}$ and let $p$ be defined
in (2.2).

Let $\alpha>0$ be a normalization parameter and let $u\in
W^{1}_{0}L^{A}(\Omega)$ be the minimizer associated to the energy $E(\alpha)$.
Let $\lambda(\alpha)$ be the eigenvalue of (1.2) associated to $u$. Then, we
have

$\frac{1}{p}\frac{E(\alpha)}{\alpha}\leq\lambda(\alpha)\leq
p\frac{E(\alpha)}{\alpha}.$

###### Proof.

First, observe that from [9, Theorem 4.23], there exists a minimizer $u\in
W^{1}_{0}L^{A}(\Omega)$ for $E(\alpha)$ and $u$ is a weak solution to (1.2).

Now, if we use $u$ as a test function in the weak formulation of (1.2), we get
that

$\lambda(\alpha)=\frac{\int_{\Omega}a(|\nabla u|)|\nabla
u|\,dx}{\int_{\Omega}a(|u|)|u|\,dx}.$ (3.2)

Next, just observe that we have the inequality

$A(t)\leq a(t)t\leq pA(t),$ (3.3)

where the first inequality follows from the monotonicity of $a(t)$ and the
second one is the $\Delta_{2}-$condition.

The result now follows just combining (3.2) and (3.3). ∎

######  Remark 3.5.

Observe that combining Lemmas 3.3 and 3.4 we immediately get that there exist
$c,C>0$ such that

$c\min\\{\alpha^{p-1},\alpha^{\frac{1}{p}-1}\\}\leq\frac{E(\alpha)}{\alpha},\lambda(\alpha)\leq
C\max\\{\alpha^{p-1},\alpha^{\frac{1}{p}-1}\\}.$

Also observe that these bounds trivialize when $\alpha\to 0$ and when
$\alpha\to\infty$.

##  4\. Continuity and differentiability of the energy function

In this section we study the regularity properties of the energy function
$E(\alpha)$ defined in (1.1). First, without assuming the $\Delta_{2}$
condition on $A$ nor $\bar{A}$, we prove the continuity of $E(\alpha)$

######  Theorem 4.1.

The energy function $E(\alpha)$ is continuous for every $\alpha>0$.

We split the proof into two lemmas.

######  Lemma 4.2.

Let $A$ be a Young function and let $E(\alpha)$ be the energy function defined
in (1.1). Then, $E$ is upper semicontinuous.

###### Proof.

Let $\alpha_{0}>0$ be fixed. Given $\varepsilon>0$, take
$u_{0}=u_{\varepsilon,\alpha_{0}}\in C^{\infty}_{c}(\Omega)$ be such that

$\int_{\Omega}A(|\nabla u_{0}|)\,dx\leq
E(\alpha_{0})+\varepsilon\quad\text{and}\quad\int_{\Omega}A(|u_{0}|)\,dx=\alpha_{0}.$

Let $\phi(r)=\phi_{A}^{u_{0}}(r)$ be the function defined in (3.1) and given
$\alpha>0$, let $r_{\alpha}>0$ be such that $\phi(r_{\alpha})=\alpha$.

Observe that
$\lim_{\alpha\to\alpha_{0}}\phi(r_{\alpha})=\lim_{\alpha\to\alpha_{0}}\alpha=\alpha_{0}=\phi(1)$.
But, since $\phi$ is strictly increasing and continuous, it follows that
$\phi^{-1}$ is also continuous and strictly increasing, from where it follows
that $r_{\alpha}\to 1$ when $\alpha\to\alpha_{0}$.

Hence $r_{\alpha}u_{0}$ is an admissible test for $E(\alpha)$ and we get that

$\limsup_{\alpha\to\alpha_{0}}E(\alpha)\leq\limsup_{\alpha\to\alpha_{0}}\int_{\Omega}A(r_{\alpha}|\nabla
u_{0}|)\,dx=\int_{\Omega}A(|\nabla u_{0}|)\,dx\leq E(\alpha_{0})+\varepsilon.$

Since $\varepsilon>0$ is arbitrary, the result follows. ∎

######  Lemma 4.3.

Let $A$ be a Young function and let $E(\alpha)$ be the energy function defined
in (1.1). Then, $E$ is lower semicontinuous.

###### Proof.

Let $\alpha_{0}>0$ be fixed and consider $\alpha_{k}\to\alpha_{0}$ as
$k\to\infty$. Let $\varepsilon>0$ and for each $k\in\mathbb{N}$ let $u_{k}\in
C^{\infty}_{c}(\Omega)$ be such that

$\int_{\Omega}A(|\nabla u_{k}|)\,dx\leq
E(\alpha_{k})+\varepsilon\quad\text{and}\quad\int_{\Omega}A(|u_{k}|)\,dx=\alpha_{k}.$

Observe that by Lemma 4.2, $E(\alpha_{k})$ is bounded, and hence
$\int_{\Omega}A(|\nabla u_{k}|)\,dx$ is bounded. Therefore, by Lemma 2.9,
$\\{u_{k}\\}_{k\in\mathbb{N}}$ is bounded in $W^{1}_{0}L^{A}(\Omega)$.

This allows to apply Proposition 2.10 to derive the existence of a function
$u\in W^{1}_{0}L^{A}(\Omega)$ such that

$u_{k}\rightharpoonup u\text{ weakly* in
}W^{1}_{0}L^{A}(\Omega)\quad\text{and}\quad u_{k}\to u\text{ strongly in
}L^{A}(\Omega).$

From this we conclude that $\int_{\Omega}A(|u|)\,dx=\alpha_{0}$ and hence, it
is an admissible competitor for $E(\alpha_{0})$, from where

$E(\alpha_{0})\leq\int_{\Omega}A(|\nabla
u|)\,dx\leq\liminf_{k\to\infty}\int_{\Omega}A(|\nabla
u_{k}|)\,dx\leq\liminf_{k\to\infty}E(\alpha_{k})+\varepsilon,$

where for the second inequality we have used [8, Lemma 2.2]. This estimate
concludes the proof. ∎

From the proof of Theorem 4.1 one immediately obtains the following
information:

######  Proposition 4.4.

Let $\alpha_{0}>0$ and let $\alpha_{k}\to\alpha_{0}$ as $k\to\infty$. Let
$\\{u_{k}\\}_{k\in\mathbb{N}}$ be a sequence of minimizers of $E(\alpha_{k})$.
Then, $\\{u_{k}\\}_{k\in\mathbb{N}}$ is weak* precompact in
$W^{1}_{0}L^{A}(\Omega)$ and each accumulation point $u\in
W^{1}_{0}L^{A}(\Omega)$ of the sequence in the weak*-topology is a minimizer
for $E(\alpha_{0})$.

Moreover, if $u_{k_{j}}\rightharpoonup u$ weakly* in $W^{1}_{0}L^{A}(\Omega)$,
then

$\int_{\Omega}A(|\nabla u|)\,dx=\lim_{j\to\infty}\int_{\Omega}A(|\nabla
u_{k_{j}}|)\,dx.$

Now, we analyze the differentiability of $E(\alpha)$. We stress that no
$\Delta_{2}-$condition is imposed on $A$ nor $\bar{A}$. In fact, we prove the
following:

######  Theorem 4.5.

Let $A$ be a Young function and let $E(\alpha)$ be the energy function defined
in (1.1). Then $E(\alpha)$ is Lipschitz continuous. Moreover,

$0\leq\liminf_{\alpha\to\alpha_{0}}\frac{E(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}\leq\limsup_{\alpha\to\alpha_{0}}\frac{E(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}\leq\lambda(\alpha).$

###### Proof.

Observe that the lower bound is immediate, since $E(\alpha)$ is increasing.

Now, Let $\alpha_{0}>0$ and let $u_{0}\in W^{1}_{0}L^{A}(\Omega)$ be a
minimizer of $E(\alpha_{0})$, that is,

$E(\alpha_{0})=\int_{\Omega}A(|\nabla
u_{0}|)\,dx\quad\text{and}\quad\int_{\Omega}A(|u_{0}|)\,dx=\alpha_{0}.$

Define $F(r,\alpha)=\int_{\Omega}A(r|u_{0}|)\,dx-\alpha$. Observe that

$\frac{\partial F}{\partial
r}(r,\alpha)=\int_{\Omega}a(r|u_{0}|)|u_{0}|\,dx,\qquad\frac{\partial
F}{\partial\alpha}(r,\alpha)=-1.$

Moreover, since $u_{0}\neq 0$, $\frac{\partial F}{\partial
r}(1,\alpha_{0})=\int_{\Omega}a(|u_{0}|)|u_{0}|\,dx>0$ . Then we can apply the
Implicit Function Theorem to obtain the existence of a function $r(\alpha)$
defined for $\alpha\in I:=(\alpha_{0}-\delta,\alpha_{0}+\delta)$ for some
$\delta>0$, for which $r(\alpha_{0})=1$, $F(r(\alpha),\alpha)=0$ for any
$\alpha\in I$. Moreover,

$r^{\prime}(\alpha_{0})=-\frac{\frac{\partial
F}{\partial\alpha}(1,\alpha_{0})}{\frac{\partial F}{\partial
r}(1,\alpha_{0})}=\frac{1}{\int_{\Omega}a(|u_{0}|)||u_{0}|\,dx}>0.$

Observe that the function $G(\alpha):=\int_{\Omega}A(r(\alpha)|\nabla
u_{0}|)\,dx$, $\alpha\in I$, near $\alpha_{0}$ behaves as

$\displaystyle G(\alpha)$
$\displaystyle=G(\alpha_{0})+G^{\prime}(\alpha_{0})(\alpha-\alpha_{0})+o(\alpha-\alpha_{0})$
$\displaystyle=E(\alpha_{0})+(\alpha-\alpha_{0})r^{\prime}(\alpha_{0})\int_{\Omega}a(|\nabla
u_{0}|)|\nabla u_{0}|\,dx+o(\alpha-\alpha_{0}).$

Since $F(r(\alpha),\alpha)=0$ for any $\alpha\in I$, the function
$r(\alpha)u_{0}$ is admissible in the minimization problem $E(\alpha)$ for any
$\alpha\in I$, giving that

$\frac{E(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}\leq\frac{G(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}=\frac{\int_{\Omega}a(|\nabla
u_{0}|)|\nabla
u_{0}|\,dx}{\int_{\Omega}a(|u_{0}|)|u_{0}|\,dx}+o(1)=\lambda(\alpha_{0})+o(1),$

as we wanted to prove. ∎

In order to obtain the differentiability of the energy function, we need to
impose some technical, but somewhat natural additional conditions on $A$.

######  Theorem 4.6.

Let $A$ be a Young function such that $t\mapsto a(t)t$ is convex.

Let $E(\alpha)$ be the energy function defined in (1.1). Then $E(\alpha)$ is
differentiable for any $\alpha>0$. Moreover,

$\displaystyle
E^{\prime}(\alpha)=\lambda(\alpha),\qquad\left(\frac{E(\alpha)}{\alpha}\right)^{\prime}=\frac{1}{\alpha}\left(\lambda(\alpha)-\frac{E(\alpha)}{\alpha}\right),$

where $\lambda(\alpha)$ is the eigenvalue defined in (1.2).

###### Proof.

In view of Theorem 4.5, it remains to see that

$\liminf_{\alpha\to\alpha_{0}}\frac{E(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}\geq\lambda(\alpha_{0}).$

Let $\\{\alpha_{k}\\}_{k\in\mathbb{N}}$ be such that $\alpha_{k}\to\alpha_{0}$
when $k\to\infty$, and

$\liminf_{\alpha\to\alpha_{0}}\frac{E(\alpha)-E(\alpha_{0})}{\alpha-\alpha_{0}}=\lim_{k\to\infty}\frac{E(\alpha_{k})-E(\alpha_{0})}{\alpha_{k}-\alpha_{0}}.$

Arguiying analogously as in the proof of Theorem 4.1, we can assume that there
exists $\\{u_{k}\\}_{k\in\mathbb{N}}\subset W^{1}_{0}L^{A}(\Omega)$ such that

$E(\alpha_{k})=\int_{\Omega}A(|\nabla
u_{k}|)\,dx,\qquad\int_{\Omega}A(|u_{k}|)\,dx=\alpha_{k}$

and that

$\begin{split}&u_{k}\rightharpoonup u\text{ weakly* in
}W^{1}_{0}L^{A}(\Omega)\quad\text{ and }\\\ &u_{k}\to u\text{ strongly in
}L^{A}(\Omega)\text{ and pointwisely}.\end{split}$ (4.1)

Moreover, by Proposition 4.4, $u\in W^{1}_{0}L^{A}(\Omega)$ is a minimizer for
$E(\alpha_{0})$.

Now, proceeding as Theorem 4.5 for each $k\in\mathbb{N}$ there exists a
function $r_{k}(\alpha)$ defined in $t\in
I_{k}:=(\alpha_{k}-\delta,\alpha_{k}+\delta)$ for some $\delta>0$ for which

$r_{k}(\alpha_{k})=1,\quad\int_{\Omega}A(r_{k}(\alpha)|u_{k}|)\,dx=\alpha,\quad
r_{k}^{\prime}(\alpha_{k})=\frac{1}{\int_{\Omega}a(|u_{k}|)|u_{k}|\,dx}>0.$

Hence, testing in $E(\alpha_{0})$ with $r_{k}(\alpha_{0})u_{k}$ and performing
a Taylor expansion around $\alpha=\alpha_{k}$ we get that

$\displaystyle E(\alpha_{0})$
$\displaystyle\leq\int_{\Omega}A(r_{k}(\alpha_{0})|\nabla u_{k}|)\,dx$
$\displaystyle=E(\alpha_{k})+(\alpha_{0}-\alpha_{k})r_{k}^{\prime}(\alpha_{k})\int_{\Omega}a(|\nabla
u_{k}|)|\nabla u_{k}|\,dx+o(\alpha_{0}-\alpha_{k}).$

Observe that, since $a(t)$ is continuous and $u_{k}$ is converging
pointwisely, it follows that $o(\alpha_{0}-\alpha_{k})$ is independent of
$k\in\mathbb{N}$.

Therefore,

$\displaystyle\lim_{k\to\infty}\frac{E(\alpha_{k})-E(\alpha_{0})}{\alpha_{k}-\alpha_{0}}$
$\displaystyle\geq\liminf_{k\to\infty}\frac{\int_{\Omega}a(|\nabla
u_{k}|)|\nabla u_{k}|\,dx}{\int_{\Omega}a(|u_{k}|)|u_{k}|\,dx}$
$\displaystyle\geq\frac{\int_{\Omega}a(|\nabla u|)|\nabla
u|\,dx}{\int_{\Omega}a(|u|)|u|\,dx}=\lambda(\alpha_{0}),$

where we have used [8, Lemma 2.2] and our convexity assumption on the function
$t\mapsto a(t)t$.

Finally, using the expression for $E^{\prime}(\alpha)$ we get that

$\displaystyle\left(\frac{E(\alpha)}{\alpha}\right)^{\prime}$
$\displaystyle=\frac{1}{\alpha}\left(E^{\prime}(\alpha)-\frac{E(\alpha)}{\alpha}\right)=\frac{1}{\alpha}\left(\lambda(\alpha)-\frac{E(\alpha)}{\alpha}\right).$

The proof is complete. ∎

######  Remark 4.7.

Observe that the eigenvalue $\lambda(\alpha)$ and the energy quotient
$E(\alpha)/\alpha$ coincide if and only if they are constants. This is exactly
the case where $A$ is homogeneous, that is $A(t)=t^{p}$ for some $p>1$.

##  5\. Sharp asymptotic behavior when $\alpha\to 0$ and $\alpha\to\infty$

In order to analyze the precise behavior at $0$ and at $\infty$, we need to
assume that $A\in\Delta_{2}^{i}$ at $i=0$ and $i=\infty$ respectively.
Moreover, we need to impose the slightly stronger assumption that the limit
defining $M_{i}(t)$ in (2.3) exists.

So, fixed a normalization parameter $\tau>0$ we define the minimization
problem for the associated Matuszewska functions,

$\displaystyle E_{i}(\tau)$
$\displaystyle=\inf\left\\{\int_{\Omega}M_{i}(|\nabla u|)\,dx\colon u\in
C_{c}^{\infty}(\Omega),\,\int_{\Omega}M_{i}(|u|)\,dx=\tau\right\\}.$ (5.1)

We then obtain the following result.

######  Theorem 5.1.

Let $A$ be a Young function satisfying the $\Delta_{2}^{i}$ condition for
$i=0$ or $i=\infty$ and assume that the corresponding limit defining
$M_{i}(t)$ in (2.3) exists. Then, it holds that:

$\limsup_{\alpha\to
i}\frac{E(\alpha)}{\alpha}=\inf_{\tau>0}\frac{E_{i}(\tau)}{\tau}.$

###### Proof.

We work in the case where $i=\infty$, the other case being identical. We will
split the proof into two steps.

_Step 1: Upper bound_.

Let $u\in C^{\infty}_{c}(\Omega)$ be a fixed function such that $u\neq 0$ and
consider the function $\phi=\phi_{A}^{u}$ defined in (3.1). Hence, given
$\alpha>0$ there exists a unique $r_{\alpha}>0$ such that
$\phi(r_{\alpha})=\alpha$. Since $\phi$ is continuous and $u$ is fixed we have
that $r_{\alpha}\to\infty$ if $\alpha\to\infty$.

Since we can write

$\frac{\alpha}{A(r_{\alpha})}=\int_{\Omega}\frac{A(r_{\alpha}|u|)}{A(r_{\alpha})}\,dx,$

according to the definition (2.3) and the Dominated Convergence Theorem, we
get that

$\lim_{\alpha\to\infty}\frac{\alpha}{A(r_{\alpha})}=\int_{\Omega}M_{\infty}(|u|)\,dx.$
(5.2)

The equality $\phi(r_{\alpha})=\alpha$ implies that the function $r_{\alpha}u$
is admissible in the minimization problem for $E(\alpha)$, therefore

$\begin{split}\limsup_{\alpha\to\infty}\frac{E(\alpha)}{\alpha}&\leq\lim_{\alpha\to\infty}\frac{1}{\alpha}\int_{\Omega}A(r_{\alpha}|\nabla
u|)\,dx\\\
&=\lim_{\alpha\to\infty}\frac{A(r_{\alpha})}{\alpha}\int_{\Omega}\frac{A(r_{\alpha}|\nabla
u|)}{A(r_{\alpha})}\,dx\\\ &=\frac{\int_{\Omega}M_{\infty}(|\nabla
u|)\,dx}{\int_{\Omega}M_{\infty}(|u|)\,dx},\end{split}$ (5.3)

where we have used (5.2) and (2.3)

Since $u\in C_{c}^{\infty}(\Omega)$ is arbitrary, from (5.3) it is deduced
that

$\limsup_{\alpha\to\infty}\frac{E(\alpha)}{\alpha}\leq\inf_{\tau>0}\frac{E_{\infty}(\tau)}{\tau}$

as needed.

_Step 2: Lower bound_. Let us derive the lower bound for the ratio
$E(\alpha)/\alpha$. For this, let $u\in C^{\infty}_{c}(\Omega)$ be fixed, and
as before, let $r_{\alpha}>0$ be such that $\phi(r_{\alpha})=\alpha$.

For $\varepsilon>0$ we have that

$\inf_{\tau>0}\frac{E_{\infty}(\tau)}{\tau}-\varepsilon\leq\frac{\int_{\Omega}M_{\infty}(|\nabla
u|)\,dx}{\int_{\Omega}M_{\infty}(|u|)\,dx}-\varepsilon.$

Then, according to (2.3), there exists
$\alpha_{\infty}=\alpha_{\infty}(\varepsilon,u)$ such that for
$\alpha>\alpha_{\infty}$,

$\frac{\int_{\Omega}M_{\infty}(|\nabla
u|)\,dx}{\int_{\Omega}M_{\infty}(|u|)\,dx}-\varepsilon\leq\frac{\int_{\Omega}\frac{A(r_{\alpha}|\nabla
u|)}{A(r_{\alpha})}\,dx}{\int_{\Omega}\frac{A(r_{\alpha}|u|)}{A(r_{\alpha})}\,dx}=\frac{1}{\alpha}\int_{\Omega}A(r_{\alpha}|\nabla
u|)\,dx.$

Taking infimum over $u\in C^{\infty}_{c}(\Omega)$, we conclude that

$\inf_{\tau>0}\frac{E_{\infty}(\tau)}{\tau}-\varepsilon\leq\frac{E(\alpha)}{\alpha}\quad\text{for
}\alpha>\alpha_{\infty},$

from where we deduce that

$\inf_{\tau>0}\frac{E_{\infty}(\tau)}{\tau}\leq\liminf_{\alpha\to\infty}\frac{E(\alpha)}{\alpha}$

and this concludes the proof. ∎

Lemma 4.2, and Theorem 5.1 lead to the uniform boundedness of the ratio
$E(\alpha)/\alpha$:

######  Corollary 5.2.

It holds that

$\sup_{\alpha>0}\frac{E(\alpha)}{\alpha}<\infty.$

As a direct consequence, we have the following.

######  Corollary 5.3.

With the same assumptions of Theorem 5.1, let $p_{i}$ be the exponent given by
Proposition 2.2. Then, for $i=0$ or $i=\infty$,

$\lim_{\alpha\to i}\frac{E(\alpha)}{\alpha}=\lambda_{p_{i}},$

where $\lambda_{p}$ denotes the first eigenvalue of the $p-$Laplacian with
Dirichlet boundary conditions in $\Omega$.

##  6\. The behavior of $E(\alpha)/\alpha$ without the $\Delta_{2}-$condition

In this section we explore the limit behavior of the energy quotient
$E(\alpha)/\alpha$ in the case where the Young function $A$ does not satisfy
the $\Delta_{2}^{i}$ condition for $i=0$ or $i=\infty$.

We will assume the stronger condition (2.4), and hence, by Proposition 2.1,
the Matuszewska function $M_{i}(t)$ is trivial. That is

$\lim_{s\to i}\frac{A(ts)}{A(s)}=M_{i}(t)=\begin{cases}0&\text{if }0\leq
t<1\\\ 1&\text{if }t=1\\\ \infty&\text{if }t>1.\end{cases}$

In this case, we have the following result.

######  Theorem 6.1.

Let $A$ be a Young function satisfying (2.4) for $i=0$ or $i=\infty$. Then, if
$\Omega\subset\mathbb{R}^{n}$ has inner radius larger that 1,

$\lim_{\alpha\to i}\frac{E(\alpha)}{\alpha}=0.$

###### Proof.

Recall that the inner radius of $\Omega$ is defined as

$R_{\Omega}=\sup\\{r>0\colon B_{r}(x)\subset\Omega\text{ for some
}x\in\Omega\\}.$

If $R_{\Omega}>1$, then there exists $r>0$, $x_{0}\in\Omega$ with
$B_{r}(x_{0})=B_{r}\subset\Omega$ and $u\in C^{\infty}_{c}(\Omega)$ such that

$u(x)=1\text{ for }x\in B_{r},\quad 0\leq u(x)\leq
1\quad\text{and}\quad|\nabla u(x)|<1\text{ for }x\in\Omega.$

Then, as usual, for $\alpha>0$ we define $r_{\alpha}$ such that
$\phi(r_{\alpha})=\phi_{A}^{u}(r_{\alpha})=\alpha$ and therefore

$\frac{E(\alpha)}{\alpha}\leq\frac{A(r_{\alpha})}{\alpha}\int_{\Omega}\frac{A(r_{\alpha}|\nabla
u|)}{A(r_{\alpha})}\,dx=\frac{\int_{\Omega}\frac{A(r_{\alpha}|\nabla
u|)}{A(r_{\alpha})}\,dx}{\int_{\Omega}\frac{A(r_{\alpha}|u|)}{A(r_{\alpha})}\,dx}.$

But, since $r_{\alpha}\to i$ as $\alpha\to i$, it follows from Proposition 2.1
that

$\int_{\Omega}\frac{A(r_{\alpha}|u|)}{A(r_{\alpha})}\,dx\to|\\{u=1\\}|=|B_{r}|>0\quad\text{and}\quad\int_{\Omega}\frac{A(r_{\alpha}|\nabla
u|)}{A(r_{\alpha})}\,dx\to 0,$

and the proof is complete.

∎

##  7\. Nonlocal case

In this section we analyze the problem of the behavior of the energy function
for the nonlocal case $E^{s}(\alpha)$ defined in (1.5).

All the program developed in Sections 3 to 5 can be carried out without any
change so we just include the results leaving the details to the interested
reader.

The results obtained in this case can be summarized in the following theorem

######  Theorem 7.1.

Let $A$ be a Young function. Then the nonlocal energy function $E^{s}(\alpha)$
defined in (1.5) is Lipschitz continuous. If in addition the function
$t\mapsto ta(t)$ is convex, then $E^{s}(\alpha)$ is differentiable and

$\frac{d}{d\alpha}E^{s}(\alpha)=\lambda^{s}(\alpha),\qquad\frac{d}{d\alpha}\left(\frac{E^{s}(\alpha)}{\alpha}\right)=\frac{1}{\alpha}\left(\lambda^{s}(\alpha)-\frac{E^{s}(\alpha)}{\alpha}\right),$

where $\lambda^{s}(\alpha)$ is the eigenvalue defined in (1.6).

Moreover, when $A\in\Delta_{s}^{i}$ for $i=0$ or $i=\infty$, and the limit in
(2.3) exists, then

$\lim_{\alpha\to
i}\frac{E^{s}(\alpha)}{\alpha}=\inf_{\tau>0}\frac{E^{s}_{i}(\tau)}{\tau}$

where

$E^{s}_{i}(\tau)=\inf\left\\{\iint_{\mathbb{R}^{2n}}M_{i}(|D^{s}u|)\,d\nu_{n}\colon
u\in C^{\infty}_{c}(\Omega),\ \int_{\Omega}M_{i}(|u|)\,dx=\tau\right\\}.$

## Acknowledgements

This work was partially supported by ANPCyT under grant PICT 2019-3837 and by
CONICET under grant PIP 11220150100032CO. Both authors are members of CONICET
and are grateful for the support.

## References

  * [1] Angela Alberico, Andrea Cianchi, Luboš Pick, and Lenka Slavíková, _Fractional Orlicz-Sobolev embeddings_ , J. Math. Pures Appl. (9) 149 (2021), 216–253. MR 4239001
  * [2] by same author, _On fractional Orlicz-Sobolev spaces_ , Harmonic analysis and partial differential equations—in honor of Vladimir Maz’ya, Birkhäuser/Springer, Cham, [2023] ©2023, pp. 45–65. MR 4592608
  * [3] Sabri Bahrouni, Hichem Ounaies, and Ariel Salort, _Variational eigenvalues of the fractional $g$-Laplacian_, Complex Var. Elliptic Equ. 68 (2023), no. 6, 1021–1044. MR 4596303
  * [4] Sabri Bahrouni and Ariel Salort, _Neumann and Robin type boundary conditions in fractional Orlicz-Sobolev spaces_ , ESAIM Control Optim. Calc. Var. 27 (2021), Paper No. S15, 23. MR 4222173
  * [5] Julián Fernández Bonder and Juan F. Spedaletti, _Ljusternik–schnirelmann eigenvalues for the fractional $m-$laplacian without the $\Delta_{2}$ condition_, ArXiv preprint (2024), 29pp.
  * [6] Julián Fernández Bonder and Ariel Salort, _Fractional order Orlicz-Sobolev spaces_ , J. Funct. Anal. 277 (2019), no. 2, 333–367. MR 3952156
  * [7] Julián Fernández Bonder, Ariel Salort, and Hernán Vivas, _Global Hölder regularity for eigenfunctions of the fractional $g$-Laplacian_, J. Math. Anal. Appl. 526 (2023), no. 1, Paper No. 127332, 15. MR 4580443
  * [8] M. García-Huidobro, V. K. Le, R. Manásevich, and K. Schmitt, _On principal eigenvalues for quasilinear elliptic differential operators: an Orlicz-Sobolev space setting_ , NoDEA Nonlinear Differential Equations Appl. 6 (1999), no. 2, 207–225. MR 1694787
  * [9] J.-P. Gossez and R. Manásevich, _On a nonlinear eigenvalue problem in Orlicz-Sobolev spaces_ , Proc. Roy. Soc. Edinburgh Sect. A 132 (2002), no. 4, 891–909. MR 1926921
  * [10] Jean-Pierre Gossez, _Nonlinear elliptic boundary value problems for equations with rapidly (or slowly) increasing coefficients_ , Trans. Amer. Math. Soc. 190 (1974), 163–205. MR 342854
  * [11] M. A. Krasnoselskiĭ and Ja. B. Rutickiĭ, _Convex functions and Orlicz spaces_ , russian ed., P. Noordhoff Ltd., Groningen, 1961. MR 126722
  * [12] Peter Lindqvist, _On the equation ${\rm div}\,(|\nabla u|^{p-2}\nabla u)+\lambda|u|^{p-2}u=0$_, Proc. Amer. Math. Soc. 109 (1990), no. 1, 157–164. MR 1007505
  * [13] by same author, _Notes on the $p$-Laplace equation_, Report. University of Jyväskylä Department of Mathematics and Statistics, vol. 102, University of Jyväskylä, Jyväskylä, 2006. MR 2242021
  * [14] by same author, _A nonlinear eigenvalue problem_ , Topics in mathematical analysis, Ser. Anal. Appl. Comput., vol. 3, World Sci. Publ., Hackensack, NJ, 2008, pp. 175–203. MR 2462954
  * [15] Lech Maligranda, _Orlicz spaces and interpolation_ , Seminários de Matemática [Seminars in Mathematics], vol. 5, Universidade Estadual de Campinas, Departamento de Matemática, Campinas, 1989. MR 2264389
  * [16] W. Matuszewska and W. Orlicz, _On certain properties of $\varphi$-functions_, Bull. Acad. Polon. Sci. Sér. Sci. Math. Astronom. Phys. 8 (1960), 439–443. MR 126158
  * [17] Marcelo Montenegro and Sebastián Lorca, _The eigenvalue problem for quasilinear elliptic operators with general growth_ , Appl. Math. Lett. 25 (2012), no. 7, 1045–1049. MR 2915122
  * [18] Vesa Mustonen and Matti Tienari, _An eigenvalue problem for generalized Laplacian in Orlicz-Sobolev spaces_ , Proc. Roy. Soc. Edinburgh Sect. A 129 (1999), no. 1, 153–163. MR 1669197
  * [19] Ariel Salort, _Eigenvalues and minimizers for a non-standard growth non-local operator_ , J. Differential Equations 268 (2020), no. 9, 5413–5439. MR 4066053
  * [20] by same author, _Lower bounds for Orlicz eigenvalues_ , Discrete Contin. Dyn. Syst. 42 (2022), no. 3, 1415–1434. MR 4385762
  * [21] Ariel Salort and Hernán Vivas, _Fractional eigenvalues in Orlicz spaces with no $\Delta_{2}$ condition_, J. Differential Equations 327 (2022), 166–188. MR 4413298
  * [22] Matti Tienari, _Ljusternik-Schnirelmann theorem for the generalized Laplacian_ , J. Differential Equations 161 (2000), no. 1, 174–190. MR 1740361


# The early Universe was dust-rich and extremely hot

Marco P. Viero,1 Guochao Sun,1 Dongwoo T. Chung,2,3 Lorenzo Moncelsi,1 Sam S.
Condon1
1California Institute of Technology, 1200 East California Boulevard, Pasadena,
CA, 91125, USA
2Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto, ON, M5S 3H8, Canada
3Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto, ON, M5S 3H4, Canada E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the dust properties and star-formation signature of galaxies in
the early Universe by stacking 111 227 objects in the recently released COSMOS
catalogue on maps at wavelengths bracketing the peak of warmed dust emission.
We find an elevated far-infrared luminosity density to redshift 8, indicating
abundant dust in the early Universe. We further find an increase of dust
temperature with redshift, reaching $100\pm 12\ \rm K$ at $z\sim 7$,
suggesting either the presence of silicate rich dust originating from
Population II stars, or sources of heating beyond simply young hot stars.
Lastly, we try to understand how these objects have been missed in previous
surveys, and how to design observations to target them. All code, links to the
data, and instructions to reproduce this research in full are located at
https://github.com/marcoviero/simstack3/.

###### keywords:

galaxies: star formation — infrared: galaxies — submillimetre: diffuse
background — early Universe

††pubyear: 2022††pagerange: The early Universe was dust-rich and extremely
hot–The early Universe was dust-rich and extremely hot

##  1 Introduction

The history of star formation at high redshift has always had a dust problem.
Star formation is well traced by ultraviolet emission, but a large fraction of
ultraviolet light is obscured by dust, thus proving it to be an incomplete
estimator. But dust also emits light, re-radiating the ultraviolet energy it
absorbed as roughly a blackbody at far-infrared (FIR) wavelengths. To measure
the blackbody is to measure star formation; to measure it across cosmic time
is to characterise cosmic star-formation history, and maybe learn something
about the evolution of the different galaxy populations.

Fundamental questions about the early Universe still remain, such as: was it
dust rich or dust poor? How much does obscured star formation contribute to
the star-formation rate density? Do evolutionary trends of galaxies with
redshift seen locally (e.g. temperature, specific star-formation rate)
continue indefinitely, or do they plateau, or even turn over? Just how active
are AGNs in this period? And how much ionizing radiation is escaping the dark
curtain of dust exactly?

The catch is that measuring the far-infrared spectrum of a galaxy at high
redshift is challenging. A modified blackbody spectrum (or spectral-energy
distribution; SED) with a temperature of 18 K peaks in the infrared at a rest-
frame of roughly 160$\,\rm\mu m$, or between 400–600$\,\rm\mu m$ in the
observer frame if emitted at the peak of star formation, $z\sim 2$–3. To fully
bracket that peak requires combining observations from space telescopes like
Spitzer and Herschel, and ground based observatories like the JCMT, LMT, and
ALMA. FIR space telescopes are wonderful except for one drawback: source
confusion. The primary beams are large enough to contain tens (or hundreds) of
galaxies, complicating optical counterpart identification. Ground-based FIR
telescopes are wonderful except for a different drawback: the atmosphere,
which is to infrared light what dust is to optical light.

Until we are able to put an actively-cooled, $50\,\rm m$ mirror in space, the
solution is stacking. Galaxies with similar physical properties estimated from
optical photometry (e.g. redshift, stellar-mass, star-forming activity) have,
on average, similar infrared properties (e.g. infrared luminosity, dust
temperature; Schreiber et al., 2018). Viero et al. (2013) developed an
algorithm to stack galaxies, grouped by their similarity, simultaneously (see
also Kurczynski & Gawiser, 2010), thus providing a way to overcome the bias
inherent to stacking in highly confused images.

The techniques in this paper are not new. The reason we are only now able to
address these questions is because of newly available catalogues released to
the public from COSMOS (Scoville et al., 2007). We combine the COSMOS2020
catalogue (Weaver et al., 2022) with infrared/submillimetre maps publicly
available through the Herschel Extragalactic Legacy Project (HELP) and SCUBA-2
Cosmology Legacy Survey (S2CLS), and a newly released, python3-compatible
version of simstack, to measure the far-infrared SEDs from redshift 0 to 10
and stellar-masses of $\log(M/\mathrm{M_{\odot}})=9.5$ to 12.

We find that the early Universe was dust-rich, and unusually hot. Effective
dust temperatures increase with redshift, following established relationships
to $z=4$, but then rise more rapidly to $z\sim 10$. Similarly, we find that
the far-infrared luminosity density ($\rho_{\mathrm{LIR}}$) tracks the
established star-formation rate density (SFRD or $\rho_{\mathrm{SFR}}$) to
$z=2$, but then decouples from it, remaining flat to $z=8$. Given the unusual
nature of the hot dust at $z>4$, we explore the mechanisms that could be
responsible for this heating, how this could have been missed previously, and
strategies for follow-up observations.

We assume a Chabrier (2003) initial-mass function (IMF) and Planck18 cosmology
(Planck Collaboration et al., 2020), with $\Omega_{\rm M}=0.315$,
$\Omega_{\Lambda}=0.685$, $H_{0}=67.4\,\rm km\,s^{-1}\,Mpc^{-1}$,
$\sigma_{8}=0.811$, and $\tau=0.054$. All code and instructions to reproduce
these results can be found at https://github.com/marcoviero/simstack3/ and at
the DOI: https://doi.org/10.5281/zenodo.6792395.

##  2 Data

###  2.1 COSMOS2020 Catalogue

Our analysis is performed on data in the
COSMOS111https://cosmos.astro.caltech.edu field (Scoville et al., 2007),
centred at $\rm 10^{h}00^{m}26^{s},2^{\circ}13^{\prime}00^{\prime\prime}$, and
covering 1.606 deg2 after masking. We use the recently released COSMOS2020
catalogue (Weaver et al., 2022), whose $K_{s}=25.2\,\rm(AB)$-selected sample
is an improvement of nearly two magnitudes over previous releases used for
stacking (Muzzin et al., 2013). This extra depth enables access to objects
reaching redshift of 10, and completeness improvements of $\sim 2\,\rm dex$ in
stellar mass.

The COSMOS2020 release comes with two options for photometry; we use
FARMER/LePhare, generally understood to be superior for faint sources at high
redshift. Following Weaver et al. (2022), we split the population into star-
forming and quiescent using the $NUV-r$ versus $R-J$ selection (Ilbert et al.,
2013).

###  2.2 Maps

Maps are selected to bracket the rest-frame dust SED at $z=0$–10. Unless
otherwise noted, they were obtained from the Herschel Extragalactic Legacy
Project (HELP222https://herschel.sussex.ac.uk).

The Spitzer/MIPS map at 24$\,\rm\mu m$ is a Spitzer Enhanced Imaging Product
(SEIP). It has an rms of $27.9\,\rm\mu Jy$, point-spread function (PSF) full-
width at half maximum (FWHM) of 6.32 arcsec, and aperture/colour corrections
equaling 1.24 (Béthermin et al., 2010). Maps are converted from MJy sr-1 to Jy
beam-1 by dividing the map by the beam solid angle, $1.55\times 10^{-9}\,\rm
sr$.

_Herschel_ /PACS maps at 100 and 160$\,\rm\mu m$ are from the PACS
Evolutionary Probe (PEP; Poglitsch et al., 2010). They have rms values of
$0.10$ and $2.11\,\rm mJy$, and FWHM of 7.5 and 11.3 arcsec. PACS maps are
also converted from MJy sr-1 to Jy beam-1, with beam solid angles of
$2.03\times 10^{-9}$ and $4.66\times 10^{-9}\,\rm sr$. Calibration uncertainty
is 4%.

_Herschel_ /SPIRE maps at 250, 350, and $500\,\rm\mu m$ are from the Herschel
Multi-tiered Extragalactic Survey (HerMES; Oliver et al., 2012), with
instrumental depths of 15.9, 13.3, and 19.1 mJy ($5\sigma$); FWHM are 17.6,
23.9, and 35.2 arcsec. Calibration uncertainty is 5%.

The $850\,\rm\mu m$ map is a product of
S2CLS333https://zenodo.org/record/57792 (Geach et al., 2017). Map depth is 1.6
mJy beam-1 ($1\,\sigma$), and FWHM 12.1 arcsec. Note, we specifically select
the non-match-filtered (NMF) maps, which retain the unresolved fluctuations on
which simstack relies.

##  3 Method

###  3.1 Unbiased Stacking with simstack

Our stacking analysis is based on simstack, an algorithm designed to overcome
the biases inherent with stacking on images dominated by confusion noise.
Briefly, simstack uses the 3D positions in a catalogue to construct and fit a
model to a map. The model is made up of layers representing bins that have
been defined by the user (e.g. by stellar-mass, redshift, etc.), and whose
best-fittings are interpreted as the average flux density of the binned
objects. The software for stacking used here is the same one detailed in Viero
et al. (2013) with some minor modifications.

The first difference, following Sun et al. (2018), is that instead of stacking
in redshift slices, now all redshift bins are also stacked simultaneously.
This change is to minimize the potential bias that may arise from objects
outside the bin boundaries, or from objects correlated with interloping
sources.

The next difference is motivated by Duivenvoorden et al. (2020), who showed
that stacked flux densities can be biased by holes in the catalogue coverage
caused by extremely bright objects like stars. Their solution is the addition
of a flat foreground layer, which is now an optional flag (default=True) in
simstack.

The last change concerns the bootstrap estimation. Previously, the bootstrap
consisted of randomly sampling the catalogue (with replacement), which is
robust in most scenarios, but in this case actually biases the fluxes low. The
modification satisfies the central principle of simstack, that all correlated
objects must be stacked at the same time to prevent biasing the estimate, by
randomly splitting each bin in two, with a 80:20 ratio, and stacking the
entire set simultaneously. The 80% bins are retained as the bootstrapped flux,
and the 20% bins are examined for outliers. The drawback with this method is
that the memory required to perform the stack is doubled, requiring either a
more powerful computer, or splitting the bootstrap into chunks.

Figure 1: Best-fitting SEDs of star-forming galaxies, binned by redshift
(increasing left to right) and stellar mass (increasing top to bottom). Data
at 24$\,\rm\mu m$ (blue circles) are from _Spitzer_ /MIPS; at 100 and
160$\,\rm\mu m$ (green circles) from _Herschel_ /PACS; at 250, 350, and
500$\,\rm\mu m$ (red circles) from _Herschel_ /SPIRE; and at 850$\,\rm\mu m$
(yellow circles) from SCUBA2. Error bars represent $1\,\sigma$ Gaussian
uncertainties, derived with the bootstrap estimator, on 150 iterations. Non-
detections are shown at $3\,\sigma$ upper limits. The SED model is a hybrid
blackbody with fixed emissivity index $\beta=1.8$ on the Rayleigh-Jeans side,
and fixed power-law approximation slope $\alpha=2$ on the Wien side. MCMC fits
(cyan) show the 50th percentile line bracketed by a shaded region marking the
25th and 75th percentiles. For comparison we also show hypothetical main-
sequence galaxy SEDs (thin magenta lines), extrapolated from the Schreiber et
al. (2015) model at $z<4$.

###  3.2 Spectral Energy Distribution Fitting

To obtain SEDs from stacked flux densities we fit them with a hybrid
blackbody, substituting the mid-infrared Wien side with a power law, $\alpha$,
with fixed slope of 2.0, and emissivity index $\beta$ fixed at 1.8 (e.g.
Casey, 2012). Fits then have two free parameters, amplitude and observed
temperature. We sample the parameter space with emcee for python (Foreman-
Mackey et al., 2013), a user-friendly, affine-invariant ensemble sampler for
Markov chain Monte Carlo (MCMC). We adopt a flat prior with generous upper and
lower limits for amplitude and temperature. For five bins at $z=1.5$–2 and
$z=2.0$–2.5, we add a prior on temperature to prevent the solution from
preferring local minima at higher temperature. Additionally, we account for
PAH emission lines contaminating the 24$\,\rm\mu m$ flux density by
artificially inflating the MIPS errors over the redshift range $0.5<z<3$. Full
covariances are derived from bootstraps to account for the highly correlated
nature of bands near to one another in wavelength. Flux density measurements
consistent with zero are treated as upper limits, whose contribution to the
log-likelihood is via a survival function (Sawicki, 2012). We use 15 000
samples, discarding the first 3000.

Infrared luminosities (LIR) are estimated as the integral of $\nu I_{\nu}$
from rest-frame 8 to 1000$\,\rm\mu m$ (Casey, 2012), and naïvely converted to
the star-formation rate (SFR) using the Kennicutt (1998) relation, with a 0.23
dex correction to convert to a Chabrier (2003) IMF. Error bars account for
additional uncertainties in the luminosity distances arising from the
photometric redshifts.

We also estimate the SFR via the rest-frame 850$\,\rm\mu m$ luminosity
$L_{850}$. As a proxy for the interstellar medium gas mass, it is insensitive
to dust temperature $T_{\mathrm{d}}$ with a relatively constant
luminosity–mass ratio, and is also a reliable SFR estimator (Scoville et al.,
2016).

###  3.3 Null Tests and Other Consistency Checks

To check that the signal we measure is not noise, we redo the stack but
randomly shuffle source positions. We find that the measured flux densities
are centred around zero, with uncertainties consistent with the rms in the
maps. We further check that the measured signal is not dominated by a few
bright sources by inspecting the distribution of the 20% group of the
bootstrap split, finding their relative variance is proportional to their
relative size.

Next, we carefully examine the impact that low-$z$ interlopers have on the
stacked SED shape at $z>4$ using simulated catalogs/maps generated with the
pysides simulation (Béthermin et al., 2010). We use the redshift probability
distributions $P(z)$ provided by the COSMOS team to i) draw a redshift for
each object; ii) look for the closest counterpart in the sides catalog; iii)
bin and average their flux densities; iv) and compare the best-fit SEDs. We
repeat this 100 times (notebooks and downloads of the simulation results can
be found in the same GitHub repository as for the main results).

Objects at higher redshift or with lower stellar masses have broader $P(z)$,
sometimes even double-peaked; provided they are reliable probabilistic
estimators of the true redshifts, this method should give a reasonable
representation of the outlier bias on the mean SED. We find that the overall
bias in rest-frame dust temperature is small, but that the scatter is large,
particularly at $z>5$. We apply these excess uncertainties to the errors in
the dust temperatures.

Finally, we check that the cumulative flux density is consistent with the
cosmic infrared background (e.g. Dole et al., 2006), and that the
contributions from different redshift and stellar-mass bins agree with Viero
et al. (2013).

##  4 Results

The full set of SEDs for star-forming galaxies are shown in Figure 1. The best
fits are in blue, with blue shaded regions outlining the 25th and 75th
quartiles of the full chain of MCMC samples. Overlaid in magenta are SEDs
drawn from the star-forming main sequence of Schreiber et al. (2015), as
described in § 3.3, showing a clear departure from existing trends at $z=0$–4.

Figure 2: Rest-frame dust temperatures from best-fittings to SEDs plotted as
open circles. A polynomial fit for $T(z)$ is shown as a dashed red line.
Closed circles are the five bins fit with priors, and are not included in the
fit. Positions on the redshift axis correspond to median redshifts of the
objects in the bins. Also plotted are $T_{\mathrm{d}}$–$z$ relationships
(Viero et al., 2013; Schreiber et al., 2018; Bouwens et al., 2020),
extrapolated to $z=9$. Temperatures of individual objects include Hashimoto et
al. (2019), Faisst et al. (2020), Sommovigo et al. (2022), Béthermin et al.
(2020), and Laporte et al. (2017)/Behrens et al. (2018), and lower limit from
Bakx et al. (2020).

###  4.1 Dust Temperature

We find dust temperatures increasing with redshift (Figure 2), in agreement
with previous relationships found to redshift 4 (Viero et al., 2013; Schreiber
et al., 2018; Bouwens et al., 2020). At $z>4$, we find temperatures diverge
from those relationships. Temperatures are incompatible with other estimates
that follow (extrapolated) existing relationships (e.g. Béthermin et al.,
2020; Faisst et al., 2020; Ferrara et al., 2022; Hashimoto et al., 2019;
Sommovigo et al., 2022). They are in better agreement with the lower limit of
Bakx et al. (2020), and estimate of Behrens et al. (2018) based on data from
Laporte et al. (2017); these findings had been considered outliers. We model
the $T_{\mathrm{d}}$–$z$ relationship with a polynomial:

$T_{\mathrm{d}}(z)=(20.9\pm 3.5)+(5.9\pm 1.9)\,z+(0.5\pm 0.2)\,z^{2}.$ (1)

Naïve explanations for these findings (after considering heating by the cosmic
microwave background; da Cunha et al., 2013) are inadequate. Ambient heating
by the elevated cosmic microwave background is subdominant in these objects
(da Cunha et al., 2013). Further, at $z\sim 7$, assuming typical, steady-state
dust grains and an interstellar radiation field with total energy density
$U\gtrsim 1500\,\rm{eV\,cm^{-3}}$, to supply enough energy by star formation
would require an SFR of greater than $1000\,\rm M_{\odot}\,\mathrm{yr^{-1}}$.
Even for a halo as massive (and rare) as $\rm 10^{13}\,M_{\odot}$, this
entails an high, order-of-unity star formation efficiency (Sun & Furlanetto,
2016).

However, the dust composition of galaxies at these redshifts is very different
from local (e.g. Behrens et al., 2018). Multiple dust components (e.g. cold
vs. warm) may be present (Strandet et al., 2017; Liang et al., 2019). As
pointed out by De Rossi et al. (2018), dust in primeval galaxies (Pop II) is
typically younger than $400\,\rm Myr$; silicate rich, being produced by low
metallicity AGB stars (Ventura et al., 2012), with poor emission efficiency at
low and high wavelengths (Koike et al., 2003); and located in environments
near to the heating source (Spilker et al., 2016). AGNs are another likely
significant source of heating (Lyu et al., 2016), whose host galaxies are very
compact (Decarli et al., 2018), with SFRs of $1000-5000\,\rm
M_{\odot}\,yr^{-1}$ (and SFR surface densities reaching $300-2000\,\rm
M_{\odot}\,yr^{-1}\,kpc^{-2}$; Lyu et al., 2016). A fuller treatment must
specifically consider the nature of galaxies in the early Universe.

Figure 3: Star-formation rate density vs. redshift, binned by stellar mass
(coloured lines), and summed (black line with 1 $\sigma$ shaded region). Star-
formation rates are converted from far-infrared luminosities – derived from
the integral of intensities from (rest-frame) 8 to 1000$\,\rm\mu m$– by
scaling with the Kennicutt (1998) relation. We find we are in broad agreement
with a selection of the latest $\rm SFRD$ from models (Béthermin et al., 2017;
Pillepich et al., 2018) and measurements (Zavala et al., 2021). Also plotted
is the rest-frame 850$\,\rm\mu m$-derived star-formation rate density as a
light blue line and shaded region ($1\,\sigma$). The LIR-derived SFRD tracks
850$\,\rm\mu m$-derived relation to $z=4$, but then they diverge, similar to
the dust temperature behaviour, suggesting that the FIR luminosity may not be
a good tracer of star-formation at high redshift.

###  4.2 The LIR and Star-Formation Rate Density

We estimate the far-infrared luminosity density ($\rho_{\mathrm{LIR}}$) as a
number-weighted sum of the infrared luminosities in different stellar-mass
bins derived above,

$\rho_{\mathrm{LIR}}=\frac{1}{V}\displaystyle\sum^{m}_{i=0}L_{i}N_{i},~{}~{}\delta\rho_{\mathrm{LIR}}^{2}=\frac{1}{V^{2}}\displaystyle\sum^{m}_{i=0}(N_{i}\delta
L_{i})^{2}+(L_{i}\delta N_{i})^{2},$ (2)

where $\delta N^{2}_{i}=N^{2}\epsilon^{2}_{CV}$, and $\epsilon_{CV}$ is the
cosmic variance term (Moster et al., 2011). Similarly, we estimate the SFR
density ($\rho_{\mathrm{SFR}}$) from SFRs derived from rest-frame 850$\,\rm\mu
m$ flux densities.

In Figure 3 we show the $\rho_{\mathrm{LIR}}$ as a grey line and shaded
region, which when converted to $\rho_{\mathrm{SFR}}$ is in good agreement
with previous measurements (e.g. Bouwens et al., 2020; Zavala et al., 2021)
and models (Béthermin et al., 2017; Pillepich et al., 2018) to $z=4$, but then
remains elevated, while the models decline. Given our finding of hot dust at
$z>4$ this excess is not unexpected, as we know that hot dust SEDs are poor
estimators of star-formation rates (Bouwens et al., 2016). On the other hand,
the 850$\,\rm\mu m$-derived $\rho_{SFR}$, shown as a blue line and shaded
region, exhibits the expected decline after peaking at $z=1-2$. Note, we
exclude the $z=8-10$ bins, whose galaxy abundances suggest a large
contribution from interlopers.

It is reasonable to ask how so much of the $\rm SFRD$ can be recovered when
the incompleteness levels of the catalogue at $z>4$ must be high. The reason
is because we are stacking on images with large PSFs, which would result in a
bias from galaxies that were not in, but were correlated with, those in the
catalogue.

The high-$z$ SFRD has significant consequences for the evolution of the early
Universe (Robertson, 2021); leveraging statistical techniques like stacking to
access the full unresolved population (see Viero et al., 2015; Cheng et al.,
2021) will complement traditional methods to provide a complete picture of
star-formation activity.

##  5 Discussion

How could so much hot dust in the early Universe have been missed? In
hindsight, all the clues were there: there was a dearth of high-redshift dusty
galaxies (Casey et al., 2018) that hot dust could explain (Sommovigo et al.,
2020). There was known tension in low IRX values observed for given $\beta$ at
high redshift (Capak et al., 2015) that could be alleviated with warmer dust
(Bouwens et al., 2016). There were unusually high dust-mass estimates that
challenged dust production time-scales (Leśniewska & Michałowski, 2019), which
is resolved with LIR coming from hotter dust (Bakx et al., 2020). And then
there were simple models showing that ALMA observations tended to miss
galaxies with high-temperature dust (Chen et al., 2022).

So why has the idea of _hot_ dust at high redshift not gained traction? First,
a handful of hot dusty galaxies at $z>6$ _had_ been observed (e.g. Behrens et
al., 2018; Bakx et al., 2020), but many others were found to be consistent
with extrapolations to previous relationships (e.g. Hashimoto et al., 2019;
Faisst et al., 2017; Sommovigo et al., 2020; Ferrara et al., 2022), so that
the behaviour of the $T_{\rm d}-z$ relationship – increasing or plateauing –
was still unclear (Faisst et al., 2017; Sommovigo et al., 2020). However,
those temperature estimates were derived from measurements with ALMA that did
not actually bracket the peak of the SED (largely staying below Band 8, so no
shorter than 600$\,\rm\mu m$). If in fact these objects were much hotter, no
adequate prior to extrapolate model SEDs even existed.

Alongside hot dust is an elevated infrared-luminosity density at $z>4$, which
appears to decouple from the cold-dust derived SFRD. The source of this excess
heating is likely some combination of young, Pop II stars and silicate rich
dust, and AGN, but in what combination, and is that enough? Such details
impact our understanding of reionization time-scales, or of the presence of
quiescent galaxies at high redshift (Glazebrook et al., 2017).

Understanding the nature of the dust, and the extreme source of heating – is
it stars, AGNs, PBHs? – ranks high on the list of fundamental questions going
forward.

##  6 Conclusion

We estimate far-infrared luminosities and dust temperatures from galaxies at
$z=0$–10, by stacking 111 227 objects from the COSMOS2020/Farmer catalogue on
maps in the far-infrared/submillimetre. We find agreement with previous
measurements to $z=4$, but beyond this find elevated temperatures and
luminosities, indicating an early Universe filled with hot dust. The
temperatures measured suggest a direct detection of Pop II galaxies – compact,
close-packed, and silicate rich – whose SEDs carry the signature of localized
hot dust.

If confirmed, this will have far-reaching consequences. It will inform future
probes into the epoch of reionization, whether directly with the _James Webb
Space Telescope (JWST)_ , or via statistics of spectral lines like CO and [C
ii] with line-intensity mapping experiments like COMAP (Cleary et al., 2022),
TIME (Sun et al., 2021), or CONCERTO (Catalano et al., 2022). In particular, a
warm dust background could present problems for [C ii] experiments by
decreasing the line-continuum contrast.

Clearly, follow-up observations targeting the peak of hot dust emission are in
order. ALMA band 10 will probe nearer to the peak than band 9 observations;
targeting the most massive and reliable subset of the $z>7$ sample should
achieve detections with reasonable signal-to-noise, but still will not fully
bracket the peak.

Long term, dedicated observatories targeting the peak of warm dust SEDs in the
submillimetre will be necessary to complement deep optical/NIR observations.
On the horizon are up-to-50 m class ground-based single-dish observatories
like CCAT-prime/FYST (CCAT-Prime collaboration et al., 2021), AtLAST (Klaassen
et al., 2020), and the Japanese 30 m Antarctic THz telescope intended for New
Dome Fuji (Chen et al., 2022), as well as space missions (e.g. _Origins Space
Telescope_ ; Wiedner et al., 2020).

## Acknowledgements

The authors are grateful to the referee for their careful reading and
constructive suggestions, thank you. DTC holds a CITA/Dunlap Institute
postdoctoral fellowship. The Dunlap Institute is funded through an endowment
established by the David Dunlap family and the University of Toronto. The
University of Toronto operates on the traditional land of the Huron-Wendat,
the Seneca, and most recently, the Mississaugas of the Credit Rive.

## Data Availability

The simstack package needed to reproduce these results, and Jupyter Notebooks
guiding the user through each step, are available at
https://github.com/marcoviero/simstack3/tree/main/viero2022/ and at the DOI:
https://doi.org/10.5281/zenodo.6792395.

## References

  * Bakx et al. (2020) Bakx T. J. L. C., et al., 2020, MNRAS, 493, 4294
  * Behrens et al. (2018) Behrens C., Pallottini A., Ferrara A., Gallerani S., Vallini L., 2018, MNRAS, 477, 552
  * Béthermin et al. (2010) Béthermin M., Dole H., Beelen A., Aussel H., 2010, A&A, 512, A78
  * Béthermin et al. (2017) Béthermin M., et al., 2017, A&A, 607, A89
  * Béthermin et al. (2020) Béthermin M., et al., 2020, A&A, 643, A2
  * Bouwens et al. (2016) Bouwens R. J., et al., 2016, ApJ, 833, 72
  * Bouwens et al. (2020) Bouwens R., et al., 2020, ApJ, 902, 112
  * CCAT-Prime collaboration et al. (2021) CCAT-Prime collaboration et al., 2021, arXiv e-prints, p. arXiv:2107.10364
  * Capak et al. (2015) Capak P. L., et al., 2015, Nature, 522, 455
  * Casey (2012) Casey C. M., 2012, MNRAS, 425, 3094
  * Casey et al. (2018) Casey C. M., Hodge J., Zavala J. A., Spilker J., da Cunha E., Staguhn J., Finkelstein S. L., Drew P., 2018, ApJ, 862, 78
  * Catalano et al. (2022) Catalano A., et al., 2022, in European Physical Journal Web of Conferences. p. 00010 (arXiv:2110.14984), doi:10.1051/epjconf/202225700010
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chen et al. (2022) Chen Y. Y., Hirashita H., Wang W.-H., Nakai N., 2022, MNRAS, 509, 2258
  * Cheng et al. (2021) Cheng Y.-T., et al., 2021, ApJ, 919, 69
  * Cleary et al. (2022) Cleary K. A., et al., 2022, ApJ, 933, 182
  * De Rossi et al. (2018) De Rossi M. E., Rieke G. H., Shivaei I., Bromm V., Lyu J., 2018, ApJ, 869, 4
  * Decarli et al. (2018) Decarli R., et al., 2018, ApJ, 854, 97
  * Dole et al. (2006) Dole H., et al., 2006, A&A, 451, 417
  * Duivenvoorden et al. (2020) Duivenvoorden S., et al., 2020, MNRAS, 491, 1355
  * Faisst et al. (2017) Faisst A. L., et al., 2017, ApJ, 847, 21
  * Faisst et al. (2020) Faisst A. L., Fudamoto Y., Oesch P. A., Scoville N., Riechers D. A., Pavesi R., Capak P., 2020, MNRAS, 498, 4192
  * Ferrara et al. (2022) Ferrara A., et al., 2022, MNRAS, 512, 58
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Geach et al. (2017) Geach J. E., et al., 2017, MNRAS, 465, 1789
  * Glazebrook et al. (2017) Glazebrook K., et al., 2017, Nature, 544, 71
  * Hashimoto et al. (2019) Hashimoto T., et al., 2019, PASJ, 71, 71
  * Ilbert et al. (2013) Ilbert O., et al., 2013, A&A, 556, A55
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ARA&A, 36, 189
  * Klaassen et al. (2020) Klaassen P. D., et al., 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 114452F (arXiv:2011.07974), doi:10.1117/12.2561315
  * Koike et al. (2003) Koike C., Chihara H., Tsuchiyama A., Suto H., Sogawa H., Okuda H., 2003, A&A, 399, 1101
  * Kurczynski & Gawiser (2010) Kurczynski P., Gawiser E., 2010, AJ, 139, 1592
  * Laporte et al. (2017) Laporte N., et al., 2017, ApJ, 837, L21
  * Leśniewska & Michałowski (2019) Leśniewska A., Michałowski M. J., 2019, A&A, 624, L13
  * Liang et al. (2019) Liang L., et al., 2019, MNRAS, 489, 1397
  * Lyu et al. (2016) Lyu J., Rieke G. H., Alberts S., 2016, ApJ, 816, 85
  * Moster et al. (2011) Moster B. P., Somerville R. S., Newman J. A., Rix H.-W., 2011, ApJ, 731, 113
  * Muzzin et al. (2013) Muzzin A., et al., 2013, ApJS, 206, 8
  * Oliver et al. (2012) Oliver S. J., et al., 2012, MNRAS, 424, 1614
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Poglitsch et al. (2010) Poglitsch A., et al., 2010, A&A, 518, L2
  * Robertson (2021) Robertson B. E., 2021, arXiv e-prints, p. arXiv:2110.13160
  * Sawicki (2012) Sawicki M., 2012, PASP, 124, 1208
  * Schreiber et al. (2015) Schreiber C., et al., 2015, A&A, 575, A74
  * Schreiber et al. (2018) Schreiber C., Elbaz D., Pannella M., Ciesla L., Wang T., Franco M., 2018, A&A, 609, A30
  * Scoville et al. (2007) Scoville N., et al., 2007, ApJS, 172, 1
  * Scoville et al. (2016) Scoville N., et al., 2016, ApJ, 820, 83
  * Sommovigo et al. (2020) Sommovigo L., Ferrara A., Pallottini A., Carniani S., Gallerani S., Decataldo D., 2020, MNRAS, 497, 956
  * Sommovigo et al. (2022) Sommovigo L., et al., 2022, MNRAS, 513, 3122
  * Spilker et al. (2016) Spilker J. S., et al., 2016, ApJ, 826, 112
  * Strandet et al. (2017) Strandet M. L., et al., 2017, ApJ, 842, L15
  * Sun & Furlanetto (2016) Sun G., Furlanetto S. R., 2016, MNRAS, 460, 417
  * Sun et al. (2018) Sun G., et al., 2018, ApJ, 856, 107
  * Sun et al. (2021) Sun G., et al., 2021, ApJ, 915, 33
  * Ventura et al. (2012) Ventura P., et al., 2012, MNRAS, 420, 1442
  * Viero et al. (2013) Viero M. P., et al., 2013, ApJ, 779, 32
  * Viero et al. (2015) Viero M. P., et al., 2015, ApJ, 809, L22
  * Weaver et al. (2022) Weaver J. R., et al., 2022, ApJS, 258, 11
  * Wiedner et al. (2020) Wiedner M. C., et al., 2020, arXiv e-prints, p. arXiv:2012.02731
  * Zavala et al. (2021) Zavala J. A., et al., 2021, ApJ, 909, 165
  * da Cunha et al. (2013) da Cunha E., et al., 2013, ApJ, 766, 13