Time-resolved Hubble Space Telescope Wide Field Camera 3 Spectrophotometry Reveals Inefficient Day-to-night Heat Redistribution in the Highly Irradiated Brown Dwarf SDSS 1557B

Brown dwarfs (BDs) in ultra-short-period orbits around white dwarfs (WDs) offer a unique opportunity to study the properties of tidally locked, fast-rotating (1–3 hr), and highly irradiated atmospheres. Here we present phase-resolved spectrophotometry of the WD–BD binary SDSS 1557, which is the fifth WD–BD binary in our six-object sample. Using the Hubble Space Telescope Wide Field Camera 3 Near-infrared G141 instrument, the 1.1–1.7 μm phase curves show rotational modulations with semiamplitudes of 10.5% ± 0.1%. We observe a wavelength-dependent amplitude, with longer wavelengths producing larger amplitudes, while no wavelength-dependent phase shifts were identified. The phase-resolved extracted BD spectra exhibit steep slopes and are nearly featureless. A simple radiative energy redistribution atmospheric model re-creates the hemisphere-integrated brightness temperatures at three distinct phases and finds evidence for weak redistribution efficiency. Our model also predicts a higher inclination than previously published. We find that SDSS 1557B, the second most irradiated BD in our sample, is likely dominated by clouds on the nightside, whereas the featureless dayside spectrum is likely dominated by H− opacity and a temperature inversion, much like the other highly irradiated BD EPIC 2122B.


INTRODUCTION
Irradiated planetary atmospheres -most prominently hot Jupiters -are generally shaped by four major processes: (1) Energy transport via radiation, convection, and large-scale circulation; (2) Atmospheric equilibrium and non-equilibrium chemistry; (3) Photochemistry, often including UV-driven reactions and dissociation; (4) Micro-and macro-physics of clouds and haze layers.With the steadily growing and increasingly high-quality datasets on hot Jupiters, significant progress has been made in understanding the combined impact of these processes.In particular, forced atmospheric circulation in the tidally locked atmospheres of hot Jupiters has been successfully explored (Showman & Guillot 2002;Heng & Showman 2015;Showman et al. 2020;Lian et al. 2022), as well as the formation of some cloud species and their spatial distribution (Wakeford & Sing 2015;Parmentier et al. 2016;Zhang & Showman 2018;Lee et al. 2018;Powell et al. 2019).
However, with any two hot Jupiters differing in multiple key parameters, it remains challenging to unequivocally separate the importance of the above processes and their impact on individual atmospheres.For example, changes in the 3D pressure-temperature distribution Figure 1.Graphic illustration of the SDSS 1557 system.Radii of the primary white dwarf and secondary brown dwarf are to scale, while orbital distances are not to scale.Colors correspond to calculated temperatures from spectral fitting.The circumbinary dust disk is believed to be accreting onto the WD (Farihi et al. 2017), although the accretion is not shown here.
and atmospheric chemistry may lead to observational changes that could also be explained by changes in the cloud properties and particle size distributions (Venot et al. 2020).
Studies of irradiated atmospheres that significantly differ in rotation and internal heat from hot Jupiters offer important opportunities to better separate the impact of these processes.As pointed out by Showman et al. (2020), highly irradiated brown dwarfs on very short-period (1-3 hr) orbits around white dwarfs (WD-BD systems) offer such targets.
While only relatively few, close-in WD-BD systems are known today (Casewell et al. 2012(Casewell et al. , 2020;;Beuermann et al. 2013;Parsons et al. 2017;Farihi et al. 2017;Owens et al. 2023), due to the range of orbital radii and WD temperatures present in the systems, they offer a broad range of irradiation levels.The irradiance in known systems range from those that correspond to cool hot Jupiters (800 K) (Samland et al. 2017) to those that correspond to ultra hot Jupiters (>2500 K) (Smith et al. 2011;Sing et al. 2013).At the same time, due to the high to very high (7700-25000 K) temperatures of the WDs, a greater fraction of the irradiation is emitted in the ultra-violet than it is the case for hot Jupiters that orbit main sequence stars.Furthermore, the mass of the BDs in WD-BD systems is typically 30-50× higher than is typical for hot Jupiters, translating to higher gravity atmospheres.Finally, while the orbital and rotational periods are synchronized (just as for hot Jupiters), the orbital periods are about an order of magnitude shorter (1-3 hr).In short, BDs in rare WD-BD systems provide an exciting opportunity to study the same processes that shape hot Jupiters, but to do so in part of the parameter space that overlaps with but also substantially different from those occupied by hot Jupiters.
As part of a larger Hubble Space Telescope (HST) program, we observed phase curves of six WD-BD systems using time-resolved G141 infrared spectrophotometry.The high photometric precision (∼0.3%) achievable with these observations allows us to probe the emission spectra of the brown dwarfs as a function of the longitude (or, equivalently, the rotational phase).Thus, the temporal resolution for these rotating targets provides spatially resolved information on how the atmosphere changes with varying irradiation from the WD host.Our target is SDSS-J155720.78+091624.7 (hereafter SDSS 1557), the fifth WD-BD system in our program, and the second-most irradiated.
This paper is organized as follows: First, in Section 2, we review the current knowledge on the WD-BD binary target of our study.Next, in Section 3 we present the new HST time-resolved observations and the data reduction.This is followed by the analysis of the broadband photometric and spectrally dispersed light curves, including model fits (Section 4).In Section 5, we analyze the rotational phase-resolved spectra, including the subtraction of the WD spectrum, extracting spectral modulations as a function of rotational phase, fitting a simple energy balance equation-based temperature distribution model, and also constraining the brightness temperatures in the day and night sides of the BD.In Section 6, we compare the extracted BD spectra to state-of-the-art irradiated atmosphere models.Finally, in Section 9, we review the main conclusions of this study.SDSS 1557 was first suggested to be a WD + debris disk system by Steele et al. (2011), who found Near-Infrared (NIR) flux excesses in the Y, J, and K photometry from the UKIRT Infrared Deep Sky Survey (UKIDSS) Large Area Survey Data Release 8. Excess K -band photometry was best-fit by 700 K blackbody, indicating a potentially cool disc.However, excesses in Y -and J -band were best matched by companion of spectral type L4, suggesting a BD companion.
Follow-up observations by Farihi et al. (2017) with Gemini+GMOS (2012) and VLT+XSHOOTER (2014)(2015)(2016) revealed that both sources of NIR excess were likely present within the system, making this the first known WD-BD binary with a metallic circumbinary debris disk.Farihi et al. (2017) used radial velocity changes in the Mg II 4482 Å absorption feature to reveal a period of 138.4 min, along with a predicted mass ratio of M 2 /M 1 = 0.14 via the Hα emission feature and assuming it originated from the companion.Further assuming that the Hα emission comes uniformly from the day side hemisphere yields a substellar companion mass of M 2 = 0.063 M ⊙ (66 M Jup ) at an orbital inclination of i ∼63 • .Note, higher orbital inclination corresponds to more equator-on viewing angles with spin-orbit aligned systems.
White dwarf effective temperature, T eff , and surface gravity, log(g), were derived by Farihi et al. (2017) from fitting atmospheric models to Balmer lines from the XSHOOTER spectra.These values, along with photometry, were cross referenced with evolutionary mod-els from Fontaine et al. (2001) to predict a cooling age of 33±5 Myr.Additionally, Farihi et al. (2017) found evidence in both the GMOS and XSHOOTER spectra that the metal polluted debris disk was actively accreting onto the WD primary, since detected metal lines such as Mg II and Ca II K should sink below the photosphere of the WD within a few weeks.Swan et al. (2020) modeled the emission of the SDSS 1557 system using a combination of two black bodies and found the model overpredicted the observed modulation amplitude at 4.5 µm.They argue that this could be due to flux dilution caused by the dust emission in the debris disk, or caused by molecular absorption and clouds in the companion's atmosphere.
Figure 1 shows a simple illustration of the SDSS 1557 binary system, along with the debris disk, at an imagined orbital inclination of 63 degrees.The three major bodies are labeled, along with their effective temperatures, or brightness temperatures in the case of the BD.Not illustrated is the process of accretion from the debris disk onto the WD.The WD's color came from (Harre & Heller 2021) and its published effective temperature of 21800 K.The color of the BD came from (Cranmer 2021) and the estimated day and night side temperatures derived in this work.Relevant properties this binary system, sourced from literature and this work, can be found in Table 1.

OBSERVATIONS & DATA REDUCTION
SDSS 1557 was observed with Hubble Space Telescope (HST ) using the WFC3/IR/G141 grism setup as part of the HST program GO-15947 (PI: Apai).Ten orbits split equally between two visits were scheduled for 2020 June.During this run, three orbits failed owing to guide star acquisition failures, one during the first visit and two during the second visit.Due to the failures, a third visit (five additional orbits) was scheduled for 2022 June with the same instrument setup as visits 1 and 2. Four out of five orbits were successful during this third run, leading to a total of 11 successful orbits for SDSS 1557 (summarized in Table 2).At the beginning of each orbit, two direct images were obtained for wavelength calibration using the F127M filter, a 256×256 subarray setup, and the GRISM256 aperture.After the direct images, eight spectroscopic exposures of 313 s each were obtained using the G141 grism, a 256×256 subarray setup, and the GRISM256 aperture.This observing sequence was successfully carried out for 4 other WD-BD systems as part of the same HST program: SDSS J141126.20+200911.1 (Lew et al. 2022), WD 0137-349, EPIC 212235321 (Zhou et al. 2022), and NLTT5306 (Amaro et al. 2023).
The spectral extraction process is largely similar to that described in Amaro et al. (2023), based on an established pipeline that combines the WFC3/IR spectroscopic software axe (Kümmel et al. 2009) and custom Python scripts.This pipeline has successfully extracted time-resolved observations of BDs (Buenzli et al. 2012;Apai et al. 2013) and WD-BD binaries (Zhou et al. 2022;Lew et al. 2022;Amaro et al. 2023).
We downloaded the CalWFC3 (version 3.5.2) pipeline product flt files for SDSS 1557 from the Barbara A. Mikulski Archive for Space Telescopes (MAST)1 .We ensured precise wavelength calibration for each orbit by using Source Extractor (Bertin & Arnouts 1996) on the median combined direct images from the beginning of each orbit.Then we took special care to manually check the positions of bad pixel flags in the G141 observations if they were near the grism spectrum.The same care was taken when checking the master sky image WFC3.IR.G141.sky.V1.0.fits.Configuration and reference files were downloaded from the aXe (hstaxe) WFC3 Extraction Example Cookbook2 .We then extracted the 1D spectra using a four pixel radius extraction window, resulting in a spectral resolution of ≈27 near 1.4 µm.
Next, we corrected the HST ramp effect -caused by charge trapping and delayed release -by fitting a RECTE model (Zhou et al. 2017) to the broadband light curve.The process for fitting the RECTE model is described in further detail in Amaro et al. (2023).To summarize, we use a custom Markov Chain Monte Carlo (MCMC) performed by emcee (Foreman-Mackey et al. 2013) to fit 6 RECTE parameters, two per visit.The results of the MCMC are presented in Table 3 and the final ramp model made with these results was divided out of each spectrum.
We then resampled our spectra to the resolving power of the observations using the 2D stamps output axecore from hstaxe.We plotted the counts of each pixel in the direction perpendicular to the dispersion axis and fit the shape of the count curves with a Moffat profile.The FWHM values were converted from pixel space to wavelength space, resulting in a spectral resolution of R∼215 at 14000 Å.The errors on the resampled spectra were calculated according to a method from Carnall (2017) (see Equation 1 in Amaro et al. 2023).
Finally, to remove unreliable data points, we performed wavelength cutoffs wherever the resampled spectral errors were greater than two times the average error value between 13,00 and 15000 Å.The resulting science spectra of SDSS 1557 is shown in Figure 2.

Broadband Light Curve
Studying the light curves of systems with tidallylocked companions reveals the presence, or absence, of any longitudinal intensity modulations that may be present within the photosphere.In the case of ultrashort period WD-BD binaries, we assume constant flux from the WD primary, with any variation on the surface being negligible compared to the BD companion.To begin characterization of the flux modulations coming from the rotation of the BD companion, we created simple phase curve models that should emulate most of the periodic behavior.
Once the RECTE ramp models were divided out and each spectrum was resampled, we determined the bestfit phase curve model to the broadband light curve with an MCMC.Similar to Amaro et al. (2023), the number of phase curve parameters varied according to the number of sine curves in each model.Specifically, we considered combinations of Fourier series from: We considered three combinations of k=1, 2, and 4 (i.e.N =1 is only k=1; N =2 is k=1 and 2; and N =4 is k=1, 2, and 4).The best fitting parameters and Bayesian information criterion (BIC) values are presented in Table 4.
Once models were created, we period-folded the data and calibrated phase = 0 to the time of the first observation (Table 2).A comparison of each period-folded phase curve for N = 1, 2, and 4 are shown in Figure 3, where the N =2 model was the preferred fit according to the BIC values.All best-fit phase curve model parameters and corresponding BIC values for the broadband light curve are presented in the top left section of Table 4.
In the N =2 broadband phase curve of SDSS 1557, the amplitudes were amp 1 = 0.105 ± 0.001 and amp 2 = 0.015 ± 0.001.We observed a narrow peak and a wider base.We believe this suggests a surface flux distribution with a compact, extremely hot substellar area, caused by the tidally-locked irradiation, and a broader area with lower temperature, caused by a low heat redistribution fraction.
An ephemeris for the orbit of SDSS 1557 was published in Farihi et al. (2017), where the authors use radial velocities from VLT+XSHOOTER observations from 20143 and 20164 .A precise ephemeris can be used to determine the presence of an absolute phase shift in the substellar hot spot, which would suggest the presence of winds strong enough to shift the hottest point in the atmosphere.However, when we extrapolate the published ephemeris to the date of our HST observations, the uncertainty grows to be of the same order of magnitude as our uncertainty in phase (0.02 phase).For future atmospheric observations of this binary system, we recommend observing precise radial velocities at or near the time of observation in order to quantify any hot spot shift.

Sub-band Light Curves
In addition to analyzing broadband data, studying light curves generated from sub-bands offers a valuable opportunity to investigate latitudinal brightness variations as a function of atmospheric pressure.This is because different wavelength ranges probe different atmospheric pressure levels.In brown dwarf studies, it is common to use J-and H-band light curve filters (e.g.Buenzli et al. 2012).Thus, we derived J-and H ′band light curves from our spectral observations.To create the sub-band J-band light curve, we convolved each spectrum with the Two-Micron All Sky Survey (2MASS) J-band transmission filter from the 2MASS All-Sky Data Release (Skrutskie et al. 2006) from 11000 Å to 13500 Å.We used the same method to create the H ′ -band light curve, except only we integrated from 15000 Å to the end of each spectrum (approximately 16700 Å), since the 2MASS H-band transmission filter extended beyond our spectra.In addition, we created a third sub-band light curve in the water-band by integrating from 13600 to 14800 Å. Figure 2 shows each filter response function and their respective wavelength ranges.
Real atmospheres are complex, and different wavelengths probe different pressures due to strongly wavelength-dependent gas and dust opacities (e.g.Buenzli et al. 2012;Yang et al. 2016).Thus, by comparing light curves from different filters, we can directly compare brightness distributions in different atmospheric layers.Similarities or differences between the light curves allow us to identify pressure-dependent structure and dynamics within the atmosphere.To quantify any similarities or differences between our three sub-band light curves (i.e., J-, water-, and H ′ -band), we created phase curve models using the same approach described in Section 4.1.
The best-fit parameters from all three models (N = 1, 2, and 4) for each sub-band light curve are presented in Table 4. Similar to the broadband light curve, the N =2 model yielded the best-fit.The best-fit amp 1 amplitudes for the Water-and H ′ -band phase curve models were 0.126 ± 0.002 and 0.133 ± 0.002, respectively.Interestingly, the J-band model had a significantly shallower amplitude of 0.084 ± 0.001.We see no evidence for relative phase shifts between the peaks of the sub-band light curves: ϕ peak (J) − ϕ peak (Water) = 0.001 ± 0.021; ϕ peak (J) − ϕ peak (H ′ ) = 0.005 ± 0.021.These phase shifts are relative to the start of the first observations, which is arbitrary.We are therefore interested in the difference between these phase values.Light curves and best-fit models are shown in Figure 4.

Spectroscopically-resolved Phase Curves
The shallower amplitude observed in the sub-band Jband phase curve suggested a potential wavelength dependence.To explore this possibility, we applied more precise wavelength bins in our analysis of the phase curve models derived from the spectra.The narrower wavelength bins had a width of ≈447 Å and offered an average signal-to-noise ratio (S/N) of 80. Building upon the success of the N =2 model for broadband and subband light curves, we focused solely on the N =2 model for these spectroscopically-resolved phase curves.
In Figure 5, the best-fitting amplitudes, shown as ∆F /F max , are presented as a function of wavelength.Within the range of 1.14 to 1.50 µm, the relative am-0.9 1.0  We observe no significant relative phase shift between the peaks or troughs of the phase curves.The Water-and H ′band phase curves both exhibit amplitudes greater than 12%, whereas the J-band phase curve exhibits a shallower amplitude around 8%.
plitude differences increased from ∼15 to ∼25 percent, with the greatest increase occurring between 1.32 and 1.44 µm.Past 1.50 µm, within the H ′ -band, the differences level-off, hovering around an average of 23.5 percent.
Spectroscopically-resolved phase curves and the analysis of wavelength dependence on their parameters, provides us a glimpse into the longitudinal-vertical circulation patterns within an atmosphere, as wavelength can be used a tracer for atmospheric pressure (e.g., Yang et al. 2016).By anchoring our interpretations using these methods, we can better understand the dynamics at play.In Figure 5, we observe that longer wavelengths generally exhibit stronger modulations.This correlation can be attributed to the BD, which is believed to be the primary contributor to the emission flux at longer wavelengths, due to relative temperature differences.

ROTATIONAL PHASE-RESOLVED SPECTRA
The ∼1.1 to 1.7 µm wavelength coverage of HST /WFC3/G141 is excellent for constraining the at- mospheric properties of irradiated brown dwarfs, since their flux contributions peak in the NIR.Further aided by the favorable flux contrast (∼10-30%) between the irradiated BD companion and WD primary, in NIR wavelengths, we are able to spectrally and spatially map irradiated BD atmospheres in great detail.Thus, we are also able to facilitate comparisons to current models of hot Jupiter and BD atmospheres.

Modeling the Debris Disk Contribution
Among our sample of six WD-BD systems, SDSS 1557 is unique due to a third component in the system, its circumbinary debris disk.Farihi et al. (2017) initially identified this disk and calculated that it should exist at an orbital radius of 3.3 R ⊙ at a temperature of T ≈1100 K.We modeled the expected flux contribution of the disk at our HST wavelengths, and determined that the disk excess emission at these wavelength is negligible compared to the brightness of the WD-BD binary, with F disk /F WD−BD ranging from 0.04 to 1.3%.Therefore, our analysis cannot yield significant insights into SDSS 1557's debris disk.For a comprehensive analysis on the disk structure, evolution, and role within the SDSS 1557 system, we refer to Farihi et al. (2017) and Swan et al. (2020).We note, that upcoming studies with the James Webb Space Telescope -with its far greater sensitivity at longer wavelengths -will likely provide powerful constraints on the disk composition and structure.

Modeling the White Dwarf Contribution
To study the atmosphere of the irradiated BD in more detail, we modeled and subtracted the WD contribution from our phase-resolved spectra.We began the modeling process by downloading a model grid for purehydrogen atmosphere WDs from Koester (2010), which features WD atmosphere models based on two parameters: T eff (effective temperature) and log(g) (surface gravity).We bilinearly interpolated along the grid to the published values for SDSS 1557 (see Table 1).Uncertainty in the adopted values was accounted for by creating 10,000 interpolated WD models, with T eff and log(g) values sampled from Gaussian distributions.For the Gaussian distributions, the means were set to 21800 K and 7.63 cm s −2 , with corresponding standard deviations of 800 K and 0.11 cm s −2 .To quantify the uncertainty, we measured the Full-Width-Half-Maximum (FWHM) of the histogram of flux values at each wavelength.These model uncertainties are represented as solid black error bars in Figure 6.
Next, we needed to scale the flux of our WD model to reproduce the expected flux at the distance and radius of our WD.The previously published value for distance was 520±35 pc (Farihi et al. 2017).Using more recent, photogeometric data from Gaia Data Release 3, we updated values for the radius and distance to R WD = 0.0162 ± 0.0012 R ⊙ and D = 500 +19.8 −18.0 pc (refer to Appendix A for further details on this calculation).The updated distance is within the 1σ uncertainty range of the previous value.We then multiplied our WD model by the scale factor R 2 WD /D 2 and propagated uncertainties associated with the WD's radius and distance.Consequently, there are two sources of uncertainty in modeling the WD contribution: the first arising during the model interpolation and the second introduced during the subsequent flux scaling.In Figure 6, the total uncertainty, considering both uncertainty sources, is represented by the gray shaded region.

Extracting Phase-resolved Spectra of SDSS 1557B
In order to understand the impact of external irradiation on an ultracool atmosphere, we set out to isolate the BD spectra from our observations of the nonspatially resolved system.To start, we focused on day vs. night hemispheres, as these phases should show the greatest differences in a tidally-locked system.Previous studies of ultra-short period WD-BD systems have shown brightness temperature differences between day and night hemispheres, which could be an indicator for differences in cloud coverage, atmospheric composition, or a combination of both.
For our SDSS 1557 data, we applied a wavelength shift to each spectrum to account for any radial velocity shifts before subtracting the WD model.We determined the magnitude of the radial velocity shift (v r ) at a certain phase using the following equation: Here, the phase (ϕ) of each spectrum was calibrated to the ephemeris and period from Farihi et al. (2017).The parameters γ 1 and K 1 represent the WD's radial velocity, while γ 2 and K 2 correspond to the brown dwarf's contribution.We adopted the values for γ 1 , K 1 , γ 2 , and K 2 from Farihi et al. (2017) as well.
Next, we performed a shift on each spectrum to its rest wavelength using the Doppler shift equation: Solving for λ rest in Equation 3, the largest wavelength shift observed in our sample was less than 20 Å; i.e. less than a third of our spectral resolution of 65.27 Å.
To ensure data quality and minimize spectral variations, we derived the BD spectra at representative phases by calculating the median of five spectra.For "Noon" and "Midnight", we used five spectra that corresponded to the five brightest and faintest broadband light curve points.For "Morning" and 'Evening", we identified the appropriate phase by adding and subtracting 0.25 from the phase with the lowest absolute flux in the broad band light curve.Then, we found the five spectra nearest to these new phases and took their median.The four phase-resolved spectra are presented in Figure 6.
We then proceeded to subtract the WD model and isolate the irradiated BD components, shown in Figure 7.The noon spectrum presents as a nearly featureless slope, likely due to the intense irradiation received by the WD primary.The midnight spectrum becomes nearly undetectable with flux levels near zero and S/Ns below 3σ for all wavelengths.The morning and evening spectra exhibit no significant features, but are both nearly half the flux of the noon spectrum.

Day and Night Brightness Temperatures of SDSS 1557B
The thermal structure in the atmosphere of a stellar companion is influenced by various factors, including internal heat flux, the profile of net absorbed stellar flux, and opacity structure (Marley & Robinson 2015).To explore the atmospheric pressure-temperature structure, we determined brightness temperatures (T B ) as a function of wavelength, i.e., the temperatures at which a blackbody of the brown dwarf's size would emit the same specific intensity as the observed flux values.Expressing emission spectra as T B is widely used in BD and exoplanet studies, as well as previous works on our targets (Casewell et al. 2015(Casewell et al. , 2018;;Lew et al. 2022;Zhou et al. 2022;Amaro et al. 2023).By combining brightness temperatures with a Pressure-Temperature (P-T) profile, we can effectively map the relative pressure regions associated with each wavelength, revealing the vertical structure of the atmosphere during each observed phase (Yang et al. 2016).
Figure 8 presents the calculated brightness temperatures versus wavelength for the four extracted spectra of SDSS 1557B.In the noon spectrum of SDSS 1557B, T B decreases nearly monotonically towards longer wavelengths, from 2600 K at 1.12 µm to 2250 K at 1.56 µm.We do observe two increases in T B that deviate from the monotonic decline: the first centered around 1.5 µm and the second starting at 1.56 µm until 1.67 µm.The morning and evening spectra exhibit similar slopes as the noon spectrum, and are qualitatively the same except for deviations at 15000 Å and 16500 Å.
The night side of SDSS 1557B is particularly faint and did not result in confident detections: At its original resolution of approximately 65 Å per data point, the spectrum primarily consisted of non-detections.We calculated upper limits for brightness temperatures by binning the spectra in wavelength, thus increasing the signal-to-noise ratio (S/N) for each bin.A S/N value greater than or equal to 3 was considered a detection, while values greater than or equal to 2 were suitable for establishing upper limits.We gradually increased the bin size in the spectrum and ultimately settled on a spectrum with a resolution of 965 Å per data point, as this was the only binning option that yielded more than one data point with an S/N of 2 or higher.The upper limits for those two data points are shown in Figure 8.
Given the temperature estimate, the day side is likely too hot to host silicate clouds (Lunine et al. 1986(Lunine et al. , 1989;;Burrows & Sharp 1999;Marley et al. 2002), but H − opacity may contribute to the spectral slope (Arcangeli et al. 2018).This scenario would argue for the absence of solid-state silicate absorption on the day side and a temperature inversion.However, in the morning and evening longitudes, the Si-O solid state (stretchingbending) absorption feature may be present -testing this will require observations at longer wavelengths, near 10 µm (e.g., Apai et al. 2005;Luna & Morley 2021;Suárez & Metchev 2022).A known signature of temperature inversion is the Paschen series atomic hydrogen n=5→3 emission feature at 1.28 µm.This signature may be present in the day.
Based on the upper limits for midnight brightness temperatures, temperatures exceeding 1500 K between 12000 and 14000 Å are not anticipated.If we assume the actual midnight brightness temperatures at a higher resolution exhibit a comparable pattern to those in other phases, it is reasonable to infer the presence of silicate cloud coverage in this hemisphere.

Irradiated Models
We compared the phase-resolved spectra of SDSS 1557B to a number of PHOENIX atmosphere model spectra.These models follow the same set-up as those of Lothringer & Casewell (2020a).In short, the companion atmosphere is irradiated by a WD whose spectrum is modeled by the Koester (2010) grid, with fundamental parameters from (Farihi et al. 2017).An iterative process is conducted to calculate the BD's thermal structure until it reaches radiative equilibrium.We run the model at a wavelength sampling of 1 Å from 10 Å to 5µm, with coarser sampling out to 1000 µm, with the irradiation spectrum defined between 0.09 and 3 µm.The models used a rainout abundance scheme, where elements that participated in condensation lower in the atmosphere were depleted in the layers above (or at lower pressure than) the condensation level of that species.
We first started with a small grid of day side heatredistribution models (T irr = 2900 K) at solar metallicity with a range of surface gravity (log(g) cgs = 4.95, 5.2, and 5.45) and internal temperatures (T int = 1000, 1500, 2000 K).To test hotter and cooler atmospheres, with the latter corresponding to other orbital phases than the day side, we tested models with irradiation temperatures of (T irr =1800, 1950, 2050, 2450, 2900, 3100, 3200 K).
In general, the models were difficult to converge because most of the irradiation was in the NUV due to the high effective temperature of the white-dwarf host, resulting in much of the irradiation being absorbed in a relatively small part of the atmosphere.This heating at a relatively discrete layer of the atmosphere results in a sharp temperature inversions at low pressures, a more extreme case of the strong inversions than those seen in ultra-hot Jupiters around early-type main sequence host stars (Lothringer & Barman 2019).At some irradiation temperatures, the temperature profile of the atmosphere also crossed important condensation boundaries, when certain short wavelength absorbers will begin to rain out of the atmosphere.This can cause numerical oscillations in the model, where the lower atmosphere will begin to warm without those important upper atmosphere absorbers, causing those same absorbers to evaporate back into the gas phase, and cooling the lower atmosphere through an anti-greenhouse effect.This cycle can continue and cause instabilities in the numerical convergence of these models towards radiative-convective equilibrium.
The morning and evening hemispheres also remaining featureless was somewhat unexpected.Based on the measured brightness temperatures, the terminator regions appeared to be between 1800-2100 K, cool enough to avoid H 2 O thermal dissociation and for H − to no longer dominate the opacity.A couple possible scenarios could explain the featureless morning and evening spectra.First, if the brightness temperatures are underestimated by 500 K, then the morning and evening hemispheres could remain warm enough for H − opacity to remain dominant.
Alternatively, the morning and evening hemispheres may become increasingly isothermal as the temperature inversion induced by the extreme irradiation on the day side hemisphere begins to wane at longitudes further from the substellar point.Emission spectra in this region would thus remain featureless.One would then expect the undetected night side to finally show evidence of H 2 O absorption as the temperature structure becomes fully non-inverted.This behavior is qualitatively similar to results from the HST/WFC3/G141 phase curve of ultra-hot Jupiter WASP-121b, which has a nearly identical equilibrium temperature as SDSS 1557 (T eq = 2350 K, Mikal-Evans et al. 2022).
Another scenario to explain the morning and evening spectra and brightness temperatures would be dilution from the still-visible portion of the day side hemisphere.The portion of SDSS 1557B visible during the morning and evening phases will include a combination of the illuminated day side and the non-illuminated night side.Thus, the flux measured at the morning and evening phases might be better represented as a combination of a representative day side and night side spectrum, or by diluted day side spectrum (Taylor et al. 2020).This latter case is what is effectively done when comparing our model spectra to the normalized observations, rather than comparing absolute fluxes.

Non-Irradiated Models
It is also illustrative to compare our extracted irradiated BD spectra to models of cloudy, non-irradiated BD atmospheres to get a baseline of what this system may have looked like without its WD primary.Thus, we conducted a comparative analysis with a set of nonirradiated models specifically designed for the study of L-type BDs (Brock et al. 2021).As with the irradiated models in Section 6.1, this grid was calculated with PHOENIX, but unlike the irradiated models, assumes local thermodynamic equilibrium.The grid covers T eff from 800 K to 2100 K and surface log(g) from 3.5 to 5.5 (cgs units).Briefly, these models include the effects of non-equilibrium chemistry (parameterized by the eddy diffusion coefficient K zz ) and the effects of cloud opacity.The cloud location, composition and particle number densities are initially determined by equilibrium chemistry.
This initial condensate distribution is then constrained by two free parameters, the pressure at the cloud top and the mean particle size (assuming a lognormal size distribution).Above the cloud top, condensate number density decreases exponentially, while the cloud base is established by equilibrium chemistry.The grid coarsely samples models with cloud top pressures (P c ) ranging from 0.1 to 10 bar and mean grain sizes of 0.25, 0.5 and 1 µm.The Brock et al. (2021) study was focused on a specific set of BD binaries and, consequently, the grid does not uniformly sample all the parameters.However, for the temperatures estimated above, the grid provides good coverage of the expected T eff and surface gravity and offers a variety of cloud scenarios to explore.with best-fit models from the Lothringer model (LM) grid (Section 6.1) and the Brock & Barman Model (BBM) grid (Section 6.2).Best-fitting blackbody curves are also shown.Ignoring Bond Albedo, the irradiation temperature at the substellar point for SDSS 1557B is 2540 K. Neither model grid is able to match the slope or features of each extracted spectrum.
In Figure 9, the best-fitting cloudy, non-irradiated models for each hemisphere integrated spectrum are labeled as BBM (for Brock and Barman Model).For each phase, Noon, Evening, and Morning, the same nonirradiation model was the closest match.As with the irradiated models, the observed data showcases a slope that is steeper than anticipated by the cloudy, non- irradiated atmosphere models.Our extracted spectra also appears more featureless, as one might expect from the constant irradiation of the primary.

LONGITUDINAL TEMPERATURE DISTRIBUTION
One of the key atmospheric processes affecting the structure of tidally-locked irradiated atmospheres is the irradiation redistribution fraction, or efficiency of the heat redistribution from day to night.The strength of this fraction is thought to be determined by a competition between radiative losses and advective heat transport (Showman & Guillot 2002;Perez-Becker & Showman 2013a;Komacek & Showman 2016;Zhang & Showman 2018).In theory, a hotter planet should be more efficient at losing heat through radiation, effectively lowering the irradiation redistribution fraction, but this might be balanced by a stronger atmospheric circulation.This trend changes for atmospheres with equilibrium temperatures higher than 2300 K, where hydro-gen dissociation/recombination should help with dayto-night heat transport (Bell & Cowan 2018;Tan & Komacek 2019).In addition, extremely rapid rotation also tends to inhibit day-to-night heat transport, as appropriate for the close-in BD around WD of our case (Tan & Showman 2020).With our phase-resolved spectroscopy of SDSS 1557B, we utilized a simple atmosphere to constrain the strength of these atmospheric processes.
Using the calculated brightness temperatures and available physical parameters of the system (Table 1), we generated a temperature distribution map of SDSS 1557, meant to constrain four key atmospheric processes that impact the temperature distribution.The construction of this map utilized a grid search involving the following parameters: Bond albedo (A B ), which determined the fraction of incident irradiation reflected away from the day side; irradiation redistribution fraction (f irr−red ), which governed how much non-reflected irradiation is taken away from the day side and evenly redistributed throughout the sphere (i.e., 0.0 = no redistribution; 0.5 = half of non-reflected irradiation flux was uniformly redistributed; 1.0 = total redistribution with uniform temperatures across all latitudes and longitudes); nonirradiated BD temperature (T non−irr ), which established the base temperature of the BD without any irradiation; and lastly, Inclination (i), which controlled the visibility of one hemisphere when the other hemisphere dominated the emission.In our set up, we assumed alignment between the BD's rotation inclination and the orbital inclination.This process was also deployed for NLTT5306B (Amaro et al. 2023).
Here, we briefly explain the steps involved in constructing this map.First, we created a sphere of uniform temperature, determined by T non−iir , meant to represent the BD temperature in the absence of external irradiation.We then introduced a source of external irradiation, in this case SDSS 1557A, where the irradiation strength was governed by T WD , orbital separation a, and R WD .Immediately, a fraction of this external irradiation was reflected away, a consequence of A B .At this point, day side temperatures were combination of T non−irr , A B , and external irradiation from the WD primary.Night side temperatures are still set at the initial T non−irr value.Quantitavely, day and night temperatures at this step are determined by: where L WD = 4πR 2 WD σ SB T 4 WD , σ SB is the Stefan-Boltzmann constant, and f rad is the fraction of the surface that reradiates (set to 1.0 in our model).
Then, we parameterized the global circulation -driven by the day-night temperature difference -with the parameter f irr−red , which represents the fraction of energy transported from the day side to the night side.In our simple model, we subtracted a percentage of the day side temperatures from each day side cell, summed them, and then redistributed this flux across the whole sphere.Essentially, this acted to decrease the day side temperatures while simultaneously increasing night side temperatures and maintaining a smooth boundary at the day-to-night terminators.Finally, we applied an inclination to the toy atmosphere.To determine the best-fit combination of parameters, the average temperature of each hemisphere (T H M ) was then compared to the respective brightness temperature (T H B ) via the following: where σ H T B is the uncertainty in the brightness temperature.
Without a detectable midnight brightness temperature, we considered phases ranging from the morning through the noon to the evening to constrain the temperature distribution.Since our current model is symmetric around the substellar point, we fit our model to the average of the morning and evening brightness temperatures in the appropriate phases.However, since morning and evening hemispheres also include half of the day side hemisphere, we decreased the weight of these phases in our fits by 50%.The final best-fit combination of parameters was thus decided by the lowest value of x noon + (0.5 × x morn ), where x morn is interchangeable with x eve in our symmetric atmosphere model.
To quantify uncertainties, we initially marginalized the distributions of each parameter into a onedimensional array.Then, to establish the 1σ range, we added 1 to the lowest value within the array and determined the two points at which this adjustment takes place.These became the upper and lower 1σ values.
Figure 10 shows the 2D longitudinal temperature distribution with best-fit values labeled in the lower left: A B = 0.18 +0.14 −0.04 , f irrred = 0.10 +0.14 −0.16 , T non−irr = 1190 +121 −120 , i = 80.0 +3.6  −3.3 .Given the high day side brightness temperatures and likely regions of molecular dissociation near the substellar point, the low Bond albedo is physically reasonable within its uncertainties.
The best-fit inclination of i = 80 +3.6 −3.3 degrees is higher than the previously published value of 62±3 degrees in Farihi et al. (2017).While the two inclination values are within 3σ uncertainty, we believe the higher value found is more accurate, given the higher-quality data and advanced phase-resolved spectroscopy.The enhanced precision afforded by our methodology leads us to have more confidence in the reliability of our result.Additionally, an inclination of 80% is still physically possible, given the orbital geometry of the system.

FUTURE OBSERVATIONS
Due to their temperatures, the spectra of cool brown dwarfs and many hot Jupiters peak in the near-infrared (typically between 1-3 µm).Therefore, our nearinfrared observations are particularly constraining for BD atmosphere models.The work described here provided insights into key atmospheric processes within ultracool atmospheres (e.g., vertical mixing, circulation, photochemistry, and condensate clouds).Additionally, through some mismatches between the state-of-the-art models and the observations, we identified some of the current limitations of brown atmospheric models.
With its superior infrared sensitivity and higher spectral resolution, the James Webb Space Telescope (JWST) could offer complementary infrared observations of these rare systems.Already, JWST has enabled orders-of-magnitude increases in the information content of atmospheres for hot Jupiters (Bell et al. 2023;Mikal-Evans et al. 2023) and BDs (Miles et al. 2023;Manjavacas et al. 2024;Lew et al. 2024).Recently, our team was awarded a Cycle 3 JWST proposal that will observe five rare WD-BD systems (JWST GO-4967; PI: Y. Zhou).These observations will reveal processes at a much greater range of pressures, allowing for a more comprehensive comparison between observations and models.In the case of SDSS 1557, phase-resolved spectroscopy from JWST will enable tighter longitudinal constraints on cloud coverage and possible thermal inversions as well as further characterization of circumbinary debris disk that is actively polluting the primary WD (Farihi et al. 2017).

CONCLUSIONS
In this study, we presented time-resolved, spectrophotometric phase-mapping of the highly irradiated BD SDSS 1557.The key findings of our study are as follows: • We presented high-quality HST /WFC3/G141 phase-resolved spectra of WD-BD binary SDSS 1557.The observations sampled 100% of the rotational phase of the BD and exhibited strong broadband photometric variations around 21% from peak-to-trough.
• We modeled the broadband light curve of SDSS 1557 to within 2.5%.Best-fit model consisted of k = 1 and 2 terms with corresponding amplitudes of amp 1 = 0.105 ± 0.001 and amp 2 = 0.015 ± 0.001.
• We reveal a strong wavelength dependence on amplitude using spectroscopically-resolved phase curves, where the amplitudes of both the k=1 and k=2 terms generally increase with wavelength.
The highest amplitude for both terms occurs around 1.50 µm.
• After modeling and subtracting the WD contribution, we extracted phase-resolved spectra of the irradiated BD SDSS 1557B.Due to inclination, we estimated up to 15% of the day side flux is visible from the night side.After subtracting 15% of the day side from the night side, the midnight phase spectrum was nearly undetectable.
• The hemisphere averaged brightness temperatures for Noon, Morning, and Evening phases exhibit wavelength dependence, with temperatures that monotonically decrease with wavelength.Upper limits for the midnight phase suggest brightness temperatures no greater than 1600 K.
• The noon-phase spectrum of SDSS 1557B is wellapproximated by the hottest models computed.However, best-fitting models for the morning and evening phases failed to capture the steeper spectral slopes and lack of spectral features.This could be explained by an upper atmosphere that becomes increasingly isothermal at longitudes farther from the equator, as well as the morning and evening phases being diluted versions of the day side hemisphere.
• Phase-resolved spectroscopy from JWST will likely capture the night side hemisphere, providing a more complete map of the atmosphere and more precise explanations for the morning and evening spectra of SDSS 1557B.WD Radius: In order to scale the flux of the WD model to the expected flux at the distance and radius of SDSS 1557A, we needed an estimate for WD radius.We solved the bolometric luminosity equation for radius: where T eff is the published effective temperature of 21800 ± 800 K and σ SB is the Stefan-Boltzmann constant.We determined L Bol from M Bol , the bolometric magnitude.like so: where X is the specific photometric filter.From a literature search, we found a bolometric correction for WDs in the Johnsons-Cousins V -band (Carrasco et al. 2014).Using data from Table 3 in Carrasco et al. (2014), we created a grid of V -band bolometric corrections, and bicubicly interpolated along T eff and log(g).Using this method, we calculate a BC V = -2.2384for SDSS 1557A.
With an estimate for bolometric correction, we then needed an absolute magnitude in the V -band, which can be calculated from apparent magnitude with: where D is the distance.For this variable, we adopted the Gaia photogeometric distance of D = 500 +19.8 −18.0 pc from Bailer-Jones et al. (2021).For V -band apparent magnitude, we use a photometric transformation5 from Gaia G, GBP , and GRP filters to Johnson-Cousins filters.Specifically, we used the equation: where c 0 = -0.02704,c 1 = 0.01424, c 2 = 0.2156, and c 3 = 0.01426.With apparent Gaia magnitudes of G = 18.609 ± 0.004, GBP = 18.519 ± 0.024, GRP = 18.810 ± 0.042, we solved for V, yielding an apparent V -band magnitude of m V = 18.659 ± 0.008 for SDSS 1557A.Thus, we calculated an absolute V -band magnitude of M V = 10.164 ± 0.241.Combining M V with BC V , we calculate a bolometric magnitude of M Bol = 7.926 ± 0.241, which translates to a bolometric luminosity of L Bol = 2.035×10 32 ± 3.057×10 28 .Finally, we solved for the radius of the WD: R WD = 0.0162 ± 0.0012 R ⊙ .
BD Radius: We estimate the radius of the irradiated BD SDSS 1557B using evolutionary models from the Sonora Bobcat model grid (Marley et al. 2021).The evolutionary portion predicts how parameters such as mass and radius evolve over time.For our system, we used the evolutionary tables with solar metallicity ([M/H]=0.0),given the compelling evidence for ongoing accretion from a debris disk that lies outside of the system (Farihi et al. 2017).
Similar to the method used in Beiler et al. (2023), we randomly draw a million pairs of mass and age values, then interpolate across the models for each pair to generate a range of possible radii.For the input masses and ages, we chose uniform distributions with the following boundary conditions: 60 < M BD < 75 M Jup and 0.03 < Age BD < 10 Gyr.We adopted the resulting range of BD radii, R BD = 0.106±0.024R ⊙ , as the radius of SDSS 1557B (see Table 1).Distance: Original value from Farihi et al. (2017) was 520±35 pc.Farihi et al. (2017) estimated WD parameters T eff and log(g) from fitting atmospheric models to Balmer lines in 60 individual XSHOOTER spectra.Published errors for those parameters were the standard deviation from all 60 measurements.Then, errors in derived stellar parameters, like distance, were calculated by propagating the uncertainties in the adopted T eff , log(g), and published photometry through WD evolutionary models (Fontaine et al. 2001).
We adopted an updated distance value, D = 500 +19.8 −18.0 pc, from the Bailer-Jones et al. ( 2021) distance catalog for WDs.For most WDs in this catalog, the authors inferred two types of distances: (1) geometric, which uses parallax with a direction-dependent prior on distance, and (2) photogeometric, which additionally uses the color and apparent magnitude of the WD.SDSS 1557 had both distances published, with D = 590.1 +74.3  −64.6 pc as the geometric distance, and D = 500.0+19.8  −18.0 pc as the photogeometric distance.We adopted the photogeometric distance in our study, as this is a more reliable estimate.

B. BEST-FIT LONGITUDINAL TEMPERATURE DISTRIBUTION
We conducted a grid search among 4 parameters: (1) Bond Albedo, A B =[0.0, 1.0] in increments of 0.025, (2) Irradiation redistribution fraction, f irr−red =[0.0, 1.0] in increments of 0.025, (3) non-irradiated brown dwarf temperature, T non−irr =[1000,1400] in increments of 10 K, and (4) Inclination,i=[40,110] in increments of 5 degrees.The method for finding the best-fit combination of parameters is largely similar to what is described in Section C of Amaro et al. (2023), with one key exception.In NLTT5306B, we compared the hemisphere integrated day and night side temperatures of the simple atmosphere model to the day and night brightness temperatures derived from the extracted BD spectra.For SDSS 1557B, the extracted midnight spectrum was just below significant detection levels.Thus, we altered the

Figure 3 .
Figure 3.Comparison of best fitting Fourier Series models (solid black lines) for orders N =1, 2, and 4. Each of the 11 HST orbits are represented in a different color.Based on BIC values, the N =2 model is marginally favored for SDSS 1557, with a model that captures light curve behavior to within 2.5%.

Figure 4 .
Figure 4. Sub-band phase curves made in the 2MASS J-(indigo triangles), Water-(coral circles), and 2MASS H ′bands (yellow squares).Best-fit models for each filter are presented as well as precise phases for each peak and valley (phase uncertainty shown as horizontal bars on arrows).We observe no significant relative phase shift between the peaks or troughs of the phase curves.The Water-and H ′band phase curves both exhibit amplitudes greater than 12%, whereas the J-band phase curve exhibits a shallower amplitude around 8%.

Figure 5 .
Figure5.Relative phase curve amplitudes as a function of wavelength for SDSS 1557.In this case, relative refers to the total change in flux relative to the maximum flux.Overall, the amplitude changes tend to increase with increasing wavelength.

Figure 6 .
Figure 6.The modeled flux contribution of the WD SDSS 1557A, shown with both sources of uncertainty.The first being in the creation of the model (black errorbars) and the second being in the flux scaling (gray shaded region).Additionally, the phase-resolved spectra of the irradiated BD SDSS 1557B are shown at 4 key flux phases: Noon, Evening, Midnight, and Morning.

Figure 7 .
Figure 7. Derived phase-resolved NIR spectra of irradiated BD SDSS 1557B.Vertical errorbars represent systematic uncertainty originating from the creation of the WD model.The size of scaling uncertainties for both hemispheres is labeled near the upper right.Solid lines and shaded regions depict our observations and uncertainties convolved with a Gaussian kernel with a width equal to the spectral resolution of the G141 grism.

Figure 8 .
Figure 8. Brightness temperatures of SDSS 1557B, calculated by solving the Planck equation for TB at each wavelength.Errorbars for the noon spectrum are shown as the shaded region, whereas errors for evening and morning are excluded for easier viewing.Errors for evening and morning are 95 K average.

Figure 9 .
Figure9.Comparison of extracted spectra of SDSS 1557B with best-fit models from the Lothringer model (LM) grid (Section 6.1) and the Brock & Barman Model (BBM) grid (Section 6.2).Best-fitting blackbody curves are also shown.Ignoring Bond Albedo, the irradiation temperature at the substellar point for SDSS 1557B is 2540 K. Neither model grid is able to match the slope or features of each extracted spectrum.

Figure 10 .
Figure 10.Top: A 2D depiction of the temperature distribution on SDSS 1557B projected onto a 3D sphere.Relative radii of WD and BD are to scale, whereas orbital separation is zoomed in for clarity.The WD primary, SDSS 1557A, is represented as a pale blue dot in the lower left.Bottom: A 2D projection of the toy model above with labeled best-fit values and 1σ uncertainties for Bond albedo, external irradiation redistribution fraction, non-irradiated BD temperature, and inclination.

Figure A1 .
Figure A1.Residual curves for each parameter explored in the simple heat redistribution atmosphere model shown in Figure 10.Best-fit values, which are the lowest points in the curves, are highlighted by the vertical line in each panel, while 1σ uncertainties are shown as the shaded regions.

Table 1 .
Properties of the SDSS 1557 binary system

Table 2 .
Log of HST WFC3 Observations for SDSS 1557

Table 4 .
Best-fit Phase Curve Parameters Note-The k=2 phase curve was always best-fit according to the BIC values.
X.T. provided guidance regarding the intellectual direction of analysis.T.B. provided existing atmosphere models for comparison with our data.All authors, including L.M., M.M., and V.P., contributed guidance and suggestions to improve the manuscript. and