# An Empirical Study of the Occurrence of Heavy-Tails in Training a ReLU Gate

Sayar Karmakar1 1 Department of Statistics, University of Florida,
Gainesville, United States<EMAIL_ADDRESS>Anirbit Mukherjee2 2
Department of Computer Science, University of Manchester, Manchester, United
Kingdom<EMAIL_ADDRESS>
###### Abstract

A particular direction of recent advance about stochastic deep-learning
algorithms has been about uncovering a rather mysterious heavy-tailed nature
of the stationary distribution of these algorithms, even when the data
distribution is not so. Moreover, the heavy-tail index is known to show
interesting dependence on the input dimension of the net, the mini-batch size
and the step size of the algorithm. In this short note, we undertake an
experimental study of this index for S.G.D. while training a
$\mathop{\mathrm{ReLU}}$ gate (in the realizable and in the binary
classification setup) and for a variant of S.G.D. that was proven in [7] for
$\mathop{\mathrm{ReLU}}$ realizable data. From our experiments we conjecture
that these two algorithms have similar heavy-tail behaviour on any data where
the later can be proven to converge. Secondly, we demonstrate that the heavy-
tail index of the late time iterates in this model scenario has strikingly
different properties than either what has been proven for linear hypothesis
classes or what has been previously demonstrated for large nets.

##  I Introduction

The ongoing artificial intelligence revolution [8, 16, 14, 18, 17] can be said
to hinged on developing a vast array of mysterious heuristics which get the
neural net to perform “human like” tasks. Most such successes can be
mathematically seen to be solving the function optimization/“risk
minimization” question, $\inf_{\mathbf{N}\in{\cal
N}}\mathbb{E}_{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
z$}}{\mbox{\boldmath$\textstyle\bf z$}}{\mbox{\boldmath$\scriptstyle\bf
z$}}{\mbox{\boldmath$\scriptscriptstyle\bf z$}}}\in{\cal
D}}[\ell(\mathbf{N},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
z$}}{\mbox{\boldmath$\textstyle\bf z$}}{\mbox{\boldmath$\scriptstyle\bf
z$}}{\mbox{\boldmath$\scriptscriptstyle\bf z$}}})]$ where members of ${\cal
N}$ are continuous piecewise linear functions representable by neural nets and
$\ell:{\cal N}\times{\rm Support}({\cal D})\rightarrow\mathbb{R}^{\geq 0}$ is
called a “loss function” and the algorithm only has sample access to the
distribution ${\cal D}$. The successful neural experiments can be seen as
suggesting that there are certain special choices of $\ell,{\cal N}$ and
${\cal D}$ for which highly accurate solutions of this seemingly extremely
difficult question can be found fast - and often by a surprisingly simple
algorithm, Stochastic Gradient Descent. And this is a profound mathematical
mystery of our times.

An increasingly popular approach to understand the learning abilities of
Stochastic Gradient Descent has been to look at S.D.E.s [15, 20, 21] which can
be motivated to be the continuous time limit of S.G.D.s. A key component of
this modelling relies on the assumption that the noise in the stochastic
gradient is Gaussian in real world learning scenarios. But this assumption has
be challenged in recent times in works like [4] and [19]. But much remains to
be understood about how the tail-index of the noise in the stochastic gradient
is at all determined. Also these S.D.E. based results derive the structure of
the corresponding stationary distribution by making assumptions about the loss
functions which are not verified for neural losses. Hence to make progress it
becomes imperative to be able to devise experiments to find out the structure
of the stationary distribution in stochastic neural training algorithms.

On the other hand in works like [6], albeit motivated by wanting to understand
factors influencing the “width” of the minima found by S.G.D., the authors
therein had isolated the important factors influencing the asymptotic
behaviour of S.G.D. to be the learning rate (say $\eta$) and the mini-batch
size (say $b$). We continue in the same spirit to try to understand the
structure of the distribution of the late time iterates of stochastic
algorithms which can train a $\mathop{\mathrm{ReLU}}$ gate. In particular, we
focus on the same parameters $\eta$ and $b$ and formulate precise conjectures
about how they affect the heavy-tail index of the distribution of asymptotic
iterates.

Via our experiments we will demonstrate that the dependence of this heavy-tail
index on $\eta$ and $b$ for training a $\mathop{\mathrm{ReLU}}$ gate has
significant differences from previously reported trends for both linear
regression as well as large nets. Thus we motivate targets for future
theoretical studies for this model system which can be seen as a sandbox for
developing the right tools of study.

###  I-A Organization

In the subsections I-B and I-C we will setup the formal definitions that we
will use for the heavy-tail index and the specific $\mathop{\mathrm{ReLU}}$
training algorithms that we will use. In Section II we will give detailed
comparison of our approach to previous investigations in the same theme. In
Section III we will explain our experimental methodology and summarize our
findings. In Section IV and V we will describe the specific plots that we will
obtain and then we conclude with suggestions for future work.

###  I-B Definitions

Suppose we consider a probability density
$p({\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}})$ with a power-law tail
decreasing as $\frac{1}{\norm{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}}^{\alpha+1}}$ where
$\alpha\in(0,2)$. Then for such distributions the $\alpha^{th}-$moment exists
for only $\alpha<2$. Random variables of such distributions satisfy a central
limit theorem too [13] and that has historically motivated the following
definition,

######  Definition 1 (Symmetric $\alpha$-Stable Distribution
$\mathcal{S}\alpha\mathcal{S}(\sigma)$)

We write an univariate $X\mathbin{\sim}\mathcal{S}\alpha\mathcal{S}(\sigma)$
if the characteristic function of $X$ is

$\phi_{X}(t)=\exp(-\sigma^{\alpha}\absolutevalue{t}^{\alpha}).$

The $\alpha=2$ special case of the above corresponds to Gaussians and the
$\alpha=1$ is the Cauchy distribution. Other values of $\alpha$ in general do
not yield distributions that have nice closed-form densities. We record the
following elementary fact about the $\mathcal{S}\alpha\mathcal{S}$
distributions.

###### Lemma I.1

If $X\mathbin{\sim}\mathcal{S}\alpha\mathcal{S}(\sigma)$ with $\alpha\neq 2$,
then the variance of $X$ is infinite.

One can further define a class of heavy-tailed distribution that contains the
class of Symmetric $\alpha$-Stable distributions.

######  Definition 2 (Strictly $\alpha$-Stable Distribution)

$\textstyle\bf w$ will be said to be a strictly $\alpha-$stably distributed
random vector if for any $m$ independent copies of it say
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{1},\ldots,{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{m}$ we have,

$\sum_{i=1}^{m}{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{i}\overset{\text{in
distribution}}{=}m^{\frac{1}{\alpha}}{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}.$

###### Lemma I.2

If $X\mathbin{\sim}\mathcal{S}\alpha\mathcal{S}(\sigma)$, then $X$ is Strictly
$\alpha$-Stable.

###### Proof:

Assume $X,X_{1},\ldots,X_{m}$ are i.i.d.
$\mathcal{S}\alpha\mathcal{S}(\sigma)$. Then the characteristic function of
$\sum_{i=1}^{m}X_{i}$ is

$\displaystyle\phi_{\sum_{i=1}^{m}X_{i}}(t)=$
$\displaystyle\prod_{i=1}^{m}\phi_{X_{i}}(t)=\exp(-\sigma^{\alpha}|m^{1/\alpha}t|^{\alpha})$
$\displaystyle=$ $\displaystyle\phi_{m^{1/\alpha}X}(t)$

and thus $\sum_{i=1}^{m}X_{i}\overset{\text{in
distribution}}{=}m^{\frac{1}{\alpha}}X$. ∎

An advantage of using the above larger class of heavy-tailed distributions is
that it directly accommodates vector-valued random variables which is the
context of this paper. Henceforth, we will only use the later definition and
study the estimation of $\alpha$. Interested reader can see [3] for different
definitions and variants of the idea of $\alpha$-stable distributions.

######  Lemma I.3 (Hill Estimator, Corollary 2.4 [12])

Let $\\{X_{i}\\}_{i=1,\ldots,K}$ be a collection of strictly $\alpha-$stably
distributed random variables and $K=K_{1}\times K_{2}$. Define
$Y_{i}=\sum_{j=1}^{K_{1}}X_{j+(i-1)K_{1}},\forall i=1,\ldots,K_{2}$. Now
define the estimator,

$\hat{\frac{1}{\alpha}}\coloneqq\frac{1}{\log
K_{1}}\cdot\left(\frac{1}{K_{2}}\sum_{i=1}^{K_{2}}\log|Y_{i}|-\frac{1}{K}\sum_{i=1}^{K}\log|X_{i}|\right)$

Then $\hat{\frac{1}{\alpha}}$ converges almost surely to $\frac{1}{\alpha}$ as
$K\rightarrow\infty$.

###  I-C The Algorithms That We Will Study

Algorithm 1 Modified mini-batch S.G.D. for training a $\mathop{\mathrm{ReLU}}$
gate with realizable data

Input: Sampling access to a distribution ${\cal
D}_{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}$ on $\mathbb{R}^{n}$
Input: Oracle access to labels $y\in\mathbb{R}$ when queried with some
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}\in\mathbb{R}^{n}$
Input: An arbitrarily chosen starting point of
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{1}\in\mathbb{R}^{n}$ and a
mini-batch size of $b$
For: $t=1,\ldots$

  * •

Sample i.i.d the data vectors
$\\{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}_{t_{1}},\ldots,{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}_{t_{b}}\\}\sim{\cal
D}_{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}}$ and query the oracle with
this set.

  * •

The Oracle replies with
$y_{t_{i}}=\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{*},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}_{t_{i}}\right\rangle)$
$\forall i=1,\ldots,b$.

  * •

Form the gradient (proxy),

${\mathchoice{\mbox{\boldmath$\displaystyle\bf
g$}}{\mbox{\boldmath$\textstyle\bf g$}}{\mbox{\boldmath$\scriptstyle\bf
g$}}{\mbox{\boldmath$\scriptscriptstyle\bf
g$}}}_{t}:=-\frac{1}{b}\sum_{i=1}^{b}\mathds{1}_{\left\\{y_{t_{i}}>0\right\\}}(y_{t_{i}}-{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t}^{\top}{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}_{t_{i}}){\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}_{t_{i}}$

  * •

${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t+1}:={\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t}-\eta{\mathchoice{\mbox{\boldmath$\displaystyle\bf
g$}}{\mbox{\boldmath$\textstyle\bf g$}}{\mbox{\boldmath$\scriptstyle\bf
g$}}{\mbox{\boldmath$\scriptscriptstyle\bf g$}}}_{t}$

For the case of the label of data $\textstyle\bf x$ being generated as
$y=\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{*},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}\right\rangle)$, in [7], it
was proven that under mild distributional assumptions, Algorithm 1 converges
to ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*}$ in linear time. Such
fast convergence for S.G.D. to the global minima on the
$\mathop{\mathrm{ReLU}}$ risk function remains unknown except for Gaussian
data.

We note that the Algorithm 1 that we study would have been the standard S.G.D.
on the risk,
$\mathbb{E}_{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}\sim{\cal
D}_{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}}}\left[\frac{1}{2}\left(\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{*},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}\right\rangle-\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}\right\rangle)\right)^{2}\right]$ if in the definition of
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
g$}}{\mbox{\boldmath$\textstyle\bf g$}}{\mbox{\boldmath$\scriptstyle\bf
g$}}{\mbox{\boldmath$\scriptscriptstyle\bf g$}}}_{t}$ in the algorithm the
indicators $\mathds{1}_{\left\\{y_{t_{i}}>0\right\\}}$ were to be replaced
with
$\mathds{1}_{\left\\{\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}_{t_{i}}\right\rangle>0\right\\}}$.

##  II Comparisons to previous literature on the emergence of heavy-tailed
behaviour in neural training

In a stochastic training algorithm, given a mini-batch gradient
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
g$}}{\mbox{\boldmath$\textstyle\bf g$}}{\mbox{\boldmath$\scriptstyle\bf
g$}}{\mbox{\boldmath$\scriptscriptstyle\bf g$}}}_{t}$, and a full gradient
$\nabla f$, suppose we define the noise vector say
$\textrm{U}={\mathchoice{\mbox{\boldmath$\displaystyle\bf
g$}}{\mbox{\boldmath$\textstyle\bf g$}}{\mbox{\boldmath$\scriptstyle\bf
g$}}{\mbox{\boldmath$\scriptscriptstyle\bf g$}}}_{t}-\nabla
f\in\mathbb{R}^{p}$. Then treating $\textrm{U}=\\{u_{1},\ldots,u_{p}\\}$ as a
set of $p$ scalars, the authors in [19] had estimated the tail index of this
set i.e., they estimated the tail index of the empirical distribution
$\frac{1}{p}\sum_{i=1}^{p}\delta_{u_{i}}(x)$. We posit that this is not
immediately comparable to the project we undertake, that is to measure the
tail index of the iterates of the stochastic algorithm.

The authors in [19] had also concluded from their experiments on fully-
connected networks at depths $2$ and $4$ that the mini-batch size does not
significantly affect the heavy-tail index that they were measuring. In
contrast we will show via our experiments that there are consistent patterns
of change w.r.t the mini-batch size of the heavy-tail index of the average
late time iterates of S.G.D.

[4] trained neural nets using S.G.D. for $10^{4}$ iterations for
$\eta\in[10^{-4},10^{-1}]$ and $b\in[1,10]$. They treat each layer as a
collection of i.i.d. $\alpha-$-stable random variables and measure the tail-
index of each individual layer separately. They observed that, while the
dependence of $\alpha$ on $\frac{\eta}{b}$ differs from layer to layer, in
each layer the measured $\alpha$ decreases with increase in the ratio
$\frac{\eta}{b}$ in both MNIST and CIFAR.

But as opposed to the above experiment we will demonstrate that for a single
$\mathop{\mathrm{ReLU}}$ gate the dependence of $\alpha$ on the mini-batch
size at a fixed step-length can have much more variations. Also in our
experiments with realizable data we scan the mini-batch dependence of $\alpha$
for $b\in[4,700]$ which is a much larger range of as opposed to the tests done
in [4] for only $b\in[1,10]$

The series of papers[9, 10, 11] observed an interesting heavy-tailed feature
of the power law of the Empirical Spectral Density (ESD) of the weight
matrices of deep neural networks. They reported that for larger well trained
nets the heavy-tailed property of this ESD becomes more prominent than for
smaller networks. It is well-known that the limiting density of eigenvalues of
$\textrm{A}^{\top}\textrm{A}$ for a rectangular Gaussian matrix A follows a
Marchenko-Pastur distribution. They report that such is the distribution of
the spectrum of the layer matrices for nets like LeNet5, MLP3 etc. However,
for a diverse array of much larger nets like AlexNet, Inception etc. they see
that the ESD of the layer matrices is more like that of a matrix
$\textrm{A}^{\top}\textrm{A}$ where the entries of A are samples from a heavy-
tailed distribution. They also suggested a multi-phase evolution of the neural
training to motivate reasons for this heavy-tail property. This series of
papers also reported that for decreasing mini-batch sizes from $500$ to $2$
the E.S.D of fully connected layers of MiniAlexNet show more heavy-tailed
behavior and while improving generalization accuracy in tandem.

##  III Summary of our experimental results

We shall denote the data as being “realizable” when the labels are exactly
generated from a $\mathop{\mathrm{ReLU}}$ gate (with unknown weights
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*}$) i.e
$y_{t_{i}}=\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{*},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}}_{t_{i}}\right\rangle)$ in
Algorithm 1. And we recall that whether in Algorithm 1 or in S.G.D these
algorithms are minimizing the $\ell_{2}-$risk. We note that the neural
experiments here as well those in [4] involve a certain leap of faith. It is
because the proof of the Hill estimator being a right measurement of the
heavy-tail index uses the assumption that the samples are coming from a
strictly $\alpha-$stable distribution - a property which we do not rigorously
know to be true for the S.G.D. iterates on a $\mathop{\mathrm{ReLU}}$ gate or
any neural net.

But it can be observed that via Theorem 4.4.15 of [1] it follows that the
iterates of the modified $\mathop{\mathrm{ReLU}}$ training algorithm given in
Algorithm 1 do have a power-law decay property of the tail for its stationary
distribution. And we will demonstrate that all experimental measurements of
$\alpha$ on this (mathematically better justified) modified-S.G.D. in
Algorithm 1 reproduces the standard S.G.D. data within $10\%$!

In the experiments in Section IV when the data is realizable, all measurements
are made by invoking the Hill estimator on the average of the last 1000
centered iterates of S.G.D.s running for 8000-8500 iterations. We do the
centering about the true weight ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$ where the $s$ index is
to denote that the global minima is randomly chosen in every sample of S.G.D.
that we use. We take $S=K_{1}\times K_{2}$ samples of S.G.D. runs (for every
choice of dimension and/or mini-batch size) as needed to match the definition
of the Hill estimator in Lemma I.3. We choose $K_{1}=K_{2}=25$ motivated by
considerations of both accuracy of the Hill estimator as well as trying to
keep the run-times manageable. By the $8000^{th}-$iterate the parameter
recovery errors i.e $\norm{{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t,s}-{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}}$ are less than
$10^{-6}$ (infact a couple of orders of magnitude lower for most experiments),
across all our settings.

Thus, combining the above experimental setup with the context of Algorithm 1
we realize that, the distribution whose heavy-tail index we are approximately
measuring is the distribution of the random vector,
“${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{\infty}({\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{*,s})-{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$”. (We use the same
starting point, the value of ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{1}$, across all samples of
the algorithms.) So the measured value of $\alpha$ is affected by the
algorithm’s parameters like $\eta$ and $b$ as well as how the
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$ was sampled. We note
that this is a feature of the experimental setup of [4] too.

As opposed to the experiments in Section IV, in Section V, we repeat some of
the experiments when the data is not realizable.

From all the experiments, there are $3$ main observations that we present
about how the heavy-tails emerge while stochastically training over a
$\mathop{\mathrm{ReLU}}$ gate in both the scenarios of non-realizable data and
the case of realizable data (with the
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$ sampled from a normal
distribution).

Firstly, we see that the heavy-tail index is low (hence the iterates have a
more heavy-tailed distribution) for low dimensions and it grows with
dimensions. In the realizable setting this growth of $\alpha$ can be
conjectured to be sub-linear in dimension.

Secondly, we note that at any dimension, which is not too low, the tail-index
drops with the mini-batch size. So larger mini-batches result in heavier tails
for the late time iterates of training the $\mathop{\mathrm{ReLU}}$ gate.

Thirdly, we note that wherever tested the properties derived from Algorithm 1
and S.G.D. are essentially the same. Given the simplicity of Algorithm 1
relative to S.G.D., this observation could potentially be a fulcrum for future
theoretical progress in this theme.

##  IV The plots from the regression experiments on realizable data

We present the following specific plots in this section,

  * •

In Figure 1 we study the input dimension dependence of the heavy-tail index
measured on the late time iterates of S.G.D. training a
$\mathop{\mathrm{ReLU}}$ gate (in the realizable) setting when the
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$s are sampled from a
normal distribution. We observe that at higher dimensions the iterates get
less heavy-tailed and the growth of $\alpha$ with dimension can be said to be
sub-linear.

  * •

In Figure 2 we repeat the experiment above for Algorithm 1 and check that the
trends continue to hold for this too.

  * •

In Figure 5 we study the heavy-tail index measured on the late time iterates
of S.G.D. training a $\mathop{\mathrm{ReLU}}$ gate (in the realizable) setting
when the ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$s are sampled from a
normal distribution. We see that beyond a threshold value of the mini-batch
size, the distribution of the iterates keep getting heavier tailed as the
mini-batch size increases. We show that this trend is the same at two
different step-lengths separated by an order of magnitude - but the curvature
of the $\alpha$ vs $b$ plot in the large mini-batch regime goes from concave
to convex as $\eta$ increases.

  * •

In Figure 6 we verify that the general trend as seen for S.G.D. above also
continues to hold for the simpler training Algorithm 1.

The above $4$ experiments with realizable data, were all done for a
$\mathop{\mathrm{ReLU}}$ gate mapping $\mathbb{R}^{100}\to\mathbb{R}$ and we
have checked that the above trends hold for other input dimensions too where
we could check within our computational resources. Also all our figures are
accompanied by a plot of how the parameter recovery error evolves with time in
the corresponding experiment and this acts as a check that our heavy-tail
index measurement happens on the iterates when the convergence is fairly
complete.

##  V The plots from the experiments in a binary classification setup

In here we implement the estimation of $\alpha$ on the average late time
iterates of S.G.D. training a $\mathop{\mathrm{ReLU}}$ gate with
$\ell_{2}-$loss when the labelled data
$({\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}},y)$ is sampled as follows :
firstly, a fair coin is tossed to sample $\textstyle\bf x$ from either of two
pre-chosen isotropic Gaussian distributions of different variances and whose
means are located symmetrically about the origin. Then the corresponding label
is assigned to be $1$ or $0$ dependent on which of the two Gaussians the data
got sampled from. Although we train on the $\ell_{2}-$loss, for this case we
also report the time evolution of the population risk in the classification
loss of the predictor ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf
x$}}}\mapsto\mathds{1}_{\left\\{\mathop{\mathrm{ReLU}}(\left\langle{\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf
w$}}}_{t},{\mathchoice{\mbox{\boldmath$\displaystyle\bf
x$}}{\mbox{\boldmath$\textstyle\bf x$}}{\mbox{\boldmath$\scriptstyle\bf
x$}}{\mbox{\boldmath$\scriptscriptstyle\bf x$}}},)\right\rangle\geq
0\right\\}}$ at the $t^{th}-$iterate. To keep the run-times reasonable, in
these experiments we are forced to work with much smaller values of input
dimension, mini-batch sizes and $K_{1}~{}\&~{}K_{2}$ than in the previous
section.

For the above setup we present the following plots,

  * •

In Figure 3 we demonstrate how the estimate of $\alpha$ increases with the
input dimension of the $\mathop{\mathrm{ReLU}}$ gate used. In these
experiments the different S.G.D. samples are all initialized from the same
point which is chosen from a normal distribution for respective dimensions,
and we choose a mini-batch size of $10$ and a step-length of $0.001$. In this
case the estimator of Lemma I.3 is invoked on samples of the average of the
last few hundred iterates when both the regression error as well as the
classification error have nearly saturated.

  * •

In Figure 4 we repeat the same experiment as above but at a fixed input
dimension of $8$ and a step-length of $0.005$ and demonstrate how the estimate
of $\alpha$ decreases with the mini-batch size used.

Note : We emphasize that the dependency of $\alpha$ on the input dimension as
seen in Figures 1, 2 and 3 and on the mini-batch size $b$ as seen in Figures
5, 6 and 4, both differ from the theorems proven for linear regression as well
as the experiments done on neural nets in [5].

$\eta=0.005$

$\eta=0.01$

Figure 1: Non-decreasing behaviour of the heavy-tail index with input
dimension for training a $\mathop{\mathrm{ReLU}}$ gate in the realizable
setting at different step-lengths $\eta$ using S.G.D.

$\eta=0.005$

$\eta=0.01$

Figure 2: Non-decreasing behaviour of the heavy-tail index with input
dimension for training a $\mathop{\mathrm{ReLU}}$ gate in the realizable
setting at different step-lengths $\eta$ using Algorithm 1. Figure 3: Change
in the heavy-tail index with input dimension for training a
$\mathop{\mathrm{ReLU}}$ gate via S.G.D. with non-realizable data.  Figure 4:
Change in the heavy-tail index with mini-batch size for training a
$\mathop{\mathrm{ReLU}}$ gate mapping $\mathbb{R}^{8}\to\mathbb{R}$ via S.G.D.
with non-realizable data.

$\eta=0.005$

$\eta=0.02$

Figure 5: Change in the heavy-tail index with mini-batch size for training a
$\mathop{\mathrm{ReLU}}$ gate in $100$ dimensions in the realizable setting at
different step-lengths $\eta$ using S.G.D. The index measurement is done via
using the Hill estimator as described in Section III In each sample of S.G.D.
used the generating parameter ${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*}$ was sampled from the
standard Gaussian distribution.

$\eta=0.005$

$\eta=0.02$

Figure 6: Change in the heavy-tail index with mini-batch size for training a
$\mathop{\mathrm{ReLU}}$ gate in $100$ dimensions in the realizable setting at
different step-lengths $\eta$ using Algorithm 1 . The index measurement is
done via using the Hill estimator as described in Section III In each sample
of the run of Algorithm used the generating parameter
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*}$ was sampled from the
standard Gaussian distribution.

##  VI Conclusion and Future work

Firstly we note that the experiments presented here are very computationally
intensive and repeating these on larger nets or too large input dimensions was
beyond our resources, particularly for the case of non-realizable data. Hence
an immediate next step would be to devise ways to increase the range of these
experiments that were presented here - and to cross-check the index
measurements via other estimators [2].

Secondly, in every experiment with realizable data and
${\mathchoice{\mbox{\boldmath$\displaystyle\bf
w$}}{\mbox{\boldmath$\textstyle\bf w$}}{\mbox{\boldmath$\scriptstyle\bf
w$}}{\mbox{\boldmath$\scriptscriptstyle\bf w$}}}_{*,s}$ being sampled from a
Gaussian, we have shown that the behaviour of the heavy-tail index for S.G.D.
while training on the $\ell_{2}-$risk on a $\mathop{\mathrm{ReLU}}$ gate
(recall that $\alpha$ is seen to be non-decreasing w.r.t data dimension &
decreasing w.r.t mini-batch size) \- is very closely reproduced while doing
the same training using Algorithm 1 \- which is simpler since it uses updates
linear in the weights. We note that the structure of Algorithm 1 is more
immediately within the ambit of existing theory of stochastic recursions as
given in say Theorem 4.4.15 of [1]. Hence towards developing a theoretical
understanding of the behaviour of the heavy-tail index, we suggest analyzing
this property on Algorithm 1 as a direction of future research.

## References

  * [1] D. Buraczewski, E. Damek, T. Mikosch, et al.  Stochastic models with power-law tails.  Springer.
  * [2] A. Clauset, C. R. Shalizi, and M. E. Newman.  Power-law distributions in empirical data.  SIAM review, 51(4):661–703, 2009.
  * [3] T. Crilly.  The Mathematical Gazette, 79(486):625–628, 1995.
  * [4] M. Gurbuzbalaban, U. Simsekli, and L. Zhu.  The heavy-tail phenomenon in sgd.  arXiv preprint arXiv:2006.04740, 2020.
  * [5] M. Gurbuzbalaban, U. Simsekli, and L. Zhu.  The heavy-tail phenomenon in sgd, 2020.
  * [6] S. Jastrzkebski, Z. Kenton, D. Arpit, N. Ballas, A. Fischer, Y. Bengio, and A. Storkey.  Three factors influencing minima in sgd.  arXiv preprint arXiv:1711.04623, 2017.
  * [7] S. Karmakar and A. Mukherjee.  Provable training of a relu gate with an iterative non-gradient algorithm.  Neural Networks, 2022.
  * [8] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [9] C. H. Martin and M. W. Mahoney.  Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning.  arXiv preprint arXiv:1810.01075, 2018.
  * [10] C. H. Martin and M. W. Mahoney.  Traditional and heavy-tailed self regularization in neural network models.  arXiv preprint arXiv:1901.08276, 2019.
  * [11] C. H. Martin and M. W. Mahoney.  Heavy-tailed universality predicts trends in test accuracies for very large pre-trained deep neural networks.  In Proceedings of the 2020 SIAM International Conference on Data Mining, pages 505–513. SIAM, 2020.
  * [12] M. Mohammadi, A. Mohammadpour, and H. Ogata.  On estimating the tail index and the spectral measure of multivariate $\alpha$-stable distributions.  Metrika, 78(5):549–561, 2015.
  * [13] J. Nolan.  Stable distributions: models for heavy-tailed data.
  * [14] G. D. Portwood, P. P. Mitra, M. D. Ribeiro, T. M. Nguyen, B. T. Nadiga, J. A. Saenz, M. Chertkov, A. Garg, A. Anandkumar, A. Dengel, et al.  Turbulence forecasting via neural ode.  arXiv preprint arXiv:1911.05180, 2019.
  * [15] M. Raginsky, A. Rakhlin, and M. Telgarsky.  Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis.  In Conference on Learning Theory, pages 1674–1703, 2017.
  * [16] T. J. Sejnowski.  The unreasonable effectiveness of deep learning in artificial intelligence.  Proceedings of the National Academy of Sciences, 117(48):30033–30038, 2020.
  * [17] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, et al.  A general reinforcement learning algorithm that masters chess, shogi, and go through self-play.  Science, 362(6419):1140–1144, 2018.
  * [18] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, et al.  Mastering the game of go without human knowledge.  nature, 550(7676):354–359, 2017.
  * [19] U. Simsekli, L. Sagun, and M. Gurbuzbalaban.  A tail-index analysis of stochastic gradient noise in deep neural networks.  In Proceedings of the 36th International Conference on Machine Learning, (ICML) 2019, 2019.
  * [20] P. Xu, J. Chen, D. Zou, and Q. Gu.  Global convergence of langevin dynamics based algorithms for nonconvex optimization.  In Advances in Neural Information Processing Systems, pages 3122–3133, 2018.
  * [21] Y. Zhang, P. Liang, and M. Charikar.  A hitting time analysis of stochastic gradient langevin dynamics.  Proceedings of Machine Learning Research vol, 65:1–43, 2017.

##  Appendix A Code source

  * •

https://colab.research.google.com/drive/1hH-p0CQo3LgvTD4Eq7f0dWeAPenijStE?usp=sharing

  * •

https://colab.research.google.com/drive/1y0S21IXcWpAmX2PzYxMTWVnAePDIQgUH?usp=sharing

  * •

https://colab.research.google.com/drive/1G2PR4k49fTTjkarBgjZ8l93886ovAv1n?usp=sharing


11institutetext:  CAS Key Laboratory of Planetary Sciences, Purple Mountain
Observatory, Chinese Academy of Sciences, No. 10 Yuanhua Road, Qixia District,
210023 Nanjing, PR China
11email<EMAIL_ADDRESS>22institutetext: School of Astronomy and Space
Science, University of Science and Technology of China, No. 96 Jinzhai Road,
Baohe District, 230026 Hefei, PR China 33institutetext: CAS Center for
Excellence in Comparative Planetology, No. 96 Jinzhai Road, Baohe District,
230026 Hefei, PR China 44institutetext: Instituto de Astrofísica de Canarias
(IAC), Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain 55institutetext:
Departamento de Astrofísica, Universidad de La Laguna (ULL), C/ Padre Herrera,
38206 La Laguna, Tenerife, Spain

# Confirmation of TiO absorption and tentative detection of MgH and CrH in the
atmosphere of HAT-P-41b

C. Jiang 11 2 2 G. Chen 11 3 3 F. Murgas 4455 E. Pallé 4455 H. Parviainen 4455
Y. Ma 11

(Received … Accpected …)

Understanding the role of optical absorbers is critical for linking the
properties of the day-side and terminator atmospheres of hot Jupiters. This
study aims to identify the signatures of optical absorbers in the atmosphere
of the hot Jupiter HAT-P-41b. We conducted five transit observations of this
planet to obtain its optical transmission spectra using the Gran Telescopio
Canarias (GTC). We performed atmospheric retrievals assuming free abundances
of 12 chemical species. Our Bayesian model comparisons revealed strong
evidence for TiO absorption ($\Delta\ln\mathcal{Z}=21.02$), modest evidence
for CrH ($\Delta\ln\mathcal{Z}=3.73$), and weak evidence for MgH
($\Delta\ln\mathcal{Z}=2.32$). When we combined the GTC transmission spectrum
with previously published Hubble Space Telescope (HST) and Spitzer data, the
retrieval results and model inferences remained consistent. In conclusion,
HAT-P-41b has a metal-rich atmosphere with no high-altitude clouds or hazes.
Further observations of its day-side atmosphere should be made to confirm the
hints of a thermal inversion in the upper atmosphere suggested by our results.

###### Key Words.:

Planets and satellites: individuals: HAT-P-41b – Planets and satellites:
atmospheres – Techniques: spectroscopic – Methods: data analysis

##  1 Introduction

Accurate characterization of exoplanet atmospheres is crucial in solving the
puzzle of planetary formation and evolution. Precise spectroscopic
measurements of exoplanet atmospheres, available with current observational
instruments, are mainly focused on close-in giant planets, including hot
Jupiters. Certain metals, metal oxides and hydrides have very strong opacity
in the optical wavelengths (Sharp & Burrows 2007), which can absorb the
stellar radiation and heat the upper atmosphere, leading to a thermal
inversion in the temperature-pressure (T-P) profile and thermal emission
features in the emission spectrum (Burrows et al. 2008; Guillot 2010). While
it is difficult to detect optical absorbers in the emission spectrum, alkali
metals (e.g., Na and K), metal oxides (e.g., TiO and VO), and hydrides (e.g.,
CrH and FeH) are expected to exhibit strong absorption features in the optical
band of transmission spectrum (Seager & Sasselov 2000; Fortney et al. 2008,
2010).

The hot Jupiter HAT-P-41b, discovered by Hartman et al. (2012), is an
intriguing target for further study due to its highly inflated nature. The
planet orbits an F-type star ($T_{\rm eff}=6390\leavevmode\nobreak\ {\rm K}$,
$M_{\star}\sim 1.4M_{\odot}$, $R_{\star}\sim 1.7\leavevmode\nobreak\
R_{\odot}$, $\rm[Fe/H]=0.21$) every 2.7 days and has a mass of
$\sim$$0.8\leavevmode\nobreak\ M_{\rm J}$, a radius of
$\sim$$1.7\leavevmode\nobreak\ R_{\rm J}$, and an equilibrium temperature of
$1941\leavevmode\nobreak\ {\rm K}$. Wakeford et al. (2020) presented a 0.2 –
0.8 $\rm\mu m$ transmission spectrum using the G280 grism of HST WFC3 UVIS,
which revealed a significant scattering slope and strong absorption features
of TiO and/or VO. Lewis et al. (2020) performed retrievals on the 0.2–5.0
$\mu$m combined spectrum, consisting of HST WFC3 (G280, G141) and Spitzer, and
found a potential combination of H-, CrH, AlO, VO, and Na to explain the
ultraviolet to optical wavelengths. However, the presence of $\rm H^{-}$ in
the planet’s terminator atmosphere remains inconclusive (Welbanks et al.
2023). Sheppard et al. (2021) conducted a comprehensive retrieval analysis of
the 0.3–5.0 $\mu$m spectrum derived from HST STIS, HST WFC3, and Spitzer
transit observations and found that the atmosphere of HAT-P-41b is cloud-free,
haze-free, and metal-rich, with distinct water vapor absorption and tentative
optical absorption. The day-side upper atmosphere of HAT-P-41b does not
exhibit strong evidence for a thermal inversion, as suggested by the emission
spectrum analyses (Mansfield et al. 2021; Fu et al. 2022; Changeat et al.
2022). The observed transmission spectrum suggest that optical absorbers such
as metal oxides and hydrides may be present in the planet’s terminator
atmosphere, but the emission spectrum did not reveal a distinct thermal
inversion in the day-side atmosphere. Further confirmation and identification
of the candidate optical absorbers in the atmosphere of HAT-P-41b is needed to
resolve this discrepancy.

Our goal is to further investigate the atmosphere of HAT-P-41b using low-
resolution transmission spectrum. Multiple transit observations with the GTC
OSIRIS instrument (Cepa et al. 2000) should enable us to obtain the optical
transmission spectrum with higher photometric precision and spectral
resolution than the HST WFC3 UVIS and STIS instruments. By performing
atmospheric retrievals and Bayesian model comparisons, we aim to obtain more
precise estimates of the abundances of the chemical species, better constrain
previously reported optical absorbers, and possibly discover new signatures in
the atmosphere of HAT-P-41b.

This paper is organized as follows. We summarize the observations and the data
reduction procedures in the next section. Section 3 illustrates the analysis
of the GTC transit light curves. The atmospheric retrievals are demonstrated
in Sect. 4. We further discuss the joint retrievals and then draw our
conclusion in Sect. 5.

##  2 Observations and data reduction

We observed five transits of HAT-P-41b using the GTC OSIRIS instrument between
2013 and 2022 (hereafter OB-1 to OB-5). The observation summary is listed in
Table 4. The instrument OSIRIS has an unvignetted field of view of
$7.4^{\prime}$ and a pixel scale of $0.254^{\prime\prime}$. Its detector
consists of a mosaic of two Marconi CCDs ($1024\times 2048$ pixels) with a
37-pixel ($9.4^{\prime\prime}$) gap between them after a $2\times 2$ pixel
binning. The observations were performed in the longslit spectroscopic mode
with the R1000R grism covering a wavelength range of 5100 – 10 000 Å. The gap
between the CCDs was parallel to the dispersion direction and thus would not
affect the spectra. We adopted a $12^{\prime\prime}$-wide slit in OB-1, OB-3,
and OB-4, and a $40^{\prime\prime}$-wide slit in OB-2 and OB-5 to avoid flux
loss due to the partial coverage of the point spread function (PSF). A
comparison star (2MASS 19493265$+$0439270 for OB-1, OB-3, OB-4, and OB-5;
2MASS 19494153$+$0437350 for OB-2; Cutri et al. 2003) was simultaneously
observed with the target star through the slit in the five observations.

The data were reduced following the methods described in our previous work
(Jiang et al. 2022) to produce the spectroscopic light curves with 5-nm
passbands from 521.8 nm to 922.4 nm, where the passband of 756.8 – 767.4 nm
covering the telluric oxygen A-band was discarded in subsequent analysis due
to the very-low signal-to-noise ratio. The broadband light curves (i.e.,
white-light curves) were obtained in the same manner but integrated from 521.8
nm to 922.4 nm excluding the oxygen A-band.

##  3 Light curve analysis

###  3.1 Transit model

We employed the Python package PyTransit (Parviainen 2015) to model the
transit light curves, assuming circular orbits. There are five transit
parameters in the model: the planet-to-star radius ratio ($R_{\rm
p}/R_{{}_{s}}$), normalized orbital semimajor axis ($a/R_{\star}$), orbital
inclination ($i$), central transit epoch ($t_{\rm c}$), and orbital period
($P$). For each transit, the light curves were fit separately.

For broadband light curve fitting, we employed a four-parameter nonlinear
stellar limb darkening law. The limb darkening coefficients ($u_{1}$ to
$u_{4}$) were derived using LDTK (Parviainen & Aigrain 2015) and kept fixed in
light curve fitting. The orbital period was determined based on the transit
ephemeris updated by Wakeford et al. (2020), with the central epoch of each
transit being a free parameter under a uniform prior. To avoid systematic
biases in the transmission spectra and to facilitate subsequent joint
atmospheric retrievals with the transmission spectra provided by Sheppard et
al. (2021), we adopted the orbital semimajor axis and inclination values
consistent with those used in Sheppard et al. (2021). Furthermore, we used the
broadband radius ratio measured with HST STIS G750L ($R_{\rm p}/R_{\rm
s}=0.10159\pm 0.00042$, Sheppard et al. 2021), which essentially matches the
wavelength coverage of OSIRIS R1000R, as the Gaussian prior for the broadband
radius ratio. This was done to ensure consistency across the five transit
observations and to obtain precise central transit epochs and broadband
systematic noise. Table 1 shows the adopted parameters of the HAT-P-41 system.

We performed separate fittings for the narrowband light curve in each
passband. The parameter settings for the narrowband fitting were kept
consistent with those of the broadband, with the exception of $R_{\rm p}$ and
$T_{\rm c}$. The radius ratio was treated as a free parameter under a uniform
prior, while the central epoch of each transit was fixed to the best-fit value
derived from the broadband light curve fitting. To mitigate the impact of
wavelength-independent systematic noise (i.e., common-mode systematics) and
enhance the precision of the derived transmission spectra, we employed the
common-mode correction method outlined in Gibson et al. (2013). This involved
dividing the best-fit broadband light curve by the best-fit transit model to
obtain the common-mode systematics. Each narrowband light curve was then
divided by the common-mode systematics to apply this correction. We note that
Gibson et al. (2013) used the GP systematics plus one to approximate the
common mode systematics, while we used the best-fit broad-band light curve
divided by the best-fit transit model. Given that the amplitudes of the GP
systematics were much smaller than one, these two approximations are virtually
indistinguishable.

Table 1: Input parameters of the HAT-P-41 system used in this work.

Symbols | Parameters | Values | References
---|---|---|---
Stellar parameters |  |  |
$T_{\rm eff}$ | Stellar effective temperature (K) | $6390\pm 100$ | [1]
$\rm[Fe/H]$ | Stellar metallicity | $0.21\pm 0.10$ | [1]
$\log_{10}g_{\star}$ | Stellar gravity (cgs) | $4.14\pm 0.02$ | [1]
$M_{\star}$ | Stellar mass ($M_{\sun}$) | $1.418\pm 0.047$ | [1]
$R_{\star}$ | Stellar radius ($R_{\sun}$) | $1.683^{+0.058}_{-0.036}$ | [1]
Planetary parameters |  |  |
$M_{\rm p}$ | Planetary mass ($M_{\rm J}$) | $0.800\pm 0.102$ | [1]
$R_{\rm p}$ | Planetary radius ($R_{\rm J}$) | $1.685\pm 0.076$ | [1]
$T_{\rm eq}$ | Planetary equilibrium temperature (K) | $1941\pm 38$ | [1]
$P$ | Orbital period (day) | 2.69404861 | [2]
$t_{0}$ | Initial transit epoch ($\rm BJD_{TDB}$) | 2456600.29325 | [2]
$a/R_{\star}$ | Orbital semi-major axis relative to the stellar radius | $5.44\pm 0.15$ | [1]
$i$ | Orbital inclination (deg) | $87.7\pm 1.0$ | [1]
$R_{\rm p}/R_{\rm s}$ | Planet-to-star radius ratio (HST STIS G750L) | $0.10159\pm 0.00042$ | [3]

[1] Hartman et al. (2012) [2] Wakeford et al. (2020) [3] Sheppard et al.
(2021)

###  3.2 Correction for the companion star contamination

Hartman et al. (2012) reported a resolved neighboring star, located
$\sim$3.56″away and 3.5 mag fainter than HAT-P-41 in the $i$ band. Wöllert et
al. (2015) further observed, through lucky imaging, that this neighboring star
has magnitude contrasts of $\Delta i^{\prime}=3.65\pm 0.05$ and $\Delta
z^{\prime}=3.40\pm 0.05$, and a spectral type between F5 and M2 assuming
physical companionship. According to the Gaia observations (Table 5; Gaia
Collaboration et al. 2016, 2023), these two stars share similar parallaxes,
radial velocities, and proper motions, suggesting their association.

To correct the flux contamination from the companion star, we performed a
supplementary observation with 10 exposures during OB-3, where the slit was
aligned to the direction from HAT-P-41 to the companion star, so that the
spectra of the two stars could be spatially resolved. The distance between the
PSF peaks of the two stars is $\sim$3.5″(14 pixels). To separate the two
stellar spectra, we first subtracted the mirrored non-contaminated side of the
target’s PSF from the other side to obtain the companion’s spectrum. Since
strong residuals existed at the peak location of the target, we fitted a
Gaussian function to target’s original PSF and made a mask for regions within
4$\sigma$ from the peak location. We then repaired the masked region by
mirroring the non-contaminated part of the companion’s PSF on the other side.
The target’s spectrum was then obtained by subtracting the extracted
companion’s spectrum from the total spectrum. The companion-to-host flux ratio
in each passband was then calculated based on the extracted stellar spectra.
The broadband flux ratio (521.8 – 922.4 nm excluding the oxygen A band) was
measured to be $0.03885\pm 0.00083$. The narrowband flux ratios are shown in
Fig. 1.

When fitting the light curve in each passband, the flux dilution effect was
corrected by

$f(t)=\frac{f^{*}(t)+\mathcal{F}}{1+\mathcal{F}},$ (1)

where $\mathcal{F}$ is the companion-to-host flux ratio and $f^{*}(t)$ is the
light curve model without flux dilution.

Figure 1: Narrowband flux ratios between the companion and HAT-P-41. The
errorbars are the observed flux ratios, while the orange line is the best-fit
curve using the PHOENIX model spectra.

We further refined the constraints on the physical parameters of the companion
star by fitting the observed flux ratios with PHOENIX model spectra (Husser et
al. 2013) based on the stellar parameters from Hartman et al. (2012) and
parallaxes measurements from the Gaia data (Gaia Collaboration et al. 2016,
2023). This method was detailed in Nortmann et al. (2016) and Jiang et al.
(2023). Figure 1 shows the fitting results of the flux ratios. The derived
effective temperature of $4498^{+28}_{-33}$ K and radius of $0.708\pm 0.022$
$R_{\odot}$ suggest a K-type companion star. However, to accurately determine
its gravity, mass, and spectral type, further high-resolution spectroscopic
measurements are necessary.

###  3.3 Gaussian process modeling

The light curve systematic noise was fitted with the multi-input Gaussian
process (GP) regression implemented by the Python package George (Ambikasaran
et al. 2015). The GP models were applied to all the broadband and narrowband
light curves individually. We adopted the 3/2-order Matérn kernel function
$k(r)$ with multi-input variables (time $t$, seeing111Measured with the full
width at half maximum of the stellar PSF fitted by a Gaussian function along
the spatial direction at the central wavelength of $\sim$750 nm. $w$, and
pixel drift222 Measured with the pixel drift of the stellar PSF fitted by a
Gaussian function along the spatial direction at the central wavelength of
$\sim$750 nm. $p$) for GP modeling:

$\begin{split}k(r)&=\left(1+\sqrt{3r^{2}}\right)\exp{\left(-\sqrt{3r^{2}}\right)},\\\
k_{\rm total}(\bm{r})&=\sigma_{k}^{2}\cdot k\left(\frac{\Delta
t}{\ell_{t}}\right)\cdot k\left(\frac{\Delta w}{\ell_{w}}\right)\cdot
k\left(\frac{\Delta p}{\ell_{p}}\right),\end{split}$ (2)

where $\sigma_{k}^{2}$ is the variance of the kernel, $\bm{r}$ is the distance
between two input vectors $(t,w,p)$ and $(t+\Delta t,w+\Delta w,p+\Delta p)$
normalized by the corresponding length scales $\ell_{t}$, $\ell_{w}$, and
$\ell_{p}$. In addition, the uncertainties of the normalized light curves were
estimated as photon-dominated noise and usually underestimated the total white
noise. To compensate for this, a jitter variance of $\sigma_{n}^{2}$ was added
to the diagonal of the covariance matrix ($\bm{r}=\bm{0}$). Therefore, the
free GP parameters in GP regression are $\sigma_{n}$, $\sigma_{k}$,
$\ell_{t}$, $\ell_{w}$, and $\ell_{p}$ with log-uniform priors.

###  3.4 Fitting results

We used the nested sampling algorithm implemented by PyMultiNest (Buchner et
al. 2014) to obtain posterior estimates and model evidence (see also Feroz &
Hobson 2008 and Feroz et al. 2009 for more information on the MULTINEST
algorithm). We used 1000 live points and a sampling efficiency of 0.3 in each
run to reach a log-evidence precision of $\sim$0.1 when fitting transit light
curves. The fitting results of the broadband light curves are shown in Fig. 6,
while those of the spectroscopic light curve can be found at
ScienceDB333https://doi.org/10.57760/sciencedb.11105. The posterior estimates
of the central transit epochs are listed in Table 2. As shown in Fig. 2, the
five transmission spectra observed with GTC OSIRIS showed good consistency in
most wavebands, while discrepancies in the bluer and the redder ends are due
to very strong systematics in the spectroscopic light curves, which were
compensated by the very large uncertainties. While the removal of the common-
mode noise helps to reduce the amplitudes of light curve systematics and
increase the precision of spectroscopic transit depths across most wavebands,
the residual wavelength-dependent systematics could still result in the
wavelength dependency of transit depth uncertainties. Given the low activity
of the host star (Hartman et al. 2012; Sheppard et al. 2021), the subsequent
retrieval analysis was performed using the weighted average of the five
transmission spectra.

Table 2: Posterior estimates of the central transit time of HAT-P-41b. # | Date | $T_{\rm c}\leavevmode\nobreak\ ({\rm BJD_{TDB}})$
---|---|---
OB-1 | 2013-06-21 | $2456465.59147^{+0.00031}_{-0.00034}$
OB-2 | 2014-07-14 | $2456853.53273^{+0.00119}_{-0.00143}$
OB-3 | 2016-07-05 | $2457575.53814^{+0.00027}_{-0.00030}$
OB-4 | 2016-07-13 | $2457583.61970^{+0.00041}_{-0.00038}$
OB-5 | 2022-08-11 | $2459803.51788^{+0.00061}_{-0.00061}$
Figure 2: Transmission spectra of HAT-P-41b derived from the five transit
observations with OSIRIS R1000R. Upper panel: Five individual transmission
spectra (colored error bars) and their weighted means (black error bars).
Lower panel: Same as the upper panel, but zoomed in for clarity.

##  4 Atmospheric retrievals

###  4.1 Model setup

We used the Python package petitRADTRANS (Mollière et al. 2019) to model the
transmission spectrum of HAT-P-41b. A one-dimensional isothermal atmosphere
was assumed with a pressure range from $10^{-6}$ to $10^{2}$ bar and 100
uniform layers in logarithmic space. We included 12 chemical species (Na, K,
MgO, AlO, TiO, VO, MgH, AlH, CaH, CrH, FeH, and $\rm H_{2}O$) that could show
significant spectral features in the OSIRIS wavelength range, and four
additional species ($\rm NH_{3}$, $\rm CH_{4}$, CO, and $\rm CO_{2}$) when
retrieving the optical-to-near-infrared (NIR) joint transmission spectrum. We
assumed free mass fractions $m_{i}$ for each chemical species $i$ and used
$\rm H_{2}$ and $\rm He$ as filler gases to keep the total mass fraction unity
with a fixed $\rm He/H_{2}$ mass ratio of 0.3. Rayleigh-like scattering and
collision-induced absorption of $\rm H_{2}$ and He were also included in the
model, with the amplitude of the scattering features multiplied by a factor of
$A_{\rm s}$. Additionally, we assumed a uniform opaque cloud cover with a
cloud top pressure of $P_{\rm c}$, below which all transmitted light would be
blocked. The forward model was originally computed with a spectral resolution
($\lambda/\Delta\lambda$) of 1000 then wavelength binned to the data
passbands.

To infer the possible species from the spectral features, we removed one
chemical species at a time from the free chemistry model and computed the
difference of the log-evidence with and without that species, which is
equivalent to the logarithmic Bayes factor. We interpret the Bayes factors
using the criteria proposed by Kass & Raftery (1995), where the evidence is
strong if $|\Delta\ln\mathcal{Z}|\geq 5$, moderate if
$3\leq|\Delta\ln\mathcal{Z}|<5$, weak if $1\leq|\Delta\ln\mathcal{Z}|<3$, and
inconclusive if $|\Delta\ln\mathcal{Z}|<1$. The Bayesian evidence and
posterior estimates were obtained using the nested sampling algorithm in the
same way as introduced in Sect. 3.4, but sampled with fewer live points
($N=500$) due to the very high computational cost, which still guaranteed
$\sim$150 000 likelihood evaluations and $\sim$3500 posterior samples.

###  4.2 OSIRIS transmission spectrum

Table 3: Bayesian model comparison of the chemical species retrieved with the
OSIRIS transmission spectrum.

Species |  $\Delta\ln\mathcal{Z}$ aa$a$Logarithmic Bayesian evidence ($\ln\mathcal{Z}$) of the full retrieval subtracting that of the hypothesis excluding one chemical species, where the uncertainty of $\ln\mathcal{Z}$ is $\sim$0.15. The values with boldface indicate —$\Delta\ln\mathcal{Z}|\geq 1$. | Inferences
---|---|---
Na | 0.11 | inconclusive
K | -0.46 | inconclusive
MgO | 0.22 | inconclusive
AlO | 0.37 | inconclusive
TiO | 21.02 | strong
VO | -0.58 | inconclusive
MgH | 2.32 | weak
AlH | -0.34 | inconclusive
CaH | 0.09 | inconclusive
CrH | 3.73 | moderate
FeH | -0.4 | inconclusive
H2O | 0.15 | inconclusive
444

Figure 3: Retrieved transmission spectrum of HAT-P-41b observed with GTC
OSIRIS. The error bars are the averaged OSIRIS transmission spectrum. The
orange lines are the median and $2\sigma$ credible intervals of the posterior
models. The black dots are the posterior model after wavelength binning to the
observational passbands. The shaded regions below the posterior model indicate
the reference models with only TiO (red), MgH (blue), or CrH (purple),
respectively.

We first computed the full atmospheric retrieval including all 12 chemical
species (Fig. 3). The posterior joint distributions of the atmospheric
parameters can be found at
ScienceDB555https://doi.org/10.57760/sciencedb.11105. The retrieved model
suggested a clear atmosphere ($\log_{10}P_{\rm c}=-0.74^{+1.64}_{-1.53}$
$\log_{10}{\rm bar}$) with a temperature of $2447^{+299}_{-310}$ K, higher
than the planetary equilibrium temperature ($1941\pm 38\leavevmode\nobreak\
{\rm K}$, Hartman et al. 2012). The Bayesian evidence of the full retrieval
($\ln\mathcal{Z}=501.29\pm 0.16$) was much stronger than that of a flat model
($\ln\mathcal{Z}=402.18\pm 0.12$) or a pure scattering model
($\ln\mathcal{Z}=467.70\pm 0.14$), suggesting significant spectral features in
the OSIRIS transmission spectrum. These spectral features can mostly be
characterized by the absorption signatures of TiO, while near the blue and red
ends some excess absorption could be characterized by the absorption of MgH
and CrH, respectively.

We then performed model comparison based on Bayesian evidence to validate the
detection of each chemical species. We removed each chemical species from the
full retrieval and estimated the Bayesian evidence for a model hypothesis
excluding that species. If the new Bayesian evidence is substantially less
than the Bayesian evidence of the full retrieval, it indicates that the
absorption features of the missing species cannot be counteracted by other
species, implying the need for that species to characterize the transmission
spectrum. According to Table 3, there is very strong evidence for the presence
of spectral features of TiO ($\Delta\ln\mathcal{Z}=21.02$), moderate evidence
for that of CrH ($\Delta\ln\mathcal{Z}=3.73$), and weak evidence for that of
MgH ($\Delta\ln\mathcal{Z}=2.32$) in the OSIRIS transmission spectrum, which
can be approximated as a detection significance of 6.8$\sigma$ for TiO,
3.2$\sigma$ for CrH, and 2.7$\sigma$ for MgH (Sellke et al. 2001; Benneke &
Seager 2013). We found no evidence for the absorption signatures of the alkali
metals or other metal oxides and hydrides.

The absorption features of TiO in the atmosphere of HAT-P-41b have been
previously detected by HST STIS measurements (Sheppard et al. 2021), while
those of CrH and MgH were discovered for the first time in this work. Sheppard
et al. (2021) compared the results of two transmission spectrum models: the
equilibrium chemistry models of PLATON and the free chemistry models of AURA.
They found suggestive evidence for the optical absorbers and decisive evidence
for water vapor in the 0.3 – 5.0 $\mu$m transmission spectrum of HAT-P-41b.
While the free chemistry retrieval from AURA supports the presence of
absorption features of Na I and AlO in the optical band, the results from
PLATON favor the absorption of TiO and VO. However, we did not detect the
absorption features of Na I and AlO in the OSIRIS transmission spectrum using
the petitRADTRANS free chemistry model. Instead, we observed significantly
higher chemical abundances of TiO compared to VO, which is different from the
predictions of the equilibrium chemical model. Additionally, we found
preliminary evidence for MgH and CrH, which are not included in the
equilibrium chemistry model.

##  5 Discussion and conclusion

Figure 4: Retrieval results of HAT-P-41b using the optical-to-NIR joint
transmission spectra. Upper panel: The error bars are the transmission spectra
measured with different instruments, where the median instrumental offsets
have been applied ($-4$ ppm for STIS, 758 ppm for G141, and 388 ppm for IRAC).
The median and $2\sigma$ credible interval of the posterior models are shown
in orange lines. The black dots are the median model after wavelength binning.
The shaded regions below the posterior model indicate the reference models
with only AlO (green), TiO (red), MgH (deep blue), CrH (purple), $\rm H_{2}O$
(sky blue), and $\rm CO_{2}$ (brown) respectively. Lower panel: Same as the
upper panel, but zoomed into the wavelength range of the OSIRIS transmission
spectrum. Figure 5: Posterior distributions of the VMRs of chemical species
constrained by the OSIRIS transmission spectrum (blue) and the optical-to-NIR
joint spectrum (red). The symbol $X$ denotes the logarithmic volume mixing
ratio $\rm\log_{10}VMR$ with the corresponding colors. The gray lines indicate
the VMRs calculated using GGChem equilibrium chemistry (Woitke et al. 2018)
assuming solar abundances at a pressure of 1 bar, where two planetary
equilibrium temperatures were considered: 1941 K (dashed, assuming uniform
heat redistribution) and 2308 K (dotted, assuming dayside heat redistribution
only).

To precisely constrain the atmospheric chemical abundances, we performed a
joint retrieval by combining the OSIRIS transmission spectrum with the 0.3 –
5.0 $\mu$m transmission spectrum presented in Sheppard et al. (2021). To
account for possible transit depth offsets due to different instruments and/or
light curve analysis methods, three arbitrary offsets were added as free
parameters with a uniform prior of $\mathcal{U}(-1000,1000)$ ppm to the HST
STIS, HST WFC3 and Spitzer IRAC transmission spectra, separately. As shown in
Fig. 4, the retrieved model from the joint transmission spectra exhibits
strong absorption signatures of TiO, $\rm H_{2}O$, and $\rm CO_{2}$, and
tentative signatures of AlO, MgH, and CrH, which are in good agreement with
the results from the OSIRIS transmission spectrum in the optical wavelength
range.

Figure 5 presents the posterior volume mixing ratios (VMRs666Volume mixing
ratio ${\rm VMR}=m_{i}\cdot\mu/\mu_{i}$, where $m_{i}$ and $\mu_{i}$ are the
mass fraction and the particle mass of the species $i$, and $\mu$ is the
atmospheric mean molecular weight.) of each chemical species obtained from the
atmospheric retrievals, along with the VMRs predicted by the equilibrium
chemistry model computed with GGChem (Woitke et al. 2018) assuming solar
abundances. The chemical abundances constrained by the OSIRIS transmission
spectrum alone are consistent with those constrained by the joint spectrum,
except for AlO and $\rm H_{2}O$. Based on the joint retrieval, we obtained
precise measurements of the abundances of AlO ($-6.09^{+0.50}_{-0.77}$ dex),
TiO ($-6.33^{+0.28}_{-0.30}$ dex), MgH ($-2.67^{+0.51}_{-0.58}$ dex), CrH
($-6.40^{+0.65}_{-1.09}$ dex), $\rm H_{2}O$ ($-3.03^{+0.33}_{-0.32}$ dex), and
CO ($-3.89^{+0.79}_{-0.92}$ dex). Sheppard et al. (2021) also reported the
detection of AlO in the optical wavelength range. We note that the main
spectral features of AlO are from 0.40 to 0.55 microns and thus the
measurement on AlO was mainly contributed by the STIS G430L spectrum rather
than the OSIRIS R1000R spectrum. In addition, the low abundances of MgO, VO,
CaH, and FeH disfavored the presence of these species in the terminator
atmosphere of HAT-P-41b. Comparing the retrieved VMRs with the equilibrium
chemistry model assuming solar abundances, we found super-solar abundances for
all well-constrained chemical species in the joint retrieval (AlO, TiO, MgH,
CrH, $\rm H_{2}O$, and $\rm CO_{2}$). Although the detections of MgH and CrH
were tentative, there is strong evidence for a super-solar TiO abundance (and
a super-solar $\rm H_{2}O$ abundance when considering the NIR data),
supporting a metal-rich atmosphere for HAT-P-41b.

The high metallicity of the host star ($\rm[Fe/H]=0.21\pm 0.10$; Hartman et
al. 2012) provides the conditions for the formation of a metal-rich atmosphere
for HAT-P-41b. Considering that the equilibrium temperature of HAT-P-41b
($T_{\rm eq}=1941\pm 38$ K) is in the empirical transition zone between hot
and ultra-hot planets around 2150 K (Stangret et al. 2022), its atmospheric
composition may contain both characteristics in the low-temperature and high-
temperature regimes from a statistical point of view. Therefore, aside from
the detected $\rm H_{2}O$, $\rm CO_{2}$, and refractory molecules, atomic
species such as H, Fe, Mg, and alkali metals may also be expected in the
atmosphere of HAT-P-41b in future high-resolution observations.

The first tentative evidence for CrH was found in the atmosphere of WASP-31b,
a hot Jupiter with an equilibrium temperature $\sim$370 K cooler than
HAT-P-41b, based on the observations of HST STIS G750L (Braam et al. 2021).
The presence of CrH in WASP-31b was recently confirmed by high-resolution
transmission spectroscopy (Flagg et al. 2023). With the addition of HAT-P-41b
that potentially bares CrH in the atmosphere, future studies could target at
the formation of Cr-bearing clouds in the planetary atmosphere (Lodders &
Fegley 2006; Morley et al. 2012).

There are other hot Jupiters with temperatures, masses, radii, and host star
parameters similar to HAT-P-41b, for example, CoRoT-1b (Barge et al. 2008),
HAT-P-7b (Pál et al. 2008), HAT-P-32b (Hartman et al. 2011), HAT-P-65b
(Hartman et al. 2016), and WASP-19b (Hebb et al. 2010). Observations of these
planets show that hot Jupiters with transition-zone equilibrium temperatures
orbiting F-G type stars can exhibit rich spectral features of both volatile
species (e.g., $\rm H_{2}O$) and refractory species (e.g., TiO) if not
obscured by high-altitude clouds or hazes (Huitson et al. 2013; Damiano et al.
2017; Sedaghati et al. 2017; Tsiaras et al. 2018; Chen et al. 2021; Sedaghati
et al. 2021; Czesla et al. 2022; Glidic et al. 2022; Changeat et al. 2022; Li
et al. 2023; Jiang et al. 2023). The day-night and morning-evening differences
could allow $\rm H_{2}$ and $\rm H_{2}O$ to be dissociated into $\rm H^{-}$ in
the hot day-side and evening-side atmospheres, while recombining into water
molecules or other metal hydrides in the cooler night-side and morning-side
atmospheres (Parmentier et al. 2018). Our preliminary detection of MgH and CrH
in the terminator atmosphere sheds lights on the presence of $\rm H^{-}$ in
the hot day-side atmosphere. We speculate that the emission features of the
day-side atmosphere might be obscured by the continuum absorption of $\rm
H^{-}$, resulting in an emission spectrum with weak features and thus tending
toward an isothermal profile, which would explain why the thermal inversion
based on the secondary eclipse measurements from HST and Spitzer was
inconclusive (Fu et al. 2022).

In this study, we analyzed the optical transmission spectrum of HAT-P-41b
using data from five transits observed by GTC OSIRIS. Our free-chemistry
atmospheric retrievals and model comparisons provided strong evidence for the
presence of TiO in the atmosphere of HAT-P-41b, and we also tentatively
detected the signatures of MgH and CrH. Our joint retrieval analysis combining
GTC, HST, and Spitzer transmission spectra yielded results consistent with the
analysis of the OSIRIS data and provided precise constraints on the abundances
of TiO, MgH, CrH, $\rm H_{2}O$, and $\rm CO_{2}$. The detection of abundant
optical absorbers in HAT-P-41b’s high-temperature atmosphere suggests the
possibility of a thermal inversion in its upper atmosphere. However, current
thermal emission spectra have not provided decisive evidence for this
inversion (Mansfield et al. 2021; Fu et al. 2022; Changeat et al. 2022), which
may be due to the presence of H- obscuring other emission features. Future
observations with broader wavelength coverage or higher spectral resolution,
such as those to be obtained by the James Webb Space Telescope (JWST) or other
large ground-based telescopes, are needed to further characterize the day-side
atmosphere of HAT-P-41b and confirm the signatures of metal hydrides.

###### Acknowledgements.

We thank the anonymous referee for their valuable comments and suggestions.
G.C. acknowledges the support by the National Natural Science Foundation of
China (grant Nos. 12122308, 42075122), the B-type Strategic Priority Program
of the Chinese Academy of Sciences (grant No. XDB41000000), Youth Innovation
Promotion Association CAS (2021315), and the Minor Planet Foundation of the
Purple Mountain Observatory. This work is based on observations made with the
Gran Telescopio Canarias installed at the Spanish Observatorio del Roque de
los Muchachos of the Instituto de Astrofísica de Canarias on the island of La
Palma, and on data from the GTC Archive at CAB (CSIC-INTA). The GTC Archive is
part of the Spanish Virtual Observatory project funded by
MCIN/AEI/10.13039/501100011033 through grant PID2020-112949GB-I00. This work
has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Barge et al. (2008) Barge, P., Baglin, A., Auvergne, M., et al. 2008, A&A, 482, L17
  * Benneke & Seager (2013) Benneke, B. & Seager, S. 2013, ApJ, 778, 153
  * Braam et al. (2021) Braam, M., van der Tak, F. F. S., Chubb, K. L., & Min, M. 2021, A&A, 646, A17
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Burrows et al. (2008) Burrows, A., Budaj, J., & Hubeny, I. 2008, ApJ, 678, 1436
  * Cepa et al. (2000) Cepa, J., Aguiar, M., Escalera, V. G., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 623–631
  * Changeat et al. (2022) Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3
  * Chen et al. (2021) Chen, G., Pallé, E., Parviainen, H., Murgas, F., & Yan, F. 2021, ApJ, 913, L16
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Czesla et al. (2022) Czesla, S., Lampón, M., Sanz-Forcada, J., et al. 2022, A&A, 657, A6
  * Damiano et al. (2017) Damiano, M., Morello, G., Tsiaras, A., Zingales, T., & Tinetti, G. 2017, AJ, 154, 39
  * Feroz & Hobson (2008) Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Flagg et al. (2023) Flagg, L., Turner, J. D., Deibert, E., et al. 2023, ApJ, 953, L19
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fortney et al. (2010) Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396
  * Fu et al. (2022) Fu, G., Sing, D. K., Deming, D., et al. 2022, AJ, 163, 190
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Gibson et al. (2013) Gibson, N. P., Aigrain, S., Barstow, J. K., et al. 2013, MNRAS, 428, 3680
  * Glidic et al. (2022) Glidic, K., Schlawin, E., Wiser, L., et al. 2022, AJ, 164, 19
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27
  * Hartman et al. (2012) Hartman, J. D., Bakos, G. Á., Béky, B., et al. 2012, AJ, 144, 139
  * Hartman et al. (2016) Hartman, J. D., Bakos, G. Á., Bhatti, W., et al. 2016, AJ, 152, 182
  * Hartman et al. (2011) Hartman, J. D., Bakos, G. Á., Torres, G., et al. 2011, ApJ, 742, 59
  * Hebb et al. (2010) Hebb, L., Collier-Cameron, A., Triaud, A. H. M. J., et al. 2010, ApJ, 708, 224
  * Huitson et al. (2013) Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Jiang et al. (2023) Jiang, C., Chen, G., Pallé, E., et al. 2023, A&A, 675, A62
  * Jiang et al. (2022) Jiang, C., Chen, G., Pallé, E., et al. 2022, A&A, 664, A50
  * Kass & Raftery (1995) Kass, R. E. & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773
  * Lewis et al. (2020) Lewis, N. K., Wakeford, H. R., MacDonald, R. J., et al. 2020, ApJ, 902, L19
  * Li et al. (2023) Li, X.-K., Chen, G., Zhao, H.-B., & Wang, H.-C. 2023, Research in Astronomy and Astrophysics, 23, 025018
  * Lodders & Fegley (2006) Lodders, K. & Fegley, B., J. 2006, in Astrophysics Update 2, ed. J. W. Mason, 1
  * Mansfield et al. (2021) Mansfield, M., Line, M. R., Bean, J. L., et al. 2021, Nature Astronomy, 5, 1224
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67
  * Morley et al. (2012) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172
  * Nortmann et al. (2016) Nortmann, L., Pallé, E., Murgas, F., et al. 2016, A&A, 594, A65
  * Pál et al. (2008) Pál, A., Bakos, G. Á., Torres, G., et al. 2008, ApJ, 680, 1450
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110
  * Parviainen (2015) Parviainen, H. 2015, MNRAS, 450, 3233
  * Parviainen & Aigrain (2015) Parviainen, H. & Aigrain, S. 2015, MNRAS, 453, 3821
  * Seager & Sasselov (2000) Seager, S. & Sasselov, D. D. 2000, ApJ, 537, 916
  * Sedaghati et al. (2017) Sedaghati, E., Boffin, H. M. J., MacDonald, R. J., et al. 2017, Nature, 549, 238
  * Sedaghati et al. (2021) Sedaghati, E., MacDonald, R. J., Casasayas-Barris, N., et al. 2021, MNRAS, 505, 435
  * Sellke et al. (2001) Sellke, T., Bayarri, M. J., & Berger, J. O. 2001, The American Statistician, 55, 62
  * Sharp & Burrows (2007) Sharp, C. M. & Burrows, A. 2007, ApJS, 168, 140
  * Sheppard et al. (2021) Sheppard, K. B., Welbanks, L., Mandell, A. M., et al. 2021, AJ, 161, 51
  * Stangret et al. (2022) Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2022, A&A, 662, A101
  * Tsiaras et al. (2018) Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156
  * Wakeford et al. (2020) Wakeford, H. R., Sing, D. K., Stevenson, K. B., et al. 2020, AJ, 159, 204
  * Welbanks et al. (2023) Welbanks, L., McGill, P., Line, M., & Madhusudhan, N. 2023, AJ, 165, 112
  * Woitke et al. (2018) Woitke, P., Helling, C., Hunter, G. H., et al. 2018, A&A, 614, A1
  * Wöllert et al. (2015) Wöllert, M., Brandner, W., Bergfors, C., & Henning, T. 2015, A&A, 575, A23

##  Appendix A Supplementary figures and tables

Table 4: Observation summary for HAT-P-41b.

# | Proposal ID | Date (UT) |  Slit (′′) |  Readout (kHz) |  Exp. time (s) | Useful exposures |  Airmass (start – min – end) |  Seeingaa$a$Seeing is measured with the full width at half maximum of the stellar PSF along the spatial direction at the central wavelength of $\sim$750 nm. (95% interval)
---|---|---|---|---|---|---|---|---
OB-1 | GTC52-13Abb$b$PI: E. Pallé. | 2013-06-21 | 12 | 200 | 8.0 | 426 | 1.71 – 1.09 – 1.23 |  $1.07^{\prime\prime}$ – $2.34^{\prime\prime}$
OB-2 | GTC9-14Abb$b$PI: E. Pallé. | 2014-07-14 | 40 | 200 | 7.0 | 668 | 2.27 – 1.10 – 1.18 |  $0.69^{\prime\prime}$ – $1.60^{\prime\prime}$
OB-3 | GTC65-16Acc$c$PI: R. Alonso. | 2016-07-05 | 12 | 500 | 7.5 | 987 | 1.74 – 1.09 – 1.33 |  $0.84^{\prime\prime}$ – $1.32^{\prime\prime}$
OB-4 | GTC65-16Acc$c$PI: R. Alonso. | 2016-07-13 | 12 | 500 | 20.0 | 707 | 1.39 – 1.09 – 2.18 |  $1.37^{\prime\prime}$ – $2.24^{\prime\prime}$
OB-5 | GTC90-22Add$d$PI: F. Murgas. | 2022-08-11 | 40 | 200 | 10.0 | 602 | 1.35 – 1.10 – 2.43 |  $0.94^{\prime\prime}$ – $1.73^{\prime\prime}$

777 The OSIRIS instrument was positioned at the Nasmyth-B focus of GTC for the
first four observations, and at the Cassegrain focus for OB-5.

Table 5: Measurements of the HAT-P-41 system and updated constraints on its
companion.

Parameters | HAT-P-41 | Companion | References
---|---|---|---
Right ascension (deg) | 297.322651868 | 297.322580691 | [1]
Declination (deg) | 4.672411376 | 4.671410231 | [1]
Parallax (mas) | $2.85\pm 0.02$ | $2.91\pm 0.03$ | [1]
Radial velocity (km/s) | $32.54\pm 0.47$ | $32.63\pm 3.81$ | [1]
Proper motion (mas/yr) | 7.30 | 7.65 | [1]
Effective temperature (K) | $6390\pm 100$ | $4498^{+28}_{-33}$ | [2], this work
Metallicity (dex) | $0.21\pm 0.10$ | $0.41^{+0.13}_{-0.17}$ | [2], this work
Gravity (cgs) | $4.14\pm 0.02$ | $4.30^{+0.19}_{-0.20}$ | [2], this work
Mass ($M_{\sun}$) | $1.418\pm 0.047$ | $0.36^{+0.19}_{-0.19}$ | [2], this work
Radius ($R_{\sun}$) | $1.683^{+0.058}_{-0.036}$ | $0.708^{+0.022}_{-0.022}$ | [2], this work

[1] Gaia Collaboration et al. (2016, 2023) [2] Hartman et al. (2012)

Figure 6: Fitting results of the broadband transit light curves of HAT-P-41b.
Each row corresponds to one transit observation. Left column: the normalized
fluxes (black error bars), the best-fit models (red), and the systematic noise
(green, also the common-mode systematics). Middle column: the detrended light
curves and the best-fit transit models. Right column: the residuals.