ACCESS: A featureless optical transmission spectrum for WASP-19b from Magellan/IMACS
Néstor Espinoza Benjamin V. \ Rackham Andrés Jordán Dániel Apai Mercedes López-Morales David J. Osip Simon L. Grimm Jens Hoeijmakers Paul A. Wilson Alex Bixel Chima McGruder Florian Rodler Ian Weaver Nikole K. Lewis Jonathan J. Fortney Jonathan Fraine


firstpage–lastpage
Exploring Sensitivity to NMSSM Signatures with Low Missing Transverse Energy at the LHC
    A. Titterton^1,2,3, U. Ellwanger^4, H.U. Flaecher^1, S. Moretti^2,3,
 C.H. Shepherd-Themistocleous^3
    December 30, 2023
========================================================================================================





The short period (0.94-day) transiting exoplanet WASP-19b is an exceptional target 
for transmission spectroscopy studies, due to its relatively large 
atmospheric scale-height (∼ 500 km) and equilibrium temperature (∼ 2100 K). 
Here we report on six precise spectroscopic Magellan/IMACS observations, five 
of which target the full optical window from 0.45-0.9μm and one targeting 
the 0.4-0.55μm blue-optical range. Five of these datasets are consistent with 
a transmission spectrum without any significant spectral features, while one shows a significant slope as a function of wavelength, 
which we interpret as arising from photospheric heterogeneities in the star. Coupled 
with HST/WFC3 infrared observations, our optical/near-infrared measurements point to 
the presence of high altitude clouds in WASP-19b's atmosphere in agreement with previous 
studies. Using a semi-analytical retrieval approach, considering both planetary and stellar 
spectral features, we find a water abundance consistent 
with solar for WASP-19b and strong evidence for sub-solar abundances for optical absorbers such as 
TiO and Na; no strong optical slope is detected, which suggests that if hazes are present, 
they are much weaker than previously suggested. In addition, two spot-crossing events are 
observed in our datasets and analyzed, including one of the first unambiguously detected bright spot-crossing events on an exoplanet host star.





techniques: spectroscopic – planets and satellites: atmospheres – 
planets and satellites: individual: WASP-19b –
stars: individual: WASP-19 – stars: activity – stars: starspots








§ INTRODUCTION




Transmission spectroscopy, the study of the variation of the planetary radius as a function of wavelength due to different opacity sources in its atmosphere <cit.>, offers  
one of the most successful approaches to date for detecting atomic and molecular absorption in exoplanet atmospheres.
From detections of water vapor <cit.>, sodium and potassium <cit.> and many signatures of aerosols, i.e., clouds and/or hazes <cit.>, transmission spectroscopy studies are paving the way for us to understand the compositions 
of distant worlds. This exploration is not only interesting in that it reveals 
the diversity of compounds present in these distant worlds and the 
thermochemical processes at play <cit.>, but it can also offer a glimpse into the formation mechanisms that generated the atmospheric composition we observe 
today <cit.>.

Among the exoplanets amenable to atmospheric characterization via transmission spectroscopy, the ultra-short period (P=0.94-day) exoplanet WASP-19b <cit.> orbiting a G8V star is a very interesting one in many ways. It is not only 
inflated <cit.> but also, due to its proximity to its host star, has a high equilibrium temperature (T = 2,100 K). This, in turn, gives 
it a large atmospheric scale-height of order H∼ 510 km, which should give rise to a signal in transmission between 200-600 ppm. Indeed, low-resolution transmission spectroscopy studies in the infrared using the Hubble Space Telescope (HST) have succeeded 
in detecting water vapor absorption in its atmosphere, with amplitude of 300 ppm <cit.>. 

What this infrared signature actually means in terms of the actual water abundance 
and what it means to, e.g., the C/O ratio of its atmosphere or its metallicity, however, is still unclear. On one hand, <cit.> identified the signature as consistent with that from a clear atmosphere, while <cit.> interpreted the signature as arising from a cloudy atmosphere. Optical signatures could, in principle, help us interpret the observed infrared water signature; however, they have 
proven to be elusive mainly due to the lower precisions in most previous optical works from both space and ground-based observatories compared to those in the infrared from HST <cit.>. 

Optical absorbers are fundamental to 
interpret infrared water absorption signatures. If no 
absorbers are detected, this would strengthen the case for clouds in the atmosphere, which could lead to a large water abundance in the exoplanet given its large observed water signature. This could hint, in turn, to a low C/O ratio in its atmosphere, 
which would make it consistent with our understanding of solid enrichment in planetary 
envelopes <cit.>. On the other hand, the actual detection of optical signatures could help us understand many underlying processes in these hot exoplanets, including cloud formation and whether cold traps in the night side could sequester 
Ti or V from the terminator region probed in transmission, preventing the formation 
of TiO/VO expected from equilibrium 
chemistry calculations <cit.>. Without the strong optical absorbers of TiO/VO, Na 
and K should be the dominant optical absorbers in a clear atmosphere; detecting any 
of those, thus, would be exciting in terms of understanding the underlying dynamics 
in a hot exoplanet such as WASP-19b. A-priori, however, one would expect the atmosphere to be cloudy, given the ubiquity of clouds in giant exoplanet atmospheres <cit.>. Data, of course, are the ultimate judge that should guide our a-posteriori beliefs.

Motivated by the exciting possibilities in the optical for WASP-19b, we obtained observations for this exoplanet as part of the Arizona-CfA-Católica-Carnegie Exoplanet Spectroscopy 
Survey (ACCESS). ACCESS is a multi-institutional effort which 
aims at obtaining a large, homogeneous library of optical spectra of exoplanet atmospheres,
providing a complement to longer-wavelength observations and key information necessary to understand 
and interpret observed molecular and/or atomic features. To this end, 
ACCESS is currently using the Inamori-Magellan Areal Camera and Spectrograph <cit.> instrument mounted at the Magellan Baade 6.5m Telescope in Las Campanas Observatory in Chile. This 
project has allowed us to study so far the atmospheres of more than a dozen exoplanets including the giant exoplanet 
WASP-6b <cit.> and the sub-Neptune GJ 1214b <cit.>.

During the writing of this manuscript, <cit.> presented an optical transmission spectrum of WASP-19b from VLT/FORS2 measurements obtained at three different epochs in three different but overlapping wavelength ranges (4,400–6,100 Å, 5,400–8,460 Å, and 7,400–10,000 Å).
Their combined spectrum shows structures which they interpret as (1) a strong scattering slope, (2) strong TiO signatures and (3) a water signature in the 
near-infrared, all of which point to a clear atmosphere in the optical. 
This finding is in contrast to previous low-resolution 
observations with HST/STIS <cit.> showing a relatively 
featureless spectrum, which was interpreted as evidence for clouds and lack of TiO absorption. 
In this work, we present six Magellan/IMACS observations obtained between 2014 and 2017, five of which cover the same wavelength 
range over which most of the strong optical signatures are seen by <cit.> (i.e., from 4,500 to 9,500 Å), and one observation which covers the bluest wavelengths (4,200–5,500 Å), where the strong scattering signature is observed in that work. We do not observe this strong scattering slope 
in our data, nor do we detect the strong TiO absorption feature reported by <cit.>.

This paper is structured as follows. In Section <ref> we present our Magellan/IMACS observations of WASP-19b along with the photometric monitoring of the host star available around the 
times of these observations. In Section <ref> we present 
the analysis of the Magellan/IMACS data, which includes the determination of the transmission spectra for each of the transit observations as well as an analysis of two observed spot-crossing events and a study of the effect of unocculted stellar 
photospheric heterogeneities on our measurements. 
In Section <ref> we discuss the implications of our findings in light of previous attempts at detecting optical absorbers in the atmosphere of this exoplanet and the HST infrared observations, and in 
Section <ref> we outline our conclusions.



§ DATA






The data presented in this work considers both spectrophotometric measurements and long-term photometric monitoring of the star. The former, described in detail in Section <ref>, consists of six transits observed with Magellan/IMACS in 2014, 2015 and 2017, which we use to extract WASP-19b's transmission spectrum in Section 
<ref>. The long-term photometric monitoring, on the other hand, 
consists of data taken during the 2014 season with the SMARTS 1.3m telescopes at 
Cerro Tololo Inter-American Observatory (CTIO) and data gathered in 2017 by 
the All-Sky Automated Survey for Supernovae <cit.>. 
These photometric data and their importance to the present work are presented in 
Section <ref>. 



 §.§ Observations with Magellan/IMACS


Transit observations of WASP-19b were carried out with the Inamori-Magellan Areal Camera and Spectrograph <cit.> instrument mounted at the Magellan Baade 6.5m Telescope 
in Las Campanas Observatory (LCO) in Chile in 2014 (March 22 and April 29), 2015 (June 3rd) and 2017 (February 22, 
April 4th and April 12). Hereafter we refer to these datasets by their observation date in the YY/MM/DD format. The f/2 camera was used in the multi-object spectrograph configuration, and a mask was cut with 10” wide slits at the positions of the target and seven comparison 
stars. The slit size was chosen to be large enough to avoid significant slit losses with typical seeing conditions at LCO. 

One of the key features of Magellan/IMACS is 
its large field of view (FOV) for the f/2 camera (27'× 27'), which allows for observations of several 
comparison stars in order to correct for common time-series variations. This is in contrast with other 
instruments, such as VLT/FORS2 
<cit.>, or 
Gemini/GMOS  <cit.>, whose relatively 
smaller FOVs (∼ 7'× 7' and 2'× 2', respectively) typically only allow 
the observations of a single comparison star. The availability of multiple comparison stars 
makes it easier to detect and confrim features in the data, including spot-crossing events, which are 
known to occur during WASP-19b transits <cit.> and which could bias our measurements if not accounted for. A 
list of the comparison stars observed simultaneously with WASP-19 in our study is shown in Table <ref>; the star J09545906-4544371 (V = 10.3), however, was not 
used in the final analysis as the observed spectra saturated the detector in the wavelength range of interest. 

The observations during the 2014 and 2015 season were performed using the  
(30s) readout mode, with the 300+17.5 (300 lines/mm and 17.5 degree blaze angle) 
grism, which provided spectra spanning the 4,100–9,350 Å range and an adequate spectral resolution for 
our purposes (∼ 1.3 Å per pixel). The observations obtained during the 2017 season were carried 
out with the  (31s) readout mode and a blocking filter (BF) in order to assess if second order 
light could be impacting our observations. The February 22 and April 4th 2017 observations were carried out 
using the same grism as the 2014 and 2015 seasons, but the GG455 blocking filter, which has a sharp cutoff 
at 4,550 Å  was used; we refer to this mode as the “red" setup. For the 17/04/12 observations, we aimed to study the blue part of the transmission spectrum, and thus we 
decided to use the same grism but with the WB3600-5700 blocking filter, which allowed us to block any light 
bluewards of 3,600 Å and redwards of 5,700 Å. We refer to this mode as the “blue" setup. All observations were acquired with a 
2×2 binning. The data were reduced with a pipeline already outlined in previous works <cit.>, and which will be fully detailed in a future paper (Espinoza et al., 2018, in prep.).   
 




 §.§ Photometric monitoring


WASP-19 is known to display photometric variability at the 5-10 mmag level <cit.>, which is most likely caused by the rotational modulation of starspots <cit.>. 
It has been long recognized that starspots (even if unocculted) can impact transiting exoplanet transmission spectra
by introducing
time- and wavelength dependent biases into the apparent size of the planet <cit.>. 
The origin of the spectral contamination is the spectral difference between the disk-averaged spectrum (obtained pre-transit) and the spectrum of the transit chord (the actual light source for the transmission measurement), which has been termed the transit light source effect (TLSE) by <cit.>. 
Several studies have recognized and dealt with this effect in the past with different 
levels of complexity <cit.>. 
However, the work of <cit.> points out that the correction of the impact of starspots and faculae, 
which can introduce not only apparent spectral slopes but also false spectral features in transmission, is not straightforward. 
In particular, they demonstrate that the commonly assumed linear relation between stellar photometric variability and starspot/facular areal covering fractions will underestimate the covering fractions in most realistic cases. 
Nevertheless, in principle if spots/faculae covering fractions and their temperatures are known, one should be able to constrain their impact on the transmission spectrum. 
Such constraints can be informed by photometric monitoring of the exoplanet host star, especially if data were acquired close to the epochs of the transmission spectroscopic observations. 

Considering this, we decided to obtain and/or use photometric observations 
when possible/available. During the 2015 season we were not able to obtain 
photometric monitoring of WASP-19. However, during the 2014 and 2017 seasons, 
photometric data covering the epochs of our Magellan/IMACS observations are 
available. We analyze those datasets in detail in Appendix <ref>. Here 
we simply state some key results from that analysis, that we will use in future sections, 
especially in Section <ref> in which we use the observed level of variability 
to constrain the possible spot covering fractions and thus predict the expected 
level of stellar contamination in WASP-19b's transmission spectrum.

For the 2014 season, WASP-19 was monitored with the SMARTS 1.3m telescope at CTIO through 
a V filter in a continuation of the long-term photometric campaign presented by <cit.>.
However, strong systematics in the data precluded any 
meaningful photometric variability measurements for WASP-19, despite the fact that we 
used the same observational setup as the long-term photometric 
monitoring presented in <cit.>. 
We re-analyzed the <cit.> photometric dataset 
in order to verify that this was not an artifact of our data reduction pipeline and 
recovered the same results for WASP-19 as the ones presented in that work. However, 
the comparison stars all showed similar levels of variability as WASP-19 in that dataset, implying that the observed variability of WASP-19 
in the work of <cit.> was mostly due to systematic effects and not due 
to intrinsic astrophysical variability (see Section <ref> 
for details). This implies that the stellar activity 
corrections made in that work to the transmission spectrum of WASP-19b are most likely 
not adequate. This is, however, unimportant for the overall shape of the transmission 
spectrum, albeit for a possible offset in transit depth between the optical HST/STIS and 
infrared HST/WFC3 data presented in <cit.>. This will be important to consider 
in Section <ref>, when we interpret the HST optical and infrared observations 
in light of our Magellan/IMACS optical data.

For the 2017 season, ASAS-SN V-band photometry was gathered from the ASAS-SN webportal[<https://asas-sn.osu.edu/>]. This data covered all of 
our Magellan/IMACS observations during the 2017 season, allowing the variability 
of WASP-19 during that observing season to be estimated (see Section 
<ref> for details). We observed a
2% peak-to-peak variability for WASP-19, similar to what has been 
observed in previous works <cit.>.



§ ANALYSIS



In this Section we analyse both our band-integrated (“white light”) and 
wavelength-dependent Magellan/IMACS light curves for WASP-19 and the observed comparison stars, 
extracted from the spectra described in 
Section <ref>. 
Our approach is to optimally use the available information in order to account for systematic trends in the data and to analyse the datasets in a uniform way in order to introduce 
the least possible amount of variation between them. 
To this end, we take into account both the comparison stars, which provide essential information 
to correct common atmospheric and instrumental effects present in our light curves, and 
external parameters (such as, e.g., the variation of the full-width at half maximum, FWHM, 
the movement of the spectral trace, etc., during the course of the observations), which can aid in 
correcting systematic effects intrinsic to our WASP-19 light curves. The modelling approach 
is first described in detail in Section <ref>, and this is later used 
in order to analyse the band-integrated (Section <ref>) and wavelength-dependent 
(Section <ref>) light curves. 



 §.§ Light curve modelling





We modeled the flux of WASP-19 as the multiplicative combination of signals and, thus, we work in logarithmic space (base 10) in order to represent components as linearly additive. 
In what follows, however, we transform all products (e.g., plots, precision of our best-fit models) to relative flux space for ease of comparison with other studies. 

We follow the approach of <cit.> and perform a Principal Component Analysis (PCA) in logarithmic space on our 6 comparison stars in order to optimally use all the information embedded in our comparison stars regarding common systematic effects (i.e., present in every star in the field).
We extract the 6 signals, S_i(t), which allow us to (linearly) reconstruct their observed flux variations. 
As explained in <cit.>, for each signal S_i(t), an eigenvalue λ_i allows us to decide how important – in terms of information content – that signal is to reconstruct the original flux variations of the comparison stars (and, we assume, to explain the flux variations in WASP-19's light curve as well). 
In what follows, we assume the indices i of the signals order the signals from the most important (larger eigenvalue, i=0 in our case) to the least important (smaller eigenvalue, i=5 in our case). 

In order to account for systematic effects which might be unique to WASP-19's light curve (due to, e.g., color differences, detector and/or instrumental artifacts), we also include external parameters that are recorded during our observations as possible linear (in logarithmic space) regressors.
These are standardized by subtracting their means and dividing by their standard deviations. 
In this study, we observe that the rotator angle as a function of time R(t), the wavelength drift of the wavelength solution as a function of time Δλ (t), the variation of the full-width at half maximum as a function of time FWHM(t), time t, and the squares of the mentioned variables (i.e., R^2(t), Δλ^2 (t), FWHM^2(t) and t^2) are useful regressors that account for the observed variations in our WASP-19 light curves. 
However, determining which of these regressors are the most important to include in our model for each of our light curves is not trivial. 
The same goes for the number of signals S_i(t). 

In order to incorporate our ignorance regarding a “best-fit" model in our retrieved transit parameters, we follow the model averaging technique outlined for transmission spectroscopy studies in <cit.>, which has already been used and tested with both space-based <cit.> and ground-based <cit.> observations. 
We define k models for the base-10 logarithm of WASP-19's flux of the form 

    M_k(t) ∼  C_k + ∑_i=0^N_kα_i,k S_i(t) + X_k(t) + log_10 T(t|θ_𝐭) + ϵ_k.

Here, C_k is a zero point in (log-)flux for the k-th model, N_k is the number of signals S_i(t) to be included in the fit weighted by α_i,k in the k-th model, X_k(t) is a linear combination of up to 3 of the 8 regressors mentioned in the previous paragraph[The decision of up to 3 terms was empirically based on our best signal-to-noise dataset, in which we observed that more terms only slow down the process without additional improvement on the fits.], T(t|θ_𝐭) is the transit model with parameters θ_𝐭, and ϵ_k is a random variable assumed to be gaussian and uncorrelated and, thus, defines a photometric jitter term Var[ϵ_k] = σ^2_k empirically determined for each fit. 
This defines a total of 552 models to be tested on our data if all signals S_i(t) are incorporated in our framework (6), and models of up to 3 terms out of 8 possible variables are generated (i.e., ∑_n=1^3 8n=92 possible models). 
This number varies slightly for the wavelength-dependent light curves, because of the varying spectral coverage of the comparison stars (i.e., some of them do not cover the full spectral range covered by WASP-19 due to the positioning of the mask slits relative to chip gaps and detector edges). 

In practice, after removing obvious outliers from our light curves using a sigma-clipping procedure, each model is first fit with a Levenberg-Marquardt algorithm using the out-of-transit data in order to estimate the linear parts of the model (C_k, the coefficients of the signals S_i(t) and X_k(t)). 
Then, the data are detrended with this fit and initial transit parameters θ_𝐭 are estimated by fitting the non-linear transit model also using a Levenberg-Marquardt algorithm.
The retrieved parameters in this two-step process are then fed to an MCMC algorithm as starting points, and using  <cit.>, the posterior distribution of all the parameters in the full model is explored (including now the jitter term, σ^2_k). 
The median of the posteriors are then defined as the best-fit parameters of the k-th model.
The corrected Akaike's Information Criterion <cit.> of each model is then computed, and the final, model-averaged parameters and uncertainties are obtained by weighing the best-fit parameters of each model by their respective AICc values as outlined in <cit.>. 

As a final note on our modelling, as noted in <cit.>, our limb-darkening assumptions can lead us to retrieve biased transit parameters if not properly accounted for in the transit light curve fitting procedure, especially for very precise transit light curves like those we present in the next sub-section. 
In addition, for active stars such as WASP-19[Here and in what follows, we refer to “activity" as stars producing large rotational modulations and/or having evidence of heterogeneities in their photospheres.], the limb-darkening coefficients can actually be modified due to the non-uniform stellar surface owing to the presence of unocculted spots <cit.>. 
We therefore chose to fit the limb-darkening coefficients in our transit light curve fitting procedure rather than to fix them to pre-computed values. 

To select the optimal limb-darkening law to use in our case, we use the method outlined in <cit.>; 
in brief, we generate light curves with similar geometric and noise properties as the ones 
observed using a non-linear limb-darkening law and then try to retrieve the transit parameters using two-parameter laws and the linear law. 
We then identify the law that retrieves the parameters closest to the input values in a mean-square error sense as the optimal law to use. 
We observed that the square-root law is, in this case, the optimal one for noise levels as low as ∼ 1000 ppm, and thus we use this limb-darkening law throughout our analysis of the band-averaged (i.e., "white-light") light curves. 
Beyond this noise level, our analysis shows that the linear law is as good as the other tested laws, and so we use it for the wavelength-dependent light curves. 
The [<https://www.cfa.harvard.edu/ lkreidberg/batman/>] <cit.> package was used to model the transit light curves, and the limb-darkening coefficients were sampled following the uninformative sampling procedure outlined in <cit.>. 
The routines used to model our transit light curves are available at GitHub[<http://www.github.com/nespinoza/exotoolbox>].



 §.§ “White-light" light curve analysis





The band-integrated (“white-light") light curves were obtained by following the scheme outlined in Section <ref>. 
In this case, 250 walkers are used on each of the k models in order to explore the parameter space. A total of 1000 steps are used, discarding the first 500 as burn-in.

In Figure <ref>, we present our systematics-removed transit light curves along with the residuals for each model.
The systematics models were obtained by averaging each of the k models fitted to the data following the procedures outlined in Section <ref>, weighted by their AICc values. 
It is evident that spot-crossing events are present in the 14/04/29 and 17/04/04 light curves and that the features are wavelength-dependent. 
For simplicity, these portions of the light curves were left out of the white-light analysis.
These portions, however, are used and modelled in the analysis of the wavelength-dependent light curves (see next sub-section). 
Following the procedures we outline in Section <ref>, we find no evidence for spots in the other light curves.

The precision of our white-light light curves is generally 4–5 times the 
photon noise in our datasets, except for the 14/03/22 dataset in which the precision is only 3 times the photon noise. 
Notably, the white-light light curve obtained on March 2014, which was taken under excellent conditions, attains a precision of ∼ 270 ppm in 33-second exposures (giving a total cadence of ∼60 s, including the 30 second readout time of the f/2 camera on  mode). 
This is one of the best photometric precisions attained from the ground. 
Our other light curves are similarly precise, with the exception of the June 2015 data (which were taken under poor photometric conditions) and the April 12th 2017 data, which were taken with our “blue" setup and, thus, on a portion of the spectrum in which (1) WASP-19 is not very bright and (2) the throughput of the instrument is not at its peak. 
We note that these observations are also more precise than both the HST ones for WASP-19, where <cit.> reports a precision of 280 ppm with 293 second exposures with HST/WFC3, and the VLT/FORS2 ones, where <cit.> reports a precision of 697 ppm with 30 second exposures for WASP-19. 
Therefore, we believe that Magellan/IMACS will be an interesting option for very precise follow-up light curves of Transiting Exoplanet Survey Satellite <cit.> targets.



Our precise white-light transit light curves allow us to refine the planetary parameters of the WASP-19b system. 
In Table <ref> we summarize the best-fit values for the transit light curve parameters for each of the observed nights. The retrieved parameters are all consistent within the errorbars between observations, which demonstrates the quality of both the precision and accuracy of our white-light transit light curves through different seasons and observing modes. Using the observed times of transit center, we 
obtain a revised ephemerides for WASP-19b (see details in Appendix 
<ref>); in addition, we combine the orbital and physical 
parameters presented in Table <ref> and compare them 
to other works performing precise photometry for WASP-19b in Table <ref>. As 
can be seen, the agreement with the work of <cit.> and 
<cit.> is, in general, excellent. We combine the precise planet-to-star radius ratios obtained in <cit.> and in 
<cit.> with ours, obtaining a revised value for this 
parameter of R_p/R_* = 0.14257 ± 0.00027. Using the stellar radius derived in <cit.> of R_s = 1.018 ± 0.021R_⊙ (where, for simplicity, we have added the systematic and random uncertainties) for WASP-19, we obtain a revised planetary radius of R_p = 1.413± 0.029 R_J. 
Using the mass for WASP-19b also derived in <cit.>, we obtain a planetary gravity of 1414± 95 cm/s^2. The importance of these revised parameters will be made evident in Section <ref>.







 §.§ Wavelength-dependent light curve analysis







For the analysis of the wavelength-dependent light curves, we followed a similar analysis as the one outlined 
for the white-light light curve analysis. However, we fix here the value of P, a/R_*, i, 
and t_0 to the values found in our white-light light curve analysis, so the only free 
parameters that define our wavelength-dependent transit light curves are the planet-to-star 
radius ratio, R_p/R_*(λ), and the limb-darkening coefficient of the linear law used to account for this effect. For the systematics, we use the same model as in for the white-light analysis. 

For the light curves that show spot-crossing events, however, we apply a special treatment. 
Although proper analysis of the spot shapes, contrasts, and positions will be performed in the next section, for the measurement of the transmission spectrum here we subtract the predicted model of the white-light light curves from the observed transit containing the spot in order to obtain an estimate of the shape of the spot feature on our transit light curves. 
We use this estimate to correct the wavelength-dependent light curves in order to analyze all the light curves from different nights in a uniform way. 
To obtain this estimate from the residuals of our best-fit white-light light curve model, we apply a median filter with a 5-point window on the residuals, which we then smooth with a Gaussian filter with a 3-point 
standard deviation in order to get a smooth version of the spot shape. The result of this procedure is shown 
in Figure <ref>. Using this spot shape, which we denote by S(t), we model the transit at each 
wavelength as T(t|θ_𝐓) - AS(t), where A is a term that modulates the amplitude of 
the spot event at each wavelength.




Figures <ref> and <ref> show the resulting light curves in the 
different wavelength bins analysed in this work.
The bin widths, detailed in Table <ref>, were selected to maximize the detectability of Na I and K I features as well as TiO/VO features (i.e., wide enough to have sufficient signal-to-noise in the light curves, but not too wide to lose spectral resolution). 
Our 2014 data span the wavelength range from 4100 to 9350 Å in bins of variable size spanning from ≈ 150 to 250 Å bins, which results in 28 bins across that wavelength range. 
Depending on the wavelength bin, the wavelength-dependent light curve jitter (i.e., σ_w in our model) varies from 
around ∼ 600-1500 ppm for the March data to around ∼ 900-2,000 ppm for the April data. 
The 2015 data, taken with the same setup but during worse weather conditions, varies between ∼ 1,500-3,000 ppm. 
For the 2017 season, the data in February and the 4th of April, taken with our “red" setup, span a wavelength range from 4,500 to 9,350 Å, in 26 wavelength bins. 
The wavelength-dependent light curves in our February and April 4th observations show jitter with values around 1000-1500 and 800-1200 ppm, respectively. 
Finally, the 2017 data taken on April 12th with our “blue" setup, spans the wavelength range from 4,100 to 5,560 Å, with jitter values close to ≈ 900 ppm for the reddest wavelengths, reaching ≈ 2500 ppm at the bluest wavelengths.





As discussed in Section <ref>, WASP-19 is an active star and as such, differences in both the mean and wavelength-dependent transit depths are expected due to unocculted stellar spots. Differences between the retrieved mean transit depths are indeed observed in our data, with typical differences on the order of ∼ 500 ppm, except for the 17/04/04 dataset, where a difference of ∼ 1000 ppm is observed. 
Because of this, we work in what follows with the mean-subtracted transit depths, although we provide a detailed analysis of each dataset in Section <ref>. 
In Table <ref> we present the transmission spectra from all our epochs.
Figure <ref> illustrates these spectra, grouped by whether we used a BF (which was the case for the 2014-2015 season) or not (which was the case for the data obtained from the 2017 season). 
Overall, almost all the datasets are consistent with a flat line within the errorbars, which suggests that whether a BF is used does not have a large impact at the level of precision of our Magellan/IMACS observations. 

The only exception to the observed featureless spectra is the February 2017 dataset, which shows an apparent slope that decreases as a function of wavelength with an amplitude of about 3000 ppm. 
Interestingly, this dataset was taken during a period of larger photometric variation for WASP-19, according to our analysis in <ref>. 
Cold unocculted spots covering a significant fraction of the star could produce such a slope, which to first order, should make the star appear dimmer and redder.
However, depending on their distribution and sizes, spots on the star can emulate a variety of photometric signatures which do not have a 1-to-1 correspondence with the spotted area, which complicates determinations of the stellar contamination signal from photometric monitoring alone <cit.>.
Alternatively, the crossing of bright areas by the exoplanet could also lead to the same effect <cit.>. 
We discuss the impact of stellar heterogeneities on our transmission spectrum in detail in Section <ref> and provide evidence that this indeed can explain the observed variation in our February 2017 dataset in Section <ref>.

In Figure <ref> we combine all of our measurements except for the February 2017 data (due to the observed slope) in order to generate our combined transmission spectrum. 
The resulting spectrum is consistent with a flat line, with a mean transit depth of -27 ± 44 ppm. 
If we assume this (null) hypothesis to be true, then the probability that we observed these data or more extreme values (i.e., the p-value) using our 28 datapoints (which imply 27 degrees of freedom, given our estimated mean) is 0.15, with which we cannot reject the hypothesis that the data actually comes from a simple flat line. 
Of course, this simple null-hypothesis testing scheme ignores the possibility of other possible hypotheses that could explain the data. 
We discuss those in the next section. 





 §.§ Spot-crossing events







The spot-crossing events mentioned in the previous section for the April 2014 and 2017 data were analyzed using  <cit.>, which allows us to 
fit for multiple spots on the stellar surface simultaneously with the transit parameters. Because the posterior 
sampling is a complicated (possibly multi-modal) function, we here use MultiNest <cit.> via the PyMultiNest package <cit.> to both (1) sample points from the posteriors and (2) compute posterior bayesian evidences, Z = ℙ(D|H), i.e., the probability of the data D given hypothesis H. This latter property of nested sampling algorithms in turn 
allows us to study how complex our models have to be in order to explain the observed light curve distortions (e.g., n spots on the surface) via the posterior odds, ℙ(H_n|D)/ℙ(H_k|D), where ℙ(H_n|D) = ℙ(D|H_n)ℙ(H_n), with ℙ(H_n) being the prior probability on the hypothesis H_n. 
Note that in the case in which the prior probability of the models is the same, 
then the posterior odds are simply the ratio of the evidences, 
ℙ(H_n|D)/ℙ(H_k|D) = Z_n/Z_k. 
In this work, we usually 
work in terms of the differences of the log-evidences, Δln Z, which are the logarithm of the odds, assuming equal prior probabilities for the hypotheses. 
A good review in terms of how the posterior odds translate to frequentist 
significance hypothesis testing can be found in <cit.>; we note from that work that absolute values of log-odds below 1 are usually considered inconclusive, whereas absolute log-odds of 2.5 can be interpreted as moderate evidence and absolute log-odds higher than 5 can be interpreted usually as highly significant. We caution, however, that this comparison with frequentist methods is usually useful for understanding but not very meaningful in practice in the 
sense that frequentist hypothesis testing has only one null hypothesis, whereas proper bayesian model comparison considers a range of possible hypotheses at play. The engine for modelling transiting exoplanet light curves using  with PyMultinest, , is available at GitHub[<http://github.com/nespinoza/spotnest>]. 

A detailed explanation of the analysis of the spot-crossing events can be found in Appendix <ref>; we here discuss the main results of that analysis.
For the April 2014 dataset, using the white-light light curve we found that the best model 
for the spot-crossing event is one with a single spot located at (x,y)=(0.21^+0.01_-0.01,0.68^+0.04_-0.06), where x goes in the direction almost perpendicular to the planet's motion and the center of the star is at (0,0).
The spot has a contrast of F_spot/F_star = 0.84^+0.01_-0.03, (i.e., a dark spot) and a considerably large radius of R_spot = 0.21^+0.01_-0.01 stellar radii. 
For the April 2017 dataset, the one-spot model is also favored; the retrieved position of the spot is 
identified with high precision to be at (x,y)=(0.41^+0.02_-0.02,0.43^+0.08_-0.08) on the stellar surface, with a size of 0.25^+0.05_-0.04 stellar radii. 
The most interesting feature of the April 2017 feature, however, is its contrast: the retrieved spot contrast is F_spot/F_star = 1.16^+0.05_-0.04, which implies the feature is actually a bright spot. 
This makes the April 2017 event, thus, one of the first unambiguously detected bright spot features, i.e. a facular or plage region, on an exoplanet host star <cit.> and the first one on WASP-19. 
A physical depiction of the star-spot-planet system is shown in Figure <ref>, along with the posterior distributions of the spot contrasts, sizes and positions retrieved from our light curve analysis.

Using the derived shapes and positions of the spots from the white-light light curves, we analysed the systematics-removed, wavelength-dependent light curves to study the spot contrast as a function wavelength. 
We followed a similar approach to the one defined in Section <ref>: we fix the orbital parameters to the joint best value found in Section <ref> and the spot position and radius to the values from the white-light analysis, and let the transit depth, the limb-darkening coefficient of the linear law and the spot contrast vary as a function of wavelength in each light curve. 
The resulting contrasts are presented in Figure <ref>. 

Using stellar models to model the contributions of the spot and the stellar surface to the observed contrasts (see Appendix <ref> 
for details), we were able to derive spot temperatures for each of the events. 
For the cold spot observed in the April 2014 dataset, we find a temperature of T_s = 5278 ± 81 K (or Δ T = 192 ± 10 K colder than the star). 
For the bright spot, we find a temperature of T_s = 5588 ± 92 K (or Δ T = 137 ± 10 K hotter than the star). 
These models are also presented in Figure <ref>.





 §.§ Stellar contamination





In addition to altering transit depths through active region crossing events, unocculted active regions can also alter observed transit depths through the transit light source effect.
We assessed the potential for stellar contamination of the transmission spectrum of WASP-19b through the transit light source effect using the method detailed by <cit.>. 
The approach consists of placing constraints on allowed spot and faculae covering fractions using a set of rotating photosphere models and then translating the covering fractions into ranges of potential stellar contamination in the transmission spectrum. 

WASP-19 displays rotational variability at the 2% level in the V band (Section <ref>). 
Our transit observations show that both cool (spotted) and hot (facular) active regions are present in the photosphere (Section <ref>). 
These active regions appear to be large, with radii of approximately 0.25 stellar radii, and display observed temperature contrasts of Δ T = -192 ± 10 K and Δ T = +137 ± 10 K, respectively. 
<cit.> also observe WASP-19b transit with a spot-crossing event that displays a larger temperature contrast (Δ T = -680 ± 50 K) and smaller radius (R_s=0.1651 ± 0.0045 stellar radii) than what we observe. 
Accordingly, we consider both smaller, high-contrast spots and larger, low-contrast spots in this analysis. 
In either case, as spots and faculae contribute approximately opposing variability signals, they can mask each other's presence in photometric monitoring <cit.>. 
However, they also contribute nearly opposing stellar contamination signals and, therefore, the net contamination of the transmission spectrum arises from the relative overabundance of one signal. 
For this reason, we considered the effects of spots and faculae separately in this analysis, deriving end-member examples to examine the range of possible contamination signals. 
In total, we examined three possible sources for the rotational variability of WASP-19:  high-contrast spots, low-contrast spots, and faculae.

For each case, we constructed 100 model photospheres, using PHOENIX model stellar spectra <cit.> of different temperatures to simulate contributions from the immaculate photosphere, cool spots, and hot faculae ignoring, for simplification, the impact 
of the position of the heterogeneities in the disk on the contrasts <cit.>. 
Utilizing constraints from the crossing events analyzed in Section <ref> and by <cit.>, we fixed the stellar parameters to the values in Table <ref> and interpolated within the PHOENIX model grid to generate component spectra with these parameters (namely, F_λ, phot, F_λ, spot, and F_λ, fac).
We successively added active regions at random coordinates drawn from a uniform distribution and recorded the corresponding rotational variability. 
We find the mean and 68% confidence intervals for the high-contrast spots, low-contrast spots, and faculae covering fractions consistent with the observed variability to be f_spot,1 = 2.0^+2.4_-0.7%, f_spot,2 = 10^+30_-5%, and f_fac = 19^+31_-10%, respectively (Figure <ref>). 
The similarity of the estimates for the low-contrast spots and faculae owes to their similar temperature contrasts. 
High-contrast spots, however, generate larger variabilities and therefore can reproduce the observed variability of WASP-19 with a smaller covering fraction.





If present in the unocculted stellar disk, these active regions will not be directly evident in transit light curves, but will nonetheless impact transmission spectra by imparting a spectral difference between the disk-integrated stellar spectrum and the resolved region illuminating the exoplanet atmosphere during the transit (the actual light source used in transmission spectroscopy). Figure <ref> shows the stellar contamination spectrum ϵ_λ produced by the spot and faculae, calculated as

    ϵ_λ, het =
    1/1 - f_het(1 - F_λ, het/F_λ, phot),

in which the subscript “het” refers to either spots or faculae <cit.>. 
The stellar contamination spectrum combines multiplicatively with the planetary transit depth to produce the observed transmission spectrum. 
In this case, low-contrast spots and faculae can have the largest effect on the observations; high-contrast spots have a comparably smaller effect given the lower allowed covering fractions. 
Low-contrast spots can increase transit depths at the shortest wavelengths of the IMACS observations (4500–4700 Å) by 1.9%, considering the mean covering fraction estimate, or as much as 8.4% for the 1σ upper estimate. 
Faculae, on the other hand, can decrease transit depths at these wavelengths by 3.1% on average and as much as 7.7% at 1σ confidence. 
Taking the transit depth of WASP-19b to be 20345 ppm, these values correspond to absolute transit depth changes of 390 ppm and -630 ppm for the mean spot and faculae coverages, respectively, and 1700 ppm and -1600 ppm for their upper limits. 
Thus, at 1σ confidence, we estimate that stellar contamination can alter the transit depth of our bluest data by at most 1700 ppm. 



Interestingly, as outlined in our analysis in Section <ref>, we do see changes in the transit depths between different epochs on the order of these effects. Even more, the slope observed in the 17/02/11 dataset can, in fact, be explained as arising from stellar contamination from unocculted (low-contrast, cold) spots with a spot coverage within 1-sigma of what is predicted in Figure <ref> (see Section <ref>). Another possibility to explain that slope is from occulted bright spots: if not taken into account in the modelling, they could produce a slope similar to the one observed. However, there is no evidence in the white-light light curve for occulted spots (either bright or cold) on the 17/02/11 night. This would mean that any occulted bright spot is below the noise level in our white-light light curves, and thus we cannot confirm this possibility with our data.  



§ DISCUSSION




 §.§ Interpreting the Magellan/IMACS optical transmission spectra



As presented in Section <ref> and Figure <ref>, our combined optical transmission spectrum from our Magellan/IMACS observations lacks spectral features at the couple of hundred parts per million level. 
This is in striking contrast with the transmission spectrum presented in <cit.>, which shows a 
strong scattering slope and some features interpreted as TiO in that work. 
In Figure <ref> we compare our mean-subtracted spectrum with those obtained by both <cit.> and <cit.>; note the caveat, however, that the data of <cit.> for wavelengths ≲ 6,000 Å, between 6,000-8,000 Å, and ≳ 8,000 Å were taken in different epochs (and, therefore, possibly under different levels of stellar activity). 
It is interesting to see that the precision of the Magellan/IMACS spectrum is similar to that of the HST spectrum <cit.>, although with better wavelength resolution. 
However, the data from <cit.> attains similar precision to ours but at an apparently better resolution; it is unclear exactly why this is so, although it might be due to our different binning strategies[For example, in <cit.>, where an early study of one of the WASP-19b datasets was presented, the authors 
binned the data in 200-Å bins, with separations between the central 
wavelengths of each bin of 100 Å. In principle, in that case 
each datapoint does not provide independent information on any possible features 
to be observed in the transmission spectrum, and for this reason we avoid performing 
such “overlapping" binning of the data.
Of course, this gives rise to larger 
errorbars, but it avoids having correlations between adjacent wavelength bins 
in the transmission spectrum.].

Figure <ref> shows that the scattering slope bluewards of ≲ 5500 Å observed in <cit.> is not evident in either our combined transmission spectrum or that of <cit.>. 
Interestingly, however, our February 2017 observations do show a decrease on the same order of magnitude as a function of wavelength of the transit depth, which could be due to stellar heterogeneities in the star (see Section <ref>). 
It is thus possible that the observed slope in the data of <cit.> is actually produced by starspots/faculae. 
This highlights the power of repeatability: 
especially in ground-based optical observations, repeated observations in the same wavelength range (and, if possible, also with different instruments/setups) are mandatory not only to increase the precision of the measurements but also their accuracy, which in turn allows us to confirm that the features observed in transmission come indeed from planetary phenomena and not from systematic effects, either of stellar or instrumental nature. 
In addition to the absence of the scattering slope, we also do not observe clear signatures of sodium and potassium. 
Regarding TiO, by simply inspecting Figure <ref> we can observe that the only possible sign of TiO absorption is with our datapoints at 7,000 Å and 7,200 Å. 
If there is a TiO feature, it is definitely smaller than (approximately half) the one presented in that work.



  §.§.§ Quantifying possible spectral signatures with atmospheric retrievals





To provide a quantitative measurement on the absence or presence of features, we performed an atmospheric retrieval on our data, simultaneously modeling the atmosphere of the planet and the spectral signal of any heterogeneities present in the star. 
For the former, we used the semi-analytical formalism for transmission spectroscopy on an isothermal, isobaric atmosphere presented in <cit.>. 
For the latter, we used the formalism of <cit.>, 
already presented in Section <ref>. 

A detailed overview of the implementation of our retrieval framework can be found in Appendix <ref>. 
In summary, as with our spot-crossing analysis detailed in Section <ref>, we implement this framework using PyMultiNest <cit.>, which allows us to compute bayesian evidences for the different possible models, allowing us to (1) quantitatively measure the evidence for different features in the atmosphere of the planet and (2) incorporate our ignorance on a “best model" for the transmission spectra in the posterior distribution of the retrieved atmospheric parameters <cit.>. 
Our full retrieval has initially 3 parameters of stellar origin: a temperature of the occulted (by the planet) stellar 
surface T_occ, a temperature of the heterogeneous, unocculted surface of the star T_het, and a fraction f_het, which defines the fraction of the projected stellar disk covered by T_het (see Section <ref>). 
The planetary atmosphere is defined by 5 + n parameters: f, a factor that scales the derived planetary radius in Section <ref> in order to find the reference radius R_0=fR_p at which the atmosphere is optically thick, P_0 the pressure at which the atmosphere is optically thick (which we could also interpret as a cloud-top pressure), two parameters that define our haze prescription <cit.>, the atmospheric temperature T, and n mixing ratios of the different elements considered in this work (Na, K, TiO, H_2O, CO_2, CO and CH_4). 
Using this framework, we now set to analyze the individual and combined transmission spectra of WASP-19b obtained with Magellan/IMACS.



  §.§.§ Analysis of the individual Magellan/IMACS transmission spectra


We first analyzed each Magellan/IMACS optical transmission spectrum obtained 
in different epochs separately in order to test them individually for possible 
stellar contamination and/or planetary spectral features using our retrieval framework. 
n_live=2000 points were used to explore the parameter 
space, and models were computed with and without heterogeneities as well as 
with and without hazes. For these retrievals, only Na, K, TiO and H_2O 
were considered as possible opacity sources in the planetary atmosphere, and 
all combinations of those elements were tried in order to see which model 
was preferred by the data. A simple flat line (i.e., a constant transit depth 
as a function of wavelength, which would be the expected signature for high 
altitude clouds at the resolutions and precision probed by our data) was also considered, as well as spectra purely dominated by stellar 
contamination (i.e., a flat planetary transmission spectrum multiplied 
by equation (<ref>)). 

As qualitatively inspected in past sections, every dataset but the 17/02/11 dataset was consistent with models that neglected stellar contamination. 
Of those datasets, every dataset but the one obtained on 14/03/22 has the simple flat line as the model with the largest evidence (i.e., as the “best-fit"), with all the other models being indistinguishable from it (i.e., all the atmospheric retrieval models have Δln Z < 2 with respect to that model). 
For the 14/03/22 dataset, there is a slight preference for models incorporating TiO (Δln Z = 2.4 in favor of the best TiO model compared to a flat line). 
However, the retrieved temperature for those models was unphysical (1,300 K) given WASP-19b's large equilibrium temperature, and thus we attribute this preference to statistical fluctuations rather than to a physical reason. 
Given our data, Na and K appear to not be present in the individual Magellan/IMACS optical spectra.

The 17/02/11 dataset is a very interesting one as the best-fit model, which is shown in Figure <ref>, includes stellar contamination, clouds and TiO, but no hazes in the planetary atmosphere. 
This model is clearly preferred over a simple flat line (Δln Z = 18.4) and is slightly preferred over models that do not include a stellar contamination component with or without hazes and with or without TiO (Δln Z > 2.7). 
The model is, however, indistinguishable over a spectrum dominated purely by stellar contamination (Δln Z = 1.3). 
The caveat with this latter model is that in order to reproduce the observed transmission spectrum it requires that both the occulted and non-occulted parts of the photosphere have temperatures < 4,500 K. On the other hand, the best-fit model shown in Figure <ref> has parameters whose posterior distributions are clearly bi-modal. 
One mode covers the same parameters as the just-described stellar contamination model: 
occulted and non-occulted temperatures having temperatures < 4,500 K and, because those can match the data perfectly well without the contribution from a planetary atmosphere, low TiO volume mixing ratios (log_10TiO<10^-10). 
The other mode, however, has a complex structure. For the stellar contamination parameters, it gives
f_het=0.58^+0.24_-0.21, 
T_het=5402^+198_-230 K and 
T_occ=5780^+148_-215 K. 
In this case, the stellar surface is consistent with WASP-19's effective temperature and, thus, it considers the scenario in which the average temperature of the occulted stellar photosphere by the planet is hotter than the average disk temperature.
Here, 58^+24_-21% of the stellar surface has a temperature T_het which is consistent with WASP-19's effective temperature. 
Interestingly, the temperature of the occulted part of the star by the planet in this case is, in fact, consistent within the errors with the bright spot temperature derived in Section <ref>; 
the areal covering fraction by this temperature of 1 - f_het = 42^+24_-21 % is, in fact, also consistent within 1-sigma with the amount of variability observed in WASP-19b, according to our analysis in Section <ref> (see Figure <ref>). 
As for the parameters on the planetary atmosphere in this latter mode, in this retrieval the planetary temperature is T=1951^+542_-564 K — very uncertain but consistent with WASP-19b's equilibrium temperature — and 
the log-mixing ratio of TiO is log_10TiO = -6.50^+2.62_-3.26; 
uncertain but consistent with a solar-composition atmosphere. 

It is intriguing that the only dataset in which TiO is apparently detected is the 
one dataset in which there is a clear stellar contamination component, 
which was detected thanks to the large wavelength coverage of Magellan/IMACS. 
It might be that our TiO detection is, in fact, a spurious detection due to 
an oversimplification of our modelling of the stellar component present 
in the observed transmission spectrum and/or occulted spots in our data 
<cit.>. We tried a more complex parameterization of the 
stellar surface in which we consider both spots and faculae with different
disk-averaged and transit-chord covering fractions, as detailed by <cit.>, but this only added complexity to 
the model without a further improvement of the fit: this more 
complex model has, in fact, an evidence Δln Z = 7 smaller than 
the simpler model mentioned in the previous paragraph. It might be that there are 
stellar effects mimicking the TiO in this dataset that we are not considering (e.g., 
distribution of spots with different sizes and temperatures, including spots 
with temperatures cooler than 4,500 K that could mimic atmospheric absorption by TiO in 
the transmission spectrum); a thorough search for 
such effects is, however, beyond the scope of the present study.





  §.§.§ A featureless combined optical transmission spectrum from Magellan/IMACS


We now turn to the study of the combined optical transmission spectrum already presented at 
the beginning of this section. We ran the same retrievals as those for the 
individual transmission spectra just discussed in order to search for possible features 
in this more precise combined dataset. A featureless optical transmission spectrum (i.e., a 
constant transit depth as a function of wavelength) is the favored model given the data, 
with a log-evidence of ln Z = -200.6; all the models including any one opacity source 
have larger evidences. In fact, the nominal model found by <cit.> in their 
VLT/FORS2 dataset (which includes a haze, water, Na and TiO) has a log-evidence of 
ln Z = -204 which means that, given our data, this latter model is 30 times less likely than the 
simpler flat transmission spectrum model in the optical. In summary, we find no evidence in our data 
for hazes and/or TiO absorption, and thus we conclude that, given our data, the transmission spectrum 
in the optical wavelength ranges probed by Magellan/IMACS is flat and consistent with high altitude clouds.
This conclusion is in agreement with predictions from the analysis of the water feature of <cit.>, the interpretation of <cit.>, and the lower-resolution HST/STIS optical transmission spectrum of <cit.>. 



 §.§ The complete optical-to-IR transmission spectrum of WASP-19b



We now join our optical transmission spectrum from Magellan/IMACS with high-precision space-based observations to see what further constraints we can place on the interpretation of the complete optical-to-IR transmission spectrum of WASP-19b.
In addition to our optical observations, we consider the data obtained with HST in the optical (through HST/STIS) and IR (through HST/WFC3) by <cit.> and <cit.> and the re-analyzed 3.6 and 4.5 μm Spitzer data for WASP-19b presented in <cit.>. We chose to use the data 
published in <cit.> and, e.g., not the analysis of the HST/WFC3 dataset presented in <cit.>, as the data in <cit.> for HST/WFC3, HST/STIS and Spitzer has been 
analyzed uniformly. We perform the same retrievals as in the previous subsection but now including the most likely infrared absorbers (H_2O, CO, CH_4 and CO_2) and three additional parameters in order to account for three possible offsets in transit depth: 
one for our optical observations, one for the optical HST/STIS observations and one for the IR observations with HST/WFC3. 
These offsets are fitted in our retrieval in order to take into account possible overall mismatches due to stellar activity, in light of our discussion in Section <ref> regarding the corrections made to the HST observations. 
To apply them, we mean-subtract those datasets, add the white-light transit depth (R_p/R_*)^2 calculated with the planetary and stellar radii value found in Section <ref>, and consider constant offsets with gaussian priors centered in zero with wide standard deviations of 1000 ppm between the datasets. 
No offset is applied to the Spitzer data, as the impact of stellar activity at those wavelengths is expected to be negligible given the precision of the data and our analysis in Section <ref>. 
Impacts of using different orbital parameters for the transit fit between these datasets and the Spitzer one are assumed to be negligible, as the ones used in <cit.> and the ones used in this work have only small differences, with any impact on the transit depths being well within the Spitzer errorbars. 

There are many tens of retrieval models that are statistically indistinguishable between one another 
(Δln Z < 2), almost all of which containing water as a common opacity source, no hazes, and 
no stellar heterogeneity component. 
The simplest model, and the one with the lowest evidence is a transmission spectrum purely dominated by clouds, no haze, no stellar contamination component and water absorption (ln Z = -304.5), which is the retrieval presented in Figure <ref>. 
It is interesting to note that this model is statistically indistinguishable from a flat line (Δln Z = 1.2 in favor of the retrieval with water absorption), which highlights the fact that more data are needed in order to claim a strong water absorption detection in the infrared given the unknown transit depth offset between 
the optical and infrared data and the relatively large errorbars of the HST/WFC3 measurements. 



In order to find the posterior distributions for the parameters of interest in our atmospheric retrievals, we perform bayesian model averaging <cit.> using the posterior distribution of all of our retrievals and our estimated evidences, assuming each model is equiprobable a-priori.
In this way, we get the marginalized (over our models) posterior distributions for all the parameters in our retrievals.

Figure <ref> presents the most important parameters for our discussion. 
As can be seen, we obtain an uncertain atmospheric temperature of 2048^+490_-540 K, whose distribution is peaked around the expected one from WASP-19b's equilibrium temperature of 2100 K (calculated using perfect day-night heat redistribution).
We also obtain an uncertain log-pressure (in bars) for the optically thick part of the atmosphere of log_10 P_0 = -0.50^+3.79_-2.81. 
The posterior distributions of both parameters are very interesting. 
On one hand, the retrieved terminator temperature is consistent with WASP-19b's day-side temperature estimated by <cit.> to be of 2372 ± 60 K. 
The retrieved distribution of the pressure P_0, on the other hand, is also interesting as it seems to be peaked towards small values (around 1 mbar), which we can interpret as the data being most consistent with high altitude clouds being present in WASP-19b. 
This interpretation is consistent with the one given by <cit.>, who also retrieve high altitude clouds for WASP-19b with a different retrieval technique. 

We also obtain a very uncertain log-water volume mixing ratio of 
logχ_H_2O = -5.83^+3.30_-4.37. 
The mixing ratio of water under chemical equilibrium at 2048 K (our median retrieved atmospheric temperature) is 7.8 × 10^-4 at pressures of 0.3 bar (our median retrieved cloud-top pressure).
This is well within our posterior distribution and hence makes our retrieved water abundance consistent with that from a solar composition atmosphere. 
The weak constraint on the water mixing ratio from our atmospheric retrieval is most likely due to three effects: (1) as already discussed, the uncertain offset due to stellar activity between the infrared and optical observations; 
(2) the fact that the water signature in the infrared is not very strong, as already discussed in the previous paragraph; and 
(3) the correlation between R_0, P_0 and the mixing ratios <cit.>. 
More data in this wavelength regime (or at any other water band in the infrared) might help put 
a tighter constraint on the water abundance in WASP-19b's terminator region. 

For TiO, we obtain a log-volume mixing ratio of logχ_TiO = -20.34^+6.83_-6.19, which is also very uncertain, but with a distribution that clearly rules out the abundance expected in a solar-composition atmosphere. 
In this case, the expected abundance of TiO is 2.2 × 10^-7 at the same pressures and 
temperatures as before. 
The probability given our data that the mixing ratio has this or larger values is ℙ(χ_TiO≥ 2.2× 10^-7|D) = 1%.
Thus, in the framework of our semi-analytical retrievals, solar or super-solar TiO abundances in the terminator region of WASP-19b are very unlikely. 

For Na, a similar scenario arises: 
the expected solar composition abundance is 5.4× 10^-6. 
The corresponding probability of the mixing ratio having this or larger values is 
ℙ(χ_Na≥ 5.4× 10^-6|D) = 2%, which also suggests that solar and super-solar Na abundances in WASP-19b are very unlikely. 

Additionally, as a final note on the retrieved parameters, we find that hazes are not clearly preferred over simpler, non-hazy models.

Overall, our retrieved water abundance seems to be consistent with that of <cit.>, 
but the abundance of TiO, the abundance of Na, and the presence of a haze are not. 
The latter discrepancies result from the fact that we do not detect their spectroscopic signatures in 
the optical in our combined Magellan/IMACS data, which is interesting and 
could perhaps be highlighting the possibility of stellar 
activity mimicking some or all of those features, as discussed in Section <ref>. The exact cause of the discrepancy between our 
combined transmission spectrum and that of <cit.>, however, remains unknown.



 §.§ Giant spots on WASP-19

We now discuss what may be gleaned regarding the spots on WASP-19 from our observations. 
In Section <ref> we reported sizes and temperatures for our observed spot-crossing events on the April 2014 and 2017 nights, with the former having a spot size of R_s = 0.21 ± 0.01 stellar radii and a temperature Δ T = 192± 10 K colder than the stellar photosphere, and the latter having a size of R_s = 0.25^+0.05_-0.04 stellar radii, and a temperature 
Δ T = 137 ± 10 K hotter than the stellar photosphere. These sizes are very large, especially compared with the Sun, whose rarest, largest spots can reach ∼ 0.1R_* (see, e.g., sunspot AR 2192, which is the largest in the last 22-year solar cycle). 

Spots this size are, however, apparently not very uncommon on WASP-19.
Spots of sizes similar to the ones we obtain in this work were observed both in 2010 <cit.> and 2012 <cit.>. 
Interestingly, however, the temperature of the detected spots apparently varies greatly. 
<cit.>, for example, report a temperature Δ T = 683 ± 103 K cooler than the stellar photosphere for the largest spot observed in that work (R_s=0.1651 ± 0.0045 stellar radii). 
This is significantly cooler than the Δ T = 192± 10 K we observe for our cool, R_s = 0.21 ± 0.01 stellar radii spot in 2014. 
This is very interesting because, on average, spots producing larger contrasts (i.e., having colder temperatures than the surrounding stellar photosphere) tend to be larger in our own Sun <cit.>, which is the opposite to what we see between the spots in this work and the one in <cit.>. 
In fact, if we compare the broadband sizes and contrasts quoted in the works of <cit.> and <cit.> with ours, it appears the two largest starspots <cit.> have the smaller contrasts, while the two smaller starspots 
<cit.> have the largest contrasts. 

There are various possibilities for explaining this behavior. 
One of them is that the relationship between size and contrast for the Sun found by <cit.> has a large dispersion.
However, if we use the same relation for WASP-19 we cannot explain the observed variation: 
for large spot areas, the dispersion in contrast seems to be of order ∼ 0.05, while the variation in contrast between our spots is of order ∼ 0.20. 
Considering this observation only, the simplest explanation is that the size-contrast relation for our Sun does not apply for more active stars like WASP-19. 
On the other hand, derived spot contrasts could also be affected by unocculted stellar spots, leading us to make direct comparisons between contrasts obtained at different epochs very difficult. 
This effect has already been pointed out for measurements of absolute transit depths in transit light curves in <cit.> and many other works, but here we argue that it could also impact  measurements of starspot contrasts, which depend on the relative dimming and brightening of the transit light curve during starspot crossing events[However, the spot size, which mostly depends on the start and end of this dimming event, is not impacted by this effect.]. 
In addition, it could also be that we are actually observing a collection of smaller spots within an active region and not giant spots themselves <cit.>; 
however, this hypothesis can only be tested with observations of spot-crossing events with better time resolution and photometric precision.



§ CONCLUSIONS


In this work, we have presented six optical transmission spectra of WASP-19b obtained in different 
epochs with the IMACS multi-object spectrograph mounted at the Magellan Baade 6.5m telescope, 
three covering the 4100-9350 Å range, two covering the 4500-9350 Å range and one covering 
the 4100-5560 Å range. All but one spectrum covering the 4500-9350 Å range are consistent 
with being flat, in striking contrast with the previous work of <cit.>, who report 
a large blue optical slope and signatures of TiO, Na and H_2O in their optical transmission spectra. 
One of our transmission spectra, however, shows a clear decrease in transit depth towards 
longer wavelengths, which we interpret as arising mainly from stellar contamination. This is the only 
dataset in which, apart from the stellar contamination signature, TiO is tentatively detected, which 
is intriguing as it is also the dataset that shows the clearest signatures of being impacted by 
stellar contamination. Whether there is a causal relation between the two (either the TiO feature 
being from stellar origin or stellar activity inducing TiO features in WASP-19b's terminator region) 
is unknown. We combine our optical measurements along with previously published HST and Spitzer data 
in order to interpret WASP-19b's panchromatic transmission spectrum and use a semi-analytical retrieval 
approach to find a water abundance consistent with that from a solar composition atmosphere, along with 
sub-solar TiO and Na abundances in WASP-19b's terminator region.

In addition to our transmission spectroscopy study, we also study spot-crossing events that occurred 
in two out of the six observed transits, one of which we interpret as being caused by the crossing 
of a bright spot, making this the first unambiguously detected bright spot-crossing 
event from an exoplanet transit light curve <cit.>. Both spots are very large (∼ 0.2R_*); however, 
this is apparently not uncommon on WASP-19b, especially given that spot-crossing events are more likely 
to be detected for these kind of giant spots. Using the wavelength-dependent light curves, we 
were able to constrain the temperatures of the spots using spectra from ATLAS model atmospheres <cit.> to model the fluxes of both the stars and the 
spots. For the bright spot, we find it to be 
Δ T = 137 ± 10 K hotter than the stellar photosphere, while for the cold spot we find 
it to be Δ T = 192 ± 10 K cooler than the stellar photosphere. We use this information in 
order to give plausible scenarios for stellar contamination in WASP-19b's transmission spectrum 
using the work of <cit.>, which in our case was fundamental to study in order to have 
a complete interpretation of the features observed by this and previous works on the transmission 
spectrum of WASP-19b.

Given our analysis and discussion, we conclude that more data are needed in the infrared in order 
to secure the HST/WFC3 water feature reported in <cit.>. Such data will be fundamental 
to understand whether WASP-19b's metal content, as well as to know if the atmosphere has a low or high C/O ratio, which in turn could give 
clues as to the possible formation pathways that led to WASP-19b present-day spectra <cit.>.




§ ACKNOWLEDGEMENTS

We would like to thank an anonymous referee for suggestions and comments that 
improved this manuscript. N.E. would like to thank J. Araya, J. Bravo, P. Jones, M. Martinez, V. Meriño, M. Navarrete, H. Nuñez, A. Pasten and G. Prieto for numerous discussions (and subsequent experiments) done at Las Campanas Observatory (LCO) which were fundamental to better understand the many sources of systematics in the instrument, and all the staff at LCO for their invaluable support which made this work possible. N.E. would also like to thank R. 
McDonald for sharing the Na and K cross-section calculation method, and 
K. Heng for fruitful discussions regarding the semi-analytical formalism 
for transmission spectroscopy. N.E. would like to thank the Gruber Foundation for its generous support to this research. N.E. and A.J. acknowledge support from the Ministry for the Economy, Development and Tourism Programa Iniciativa Científica Milenio through grant 
IC 120009, awarded to the Millennium Institute of Astrophysics.  A.J. acknowledges support by Fondecyt grant 1171208 and by CATA-Basal 
(PB06, CONICYT). N.E. acknowledges support from Financiamiento Basal PFB06 and from the CONICYT-PCHA/Doctorado Nacional graduate fellowship. 
B.R. acknowledges support from the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1143953. 
D.A. acknowledges support from the Max Planck Institute for Astronomy, Heidelberg, for a sabbatical visit. We 
thank the anonymous referee for useful comments that improved this manuscript. The results reported herein benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory (LCO), Chile.








mnras
 



§ PHOTOMETRIC MONITORING OF WASP-19




 §.§ 2014 season


For the 2014 season, WASP-19 was monitored with the SMARTS 1.3m telescope at CTIO with a V filter. This project was a continuation of the CTIO photometric monitoring of WASP-19 presented in the 
work of <cit.>. The data were downloaded from the NOAO archive[<http://archive.noao.edu/>], 
and reduced with a standard photometric pipeline that performs aperture photometry on the stars using the  package <cit.>. The five stars closest in brightness to WASP-19 
were used to obtain the relative photometry of WASP-19. The resulting photometry is shown in Figure <ref> along 
with that obtained for a comparison star close in brightness to WASP-19, in order to illustrate the long-term photometric precision and stability achieved by these observations.



Although the formal precision (∼ 0.3%) should have allowed us to 
detect variability at levels ≲ 1%, the photometry for both WASP-19 and the comparison star vary significantly 
more than the formal errors computed assuming photon noise (up to 2% in some portions of the light curve) and thus are of insufficient quality to observe the rotational modulations of the star. The conditions in this 
season were similar to those of the 2013 season shown by <cit.>, in which precisions of the same order of magnitude 
were reported. Because of this, we decided to re-reduce the data published in that work which is also available from the 
NOAO archive in order to test for any problems with our reductions and/or methods. Doing this we recovered the same variability 
and precision reported in that work. However, we observed the same qualitative behaviour as that shown in Figure <ref>: 
the variability is present not only in WASP-19 but also in the comparison stars. This implies that the observed variability 
reported in <cit.> is not entirely of astrophysical origin but, rather, mostly due to a strong instrumental component 
like that observed in Figure <ref>. This large instrumental scatter is most likely due to 
the fact that the nightly observations were obtained by placing WASP-19 and the comparison stars at different positions in the detector every night. 
Different methods were implemented in order to try improve the photometry, but none of them gave good results.

We describe this dataset in detail because 
<cit.> corrected the HST transmission spectrum assuming the 
observed variability was only astrophysical in origin, i.e., due to the rotational modulation caused by stellar spots. However, as we 
have shown, the variability in the CTIO dataset is mostly spurious.
Therefore, the correction to the transmission spectrum done in <cit.> is not appropriate. This correction mostly adds an overall shift in the transit depth, because the effect as a function of 
wavelength is smaller than the optical errorbars published in that work and it has a negligible impact on the infrared. 
We discuss this in the following sections when deriving and comparing the optical transmission spectrum of WASP-19b.



 §.§ 2017 season


For the 2017 season, we used the photometry published by the ASAS-SN, whose precision suffices for the purpose of detecting stellar activity of 
relatively bright stars such as WASP-19. 
In Figure <ref> we show the photometry and we also indicate the epochs in which the observations with 
Magellan/IMACS where performed. In order to analyze this photometry, we fitted a Gaussian process regression using 
 <cit.>. In particular, we use a kernel defined by

    k (τ) = B/2+Ce^-τ/L[cos(2πτ/P_rot) + (1+C)],

where τ = t_i - t_k, with i,k ∈ [1,2,...N], where N is the number of datapoints, and B, C, L and P_rot 
are the hyperparameters of the model, with the latter corresponding to the period of the quasi-periodic oscillations defined by 
this kernel. The usage of this kernel and its application to the recovery and modelling of stellar activity has been discussed already
in <cit.>; we refer the reader to that work for the details. However, we briefy mention here that this kernel has the key 
properties of any kernel used to model stellar activity: it is quasi-periodic (due to the cosine term), it captures additional long-term trends in 
the data not captured by quasi-periodic modulations (due to the multiplicative exponential term), and it has the key property of being positive 
for all values of τ. 



In order to fit this Gaussian process to the data, we first obtained the maximum-likelihood estimates of the parameters through a simple optimization using 's  function. With these best-fit parameters as input, we ran a 
Markov Chain Monte Carlo (MCMC) using the  ensemble sampler <cit.>, by using 100 walkers around this 
best-fit value, and letting them explore the parameter space for 200 steps, with 100 steps as burn-in, and wide priors for the parameters being optimized. This resulted in 100,000 samples from the posterior density. The posterior probability 
distribution obtained for the rotation period (P_rot) is shown on the lower-right panel of Figure <ref>. 
The posterior value is P_rot = 10.44± 0.07 days, which is in excellent agreement with 
the 10.5-day period reported in <cit.> and in <cit.>.

The flux obtained and modelled for WASP-19 with the ASAS-SN photometry is relative. However, it is desirable to find an estimate of 
the absolute flux, taking into account 
the fact that the maximum brightness does not correspond to an unspotted disk in general. In order to take this into account, we obtained an estimate of the absolute flux level of the star in the observed band using the approximate method outlined 
in <cit.>, which divides the relative flux measurements by F_* = F_max + σ_F, where F_max is the maximum observed brightness 
of the star and σ_F is the standard deviation of the light curve. We used the retrieved Gaussian process conditional on our data in order to obtain F_*; this predicted photometry 
is shown in the upper and lower-left panels of Figure <ref>. Our Gaussian process regression returns an estimate for the absolute 
flux in the V-band for WASP-19 relative to an immaculate star of 
0.972± 0.003 for the 17/02/11 observations, of 0.964 ± 0.003 for the 17/04/04 observations and 
0.966± 0.003 for the 17/04/12 observations, i.e., a fairly similar flux level for our observations made in April, but a slightly brighter 
flux level on our February observations. It is interesting to note, however, that during the February observations, the ASAS-SN photometry seems to vary significantly more around the time of the observations.



§ UPDATING THE TRANSIT EPHEMERIDES FOR WASP-19B



In order to update the ephemerides of WASP-19b using our measured transit mid-times 
we first transform our times in Table <ref> which are in BJD UTC to BJD TDB by adding 32.184 +N seconds to our BJD UTC times, 
where N is the number of leap seconds elapsed since 
1961[Obtained from the  column in <http://maia.usno.navy.mil/ser7/tai-utc.dat>], 
which were 35 for our 2014-2015 season and 37 for the 2017 season <cit.>. We couple those BJD 
TDB times with the ones in <cit.>, and derive a new ephemerides (in BJD TDB 
units) of

    t_0 = 2454775.337777(42) + 0.788839316(17)E,

where, as in <cit.>, E is the number of cycles elapsed since 
the first transit observed by <cit.>. This new ephemerides is both 
consistent with and more precise than the one derived in <cit.>. 
We observe, however, that there is a systematic difference of around 40 
seconds between the timings on our 2014-2015 observations and the 2017 
observations (see Figure <ref>). Given that these have errors between 4–10 seconds, the 
timing difference seems to be significant. A full analysis of the entire 
transit timing dataset, which is beyond the scope of this work, could elucidate
what could be producing these changes in transit time, which is also seen in 
previous epochs <cit.>.



§ SPOT-CROSSING EVENTS ANALYSIS


The analysis of the spot-crossing events described in this work were separated in two 
steps. On the first step, we analyzed the band-integrated (“white-light”) light curves 
in order to obtain the size of the spots, and then this parameter was used in order to 
fit the wavelength-dependent light curves with which we were able to obtain wavelength-dependent 
constrasts for the spots.

For the analysis of the band-integrated transit distortions, we gave as inputs to  the systematics-corrected light curves along with the priors on the transit parameters obtained in Section <ref>, and let n_live=500 points explore the parameter space, with which the log-evidences were calculated with errors below 0.1% for spotless, one-spot and two-spot models, and with which we were able to generate 3500 samples from the posterior distribution for each model. The results for our “best" models (in terms of the calculated evidences) are presented in Figure <ref>. 

For the April 2014 dataset, the spot models are favored with very high significance, as expected from the observed light curve distortions (Δln Z = 163 in favor of the one-spot model when compared to the spotless model which, in the case of both models being equally likely a priori, would imply that the model with spots is (Z_spot)/(Z_no-spot)= 6 × 10^70 times 
more likely than the non-spotted model). However, the two-spot model is not favored by our data: the difference in log-evidence is actually higher for the one-spot model, with Δln Z = 2 in favor of this simpler model. The resulting one-spot model led to a very precise measurement on the position of the spot on the (projected) stellar surface of (x,y)=(0.21^+0.01_-0.01,0.68^+0.04_-0.06), where x goes in the direction almost perpendicular to the planet's motion and the center of the star is at (0,0), with contrast of F_spot/F_star = 0.84^+0.01_-0.03, i.e., a dark spot, and a considerably large radius of R_spot = 0.21^+0.01_-0.01 stellar radii. For the April 2017 dataset, the spot model is also favored with high significance, with a Δln Z = 58 between the spotless and one-spot models in favor of the latter. Again, the two-spot model is not favored by our data, with a Δln Z = 2 in favor of the simpler, one-spot model. The position of the spot is also identified with high precision to be at (x,y)=(0.41^+0.02_-0.02,0.43^+0.08_-0.08) in the stellar surface, with a size of 0.25^+0.05_-0.04 stellar radii. The most interesting feature of the April 2017 feature, however, is its contrast: the retrieved spot contrast is F_spot/F_star = 1.16^+0.05_-0.04, which implies the feature is actually a bright spot. This makes the April 2017 event, thus, one of the first unambiguously detected bright spot features, i.e. a facular or plage region, on an exoplanet host star <cit.>, and the first one on WASP-19.

In order to extract the spot contrasts as a function of wavelength, we performed the same analysis on the systematics-removed, wavelength-dependent light curves, following a similar approach to the one defined in Section <ref>: we fix the orbital parameters to the joint best value found in Section <ref>, the position of the spot on the stellar surface and the spot radius, and let the transit depth, the limb-darkening coefficient of the linear law and the spot contrast to vary as a function of wavelength on each light curve. Once all the contrasts were obtained for both datasets, we used ATLAS stellar model atmospheres <cit.> in order to constrain the temperature of the spots using the wavelength-dependent spot contrasts, by assuming that the spots produce the same flux as a star with the same properties except for the temperature. In order to perform this fit, model spectra for stars from 4000 K to 6500 K were used, with gravities and metallicities similar to WASP-19 according to the works of <cit.> and <cit.> (i.e., log g=4.5 and [Fe/H]=0.1). We interpolated these spectra in order to have a function that returned the shape of the spectrum for a given temperature. With this function, we fitted the contrasts by a two-parameter function F(T_s)/F(T_W19), where F(T) is the (integrated, in the respective wavelength bins) flux at a given temperature, T_s is the temperature of the spot and T_W19 is the temperature of WASP-19. For T_s a uniform prior is used between 4,000 K and 6,500 K, while for T_W19 we put a normal prior with the temperature measured by <cit.> for WASP-19, which is 5,460 ± 90 K. The fit then was free to explore the parameter space, with the constraint that for the cold spot T_s<T_W19 and for the bright spot T_s>T_W19. The resulting contrasts and their best-fit variations as a function of wavelength with this simple model are shown in Figure 
<ref>.   

The fits to the observed contrasts are excellent. For the cold spot observed in the April 2014 dataset, we find a temperature of T_s = 5278 ± 81 K (or Δ T = 192 ± 10 K colder 
than the star). The variation as a function of wavelength is highly significant in this dataset, with a Δln Z = 20 in favor of a variation in wavelength versus no variation (i.e., constant contrasts as a function of wavelength). For the bright spot, we find a temperature of T_s = 5588 ± 92 K (or Δ T = 137 ± 10 K hotter than the star), and find moderate evidence in favor of a variation as a function of wavelength with a log-evidence difference of Δln Z = 3 when compared with a constant contrast as a function of wavelength, which in the equiprobable case, would imply the variation as a function of wavelength is 20 times more likely than a contrast which is constant as a function of wavelength. 



§ SEMI-ANALYTICAL ATMOSPHERIC RETRIEVALS


Our implementation of the semi-analytical framework for transmission spectroscopy of 
<cit.> follows the equation that describes the radius as a function of wavelength in 
transmission of a transiting exoplanet, R(λ), which is given by

    R(λ) = R_0 + H[γ + E_1 +  ln(τ_0) ],

where

    E_1 = -γ - lnτ_0 + ∫_0^τ_01-e^-τ/τdτ.

Here H=k_BT/μ g is the atmospheric scale-height given the mean molecular weight, μ, the temperature, T, the planetary gravity, g, and Boltzmann's constant, 
k_B = 1.38×10^-16 ergs/K, τ_0 = (P_0 σ/kT) √(2π HR_0) is 
a reference optical depth, that is defined by the (temperature and 
wavelength-dependent[Note that the cross-section is also pressure dependent. However, due to the low ranges of pressures usually probed in transmission, and the high expected temperature of WASP-19b, we omit the pressure-dependence of the cross-sections in this work.]) cross section σ, a reference radius R_0 and a reference pressure P_0. γ = 0.57721 is the Euler-Mascheroni constant. For the cross section, we use

    σ(λ,T) = σ_haze + ∑_i=0^N χ_i σ_i(λ,T),

where χ_i is the volume mixing ratio of species i which has a cross-section 
σ_i(λ,T) and where σ_haze(λ) = aσ_0(λ/λ_0)^γ_haze is used to emulate 
hazes in the atmosphere; this is the same haze prescription used in <cit.>, with σ_0 = 5.31 × 10^-27 cm^2, λ_0 = 3500 Å and free parameters a and γ_haze. In this formalism, R_0 can be interpreted as a “surface" where the atmosphere is optically thick (i.e., where the opacity τ→∞), which can be produced either by the increased atomic/molecular absoportion deep in the atmosphere and/or 
to a cloud deck below which no radiation passes through; as such, clouds are naturally considered here. Similarly, P_0 can be interpreted as the 
pressure where this optically thick portion of the atmosphere lies, and so the transmission spectrum is actually generated 
by interaction of species at pressures P>P_0. We let both of these as free parameters in our retrievals; the pressure is explicitly fit, while for R_0 we 
fit for a factor, f, and set R_0 = fR_p, where R_p is the 
planetary radius found in Section 
<ref>. For any additional grey or close-to-grey opacity source, the parametrization 
of σ_haze(λ) covers various ranges of scattering sources, from 
grey opacities (γ_Haze→ 0) to Rayleigh scattering-like 
opacities (γ_Haze→ -4).

To model heterogeneities in the 
star we used the framework outlined in Section <ref>, i.e., we 
use equation (<ref>) to model any stellar contamination simultaneously with 
our transmission spectrum retrieval. This is parametrized in our model by 
three terms: a temperature of the occulted (by the planet) stellar 
surface T_occ, a temperature of the heterogeneous, unocculted 
surface of the star T_het, and a fraction f_het, 
which defines the fraction of the projected stellar disk covered by T_het. 
The variation as a function of wavelength of 
those portions of the stellar surface is modelled with PHOENIX model atmospheres, 
as described in Section <ref>.

We fix the planetary gravity to 1414 cm/s^2 which was the one derived in Section 
<ref> and as possible opacity sources in the atmosphere of WASP-19b, 
we consider Na, K, TiO, H_2O, CO_2, CO and CH_4. All cross-sections but that of Na and K are calculated using  <cit.>. The 
line-lists for H_2O, CO_2 and CO are obtained from the HITEMP database 
<cit.>, while those of CH_4 are obtained from the ExoMol database <cit.>. TiO is synthesized using the latest (2012 version) 
line lists calculated by B. Plez, obtained from the Vienna Atomic Line 
Database <cit.>. For Na and K, we use the analytical lorentzian-profiled doublet shapes used in 
<cit.>. All the cross sections are calculated at 1 mbar — low 
enough as for pressure broadening effects to be negligible, effects that are small and that we would not be able to detect in this work. We consider the atmosphere 
to be dominated by molecular hydrogen, with a mixing ratio of H_2 of 0.83 
and He of 0.17 (which is based on chemical equilibrium calculations). 

As explained in Section <ref>, our full retrieval has initially 3 parameters of stellar origin (T_occ, T_het, f_het), and 
5 + n parameters from the planetary atmosphere: f, P_0, the haze 
parameters a and γ_haze, the atmospheric temperature T 
and n mixing ratios of the different elements considered. As with our 
spot modelling, we explore the parameter space using MultiNest <cit.> 
with the PyMultinest Python wrapper <cit.>, which allows 
us to properly compare the probability of different models given the data. For 
the parameters of stellar origin we consider uniform priors between 
2,500 and 6,500 K for the temperatures and a uniform prior for the fraction 
f_het between 0 and 1. For the planetary atmosphere parameters 
we set a wide uniform prior between 0.8 and 1.2 for f, a wide log-prior 
between 10^-6 and 10^6 bars for P_0, a wide log-prior between 10^-30 
and 10^30 for a, a wide uniform prior from -50 to 10 for 
γ_haze and a uniform prior between 1,000 and 3,000 K for 
the atmospheric temperature T. For the volume mixing ratios of the elements 
under consideration we set log-uniform priors between 10^-30 and 1.





