ACCESS, LRG-BEASTS, and MOPSS: Featureless Optical Transmission Spectra of WASP-25b and WASP-124b

We present new optical transmission spectra for two hot Jupiters: WASP-25b (M = 0.56 M J ; R = 1.23 R J ; P = 3.76 days) and WASP-124b (M = 0.58 M J ; R = 1.34 R J ; P = 3.37 days), with wavelength coverages of 4200–9100 Å and 4570–9940 Å, respectively. These spectra are from the ESO Faint Object Spectrograph and Camera (v.2) mounted on the New Technology Telescope and Inamori-Magellan Areal Camera & Spectrograph on Magellan Baade. No strong spectral features were found in either spectra, with the data probing 4 and 6 scale heights, respectively. Exoretrievals and PLATON retrievals favor stellar activity for WASP-25b, while the data for WASP-124b did not favor one model over another. For both planets the retrievals found a wide range in the depths where the atmosphere could be optically thick (∼0.4 μ–0.2 bars for WASP-25b and 1.6 μ–32 bars for WASP-124b) and recovered a temperature that is consistent with the planets’ equilibrium temperatures, but with wide uncertainties (up to ±430 K). For WASP-25b, the models also favor stellar spots that are ∼500–3000 K cooler than the surrounding photosphere. The fairly weak constraints on parameters are owing to the relatively low precision of the data, with an average precision of 840 and 1240 ppm per bin for WASP-25b and WASP-124b, respectively. However, some contribution might still be due to an inherent absence of absorption or scattering in the planets’ upper atmospheres, possibly because of aerosols. We attempt to fit the strength of the sodium signals to the aerosol–metallicity trend proposed by McGruder et al., and find WASP-25b and WASP-124b are consistent with the prediction, though their uncertainties are too large to confidently confirm the trend.


Introduction
Exoplanet science is at a new frontier, where the need to repurpose telescopes to observe planetary atmospheres is being replaced with telescopes that are designed with the goal of exoplanet atmosphere characterization in mind. Exoplanet scientists have made plenty of advancements with the current generation of telescopes. Specifically with low-resolution transmission spectroscopy, strives have been made with ground-based telescopes (e.g., Sing et al. 2012;Nikolov et al. 2016;Diamond-Lowe et al. 2018;Todorov et al. 2019;Spyratos et al. 2021), the Hubble Space Telescope (HST; e.g., Charbonneau et al. 2002;Kulow et al. 2014;Sing et al. 2016;Tsiaras et al. 2016;Wakeford et al. 2020;Rathcke et al. 2021), and Spitzer (e.g., Knutson et al. 2011;Pont et al. 2013;Alam et al. 2020;Alderson et al. 2022). However, the advancements with the next generation of telescopes will be revolutionary. This has already begun with the launch and utilization of the JWST, where novel science has been conducted with outstanding quality of data and newly discovered molecular features (e.g., Tsai et al. 2023;Ahrer et al. 2023b;Alderson et al. 2023;Feinstein et al. 2023;Rustamkulov et al. 2023). Furthermore, soon-to-belaunched telescopes like Pandora (Quintana et al. 2021;Hoffman Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. et al. 2022), the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (Tinetti et al. 2018), and the next generation of ground-based telescopes: the Extremely Large Telescope (ELT), 20 Thirty Meter Telescope (TMT), 21 and Giant Magellan Telescope (GMT) 22 will have designs and instruments specific for exoplanet atmospheric studies. These telescopes will undoubtedly be cornerstones in advancing our understanding of exoplanet atmospheres.
There is still much about exoplanet atmospheres that alludes us. For example, we have no direct link with physical conditions that cause high-altitude aerosols to form in the upper atmosphere of some observed planets (e.g., Alam et al. 2018;Chachan et al. 2019;Estrela et al. 2021) but not others (e.g., Sing et al. 2016;Kirk et al. 2019;Alam et al. 2021;Ahrer et al. 2022;McGruder et al. 2022). Here we refer to aerosols as clouds-condensation material due to specific atmospheric conditions-or hazes-material formed from photochemical reactions. The formation of aerosols likely occurs in most (if not all) atmospheres, just like in our solar system; however, high-altitude aerosols are normally the main concern when probing atmospheres. This is because the most used method to probe exoplanet atmospheres is transmission spectroscopy, which more easily probes the upper atmospheric limbs of planets, due to geometry and opacities (Lecavelier Des Etangs et al. 2008;Kreidberg 2018;Sing 2018). There are a number of studies aimed at understanding the composition and formation of high-altitude hazes (e.g., Moses et al. 2011Moses et al. , 2013Fleury et al. 2019) and clouds (e.g., Helling 2019; Gao et al. 2020;Estrela et al. 2022), and though there has been a lot of headway toward this, there is little observational support for leading theories.
Additionally, finding observational trends in aerosol formation has proven illusive, with many contradicting or inconclusive studies (Heng 2016;Stevenson 2016;Fu et al. 2017;Fisher & Heng 2018;Tsiaras et al. 2018;Dymont et al. 2022;Estrela et al. 2022). A possible reason why no correlation has been clearly identified is the limited number of observed planets relative to the parameter space, where tens of planetary atmospheres have been used for studies, but tens of parameters (host star, orbital, and planetary parameters) could be correlated to aerosol formation rates. Possible solutions are to either increase the number of planetary atmospheres observed or to reduce the parameter space by observing selected targets with many parameters that nearly match. The observations presented here aim to address both methods.
We observed the atmospheres of WASP-25b (M = 0.6 M J , R = 1.2 R J , P = 3.764 days, host star = G4,V mag = 11.9; Enoch et al. 2011;Brown et al. 2012;Southworth et al. 2014), and WASP-124b (M = 0.6 M J , R = 1.3 R J , P = 3.373 d, host star = F9,V mag = 12.7; Maxted et al. 2016). We obtained three spectroscopic transits of WASP-25b and five of WASP-124b as part of ACCESS. 23 We obtained one additional transit of WASP-25b with the ESO Faint Object Spectrograph and Camera (v.2; EFOSC2) instrument on the ESO New Technology Telescope (NTT) as part of LRG-BEASTS. 24 There was also a full and partial transit of WASP-124b obtained by the MOPSS team 25 that we added to our data set. Neither of these planets have atmospheric observations published. Furthermore, they have very similar parameters to one another and are part of a sample of seven planets proposed by McGruder et al. (2023) that could be key for identifying correlations with high-altitude aerosols and observed parameters.
The observation and reduction of all data are described in Section 2, followed by their light curve analysis in Section 3. In Section 4 we introduce the combined transmission spectra of both targets, discuss our retrieval analysis of the data (Section 4.1), and interpret the retrieval results (Section 4.2). We then compare our results with expectations from the tentative aerosol-metallicity trend proposed by McGruder et al. (2023) in Section 5. Finally, in Section 6 we recapitulate and provide conclusions.

Magellan/IMACS Transits
We observed three transits of WASP-25b (UTYYMMDD: UT180620, UT210306, UT220325) and five transits of WASP-124b (UT190826, UT210809, UT210905, UT211002, UT220605) with the Inamori-Magellan Areal Camera & Spectrograph (IMACS; Dressler et al. 2011) mounted on Baade, one of the twin 6.5 m Magellan telescopes. Those transits were observed as part of ACCESS and used a setup similar to previous ACCESS observations (i.e., Kirk et al. 2021;Weaver et al. 2021;Allen et al. 2022;McGruder et al. 2022). Our general setup uses the 8 × 8 K CCD mosaic camera at the f/2 focus, a 300 line mm −1 grating at a blaze angle of 17°.5, and a GG455 filter. This gave a wavelength coverage of (4550-9900 Å) and dispersed the spectra across two chips, but we managed to fit the spectrum of the target from 4550 to 9100 Å on one CCD, preventing gaps from occurring at wavelengths of particular interest (see Figure 1). We used 2 × 2 binning and the FAST readout mode to reduce the readout time and improve the observational duty cycle. We used 10″ by 90″ slits (0 5 by 90″ for HeNeAr wavelength calibration lamps), putting the observations in a seeing-limited regime, with an average resolving power of R = 1350. The number of exposures, range of airmasses, instrument setup, and resolution per night can be found in Table 1.
We combined our ACCESS observations of WASP-124b with two observations obtained by the MOPSS team on UT180915 and UT190615 (UT190615 was a partial transit). We used Magellan/ IMACS with a similar observational setup to that of ACCESS, but with the 300 line mm −1 grating at a blaze angle of 26°.7 and a slit size of 15″ by 20″ (1″ by 1″ for HeNeAr wavelength calibration lamps). We also had different comparison stars, which are highlighted in Table 1.
All IMACS observations utilize the multiobject spectrograph (MOS) mode to observe multiple comparison stars simultaneously. The best comparison stars were selected based on the procedure outlined in Rackham et al. (2017), where we consider a nearby star suitable if it has a color difference of D < 1 with the target. D is defined as where the uppercase letters correspond to the Johnson-Cousin apparent magnitudes of the stars, and the subscripts t and c  (Jordán et al. 2013;Rackham et al. 2017;Bixel et al. 2019;Espinoza et al. 2019;McGruder et al. 2020;Weaver et al. 2020;Kirk et al. 2021;Weaver et al. 2021;Allen et al. 2022;McGruder et al. 2022). 24 The Low Resolution Ground-Based Exoplanet Atmosphere Survey using Transmission Spectroscopy Louden et al. 2017;Kirk et al. 2018 indicate the target and potential comparison, respectively. The sky coordinates of each comparison and their D relative to the target can be found in Table 1.

IMACS Reduction
The reduction process for all IMACS data was the same and uses a custom ACCESS pipeline, which has been described in detail by Jordán et al. (2013), , Rackham et al. (2017), andBixel et al. (2019). In general, this includes wavelength calibration using the HeNeAr arc lamp measurements, bias subtraction with the overscan region, pixel tracing, and sky background subtraction utilizing the median counts outside the science aperture. The radius of the science aperture was determined by taking the average FWHM over time and wavelength. This value was then added to 3 times the standard deviation (STD) of all FWHM values (all calculated wavelengthand time-dependent FWHM values), i.e., aperture = 〈FWHM〉 + 3 × STD FWHM . Allen et al. (2022) found that optimal extraction (Marsh 1989) has the potential to be a more effective way of identifying and correcting for bad pixels and cosmic rays in ACCESS data. When testing the effectiveness of this reduction step versus the traditional pipeline steps, we see slightly less scatter in the resulting white light curve with optimal extraction, so we adopt the results of this method. The final reduced spectra for the targets from each night are shown in Figure 1.

NTT/EFOSC2 Transit
WASP-25b had an additional transit observation with the European Faint Object Spectrograph and Camera 2 (EFOSC2; Buzzoni et al. 1984) mounted on the 3.6 m New Technology Telescope (NTT) as part of LRG-BEASTS and ESO program 0100.C-0822(A) (PI: Kirk). The observation was taken on the night of UT180329 with the same instrument and setup used to detect Na in the atmosphere of WASP-94Ab (Ahrer et al. 2022) and clouds in the atmosphere of HATS-46b (Ahrer et al. 2023a). This was a 27″× 3 53 slit and Grism #11, which provided spectra from 3960 to 7150 Å and an average seeinglimited resolution of R = 150, which was dispersed on a 2048 × 2048 pixel Loral/Lesser CCD. Our long slit allowed us to simultaneously observe the comparison star 2MASS J13011275-2730485 with a V = 12.5 and B − V = 1.5, whereas WASP-25 has a V = 11.9 and B − V = 0.7. We used 2 × 2 binning and the fast readout mode.
We also obtained 54 biases, 7 lamp flats, 17 morning twilight sky flats, 26 and 3 arc lamps at the beginning and end of the . Each spectra is plotted in the same order in which it is printed in the legend (top to bottom). The shaded regions of the same color extending past the median lines are the 1σ range of counts extracted for that night. Each spectroscopic bin used is demarcated by dotted vertical lines, where the only gaps in the binning scheme are the strong telluric region from 7594 to 7638 Å (lightly red shaded region), the CCD gap for the first two WASP-124 observations (bottom left) at 6317-6424 Å, and the CCD gap for the other five WASP-124 observations (bottom right) at 9100-9225 Å (gray shaded region). The specific instrument and setup used for each spectrum are printed, where EFOSC2, IMACS-17°. 5, and IMACS-26°. 7 refer to the LRG-BEASTS (see Section 2.3), ACCESS, and MOPSS (see Section 2.1) setups, respectively. The differing throughputs between the grism used for the ACCESS and MOPSS data explain the different spectral shapes between the two WASP-124 spectra on the bottom left and the five spectra on the bottom right. Note the plotted spectrum of the NTT/EFOSC2 data (top right) was created without the sky flat, in order to visually compare with the other spectra. However, the final spectrum used for data analysis did use the sky flat as discussed in Section 2.4.
observations. The flats were taken with the same 27″-wide slit as the science observations but the arc lamps were taken with a 1″ slit to avoid saturation and ensure narrow lines for calibration. This meant that the arc lamps were only used for an initial wavelength calibration with the final wavelength calibration performed using absorption lines in the stellar spectra.

EFOSC2 Reduction
We reduced and processed the data using the LRG-BEASTS pipeline, which is described in more detail in Kirk et al. (2017Kirk et al. ( , 2018Kirk et al. ( , 2021. For the flats, we created two sets of reductions, one without a flat-field correction, as is standard for LRG-BEASTS observations (e.g., Alderson et al. 2020;Kirk et al. 2021;Ahrer et al. 2022), and one with a flat-field correction. This flat-field correction was performed in a novel way whereby we used the master sky flat without removing the sky spectrum from the sky flat. The motivation behind this approach was to avoid uncertainties associated with fitting out the sky spectrum (with a running median for example) while also capitalizing on the higher number of blue photons from the sky flat compared to a lamp flat. As we did not remove the sky spectrum from the sky flat, this meant that the stellar spectra (F 1 , F 2 ) we extracted were contaminated by the sky background imprinted into the sky flat (F sky ). Therefore, the stellar spectra we extracted were actually F 1 /F sky,1 and F 2 /F sky,2 . The target and reference stars drifted across the course of the observations, leading to changes in F sky,1 and F sky,2 . This meant that the F sky terms did not fully cancel out after dividing the target's light curve (Σ(F 1 /F sky,1 )) by the comparison's light curve (Σ(F 2 /F sky,2 )). However, we found that using the sky flat led to light curves and transmission spectra that deviated by < < 1σ from the same reduction without the flat field, while also decreasing the white noise in the light curves by 6% and leading to a small improvement the precision in the transmission spectrum (∼3%). This demonstrates that impacts on the spectrophotometry due to Δ(F sky,1 ) and Δ(F sky,2 ) are insignificant, likely due to the fact that the stellar traces drift by <3 pixels, which is corrected for by cross correlation, and this constitutes only 8% of our average bin width used to make the transmission spectrum. Due to this test, we adopted the reduction using the sky flat for the rest of our analysis.
To extract the stellar spectra we experimented with different aperture widths and background widths and compared the noise of the resulting white light curve in each case. We found that the lowest white light noise resulted from a combination of a 22-pixel-wide aperture, with two 15-pixel-wide background regions on either side of the aperture, each separated from the aperture by 15 pixels. We fit a linear polynomial across these two background regions to remove the sky background from our spectra. Following the extraction of the stellar spectra, we Note. The resolutions were calculated from the FWHM near the peak of the spectrum. The magnitudes used to calculate D were obtained from the UCAC4 Catalog (Zacharias et al. 2013). The numbers in parentheses represent the comparison labels that refer to a specific star.
clipped cosmic-ray hits by identifying 5σ outliers in the spectral time series and replaced these with a linear interpolation between the nearest two neighboring pixels. We then corrected for shifts in the stellar spectra of ±2 pixels across the night. This was done by cross correlating each spectrum with a spectrum taken in the middle of the observations and then performing a flux-conserving resampling of each spectrum onto the reference wavelength grid. Figure 1 (top right) shows the final spectrum extracted for WASP-25 with this pipeline.

Light Curve Analysis
The general steps in our light curve analysis are the same as what has been implemented in previous ACCESS papers (e.g., Allen et al. 2022;McGruder et al. 2022). This includes first creating a photometric (white) light curve by combining all counts of the entire spectrum for each exposure. The systematics are removed from the white light curve, and the transit is fit to obtain the transit parameters. These parameters are then used to constrain the priors for fitting the spectrophotometric (binned) light curves. The binned light curves are produced by summing the light within a band of wavelengths. The appropriate widths and centering of bins were determined by considering spectrophotometric precision, the overlap of spectral bands from different observations, high telluric absorption regions, and the desire to properly probe for atmospheric features. The average bin widths were ∼150 Å for WASP-25b and 160 Å for WASP-124b, with 90 Å bins centered on the Na and K doublets and wider bins on lowthroughput edges of the spectra. The final binning scheme used is demarcated with dotted vertical lines in Figure 1 and written for each bin in the first column of the Figures in Appendix A.

White Light Curve Fitting
We detrended the white light curve (WLC) using a combination of principal component analysis (PCA) and Gaussian processes (GPs), which we refer to as PCA+GP. 27 This routine is identical to what was used by McGruder et al. light curves in magnitude space, which yield eigenvectors and values that allow us to identify principal components that capture features common to all comparison light curves. However, we still needed to reduce systematics unique to the target star, for which we use GPs. We used george (Ambikasaran et al. 2015) to construct and evaluate the likelihoods of a multidimensional squared-exponential kernel dependent on the auxiliary observables of airmass, FWHM, sky flux, trace, and wavelength solution drift. Details of the GP hyperparameter priors are shown in Table 2 Table 1). The priors for the radius of the planet relative to the star (R p /R s ) and the parameterized quadratic limb darkening (LD) parameters q 1 and q 2 (Kipping 2013) were set wider. Table 2 also provides information on those and the other transit priors.
After the transit parameters of the WLC were acquired for each night, we weight-averaged the P, a/R s , and b values from each night to obtain more constrained terms for each. These means were held fixed for an additional pass of the PCA+GP run, to obtain final values for t 0 , R p /R s , q 1 , and q 2 . 28 Figure 2 shows the final detrended light curves of all transits, and Table 3 shows the values obtained when all parameters were fitted for (first four columns) and when the common parameters were fixed (last four columns). From the table, one can see that only the LRG-BEASTS transit depth and the MOPSS partial transit depth differ by more than 2σ between each other transit of a specific target. However, those two outliers are likely due to the LRG-BEASTS observations being bluer than those from Note. The priors for the GP hyperparameters are amplitude (α), jitter (ξ), and inverse squared length scale (1/λ). For 1/λ, when the a parameter of a gamma function is set to 1, it becomes an exponential function (i.e., e − x ). The mean and standard deviation (σ n ) values of the transit parameters (variables defined in Section 3.1) were obtained directly from McGruder et al. (2023). With the exception of t 0 , which was deduced for a given night from the period and t 0 obtained from McGruder et al. (2023). The σ n for this parameter was set to 30 minutes for each night. 27 Prior to the detrending process, we removed outliers via visual inspection after dividing the target light curve by the sum of the comparisons' light curves. 28 We did not include a transit fitting pass with the common parameters free for the partial transit of UT190615 and just used the weighted means of the other six WASP-124b transits for fitting t 0 , R p /R s , q 1 , and q 2 .  Note. P, b, a/R s , and i are the values obtained when leaving the transit parameters free, and R p /R s , t 0 , q 1 , and q 2 are values obtained from the second pass where only those transit parameters were allowed to be free and all others were fixed on the weighted mean values. Each printed value of t 0 is subtracted by 2,450,000 days. Transit UT190615 does not have the first four parameters because we did not allow those parameters to be free for the partial transit and just used the weighted means of the other six WASP-124b transits for fitting t 0 , R p /R s , q 1 , and q 2 .
ACCESS and the lack of full transit coverage from MOPSS hindering the detrending process.

Spectrophotometric Light Curves
We used two separate detrending routines to reduce the binned light curves (BLC). Both routines were discussed and tested by McGruder et al. (2022). For all WASP-25b transits and WASP-124b transits on UT190826, UT210809, and UT211002 we used common-mode correction (CMC) followed by polynomial correction (Poly), CMC+Poly. This method assumes that most of the systematics are captured when detrending the WLC (following the procedures of Section 3.1). As such it uses the quotient of the final detrended WLC and the "raw" light curve, produced from the normalized target divided by the sum of comparison starlight curves, as a common-mode term. The raw light curve also had a moving average 4σ clipping, similar to what was done in McGruder et al. (2022). Each BLC is then divided by this common-mode term. We then apply polynomial regression models dependent on the auxiliary observables (i.e., airmass, FWHM; see Section 3.1), to remove any additional chromatic systematics unique to each bin. For each auxiliary observable we allowed the polynomial to go up to the fourth order, aside from the airmass, which only went up to the second order because of the smooth correlation with it and the light curves. We tested all combinations of each auxiliary parameter and order polynomials 29 using scipy. optimize.minimize (Virtanen et al. 2020). We then ranked each combination by polynomial corrections based on a standardized sum of χ 2 , the Bayesian information criterion (BIC), the Akaike information criterion (AIC), and the rms of the residuals (rRMS), and took the highest-ranked 100 models (lower sums) to do a final pass of fitting with PyMultiNest. The priors on polynomial coefficients were normal with a mean set by the scipy.optimize.minimize fits and standard deviation of one. We used BMA to combine posteriors and arrive at our final binned transit depths.
The detrending method used for transits UT180915, UT190615, and UT210905 (three WASP-124b transits) was PCA+GP, which is the same algorithm used for our WLC analysis (see Section 3.1). The reasoning for using this routine instead of the CMC+Poly is that the systematics found in each bin are significantly different from one another (see the corresponding first columns of figures in Appendix A), so the CMC term poorly corrects the bins and likely introduces more systematics, which simple polynomials have trouble modeling. In fact, the transits that required PCA+GP had obvious issues with their observations such as missing ingress, scattered cloud coverage, and poor seeing. However, we also tested the performances of both methods by examining the binned spectra each produced; in particular we compared the variance of their transmission spectra and the average residual red noise of each bin β (McGruder et al. 2022; see appendix D). We found that for the three aforementioned WASP-124b transits, β and the variance were substantially worse. This supports a finding of McGruder et al. (2022), in which the appropriate detrending method should be considered by testing multiple approaches for each data set.
For both fitting methods all transit parameters were fixed to the WLC best fit (columns 2-5 of Table 3), but the LD parameters were uniform from 0 to 1 and R p /R s had a normal prior with a mean set by the WLC best fits and a standard deviation of 0.02.

Transmission Spectra
We produced our transmission spectra by comparing the bestfit R p /R s found for each detrended wavelength bin. We then combined the transmission spectra from each night to form a global transmission spectrum for each target. The spectra from each night had an offset applied so the mean depths across overlapping spectral bins were the same, as was done by McGruder et al. (2022). The spectra were then combined by weight averaging each overlapping bin. The new points were obtained using the python numpy.average (Harris et al. 2020) function where the weight of each R p /R s was given by the inverse squared R p /R s errors ("inverse variance weighting"). The errors of the resulting weighted R p /R s were calculated as the square root of the weighted variance (again using inverse squared R p /R s error as weights). For WASP-25b the overlapping bins were from 4570-7113 Å and had a mean depth of 0.14040 R p /R s , which corresponded to weighted offsets of 0.00584, 0.00016, 0.00265, and −0.00552 R p /R s for each night in chronological order. For WASP-124b the mean depth was 0.12811 R p /R s where all wavelengths overlapped aside for 6317-6424 and 9100-9225 Å. This yielded offsets of −0.00047, 0.00122, −0.00494, 0.00049, 0.00093, and 0.00295, respectively.
For WASP-124b, the partial transit on UT190615 was not included because the scatter of that transmission spectrum was too high to positively contribute to the combined spectrum. This is likely because the lack of ingress prevented both detrending methods (PCA+GP and CMC+Poly) from accurately constraining the systematics. Figure 4 shows the transmission spectra of each night (except transit UT190615) and the combined transmission spectra.
The average precision (68% confidence interval) of the combined spectra per bin 30 in R p /R s for the WASP-25b data was 0.00301 (841 ppm in depth) and 0.00485 (1238 ppm in depth) for WASP-124b. With these precisions, we can only probe as low as 4.03 and 5.71 atmospheric pressure scale heights for WASP-25b and WASP-124b, respectively, for which the scale heights are 453.2 and 634.9 km. Thus, it is difficult to determine if the relatively flat spectra seen in both targets (see Figure 4) are due to a true dearth of planetary features or a precision limitation.
It should also be noted that even though there are more transits for WASP-124b, its transmission spectrum is less precise than the four transits of WASP-25b because WASP-124 (V mag = 12.7) is dimmer than WASP-25 (V mag = 11.9) and the overall quality of the WASP-124b observations were not as high, as implied in Section 3.2 where we apply the PCA+GP detrending routine for two out of the six used transits because of aggressive systematics.

Retrieval Analysis
We ran a series of retrieval models with PLATON (Zhang et al. 2019) and Exoretrievals (Espinoza et al. 2019) on our final combined transmission spectra of WASP-25b and WASP-124b. Our analysis process was to run models including stellar activity, scattering features, or common atomic/ molecular species observed in planets of this type (i.e., H 2 O, Na, K) and the different combinations of each with 29 1875 combinations for five observables and up to fourth order for everything but the airmass, which is up to second order. 30 The average bin size was ∼150 and 160 Å for WASP-25b and WASP-124b, respectively.
Exoretrievals. Concurrently, we run models including scatters, stellar activity, neither, and both with PLATON, which assumes equilibrium chemistry and fits for the C/O ratio and planetary metallicity to extrapolate the abundances. This retrieval analysis workflow is the same as was done by McGruder et al. (2020) for WASP-31b, McGruder et al. (2022) for WASP-96b, and McGruder et al. (2023) for WASP-6b and WASP-110b. To determine which models are preferred over the others we used the Bayesian evidences (Z), given for each model because both retrieval routines use nested sampling to explore the posterior space (PyMultiNest with Exoretrievals and dynesty (Speagle 2020) with PLATON). Following the reasoning of Trotta (2008) and Benneke & Seager (2013), we considered a ΔlnZ less than 2.5, not significantly favoring one model over another; ΔlnZ between 2.5 and 5, moderately favoring the higher evidence model; and a ΔlnZ greater than 5 strongly supports the higher-lnZ model. 31 . A table of the difference in natural-log evidences (ΔlnZ) of a given model relative to the least complex model for the given retrieval is shown in Table 4. The term we use for when there are no scatters in the planetary atmosphere or activity in the stellar photosphere added to the model is plain. Given that PLATON assumes equilibrium chemistry and atomic/molecular abundances are inherently determined through this, a plain model is the least complex model used for PLATON. In the case of Exoretrievals, its least complex model is when the presence of species is turned off, i.e. a plain model (no scatters or activity), but also without atomic/molecular species included. In that case the spectrum is completely. In that case the spectrum is completely flat.
When determining the priors for the stellar activity parameter in our retrieval analysis, we used log 10 ( ¢ R HK ) and stellar rotational periods obtained from high-resolution spectra and photometric observations in McGruder et al. (2023 ; Table 1). WASP-124 has a log 10 ( ¢ R HK ) of −4.765 ± 0.056, consistent with WASP-96b's (log 10 ( ¢ R HK = −4.781 ± 0.028), which has been established to be quiet (McGruder et al. 2022;Nikolov et al. 2022). However, WASP-124 has a faster rotational period (10.65 -+ 3.01 3.27 days), which has been found to be correlated to  Note. For WASP-124b the heterogeneity parameters with Exoretrievals were obtained using the model that only included activity and the other parameters were obtained using the plain model that included K and Na. According to the evidence, both of those models were indistinguishable from each other (see Table 4). We used the model with activity, sodium, and potassium to obtain the best-fit parameters for WASP-25b. This model was used because all models including activity were indistinguishable from each other, and this model obtained elemental mixing ratios. The difference in obtained overlapping parameters (i.e., T p , ( ) P log 10 0 , ΔT het , f het ) with that model and the one with just activity and water, were well within their uncertainties. The obtained water abundance was not constrained given there are no water features in the data, so we do not show its best fit here. Given that there are no carbon or oxygen-bearing species in the wavelength coverage of our data, we do not report the C/O ratios retrieved by PLATON. The pressure in log 10 (P 0 ) is given in bars for both retrievals. activity levels (e.g., Pizzolato et al. 2003;Wright et al. 2011Wright et al. , 2013. With this in consideration we allow the covering fraction of unocculted inhomogeneities to vary uniformly from 0 to 6.8%, which is consistent with the 2σ upper level of activity for stars of this type found by Rackham et al. (2019; see their Tables 2 and 3 and Equation 2). For WASP-25, the log 10 ( ¢ R HK ) of −4.507 ± 0.119 and faster stellar rotational period of -+ 16.93 1.55 2.02 days suggests it is a somewhat active star; as such, we do not limit its stellar inhomogeneity coverage and set the covering fraction priors to be uniform from 0% to 50%. The priors used for each retrieval run are given in Appendix B.

Retrieval Interpretation
The best-fit parameters from the retrieval analysis can be seen in Table 5. Overall, the PLATON and Exoretrievals results agree with each other well for each target. The lack of prominent features in either spectrum makes it difficult for planetary atmospheric properties to be constrained, which is outlined by the wide 1σ range given for every parameter in Table 5. For WASP-25b, the best-fit models were plain atmospheres (i.e., no scatterers or atomic/molecular features) for the exoplanet and an inhomogeneous photosphere for the stellar host, with ∼20% coverage of cold spots at a temperature contrast of ΔT ∼ −2000 with respect to the quiescent photosphere. However, the uncertainty of these inhomogeneity parameters is large (see Table 5), with the most extreme case being the retrieved PLATON inhomogeneity temperature of −2001 -+ 531 1473 K. For WASP-124b, using Exoretrievals, the highest-evidence models were those with low levels of stellar activity or a plain planetary atmosphere. Even still, those evidences were indistinguishable from a flat-line model. The PLATON retrievals had the same issue where no model was significantly favored over another. This emphasizes the difficulty of constraining the atmosphere of WASP-124b with the data at hand.
The limb temperatures obtained for WASP-25b and WASP-124b using both retrieval methods (see Table 5) are in agreement with their corresponding effective temperatures of 1217 ± 101K and 1481 ± 123K, respectively (McGruder et al. 2023; Table 1). However, again the uncertainties in the retrieved values are quite large. The pressures where the atmospheres are optically thick were also poorly constrained, with values ranging from log 10 [bars] of −3.0 -+ 3.4 3.8 and −2.8 - for WASP-25b and from −3.1 - 2.2 for WASP-124b, with Exoretrievals and PLATON, respectively. Thus pressures from ∼0.4 μ to 0.2 bars are all within a 1σ interval for WASP-25b and from 1.6 μ to 32 bars for WASP-124b. Figure 3 shows the corner plot obtained for the PLATON best fit of WASP-25. It also highlights the difficulty in retrieving precise atmospheric parameters.
We compared these results to the analysis of WASP-31b . These planets were chosen because they underwent the same retrieval analysis, minimizing differences that may arise from varying model assumptions and priors (e.g., Kirk et al. 2019;Barstow et al. 2020). 32 The upper bounds of the pressure ranges for both targets are consistent with the 68% interval of the WASP-96b fit, which strongly indicates the absence of aerosols in the observed wavelength range of 0.4-1.24 μm. However, their lower bounds are also consistent with the cloud top pressure found for WASP-110b, which has the highest retrieved cloud top altitude (i.e., the largest amount of high-altitude aerosols) of the four planets. Thus, we reaffirm that further observations are needed to constrain the atmospheres of WASP-25b and WASP-124b.
When attempting to interpret the relatively featureless spectra, we can deduce that it is unlikely that hazes are prominently present in the upper atmospheres of the planets, because there is no scattering slope observed. Strong scattering slopes in the optical have been suggested to signify hazes and could have signals as high as 15 scale heights (Ohno & Kawashima 2020), well within the precision of the data. Therefore, the more likely cause for the observed features is either high-altitude clouds or observational limitations due to the lower precision, with one scenario not necessarily explaining both atmospheres. To obtain a better understanding of these atmospheres, higher precision optical observations with the HST and longer wavelength observations, ideally with the JWST, are required. McGruder et al. (2023) proposed that there is a tentative trend with observed high-altitude aerosols and the host-star metallicity for a group of seven planets with very similar system properties, which includes WASP-25b and WASP-124b. They claim that, if this trend is true, then WASP-25b and WASP-124b would sit on opposite ends, where WASP-25b would be obscured by aerosols (like WASP-6b and WASP-110b), and WASP-124b would be relatively clear (like WASP-96b).

Similar Seven
We explore how well this trend holds here, using the sodium signal as a proxy for aerosol levels, as was done in McGruder et al. (2023). We find a small hint of Na in the spectrum of WASP-124b, which is stronger than for WASP-25b, even though the data precision of WASP-25b probes deeper. However, when looking at the retrieval analysis results of the WASP-124b data (see Table 4), we see no strong favoring of the Exoretrievals model, which includes Na. Furthermore, the log mixing ratio of Na found with the Exoretrievals fit that included it is not well constrained and suggests marginal amounts of Na (−9.7 -+ 14.1 6.8 ). As such, we have no detection of Na in WASP-124b.
This seems inconsistent with the proposed trend, but the scale height probed with the WASP-124b data is 5.71. This is 3 times higher than what was able to be probed with WASP-6b, 1.998,2.212,respectively), which were used to identify the tentative trend. The precision of the other targets is likely higher because they include observations from larger (VLT) or space-based (HST) telescopes. In Figure 5 we plot a linear aerosol-metallicity trend, where the Na feature was used as a proxy for aerosols, similar to what was done by McGruder et al. (2023). However, here we divide the Na amplitude values by their theoretical Na signal when no aerosols are present in the atmosphere, ΔR p /R s (see Equation (10) of Heng 2016). This was done, because even though the planets are twin-like, they are not exactly the same and would have slightly different maximum possible signals. Yet, because of the planets' similarity, this modification had little effect on the previous trend found by McGruder et al. (2023), as shown in Figure 5.
In Figure 5, we plot a linear fit with and without the WASP-25b and WASP-124b Na signals and find that given the errors both WASP-25b and WASP-124b are consistent with the trend found using just WASP-6b, WASP-96b, and WASP-110b. In Figure 5 our linear trend is fit with scipy.odr (Virtanen et al. 2020) where we used the inverse of each Na signals' errors as weights. The regression score, R 2 , of the weighted linear fit with the Na signals from WASP-6b, WASP-96b, and WASP-110b is 0.755. When including the Na signals from WASP-25b and WASP-124b, the score was 0.754, showing that the trend continues to hold. We also compare the data to a flat-line fit, i.e., no trend, and find a mean Na amplitude of 0.083 and R 2 of −0.051, emphasizing that a linear trend is much more appropriate given the data.
Though there is no strong support for WASP-124b having high-altitude aerosols, if it does and is inconsistent with the tentative trend found by McGruder et al. (2023), a possible explanation might be due to its difference in equilibrium temperature. All the planets' equilibrium temperatures lay around ∼1250 K, and though the equilibrium temperature of WASP-124b (1481 ± 123 K) is less than 2σ from the coolest planet in our sample (WASP-6b, T = 1167 ± 96 K), it is possible that this temperature difference is important. Many studies find that equilibrium temperature is important in aerosol formation (e.g., Stevenson 2016;Fu et al. 2017;Gao et al. 2020;Estrela et al. 2022), but this literature is not consistent on whether an increase in temperature for this class of planets would produce more or less aerosols. Thus, there is no strong support in the literature that the ∼300 K difference in equilibrium temperature of WASP-124b is an important factor in the tentative trend found. Still, if metallicity does not strongly correlate with aerosol formation rates, then perhaps an unobserved parameter, such as high-energy emissions, is correlated to the higher aerosol rates observed in some of the Similar Seven planets. Alternatively, the complex nature of aerosol formations in their extreme environments might make it difficult to correlate one or two parameters to the observed aerosol rates.
To have a better grasp of whether the aerosol-metallicity trend truly holds, at minimum, more optical observations of The final weighted-average spectra are shown in violet for WASP-25 and blue for WASP-124, with their individual transmission spectra used for the combined spectra plotted in transparent colors. The best-fit PLATON retrieval models are plotted as a black line with the 1σ confidence interval highlighted in light blue. For both targets the best-fit models are the ones that just include activity, but in the case of WASP-124b, this model is not significantly preferred over any other model.  2023), "Na I Amplitude" was the depth of the bin centered around the Na feature minus the average depth of all bins in the surrounding continuum. Here we take that value and divide it by the theoretical difference of the peak depth from Na and the continuum, in order to incorporate slight differences in the planets' scale heights and temperature. scipy. odr was used to obtain a weighted linear fit with these three planetary signals and is shown as a dashed blue line. Its regression score, R 2 , is 0.75. The 1σ interval of the fit is shaded around the dashed line and was calculated including the uncertainties in both metallicity and Na signal. Using the same method for WASP-6b and WASP-110b (i.e., average R p /R s centered at 5892.9 Å minus average R p /R s within 5340-5820 Å and 5960-6440 Å), we calculated the Na signals of WASP-25b (purple squares) and WASP-124b (purple triangles). A weighted linear fit with all five planet signals is shown as a purple dashed-dotted line, with an R 2 of 0.71. Although our new measurements have larger uncertainties they are consistent with the trend identified in McGruder et al. (2023). The metallicity range of WASP-55b and HATS-29b are plotted as yellow shaded regions.
WASP-124b are needed to precisely probe its atmosphere. This is achievable with the HST and/or larger ground-based telescope observations. Further Magellan/IMACS observations would also improve precision, especially if such observations were of good quality (i.e., the full transit, sufficient baselines, and good night conditions).

Summary and Conclusions
We observed four transits of WASP-25b, one with NTT/ EFOSC2 and three with Magellan/IMACS, and seven transits of WASP-124b with Magellan/IMACS. We combined the transmission spectra from each night for each target (excluding the partial transit of WASP-124b on 2019 June 15) to produce near-continuous final transmission spectra from 4200 to 9100 Å for WASP-25b and from 4570 to 9940 Å for WASP-124b. Our transmission spectra have an average precision in a depth of 841 ppm for WASP-25b and 1238 ppm for WASP-124b, corresponding to 4.03 and 5.71 scale heights, respectively. The spectrum of both targets lacked significant features.
Nevertheless, we ran a set of retrieval models utilizing Exoretrievals and PLATON on each final transmission spectrum. In doing so, we found that the most favored model for the retrievals (with D < Z ln 5) for WASP-25b is one that included ∼20% covering a fraction of unocculted cold spots at ∼2000 K cooler than the surrounding photosphere, but no molecular/atomic features. For WASP-124b there is no model strongly favored over another, but there are marginal hints of Na and K features and low levels of activity in the host star. The retrieved atmospheric parameters from the best-fit models have wide uncertainties for both planets, but the retrieved limb temperatures are consistent with the calculated equilibrium temperatures. Given that there are no strong atomic or molecular features in either spectrum, the pressure levels where the atmosphere is optically thick and the atmospheric metallicities are poorly constrained. The lack of features is possibly due to low precision being unable to probe the required depths for feature detection and/or high-altitude clouds obscuring the spectra.
We then put these planets' atmospheres in context with the aerosol-metallicity trend proposed by McGruder et al. (2023), and plot the sodium signal of these targets relative to their hoststar metallicities. We find that the uncertainties of the Na signal, caused by lower data precision, make it difficult to provide clear insight into the existence of such a trend. We believe that further observations with higher-quality data in the optical are necessary to confirm a trend. This could be done with the HST and/or further ground-based telescopes. JWST observations would provide a broader context of the nature of these planets' atmospheres. Proving this trend has the potential to drastically direct theoretical understandings of aerosol formation and could yield more efficient target selection criteria. Figure 6. Light curves (LCs) in each step of the detrending process when utilizing the CMC+Poly method (see Section 3.2), for each bin of transit UT180329. Left: LC produced from the quotient of the WASP-25 and comparison LCs. For all other transits there are multiple comparisons; in that situation the sum of comparison light curves is used instead. The wavelength range [nm] used is printed in the bottom left regions. Middle: quotient of the LCs on the left and common-mode correction term. Overplotted in violet is the best-fit polynomial correction term. If there were any frames clipped for a given bin, it is marked with a black x and the total number of clipped values is in the lower left regions. The precision, assuming only photon noise, is printed as σ w . Right: final detrended LC with the best-fit transit model overplotted in violet. The standard deviations of the residuals (σ r ) and σ r /σ w (β) are printed. The error bars shown are determined by multiplying the photon noise precision of each frame by the β for that bin.  Figure 1, in the middle the binned transit LCs of WASP-124b corrected with PCA (the best GP systematic fits are overplotted in violet), and on the right the final detrended LC with the best-fit transit model overplotted in violet. The best-fit GP jitter term was used to approximate the error bars shown in the middle and right columns. The wavelength ranges used, clipped frames, photon noise precisions, standard deviations of residuals, and β are printed on the figures just as in Figure 6. Right: same as Figure 6 but for transit UT190826.  Note. Each spectrum was produced by weight averaging (with an offset) each of the transits for the given target. The transmission spectra for individual nights are uploaded on zenodo.org/record/8047731.