Measuring the variability of directly imaged exoplanets using vector Apodizing Phase Plates combined with ground-based differential spectrophotometry

Clouds and other features in exoplanet and brown dwarf atmospheres cause variations in brightness as they rotate in and out of view. Ground-based instruments reach the high contrasts and small inner working angles needed to monitor these faint companions, but their small fields-of-view lack simultaneous photometric references to correct for non-astrophysical variations. We present a novel approach for making ground-based light curves of directly imaged companions using high-cadence differential spectrophotometric monitoring, where the simultaneous reference is provided by a double-grating 360{\deg} vector Apodizing Phase Plate (dgvAPP360) coronagraph. The dgvAPP360 enables high-contrast companion detections without blocking the host star, allowing it to be used as a simultaneous reference. To further reduce systematic noise, we emulate exoplanet transmission spectroscopy, where the light is spectrally-dispersed and then recombined into white-light flux. We do this by combining the dgvAPP360 with the infrared ALES integral field spectrograph on the Large Binocular Telescope Interferometer. To demonstrate, we observed the red companion HD 1160 B (separation ~780 mas) for one night, and detect $8.8\%$ semi-amplitude sinusoidal variability with a ~3.24 h period in its detrended white-light curve. We achieve the greatest precision in ground-based high-contrast imaging light curves of sub-arcsecond companions to date, reaching $3.7\%$ precision per 18-minute bin. Individual wavelength channels spanning 3.59-3.99 $\mu$m further show tentative evidence of increasing variability with wavelength. We find no evidence yet of a systematic noise floor, hence additional observations can further improve the precision. This is therefore a promising avenue for future work aiming to map storms or find transiting exomoons around giant exoplanets.

. An on-sky example image of a ∼ 7 mag target observed with the dgvAPP360 coronagraph on the Large Binocular Telescope (LBT), copied and annotated on the right. The dgvAPP360 produces a dark hole of deep flux suppression, seen here surrounding the target star, allowing high-contrast companions to be detected in this region. The star itself remains unsaturated in the middle, allowing it to be used as a simultaneous reference PSF when making a differential light curve for a companion. The outer edge of the dark hole lies beyond the field of view in the LBT observations used in this work. and giant storms that show great diversity in size, shape, lifetime, and brightness (e.g. Simon et al. 2016Simon et al. , 2021Stauffer et al. 2016;Ge et al. 2019;Coulter et al. 2022). These features rotate in and out of view throughout the planet's rotation period, modulating its overall brightness and thus allowing us to map out its atmosphere (e.g. Kostov & Apai 2013;Karalidi et al. 2015Karalidi et al. , 2016Fletcher et al. 2016;Apai et al. 2017;Plummer & Wang 2022). Beyond the solar system, variations in the light curves of stars deliver information on the distribution of features such as star spots (Barnes et al. 2002;Jeffers et al. 2007;Frasca et al. 2009;Strassmeier 2009;Morales et al. 2010;Goulding et al. 2012;Herbst 2012;Nielsen et al. 2013;Park et al. 2021;Thiemann et al. 2021). By measuring the photometric variability of exoplanets and brown dwarfs in the same way, we can gain not only an insight into their visual appearance, but also key information on the distribution of condensate clouds that strongly affect the infrared spectra of directly imaged companions, allowing degeneracies between atmospheric models to be broken (e.g. Yang et al. 2016;Rajan et al. 2017;Charnay et al. 2018;Zhou et al. 2019;Zhang 2020;Tan & Showman 2021;Ward-Duong et al. 2021). Space-based photometric monitoring with the Hubble Space Telescope (HST) has already shown that giant planetary-mass and brown dwarf companions do exhibit such variability, at a range of amplitudes and periods Buenzli et al. 2014Buenzli et al. , 2015aZhou et al. 2016Zhou et al. , 2020aManjavacas et al. 2018Manjavacas et al. , 2019bMiles-Páez et al. 2019;Bowler et al. 2020b;Lew et al. 2020a). These results are in good agreement with observations of isolated brown dwarfs and giant exoplanet analogues Vos et al. 2018Vos et al. , 2020Lew et al. 2020b;Ashraf et al. 2022), including a large Spitzer survey by Metchev et al. (2015) who found that photospheric spots causing ≥0.2% variability at 3-5 µm are ubiquitous. Several studies have identified objects with much stronger variability, at the >10% level, with some even varying with peak-to-peak amplitudes as high as 26% (e.g. Radigan et al. 2012;Wilson et al. 2014;Biller et al. 2015;Eriksson et al. 2019;Bowler et al. 2020b). Vos et al. (2022) further found that young, low-mass brown dwarfs with similar colours and spectra to directly imaged exoplanetary companions are highly likely to display variability in the L2-T4 spectral type range, with an enhancement in maximum amplitudes compared to field dwarfs.
The rotation periods of brown dwarf and planetary-mass companions are consistent with those of the isolated low-mass brown dwarf population, suggesting that they may share similar angular momentum histories (Bryan et al. 2018;Vos et al. 2022). These periods are generally short, ranging from ∼1 hour to 20 hours (e.g. Radigan et al. 2014;Metchev et al. 2015;Apai et al. 2021;Tannock et al. 2021), within the range expected when evolutionary models and the age-and mass-dependent breakup velocities are considered (Leggett et al. 2016;Vos et al. 2020). These periods, derived from photometric measurements, are complementary to measurements of companion spin obtained from their spectra Schwarz et al. 2016;Bryan et al. 2018Bryan et al. , 2020bXuan et al. 2020;Wang et al. 2021b). When combined, rotation period and spin measurements can be used to constrain companion obliquities (Bryan et al. 2020a(Bryan et al. , 2021. However, the population of directly imaged companions accessible to space-based facilities such as HST remains small as most companions lie at close angular separations within the inner working angles of these facilities. Equipped with coronagraphs and extreme adaptive optics (AO) systems operating in the infrared, large ground-based observatories have the resolution and photon collecting power needed to overcome the glare of the host star and reach the high contrasts and close angular separations of substellar companions currently inaccessible to space telescopes (Bowler 2016;Hinkley et al. 2021;Currie et al. 2022). Although this provides the opportunity for variability studies of such companions, precise photometric monitoring is difficult as the companion light curve is contaminated by variability arising from Earth's atmosphere and other systematics. Therefore, a simultaneous, unsaturated photometric reference is required to remove this contaminant variability from the companion light curve. For non-coronagraphic, ground-based observations of isolated brown dwarfs and planetary-mass objects, nonvariable stars present in the field of view have often been used as photometric references to enable many successful measurements of variability (e.g. Gelino et al. 2002;Biller et al. 2013Biller et al. , 2015Girardin et al. 2013;Radigan et al. 2014; Wilson et al. 2014;Naud et al. 2017;Eriksson et al. 2019;Vos et al. 2019;Manjavacas et al. 2021Manjavacas et al. , 2022. However, the typically narrow field of view of ground-based coronagraphic imagers generally precludes the use of field stars as photometric references for observations of companions, and widely used focal-plane coronagraphs block the host star to enable the detection of the companion (Soummer 2005;Mawet et al. 2012;Ruane et al. 2018).
One solution to this problem is to use off-axis satellite Point Spread Functions (PSFs), or satellite spots, which can act as simultaneous photometric references even when a host star is blocked by a coronagraph (Marois et al. 2006b;Sivaramakrishnan & Oppenheimer 2006). Satellite spots can be created by adding a periodic modulation to the deformable mirror of an AO-equipped telescope or by placing a square grid in the pupil plane to produce spots through diffraction of starlight (Langlois et al. 2013;Wang et al. 2014;Jovanovic et al. 2015b). The former approach has been used by Apai et al. (2016), Biller et al. (2021), and Wang et al. (2022) for observations of the multi-planet HR 8799 system (Marois et al. 2008(Marois et al. , 2010. The first two of these studies observed HR 8799 with the Spectro-Polarimetric High-contrast imager for Exoplanets REsearch (SPHERE; Beuzit et al. 2019) instrument at the Very Large Telescope (VLT), and the latter used the CHARIS integral field spectrograph (IFS) in combination with the Subaru Coronagraphic Extreme Adaptive Optics instrument at the Subaru Telescope (Jovanovic et al. 2015a;Groff et al. 2017). Biller et al. (2021) used the satellite spots with a broadband-H filter to successfully constrain their sensitivity to variability to amplitudes >5% for HR 8799b for periods <10 hours, and amplitudes >25% for HR 8799c for similar periods, noting that the observed amplitude of any variability would be muted by the likely pole-on viewing angle of these planets (Vos et al. 2017;Wang et al. 2018;Ruffio et al. 2019). They also rule out non-shared variability between HR 8799b and HR 8799c at the <10-20% level over a 4-5 hour timescale by using one planet as a photometric reference for the other. Using a spectrophotometric approach, Wang et al. (2022) further constrained the variability amplitudes of HR 8799c to the 10% level, and HR 8799d to the 30% level, and found that there was no significant variability in the planet's colours. However, all three studies found that satellite spots are anti-correlated with each other and can demonstrate individual variations of their own, potentially setting a limit to the precision that can be achieved with this technique (although Biller et al. (2021) and Wang et al. (2022) note that the satellite spot light curves can be flat in their most stable epochs).

Ground-based differential spectrophotometry
In this paper we present a novel, alternative ground-based approach for constructing light curves of high-contrast companions directly through the technique of differential spectrophotometric monitoring, akin to that used highly successfully to study exoplanet transmission spectra (e.g. Kreidberg et al. 2014;Stevenson et al. 2014;Wilson et al. 2020;Arcangeli et al. 2021;Panwar et al. 2022a,b). We use a double-grating 360°vector Apodizing Phase Plate (dg-vAPP360) coronagraph (Doelman et al. 2017(Doelman et al. , 2021Wagner et al. 2020), which enables high-contrast companions to be detected without blocking the host star, hence leaving an unsaturated image of the host star available for use as a simultaneous reference. The more widely used grating vector Apodizing Phase Plate (gvAPP) coronagraph adjusts the phase of the incoming wavefront to modify the PSFs of all objects in the field of view, creating two images of the target star each with a 180°D-shaped 'dark hole', a region of deep flux suppression in which high-contrast companions can be observed, on opposing sides (Snik et al. 2012;Otten et al. 2014aOtten et al. , 2017Bos et al. 2020;Doelman et al. 2021;Sutlieff et al. 2021). The dgvAPP360 instead creates a 360°dark hole surrounding each of the two images of the target star, and then uses an additional grating to overlap the images to produce a single image of the star (Doelman et al. 2022). An example image obtained with the dgvAPP360 is shown in Figure 1, with the target star in the centre and the dark hole surrounding it.
In addition, we combine the dgvAPP360 with an IFS, enabling us to use differential spectrophotometry for high-contrast directly imaged companions. The incoming light is first dispersed into individual spectra, and then recombined into a single 'white-light' data point. This has the advantage of smoothing out wavelengthdependent flat-fielding errors and allows wavelength regions with instrumental absorption or highly variable telluric bands to be excluded, meaning systematic effects can be significantly reduced, thus yielding greater stability and precision in the final white-light curve.

The HD 1160 system
We test this approach using observations of the HD 1160 system, which is located at a distance of 120.4±0.6 pc (Gaia Collaboration et al. 2016Bailer-Jones et al. 2021) and consists of host star HD 1160 A (spectral type A0V; Houk & Swift 1999) and two high-contrast companions, HD 1160 B and C (Nielsen et al. 2012). HD 1160 B and C lie at separations of ∼80 au (∼0.78 ) and ∼530 au (∼5.1 ), respectively. Several key properties of HD 1160 A are listed in Table 1. HD 1160 A is bright (K = 7.040±0.029 mag, Cutri et al. 2003), and the contrast ratio between HD 1160 A and HD 1160 B is Δ = 6.35 ± 0.12 mag (Nielsen et al. 2012). This makes it an ideal target for demonstrating our technique, as a high signal-to-noise detection of the companion allows high cadence monitoring and a deep investigation into any residual systematic effects. A-type stars such as HD 1160 A generally vary below the millimagnitude level, which corresponds to variability amplitudes comfortably below the ∼1% level (Ciardi et al. 2011). We assess the variability of the host star in Section 4.2 using observations from the Transiting Exoplanet Survey Satellite (TESS) mission. The age of the system is poorly constrained, with estimates ranging from 50 +50 −40 Myr (Nielsen et al. 2012) to 100 +200 −70 Myr (Maire et al. 2016). The HD 1160 system may be a member of the Pisces-Eridanus stellar stream, which would place its age at ∼120 Myr, but this has yet to be confirmed (Curtis et al. 2019).
HD 1160 B has a spectral type close to the brown dwarf/stellar boundary; Nielsen et al. (2012) found a spectral type of ∼L0, but more recent papers suggest that it lies between M5-M7 (Maire et al. 2016;Garcia et al. 2017;Mesa et al. 2020). However, Mesa et al. (2020) found its spectrum to be highly peculiar and that no spectral model or template in current libraries can produce a satisfactory fit. Although the cause of this peculiarity has not yet been explained, Mesa et al. (2020) hypothesise that possible causes could include a young system age, dust in the photosphere of HD 1160 B, or ongoing evolutionary processes. The mass of HD 1160 B also remains unclear, primarily due to the poorly constrained age of the system, with estimates ranging from that of a low mass brown dwarf (∼20 M Jup , Mesa et al. 2020) to decisively in the stellar mass regime (0.12±0.01M ≈ 123M Jup , Curtis et al. 2019) if the system is indeed a member of the Pisces-Eridanus stellar stream. HD 1160 C is a low-mass star (spectral type M3.5), and its separation of ∼530 au (∼5.1 ) places it beyond the field of view of most vAPPs currently in use (Maire et al. 2016;Doelman et al. 2021).
The observations carried out on the HD 1160 system are References: (1) Houk & Swift (1999); (2) (Cutri et al. 2003;Skrutskie et al. 2006) described in Section 2, and in Section 3 we describe the spectral extraction and data reduction processes. In Section 4 we produce and present our differential spectrophotometric light curves of HD 1160 B. We then examine various factors that may be correlated with the light curves and detrend them in Section 5. These results and their implications are then discussed in Sections 6 and 7. Lastly, the conclusions of the work are summarised in Section 8.

OBSERVATIONS
We observed the HD 1160 system on the night of 2020 September 25 (03:27:31 -11:16:14 UT) with the double-grating 360°vector Apodizing Phase Plate (dgvAPP360) coronagraph (see Section 1.1) and the Arizona Lenslets for Exoplanet Spectroscopy (ALES) IFS (Skemer et al. 2015;Hinz et al. 2018;Stone et al. 2018). ALES is integrated inside the Large Binocular Telescope Interferometer (LBTI) (Hinz et al. 2016;Ertel et al. 2020), which works in conjunction with the LBT mid-infrared camera (LMIRcam), on the 2 x 8.4-m LBT in Arizona (Skrutskie et al. 2010;Leisenring et al. 2012). For these observations, ALES was in single-sided mode and was therefore fed only by the left-side aperture of LBT. Atmospheric turbulence was corrected for by the LBTI adaptive optics (AO) system (Bailey et al. 2014;Pinna et al. 2016Pinna et al. , 2021. We used the ALES L-band prism, which covers a simultaneous wavelength range of 2.8-4.2 µm with a spectral resolution of R∼40 . The plate scale is ∼35 mas spaxel −1 . The other LBT aperture was used to feed the Potsdam Echelle Polarimetric and Spectroscopic Instrument (PEPSI), which obtained R = 50,000 combined optical spectra of HD 1160 A and B in the 383-907nm wavelength range (Strassmeier et al. 2015(Strassmeier et al. , 2018, which is subject to analysis in other forthcoming works. Conditions were exceptionally clear and stable throughout the night, with no time lost to weather, and seeing ranged from 0.7-1.4 . We acquired 2210 on-target ALES frames, with an integration time of 5.4 s per frame, giving a total on-target integration time of 11934.0 s (∼3.32 h) spread over ∼7.81 h once readout time, nodding, and wavelength calibrations are included. The integration time was chosen such that the stellar PSF remained unsaturated in the core so that it can be used as the photometric reference for the companion. We used an on/off nodding pattern to enable background subtraction, nodding to a position 5 away for the off-source nod position. As HD 1160 C is located at a similar separation (∼5.1 ), we nodded in a direction away from this companion to prevent it from contaminating the frames obtained in the off-source nod position. Beam-switching in this way is possible because of the intrinsic stability provided by the dgvAPP360 coronagraph's placement in the pupil plane. We obtained dark frames with the same exposure time at the end of the night, and 6 wavelength calibrations were acquired at irregular intervals during the night. LBTI operates in pupil-stabilized mode, such that the field of view was rotating throughout the observations. The total field rotation across the observing sequence was 109.7°. The HD 1160 system was observed from an elevation of 29.4°at the start of the night to a maximum elevation of 61.7°, and then back down to an elevation of 27.5°. The dgvAPP360 creates an annular dark hole around the target PSF, with an inner radius close to the PSF core and an outer radius at the edge of the 2.2 field of view of the detector (2.7 -15 /D in ALES mode) (Doelman et al. , 2021. For these observations, this meant that HD 1160 B was located in the dark hole of HD 1160 A across the entire wavelength range covered by the ALES L-band prism. HD 1160 C (separation ∼5.1 ) remained beyond the ALES field of view.

DATA REDUCTION
The raw ALES images contain the spectra that have been projected onto the detector. Ultimately, we are aiming to produce a light curve for the companion that is made from the 'white-light', i.e. combined in the wavelength dimension. To do this, we must first extract the spectra from the raw data along with bad pixel correction and flatfielding. We can then extract the photometry of the star and the companion at each wavelength before collapsing the data in the wavelength dimension to obtain white-light fluxes. A light curve for the companion is then obtained by dividing the companion flux by the stellar flux to remove systematic trends shared by both. In the following subsections we describe the methods used to carry out each of these steps and obtain the white-light curve of the companion.

Spectral data cube extraction
Raw ALES data consists of a two-dimensional grid of 63×67 microspectra over a 2.2 x 2.2 field of view, which must be extracted into three-dimensional data cubes of -position, -position, and wavelength (Stone et al. 2022). To do this, we first performed a background subtraction using the sky frames obtained in the offsource nod position. For each ALES image, we subtracted the median combination of the 100 sky frames closest in time. We then extracted the micro-spectra into cubes using optimal extraction, which is an inverse variance and spatial profile weighted extraction approach (Horne 1986;Briesemeister et al. 2018;Stone et al. 2020). The extraction weights were obtained by measuring the spatial profile of each micro-spectrum in the dark-subtracted sky frames. As there is no significant change in the spatial profile as a function of wavelength, we average the spatial profile of each micro-spectrum over wavelength to obtain higher signal-to-noise (S/N). The wavelength calibration of the raw ALES data was then obtained using four narrow-band photometric filters at 2.9, 3.3, 3.5, and 3.9 µm, positioned upstream of the ALES optics. These filters are each of a higher spectral resolution Δ ∼100 than ALES, so are therefore unresolved and provide four single-wavelength fiducial spots with which each micro-spectra can be calibrated (Stone et al. , 2022. For each micro-spectrum, we performed this calibration by fitting a second-order polynomial to the calculated pixel positions of these four spots, therefore mapping pixel position to wavelength. Each of the 63×67 micro-spectra was thereby converted into a corresponding spaxel in the three-dimensional data cube (Briesemeister et al. 2019;Doelman et al. 2022). The resulting data cube contained 100 channels in the wavelength dimension ranging from 2.8-4.2 µm.
The primary wavelength calibration used to process this data set was obtained at 08:25:00 UT on the night of observations, i.e. 4 hours 57 minutes into the observing sequence. To test whether a wavelength calibration obtained at a different point in the night has a significant effect on the photometry of the target star and the companion, we also separately processed the data using an alternative wavelength calibration obtained at 04:54:00 UT, 1 hour 26 minutes into the observing sequence. We compare and discuss the results of the two wavelength calibrations in Section 4.3.

Data processing
Once the spectra were extracted into background-subtracted threedimensional cubes of images for each exposure, we applied several data reduction steps to remove systematics and improve the S/N at the location of the companion. We first removed 8 time frames from the data in which the AO loop opened while the data was being collected. Next, we identified bad pixels using a 6 filter and replaced them with the mean of the neighbouring pixels. We then applied a flat-field correction to calibrate the data against the response of the detector. For each wavelength channel, the corresponding flat was a time-average of the frames obtained in the 'off' nod position, which had then been corrected for bad pixels in the same way as the science frames and smoothed over using a Gaussian filter. These flats were then divided by the maximum value in the frame such that the value of every pixel was between zero and one. We then divided the science frames by these smoothed, normalised sky flats. This flat-fielding process was also repeated separately using a median filter instead of a Gaussian filter as a means to test the robustness of this step in our method. We proceed with the Gaussian filter and discuss the impact of the choice of flat frame on the photometry of the star and companion in Section 4.3. Doelman et al. (2022) previously identified that backgroundsubtracted, flat-fielded ALES images contain residual structure that cannot be described by purely Gaussian noise, in the form of timevarying row and column discontinuities (faintly visible in the top panel of Figure 2). Such discontinuities are expected and arise from the way in which the micro-spectra lie across multiple channels of the LMIRcam detector (Doelman et al. 2022). We followed the method of Doelman et al. (2022) to characterise and remove these discontinuities by fitting a third-order polynomial to each row and column in each frame ( Figure 2). Removing these systematics is important as they could impact the precision of our differential light curve, or even generate a false variability signal if the target moves over them throughout the observing sequence. Prior to fitting, we applied circular masks at the locations of the star and the companion in each frame such that their flux did not contaminate the fit. To find the position of the star in each time frame, we selected a wavelength channel with a high stellar flux per frame (channel 52, ≈ 3.69 µm) and fit the PSF core with a 2D Gaussian. The position of the companion in each time frame was then identified using its separation and position angle relative to the star and accounting for the effect of the field rotation over time. We then masked the star and the companion across all wavelength channels using circular masks with diameters of 18 pixels and 5 pixels, respectively, before fitting the third-order polynomials to each column. The resulting values were then subtracted from the data to remove the column discontinuities. This process was then repeated for each row in the resulting image to remove the row discontinuities.
In addition to removing these systematic discontinuities, this process has the effect of removing residual background flux not eliminated by earlier processing steps. This is indicated by the histograms in Figure 2, which show that the noise distribution of the data was offset from zero prior to the removal of the discontinuities (in blue) but is approximately consistent with zero after this process has been applied (in orange). We then used the position of the star in each frame, found when applying the masks in the previous step, to spatially align the data such that the star was in the centre of each frame. Finally, we rotationally aligned the images by applying an anticlockwise rotation corresponding to their parallactic angles.
We did not use any further post-processing methods that reduce quasistatic speckle noise through the subtraction of reference PSFs, such as Angular Differential Imaging (ADI, Marois et al. 2006a). Although the field rotation over the night of observations was sufficient enough to use these to remove noise with minimal companion self-subtraction, doing so would also remove the unsaturated stellar reference PSF, which is required to eliminate systematics in the companion photometry. It would not be possible to use the host star PSF prior to ADI subtraction as a photometric reference for the companion PSF after ADI subtraction, as the two would no longer share the same systematic trends. Furthermore, HD 1160 B is sufficiently bright that it can be detected at ample S/N for our purposes without further noise reduction.

Wavelength channel selection
Although the final processed cubes contain data from 100 wavelength channels across the observed wavelength range of 2.8-4.2 µm, not all of these channels are suitable for further analysis. The first 3 and final 10 channels contain flux from the adjacent spaxel in the dispersion direction as an oversized spectral length is required to extract the spectrum at each position, so the extracted data overlaps slightly in the wavelength dimension (Stone et al. 2022). Furthermore, the dgvAPP360 contains a glue layer that causes up to 100% absorption between 3.15 and 3.55 µm (Otten et al. 2017;Doelman et al. 2021). As described in Section 1.1, one of the key advantages of a spectrophotometric approach is the option to exclude channels that are known to cause systematic variability in the 'white-light' curve, hence improving the companion S/N and stability in the combined image when compared to combining all wavelength channels without any selection. This is key to reducing large systematic effects that may otherwise dominate the variability signals that we are aiming to measure. For the purposes of this technique demonstration, we proceed using 30 sequential wavelength channels (45-74, spanning 3.59-3.99 µm) which all have a high throughput and do not lie in the regions affected by the dgvAPP360 glue absorption, significant telluric absorption, or the overlapping spectral traces. In the left panel of Figure 3 we show the median combination of these channels in both wavelength and time, while the centre and right panels respectively show example images from the median com- Demonstration of the process for removing systematic row and column discontinuities present in background-subtracted ALES images, as applied to a single frame of data in the 52nd ALES wavelength channel ( ≈ 3.69 µm). Top panels: input frame prior to the removal of the discontinuities, which are faintly visible as a chequered pattern. All three panels are the same. Both HD 1160 A and B are masked. Second row: results of the third-order polynomial fits individually to the columns (left panel) and rows (centre panel), and the combination of both (right panel). The combination of both was produced by first fitting and removing the column discontinuities, and then repeating the process on the resulting image to fit and remove the rows. We show row and column fits to the input frame separately here to highlight their individual contribution to the original systematics. The star and the companion were masked during this process, and the values to be removed at their locations were found through interpolation of the fits. Third row: data frame with the discontinuities removed by subtracting the fits. The histograms in the bottom row show the distributions of the counts in the unmasked regions of the third row images, with the original noise distribution in blue and the noise distribution after the discontinuities were removed in orange. The bottom right panel shows the version used in the analysis.
bined cubes in the time and wavelength dimensions only. The images shown are those processed using the flat frame that was smoothed using a Gaussian filter; the equivalent images as processed using the median-smoothed flat frame are visually indistinguishable from these.

GENERATING DIFFERENTIAL SPECTROPHOTOMETRIC LIGHT CURVES
Variability arising from instrumental systematics and the effects of Earth's atmosphere, such as airmass, seeing, and tellurics, con-taminate the raw flux of the companion. Simultaneous flux measurements of a photometric reference are required to eliminate this contaminant variability and produce a differential light curve of the companion, relative to the photometric reference. Although suitable photometric references are generally absent when using coronagraphs (Section 1), the dgvAPP360 coronagraph uniquely provides an image of the host star simultaneously to the companion, allowing the star to be used as the photometric reference when it is not saturated. Its placement in the pupil plane also makes it inherently stable and insensitive to tip/tilt instabilities (Otten et al. 2017;Doelman et al. 2022).

RA [arcsec]
Figure 3. Reduced images of the HD 1160 system obtained with LBT/ALES and the dgvAPP360 coronagraph after all data processing steps and wavelength channel selection. Left: the median combination, in both wavelength and time, of all wavelength channels in the 3.59-3.99 µm range, covering a total integration time of 11891 s (∼3.31 h). Centre: example image from the same data cube but median combined only in the time dimension, along the 3.64 µm wavelength channel. Right: example time frame (integration time = 5.4 seconds) resulting from a median combination in wavelength over the 3.59-3.99 µm wavelength range. All three images use the same arbitrary logarithmic colour scale, and are aligned to north. North is up, and east is to the left.

Aperture photometry
We used version 1.4.0 of the Photutils Python package (Bradley et al. 2022) to simultaneously extract aperture photometry of HD 1160 A and B. We carried out this process for every individual frame in each of the 30 wavelength channels in the 3.59-3.99 µm range chosen in Section 3.3, with the aim of then combining these in the wavelength dimension to produce the white-light flux for each object. Circular apertures with radii of 9 pixels (3.1 /D) and 2.5 pixels (0.9 /D) were used for the host star and the companion, respectively. To estimate the background flux at the position of the star, we also extracted photometry in a circular annulus centred on the stellar location with inner and outer radii of 11 and 16 pixels, respectively. It was not possible to use this method to estimate the background flux at the position of the companion as the companion lies close to the edge of the field of view in some frames, limiting the space available to place an annulus that would be statistically wide enough. We therefore instead followed an approach used by Biller et al. (2021), estimating the background at the companion location by masking the companion and extracting flux in a circular annulus centred on the host star, with a width of 6 pixels at the radial separation of the companion. As most of the residual background in each frame was eliminated by the data processing steps in Section 3.2, these background values are close to zero. These apertures and annuli are shown in Figure 4, overlaid on a single processed time frame of data in the 52nd ALES wavelength channel ( ≈ 3.69 µm). We then removed the residual background from our stellar and companion flux measurements by multiplying the mean flux per pixel in the background annuli by the area of the corresponding apertures and subtracting the resulting values from the aperture photometry. We then produced single white-light measurements for both the companion and the star at each time frame by taking the median combination of the photometric measurements across the 30 wavelength channels. These raw time series, uncorrected for shared variations introduced by Earth's atmosphere (i.e. before division), are shown in grey in the top two panels of Figure 5. The discrete gaps in integration reflect time spent off-target due to the two-point on/off nodding pattern used to enable background subtraction. We also plot the data binned to 18 minutes of integration time. We binned the data by taking the median value in each time bin. The error on the binned fluxes are the Gaussian approximation of the root mean square (RMS) i.e. median absolute deviation (MAD) × 1.48 of the flux measurements inside each bin divided by √ − 1, where is the number of frames per bin. Next, we removed variability due to Earth's atmosphere and other systematics from the unbinned raw flux of the companion using the unbinned raw flux of the host star, which acts as a simultaneous photometric reference. By dividing the unbinned companion flux by the unbinned stellar flux, we eliminate trends common to both and produce a differential light curve that only contains non-shared variations. Assuming that the host star is not itself varying (see Section 4.2), the resulting differential light curve reflects the intrinsic variability of the companion plus any contamination arising from non-shared systematics. We show this raw differential light curve in the third panel of Figure 5. The bottom panel shows a zoomed-in view of the binned version, with tighter limits on the y-axis.
In the following sections we examine a number of physical, instrumental, and processing factors that may be correlated with non-astrophysical features in the differential white-light curve, and in Section 5 we model and remove non-shared variations arising from some of these factors.

TESS light curves of host star HD 1160 A
Although the vast majority of A-type stars generally vary well below the ∼1% level, a small fraction vary at a much higher level, showing up to ∼15% variations (Ciardi et al. 2011). To test our assumption that the host star HD 1160 A is not varying at a level that impacts our differential light curve, we used data from the TESS mission, which is publicly available on the Mikulski Archive for Space Telescopes . The image is a single time frame from the 52nd ALES wavelength channel ( ≈ 3.69 µm). North is up, and east is left. The orange aperture is placed at the location of HD 1160 B, which is too faint to be visible in a single frame. The companion was masked when extracting photometry in the annulus for the companion background.
October 11), 51 days in total, with 2 minute cadence. The TESS detector bandpass covers a broad-band wavelength range of 0.6-1.0 µm (Ricker et al. 2015), which does not overlap with our LBT/ALES observations in the 2.8-4.2 µm range. However, as stars are generally less variable in the infrared than in the optical regime, any variations in the TESS light curve of HD 1160 A should represent an upper limit for its variability at the wavelengths covered by ALES (e.g. Solanki & Unruh 1998;Unruh et al. 1999;Fröhlich & Lean 2004;Davenport et al. 2012;Goulding et al. 2012;Ermolli et al. 2013;Rackham et al. 2022). Each TESS pixel covers 21 on sky. This means that HD 1160 A, HD 1160 B (at a separation of ∼0.78 ), and HD 1160 C (at a separation of ∼5.1 ) are not resolved separately in the TESS images and appear as a single object. However, as both HD 1160 B (Δ = 8.85 ± 0.10 mag, Δ = 6.35 ± 0.12 mag) and C (Δ = 6.33 ± 0.04 mag, Δ = 4.803 ± 0.005 mag) are far fainter than HD 1160 A, especially at shorter wavelengths, it is reasonable to assume that flux of HD 1160 A will dominate in the TESS data (Nielsen et al. 2012).
We first masked out any bad quality exposures using the onehot encoded quality mask in the 'QUALITY' keyword in the header of the light curve files provided by the TESS Science Processing Operations Center (SPOC, Jenkins et al. 2016) on MAST. We then used the 'CROWDSAP' keyword in the header to get an estimate of the ratio of target flux to total flux in the optimal aperture used for the PDC SAP (Pre-search Data Conditioning Simple Aperture Photometry) flux (e.g. Panwar et al. 2022b). The 'CROWDSAP' value for sector 42 and 43 indicates that 0.17% and 0.2% flux is from dilution by nearby sources. We subtract the estimated diluted flux from each exposure in both the sectors. The resultant light curves for both the sectors are shown in Figure 6. Although the TESS observations are not contemporaneous with our LBT/ALES observations, we do not see variations above 0.03% in the light curve of HD 1160 A over the timescale covered by the two TESS sectors (51 days). As this is far smaller than the precision of our differential light curve, we proceed with the assumption the host star HD 1160 A is non-varying within the flux precision of our analysis of the variations in the light curve.

Impact of wavelength calibration and flat-field smoothing
In Section 3.1, we described the process used to perform the wavelength calibration of the raw data and to extract the micro-spectra into a three-dimensional image cube. This step was repeated separately using the wavelength calibration that was the most divergent of the 6 obtained throughout the observing sequence, i.e. the 3.9 µm fiducial spots for this wavelength calibration were the most significantly offset compared to the one that was originally used. Over the course of the night, the projection of the micro-spectra onto the detector drifts slightly. If this drift is significant then a particular wavelength calibration may not remain accurate for the entire observing sequence, potentially producing a false variability signal when the wrong part of the spectrum is assigned to a given channel. Repeating our spectral extraction using a wavelength calibration obtained at a different point during the observations allows us to test whether this effect has a significant impact on the photometry of the target star and the companion. After extracting the micro-spectra using the alternative wavelength calibration, we then processed the data again in full to produce an alternative differential light curve. The original and alternative wavelength calibrations were obtained at 4 hours 57 minutes (08:25:00 UT) and 1 hour 26 minutes (04:54:00 UT) after the beginning of the observing sequence, respectively. We plot the resulting alternative stellar and companion fluxes, and the differential white-light curve (binned to 18 minutes) in Figure 7, alongside the originals from Figure 5 for comparison. The differential light curves in each case are consistent within 1 , indicating that the extracted photometry is sufficiently robust to changes in the wavelength calibration and sub-pixel mis-registration of the spatial profiles for each microspectrum.
In Section 3.2, we described our method for applying a flatfield correction to the data to calibrate for the non-uniform response of the detector. Incorrect flat-fielding can lead to a false variability signal if the companion moves over regions of the detector with a non-uniform response that has not been properly calibrated. To test the robustness of our differential light curve to differences in the flat used we processed the data in two separate streams, using a flat that had been smoothed over using a Gaussian filter and a median filter, respectively. We plot the resulting differential light curves in Figure 8 for comparison. The two differential light curves are in close agreement and every binned data point lies well within their 1 error bars, indicating that the method for producing the flat is robust and does not significantly affect the final images or extracted photometry of the star or companion.

DETRENDING THROUGH LINEAR REGRESSION
Trends shared by the star and companion fluxes are removed in the differential light curve (see bottom panel of Figure 5), but any non-shared trends will still be present, including the intrinsic companion variability signal that we aim to measure. To improve our sensitivity to the companion's variability, we needed to remove the non-astrophysical residual trends in the differential light curve, which can arise from both telluric and instrumental sources. In  . Unlike traditional transmission spectroscopy observations, our target is significantly fainter than the simultaneous reference that we use for detrending. Furthermore, the target was not pixel-stabilised for these observations and moved across the detector throughout the night, so we might predict that the measured light curves could be significantly correlated with the change in position of the companion and the star on the detector over time. Knowing how to remove these systematics is key to obtaining high-precision light curves in future observations of directly imaged exoplanets. For space-based observations the instrumentation is generally sufficiently stable such that the systematics are repeatable over time, allowing them to be well characterised. This is more challenging for ground-based observations, like those here, which are inherently less stable as Earth's atmosphere introduces systematics that can vary night by night. As a basic demonstration of how to remove such residual trends from ground-based differential light curves of directly imaged planets, we here used a multiple linear regression approach to simultaneously fit several possible sources of systematics. This is not intended as a strictly rigorous statistical analysis of the trends in the light curve, but is done to perform an initial investigation into which parameters have the greatest impact and to illustrate an example approach of how to do this for future observations. As studies move towards increasing precision to measure smaller amplitude variability, one might instead consider approaches using Gaussian processes or similar.
An investigation of the LBT telemetry and white light images revealed eight physical and instrumental factors, shown plotted against time in Figure 9, that varied notably during the observing sequence and may be correlated with the residual trends in the dif-ferential light curve. We therefore included these parameters in the linear regression. The first three of these were air temperature, wind speed, and wind direction, shown in the three panels on the left-hand side of Figure 9. We also considered airmass, which is shown in the top-right panel. While the light from the companion and its host star pass through almost identical airmass (maximum difference ∼10 −5 ), their significantly different colours mean that atmospheric extinction due to absorption and e.g. Rayleigh scattering can result in a differing airmass dependence, even when such scattering effects are reduced at our longer 2.8-4.2 µm wavelength range (Allen 1955;Broeg et al. 2005;Croll et al. 2015;Panwar et al. 2022a).
The remaining four parameters included in the linear regression were the x-and y-pixel positions of the star and the companion in the images prior to spatial and rotational alignment, shown in the centre-right and bottom-right panels of Figure 9. The dgvAPP360 is located in the pupil plane, meaning that drifts in the locations of the target PSFs on the detector do not affect its response and performance as it applies the phase modification to every source in the field of view (Otten et al. 2017;Doelman et al. 2022). However, systematics could still be introduced by such drifts if there are variations in the instrumentation or detector response. A number of sharp discontinuities in the y-position, and a singular discontinuity in the x-direction, can be seen in Figure 9. These discontinuities are the result of manual positional offsets applied during the observations to keep the star close to the centre of the small (∼2.2 x 2.2 ) ALES field of view. These offsets were always applied while in the off-source nod position, and always along one axis at a time. The largest discontinuities in the x-and y-positions correspond to shifts of 0.1 along the given axis. Neglecting the discontinuities, the drift of the stellar PSF in both the x-and y-directions follow arcs with turn-overs approximately 4.5 hours into the observing sequence. This slow drift is correlated with the pointing altitude of the telescope and arises from flexure of the ALES lenslet array as the telescope rotates. As the observations were obtained in pupilstabilized mode such that the field of view was rotating over time, the change in position of the companion throughout the data has an additional rotational component compared to the star. The drift arising from the flexure of the lenslet array therefore instead produces an inflection point ∼4.5 hours in the case of the companion.
We used the linear regression tools in version 1.0.2 of the scikit-learn Python package (Pedregosa et al. 2011) to simultaneously fit these input parameters and produce a model fit to our differential light curve. The coefficients of the linear model produced by the linear regression are shown in Table 2, and the model itself is shown in green in the top panel of Figure 10, relative to the raw differential white-light curve in grey. We then divided the raw differential light curve by this model to produce a detrended version, shown in red in the bottom panel of Figure 10. We also overplot the raw differential light curve from the bottom panel of Figure 5, prior to detrending, in purple to allow the two to be compared. We also repeated this detrending process to produce detrended differential light curves for each of the 30 individual wavelength channels over the 3.59-3.99 µm wavelength range that comprise the white-light curve. These are shown in Figure 11, again binned to 18 minutes. We note that while ALES has a resolution of R∼40, the raw data were spectrally extracted into 100 wavelength channels and so there is some correlation between wavelength channels.

RESULTS
The detrended differential white-light curve ( Figure 10) shows sinusoidal-like variability over a timescale of a few hours. In the binned light curve, the normalised flux ranges from a minimum of 0.91 at 0.874 hours to a maximum of 1.13 at 2.961 hours. To better allow us to estimate the differential precision that we achieve in our light curve, we fitted the variability and removed it from our light curve. As the variability signal appears periodic, we used the NASA Exoplanet Archive periodogram service 2 to apply the Lomb-Scargle algorithm to the unbinned detrended differential white-light curve and thereby search for sinusoidal periodic signals (Lomb 1976;Scargle 1982). The strongest peak in the resulting periodogram (top panel, Figure 12) is at a period of 3.242 hours. We then fit a sinusoid to the light curve using this period as an initial guess, which returned a function with the same period (3.239), a semi-amplitude of 0.088, a phase shift of 0.228, and a -offset of 0.993. This sinusoid is shown overplotted on the detrended light curve in the centre panel of Figure 12. The differential light curve was then divided by the fitted sinusoid to remove the variability signal to the first order (centre panel residuals, Figure 12). The bottom panel of this figure is the same as the panel above but with the data phase-folded to the period of the fitted sinusoid. Next, we followed the method of Kipping & Bakos (2011) to assess the degree of "red" (correlated) noise in our light curve. We binned our detrended differential white-light curve, Normalised Flux Figure 8. The raw differential companion/star white-light curves that are produced when the flat-field correction uses a flat that was smoothed using a median filter and a Gaussian filter, in turquoise and purple, respectively. The latter light curve is the same as that in the bottom panel of Figure 5, reproduced here for comparison. with the sinusoid removed, to a range of bin sizes, before normalising and subtracting one to centre around zero. We then measured the RMS of each resulting binned light curve. These RMS values are plotted against bin size in Figure 13, alongside the expectation of independent random numbers as a function of bin size i.e. the white noise. For our chosen binning of 18 minutes of integration time per bin, which has a bin size of 200 frames per bin, we find an RMS of 0.037. For comparison, the RMS of the detrended differential white-light curve prior to the removal of the sinusoid, but with the same binning, is 0.073. We therefore conclude that the light curve of HD 1160 B shows variations with a semi-amplitude of ∼8.8% or peak-to-peak amplitude of ∼17.6%, and that the differential precision achieved in the binned light curve is at the 3.7% level. The amplitude of the variations is therefore above the measured precision. Furthermore, this estimate of the precision is likely a conservative one; the variability signal is unlikely to have been perfectly removed by the sine fit and so the measured RMS values may be higher than the true limiting precision. A caveat of this result is that the baseline of our observations is only ∼7.81 hours, so we only cover ∼2.4 periods. Additional data is therefore needed to confirm the periodicity and amplitude of the variability of HD 1160 B.
We discuss these results further and compare the precision achieved to similar studies in the literature in Section 7.2.

DISCUSSION
In this section we discuss the relative impact of the decorrelation parameters on our results, and the physical explanations of the systematics they introduce. We also compare the precision that we achieve to other variability studies in the literature that use different techniques, and discuss the potential application of differential spectrophotometry in future work.

Impact of decorrelation parameters
Using the linear regression coefficients of each parameter from Table 2, we can assess which parameters have the greatest impact on the light curve of HD 1160 B. The small angular separation of the companion and star might suggest that the effect of airmass would be small. Airmass can be an important systematic for differential spectrophotometric observations where other stars are used as the simultaneous photometric reference as the angular separation between the target and the reference can be large, causing the light from the two objects to pass through different atmospheric columns (e.g. Broeg et al. 2005;Panwar et al. 2022a). However, in our case we use the companion's host star as the photometric reference, so the angular separation between the two is much smaller than is generally the case for observations using reference stars in the field. However, there is a significant colour difference between the star and the companion, which leads to different degrees of extinction at a given airmass. Therefore, we expect an airmass dependence and indeed it has the largest coefficient in the detrending model. Similar extinction effects are often seen in studies of transiting exoplanets, and can also exhibit a non-linear wavelength dependence (e.g. Panwar et al. 2022a). We also find the air temperature at LBT to be one of the parameters that is most correlated with our differential light curve. As the air temperature changes, this can potentially cause slight changes in the optical path of the telescope or instrument that lead to this correlation.
We further predicted that the positions of the companion and the star could introduce significant non-shared systematics to our measurements; as HD 1160 A was not pixel-stabilised during our observations, we are sensitive to both intra-and interpixel variations as the star and the companion move across the detector. Both the star and the companion changed position on the detector due to flexure of the ALES lenslet array as the instrument moved and positional offsets applied intentionally to keep the target close to the centre of the field of view. We also nodded on and off of the target to enable background subtraction. While the process of nodding is itself relatively accurate, it is not repeatable at a sub-pixel pointing precision, introducing a slight offset error between nods. Furthermore, LBTI data is always pupil-stabilized, so the field of view was rotating throughout the night. Although it may have been optimal to fix the star and companion positions to the same detector pixels for the duration of the observations, this was not possible due to the effect of the lenslet array flexure and the lack of instrument derotator in the LBTI architecture (Doelman et al. 2022). However, we find that these positional changes are not the most correlated with our differential light curve. It is possible that this is because the light from the star and companion is spread out across multiple detector pixels when spectrally dispersed, reducing the impact of . The decorrelation parameters used when detrending the differential light curve via linear regression. From top to bottom, the left panels show the air temperature (in ℃), wind speed (in m s −1 ), and wind direction (in degrees east of north) at the observatory as a function of time, as extracted from the FITS headers of the raw data. The panels on the right show airmass, the x-and y-positions of the star in pixels, and the x-and y-positions of the companion in pixels, as a function of time. The gaps in time reflect the on/off nodding pattern of the observations. The companion and star positions are those in the images prior to spatial and rotational alignment, and the sharp discontinuities in pixel position are due to manual offsets applied to keep the star close to the centre of the small field of view. The stellar positions follow arc-shaped trends aside from these discontinuities, which correlate with the pointing altitude of the telescope and arise from flexure of the ALES lenslet array as the telescope rotates throughout the night. The change in the companion position has an additional trend due to the 109.7°rotation of the field of view over the observing sequence. Normalised Flux Figure 10. The model produced by the linear regression using the decorrelation parameter coefficients (in Table 2) is shown in the top panel in green alongside the raw differential light curve in grey. The bottom panel then shows in red the light curve produced when the raw differential light curve is divided by the linear regression model to remove the modelled trends that are not shared by the stellar and companion fluxes, binned to 18 minutes of integration time per bin. The raw differential light curve (i.e. prior to detrending) is also shown faintly in purple, reproduced for comparison from the bottom panel of Figure 5.
any systematic issues arising from any single pixel. The light of our target is dispersed across wavelength, similar to a technique often used for high-precision differential photometry whereby a telescope is intentionally defocused to disperse starlight over the detector (e.g. de Mooĳ et al. 2011Mooĳ et al. , 2013Crossfield et al. 2012;Croll et al. 2015). This has the effect of reducing systematics due to intrapixel variations and minimises the impact of any residual flat-field errors. The step of recombining our data into a white-light curve may therefore have helped to reduce systematic trends that would otherwise have been introduced by the positional movement of the targets throughout the night. The use of a spectrograph also has the additional benefit of allowing us to remove individual wavelength channels that are found to contain defects. For this dataset, this allowed us to leave out channels affected by overlapping spectral traces, significant telluric absorption, and absorption by the glue layer of the dgvAPP360, as described in Section 3.3.
Lastly, we find that the wind speed and direction are the least correlated with our differential light curve, but note that the conditions on the night were exceptionally stable and so cannot rule out that these factors could have an impact in less optimal conditions. For future observations, residual systematics in differential light curves of directly imaged companions could be removed using more advanced methods from the exoplanet transmission spectroscopy and high-precision secondary eclipse literature such as fitting the data using a Gaussian process regression (e.g. Gibson et al. 2012;Evans et al. 2017;Nikolov et al. 2018a,b;Diamond-Lowe et al. 2020a,b;Wilson et al. 2021). In general, methods used to identify trends in transmission spectroscopy data also include a model of the exoplanet transit itself, allowing the transit to be detected even when the strength of the signal is very low. However, this is possible because the expected shape of the signal is well understood. This is more difficult for searches for variability in directly imaged companions, where the expected shape of the variability signal is not necessarily well-known in advance. Furthermore, many of these methods assume a linear relationship between systematics and trends in the light curve, while the telluric and instrumental systematics present in time-series data can be complex and non-linear. In the future, an optimal approach should account for the functional form of the correlated parameters (Panwar et al. 2022a,b). Normalised Flux Figure 11. In red, we show the detrended versions of the individual differential light curves in each of the 30 wavelength channels that were combined to produce the white-light curve. These wavelength channels cover =3.59-3.99 µm, and are binned to 18 minutes of integration time per bin. An offset factor of 2 has been applied between each light curve to separate them from each other. Overall variability appears to increase with longer wavelength.

Variability interpretation
In Section 6 we found that HD 1160 B shows variations with a semi-amplitude of ∼8.8% (or peak-to-peak amplitude of ∼17.6%).
To compare this result to literature observations of similar objects, we must consider both spectral type and the wavelengths at which variability has been observed. Vos et al. (2022) found that virtually all L dwarfs are likely to be variable at the 0.05-3% range, and several studies have measured higher variability, up to 26% (e.g. Radigan et al. 2012;Radigan 2014;Lew et al. 2016;Biller et al. 2018;Bowler et al. 2020b). However, there is evidence that brown dwarf variability amplitude may have a strong wavelength dependence. For example, HST observations of highly variable L dwarf companion VHS 1256-1257 b identified a large variability amplitude of 24.7% at 1.27 µm, while Spitzer observations at 4.5 µm found a far lower amplitude of 5.76±0.04% Zhou et al. 2020b;Miles et al. 2022). Zhou et al. (2020b) do note, however, that the HST and Spitzer observations were not obtained contemporaneously and so the atmosphere of VHS 1256-1257 b and hence its variability properties are likely to have changed substantially over the intervening timescale. Comparisons of large surveys also suggest that variability amplitudes are lower in the mid-infrared than in the near-infrared, although there is evidence for a weaker wavelength dependence and enhanced mid-infrared variability amplitudes for the young isolated brown dwarfs most similar to substellar companions (e.g. Radigan et al. 2014;Metchev et al. 2015;Biller et al. 2018;Vos et al. 2022).
The amplitude of the L-band variability we measure for HD 1160 B is quite extreme compared to literature results, but it is important to note that there have been very few variability studies for substellar companions in the mid-infrared, making direct comparison difficult. High-amplitude variability in brown dwarfs is generally attributed to heterogeneous surface features, such as spots or clouds of varying thickness, rotating in and out of view as the object rotates (e.g. Apai et al. 2013;Biller 2017;Artigau 2018). Some light curves show more complex features that cannot be modelled with a single atmospheric feature, or features that evolve over short or long timescales (e.g. Artigau et al. 2009;Metchev et al. 2015). These features and time evolution may arise from changing weather systems, or bands of clouds which rotate within the target's atmosphere and generate waves on a global scale (e.g. Apai  a longer baseline are required to be able to characterise any time evolution in the variability signal. However, the spectral type of HD 1160 B is unclear as it has a highly peculiar spectrum that cannot be satisfactorily fit with spectral models or templates in current libraries, and some studies suggest that it could instead be a late-M dwarf (Maire et al. 2016;Garcia et al. 2017;Mesa et al. 2020). If HD 1160 B is an M dwarf, its variability would most likely arise from cool star spots caused by magnetic activity in its photosphere, which are common in M- dwarfs and would rotate in and out of view in much the same way as the cloud features of lower mass objects (e.g. Barnes et al. 2002;Frasca et al. 2009;Scholz et al. 2009;Goulding et al. 2012;Ghosh et al. 2021;Johnson et al. 2021). The properties of such variability can be highly dependent on the spot distribution and fractional spot coverage of a given object; some M-dwarfs have a very high coverage with multiple starspots covering as much as 20-50% of their fractional surface area inhomogeneously (e.g. O'Neal et al. 2004;Morales et al. 2010;Irwin et al. 2011;Goulding et al. 2012;Jackson & Jeffries 2013). Spot-induced variability amplitudes for M-dwarfs generally range from the subpercent level up to around ∼5% (e.g. Rockenfeller et al. 2006;Charbonneau et al. 2009;Birkby et al. 2012;de Mooĳ et al. 2012;Nefs et al. 2013). Although observed far less often in the infrared compared to the optical, flaring events can induce far stronger variability in M-dwarfs at amplitudes ranging from the subpercent level up to tens of percent (e.g. Goulding et al. 2012;Tofflemire et al. 2012). Our measured variability amplitude for HD 1160 B is therefore also on the higher end of what has been observed for earlier spectral types such as late M-and early Ldwarfs, barring flares, although we again note the lack of literature studies of similar objects in the L-band. While the variability observed for HD 1160 B appears high, another point to consider is its orbital inclination, which the latest orbital fits suggest is close to edge-on as viewed from Earth (92.0 +8.7 −9.3°; Bowler et al. 2020a). If the obliquity of HD 1160 B is aligned with its orbit such that we are viewing its rotation close to edge-on, the observed variability amplitude is likely to compose a much larger fraction of its true variability compared to if it were viewed face-on. Indeed, Vos et al. (2017) demonstrated that the highest variability amplitudes are seen for targets with close to edge-on viewing angles.
If we interpret the 3.24 h sinusoidal variation we observed as the true rotation period of HD 1160 B, we can further consider this within physical limitations. The breakup period of a rotating object is dependent on its radius, which is itself age-dependent, and mass. Both age and mass are poorly constrained for HD 1160 B: literature results place the system's age in the 10-300 Myr range (Nielsen et al. 2012;Maire et al. 2016;Curtis et al. 2019), and mass estimates range from ∼20 M Jup to 123 M Jup (Curtis et al. 2019;Mesa et al. 2020). Vos et al. (2020) calculated the breakup periods of brown dwarfs as a function of age by equating equatorial velocity with the escape velocity, accounting for radial contraction over time. When we compare our measured HD 1160 B variability period of 3.24 hours to their results (their Figure 13), we find that this is a physically feasible rotation period for most possible combinations of mass and age from the literature, albeit very close to the breakup period in many cases. An alternative explanation is that the 3.24 hour variability signal that we see is produced by multiple features in the atmosphere of HD 1160 B, and that its rotation period is actually longer (Leggett et al. 2016). Additional observations of this variability over a longer baseline will help to further characterise its origin and confirm whether its periodicity reflects that of the companion's rotation.
The detrended differential light curves for each of the 30 individual wavelength channels over the 3.59-3.99 µm wavelength range that comprise the white-light curve (Figure 11) show increasing statistical errors at longer wavelengths, as expected as our S/N is lower here. Using the RMS of each light curve as a metric to compare the scatter in each channel, we do see tentative evidence of increasing variability towards longer wavelengths beyond the increase of RMS expected from the S/N. However, as the total baseline of our observations is only a single night, we see too few repetitions to be confident of variability patterns in individual wavelength channels. Additional spectrophotometric data will therefore be required to confirm this. Although we modelled the overall variability in our white-light curve with a single sinusoid, it is also possible that the phase and amplitude is different per wavelength channel as distinct atmospheric features at separate locations in the atmosphere of the companion may produce variability with a different wavelength dependence. Furthermore, although the overall wavelength range of these 30 channels is relatively small, different wavelengths probe different pressure levels in the companion's atmosphere, and hence different layers of the atmosphere (e.g. Buenzli et al. 2012;Biller et al. 2013;Apai et al. 2017;Ge et al. 2019).

Light curve precision
In Section 6, we found that after removing a single sinusoid from the data, we achieved a precision of 3.7% in the detrended differential light curve when it is binned to a bin size of 200 data points, corresponding to 11 bins of 18 minutes of integration time. There have been three previous studies searching for variability in substellar companion from the ground; Apai et al. (2016, and Wang et al. (2022) each conducted variability searches on the HR 8799 planets using satellite spots as photometric references. In a pilot variability study, Apai et al. (2016) reach a ∼10% planet-to-planet photometric accuracy for SPHERE observations of 25 minute cadence when data from different nights are combined for a total telescope time of 3.5 hours. Biller et al. (2021) goes further with SPHERE to conduct a longer (>4 hours) search, successfully constraining the sensitivity to variability to amplitudes >5% for HR 8799b and >25% for HR 8799c. More recently, Wang et al. (2022) used SCExAO/CHARIS to improve the variability constraints of HR 8799c to the 10% level, and HR 8799d to the 30% level. They did this by combining the use of satellite spots with a spectrophotometric approach similar to the one we present in this paper, using the CHARIS IFS to disperse the light into individual spectra before recombining the channels into wider bands. At first glance, the sensitivity that we achieve for HD 1160 B here may appear to compare favourably with these results. However, it is important to consider a number of caveats that make direct comparison unjustified. All three of the HR 8799 studies were conducted at near-infrared wavelengths with 8.2-m telescopes, while ours was in the mid-infrared with an 8.4-m telescope. More significantly, the HR 8799 planets are fainter than HD 1160 B, with contrasts of Δ = 8 − 10 mag compared to their host star, which has a H-band magnitude of 5.28 mag (Marois et al. 2008(Marois et al. , 2010. HD 1160 B is brighter, with a contrast of Δ = 6.35 mag compared to a host star with an -band magnitude of 7.06 mag (Nielsen et al. 2012). The lower sensitivity to variability achieved by these studies is therefore partially a reflection of the intrinsically lower fluxes of their targets, which leads to higher errors on their photometry.
However, each of the HR 8799 variability studies also found that the satellite spots can demonstrate individual variations of their own and are often anti-correlated with each other. This means that they may not always serve as appropriate photometric references with which to detrend the light curve of a companion. Wang et al. (2022) found that the flux ratio of the SCExAO satellite spots shows time variation with a scatter of ∼3% across a night, and can show even larger variations on a shorter timescale, up to 10%. This potentially sets a limit to the precision that can be achieved using satellite spots, particularly on nights where observing conditions are less stable.
A key advantage of the dgvAPP360 compared to satellite spots is its simplicity; the photometric reference it provides is simply an image of the host star, and so it does not suffer from the same correlated systematics as the satellite spots. Differential photometry between the companion and the star can be carried out directly. It may be possible to reach an even deeper precision through differential spectrophotometry with a dgvAPP360 than the 3.7% level that we achieve here. Indeed, if we compare the detrended light curve RMS as a function of bin size to the white noise expectation ( Figure  13), it continues to follow the trend of the white noise and does not plateau implying that we have not yet reached any noise floor. This means that in principle the precision of the differential light curve would improve further if more data from additional epochs was added. This also indicates that this technique should remain usable for companions with less favourable contrasts than HD 1160 B, such as those in the planetary-mass regime. However, more data per bin will be required to achieve the same precision for a fainter companion, so the time-sampling in the binned light curves may be less fine in these cases.
Many transiting exoplanet studies make use of a region of the target light curve that is expected to be flat (i.e. an out-oftransit baseline) to test the degree to which systematics have been corrected. While we have shown here that key systematic trends are successfully removed in the detrended differential white-light curve of HD 1160 B, we do not have such a baseline to verify the level of impact of any remaining systematics. The possibility therefore remains that an unknown systematic could be present that has not been accounted for by any of the processes that we have applied here, and could be responsible for the variability that we see in the light curve of HD 1160 B. However, this is also inherently the case for any study that explores the variability of isolated brown dwarfs and planetary-mass objects, stellar variability due to star spots or other sources, or transiting exoplanet studies where exoplanets transit variable stars (e.g. Gelino et al. 2002;Rockenfeller et al. 2006;Biller et al. 2013Biller et al. , 2015Biller et al. , 2021Girardin et al. 2013;Radigan et al. 2014;Wilson et al. 2014;Naud et al. 2017;Eriksson et al. 2019;Vos et al. 2019;Manjavacas et al. 2021Manjavacas et al. , 2022. A subsequent study is forthcoming in which we further investigate the precision that can be reached with this technique, and use injection-recovery tests to assess the extent to which known, simulated variability signals can be recovered (Sutlieff et al., in preparation).
Nonetheless, ground-based differential spectrophotometry with the vAPP is highly complementary and advantageous to spacebased approaches for measuring the variability of high-contrast companions. There have been many successful space-based measurements of companion variability using HST, detecting variability with amplitudes down to the 1-2% level in some cases (e.g. Manjavacas et al. 2018Manjavacas et al. , 2019bBowler et al. 2020b;Zhou et al. 2020a,b). Zhou et al. (2016) was further able to detect sub-percent variability using HST observations of planetary-mass companion 2M1207b, which lies at roughly the same angular separation as HD 1160 B, albeit with a more favourable contrast. Furthermore, the first variability monitoring with JWST, which should reach an even greater precision, is currently underway as part of the Early Release Science Program (Hinkley et al. 2022). However, while JWST has the sensitivity to image fainter, lower mass companions and measure their variability with great precision, its ∼6.5 m mirror is smaller than those of the largest ground-based telescopes, and thus is cannot outperform large ground-based telescopes with extreme AO at small separations 0.5 at ∼3.5 µm (Girard et al. 2022). This means that companions at the closest angular separations such as Jupiter analogues are for now likely only accessible with groundbased monitoring techniques, for all but the nearest stars (Carter et al. 2021;Kammerer et al. 2022). Ground-based telescopes also uniquely provide access to higher resolution spectrographs, such that line profile variability could be used in Doppler imaging to create 2D global maps of features in exoplanet atmospheres such as storms similar to Jupiter's Great Red Spot (e.g. Crossfield 2014).

Observing strategy
During the observing sequence, we obtained 6 wavelength calibrations at intervals throughout the night to allow us to test whether differences in the wavelength calibration used would lead to differences in the differential light curve. In Figure 7 in Section 4.3, we found that the differential light curve is robust to changes in the wavelength calibration. It is therefore preferable to acquire wavelength calibrations at the start or end of future observations, perhaps along with a single precautionary wavelength calibration at high elevation, and instead obtain additional data on target and minimise pixel offsets.
We also used an on/off nodding pattern to enable background subtraction. However, future studies may wish to consider alternative methods to remove the thermal background such that the amount of time spent on target can be doubled and the entire system stabilised. For example, Doelman et al. (2022) developed an approach whereby 93% of the frames obtained were on-target and the thermal background was modelled and removed using the science frames themselves and a small number of background frames obtained before and after the observing sequence. As Figure 13 shows that we approach the photon noise limit, increased on-target time would therefore allow a greater differential light curve precision to be obtained in a single night of observations. Lastly, the absence of an instrument derotator in the LBTI architecture meant that the field of view was rotating throughout the observing sequence. In addition to the drift due to lenslet flexure and the manual offsets applied to keep HD 1160 A centred in the field of view, this meant that HD 1160 A and B were not pixel-stabilised during our observations. However, the linear regression correlation coefficients of the positions of the star and the companion are small relative to those of airmass and air temperature. This suggests that that, when present, these factors dominate over any effect from HD 1160 A and B not being pixel-stabilised.
Understanding whether the host star of a given target is itself varying is important when interpreting the trends in a differential light curve. Most, if not all, potential targets for differential spectrophotometry will be present in at least the TESS full frame images available on the MAST archive. Even though this data will most likely not be contemporaneous with a particular set of observations, the total baseline of the coverage should be relatively long and therefore sufficient to check a host star for variability at the required precision, especially if the target appears in multiple TESS sectors. We therefore recommend this method as a good way to verify the level of variation shown by the host star of a target for differential spectrophotometry, and hence whether it is stable enough to act as a simultaneous photometric reference without requiring further analysis to account for stellar variability.

Outlook
In principle, the technique presented in this paper can be applied to any vAPP coronagraph used in combination with an IFS. Although ALES is currently the only IFS operating over the L and M bands (Doelman et al. 2021), a vAPP coronagraph is also available on SCExAO/CHARIS on the 8.2-m Subaru Telescope, offering R∼19 spectrographic coverage over the J, H, and K bands (1.13-2.39 µm) (Groff et al. 2016;Doelman et al. 2017;Bos et al. 2019;Lozi et al. 2020;Miller et al. 2021). There will also be two different vAPPs on the Mid-infrared ELT Imager and Spectrograph (METIS) instrument on the upcoming 39-m Extremely Large Telescope (ELT), for which this work is a pathfinder. METIS will provide high spectral resolution spectroscopy (R∼100,000) over the L and M bands (Carlomagno et al. 2016;Brandl et al. 2018Brandl et al. , 2021Kenworthy et al. 2018b). Variability measurements using a vAPP may even be possible for broad-band imaging data where an IFS is unavailable, although sensitivity will be inherently more limited without the benefits of using differential spectrophotometry to reduce the effects of systematics. There are several vAPPs currently available on such coronagraphic imagers, such as MagAO (Morzinski et al. 2016;Otten et al. 2017;Sutlieff et al. 2021) and MagAO-X (Miller et al. 2019;Close et al. 2020b) on the 6.5-m Magellan Clay Telescope, and the recently-commissioned Enhanced Resolution Imager and Spectrograph (ERIS) instrument on the VLT (Boehle et al. 2018(Boehle et al. , 2021Kenworthy et al. 2018a;Dubber et al. 2022), with more planned, including for GMagAO-X on the GMT (Close et al. 2020a) and MICADO on the ELT Perrot et al. 2018). Using the vAPPs that will be available on larger telescopes, variability monitoring through differential spectrophotometry will be possible for fainter companions at closer angular separations, including those in the exoplanet mass regime. While this will be inherently more challenging for such companions at greater contrasts, these will remain accessible to this technique through the addition of data from multiple epochs as long as the systematic noise floor is not reached, albeit with a trade-off between light curve precision and time-sampling. In the era of extremely large telescopes, high-contrast imaging combined with high resolution spectroscopy will provide access to fainter companions at lower masses and older ages and allow their orbital velocities and spin to be measured (e.g. Snellen et al. 2014Snellen et al. , 2015Schwarz et al. 2016;Birkby 2018;Wang et al. 2021b;Xuan et al. 2022). Measurements of how individual absorption lines change in depth and width as an exoplanet rotates will allow twodimensional surface maps of exoplanet atmospheres and weather to be produced, through techniques such as Doppler imaging (e.g. Luger et al. 2021;Plummer & Wang 2022). Further in the future, multi-wavelength variabil-ity measurements obtained in reflected light may even enable exocartography of directly imaged Earth-like exoplanets (e.g. Luger et al. 2019Luger et al. , 2022Kawahara 2020;Kuwata et al. 2022;Teinturier et al. 2022).
A limitation of using the dgvAPP360 for variability measurements is that post-processing algorithms relying on angular diversity, such as ADI and PCA, cannot be used without also removing the central PSF of the star that we use as the simultaneous photometric reference. Furthermore, the stellar PSF must remain unsaturated throughout the observing sequence. This potentially limits the sample of targets with bright enough companions. Although HD 1160 B is bright enough that additional noise reduction techniques were not necessary to produce a detection of ample S/N, this may not be the case for fainter directly imaged companions in the exoplanet mass regime. However, it may be possible to use novel alternative algorithms to reach deeper contrasts.
For example, the Temporal Reference Analysis of Planets (TRAP; Samland et al. 2021, Liu et al., submitted) algorithm instead relies on temporal diversity. TRAP reconstructs the systematics in a given region in the data using reference pixels that share the same underlying noise sources. By simultaneously fitting the model of a companion signal 'transiting' over detector pixels and the light curves of the reference pixels, TRAP can then remove these systematics. It may be possible to leverage the information provided by TRAP to improve the companion S/N without removing the stellar PSF, or even to extract detrended light curves directly. Another option would be to use the gvAPP coronagraph, which is different from the dgvAPP360 in that it creates two images of the target star each with a 180°D-shaped dark hole on opposing sides, as well as an additional 'leakage term' positioned between the two (Snik et al. 2012;Otten et al. 2014b;Sutlieff et al. 2021). The leakage term is an entirely separate PSF of the star that appears at a fraction of its full brightness, making it ideal as a simultaneous photometric reference and enabling observations of systems with host stars that would otherwise be too bright. Post-processing algorithms can be applied to the main PSFs to reach deeper contrasts without impacting the leakage term, enabling differential variability measurements provided that the impact of the algorithms on the photometry of the companion can be characterised precisely. However, the gvAPP coronagraph can suffer from wavelength-dependent smearing, which would make such measurements more complex than those obtained with a dgvAPP360 (Otten et al. 2017). In addition to the leakage term, some vAPPs (such as the VLT/ERIS gvAPP) produce other faint reference spots specifically for use as photometric references in situations where the core of the target star PSF is saturated (Doelman et al. 2021;Kravchenko et al. 2022).
In addition to probing the intrinsic variability of high-contrast companions, differential spectrophotometry could also be used to observe the transits of satellites such as exomoons or binary planets passing in front of these companions (e.g. Heller 2016; Lazzoni et al. 2022). Candidate satellites have been identified around transiting exoplanets, directly imaged companions, and isolated planetarymass objects using a range of techniques, but none have yet been definitively detected (e.g. Teachey & Kipping 2018;Rodenbeck et al. 2018;Heller et al. 2019;Kreidberg et al. 2019;Lazzoni et al. 2020;Teachey et al. 2020;Limbach et al. 2021;Vanderburg & Rodriguez 2021;Kipping et al. 2022). For directly imaged companions, variability arising from transit events could be distinguished from that caused by inhomogeneous atmospheric features in similar ways to transiting exoplanets and star spots, by considering the companion light curves across the different wavelength channels. Transit signals are expected to be almost achromatic, while intrin-sic variability is generally wavelength-dependent (Manjavacas et al. 2019a;Limbach et al. 2021Limbach et al. , 2023Lazzoni et al. 2022). Lazzoni et al. (2022) found using simulations that although the probability of successfully detecting smaller exomoons around a directly imaged companion is very low with current instrumentation and techniques, detections of larger binary planets are already within reach. New techniques to improve differential light curve precision for directly imaged companions, including differential spectrophotometry with the dgvAPP, will help to increase these probabilities further and potentially enable the first definitive detections of satellites around substellar companions.

CONCLUSIONS
We present a novel, ground-based approach for constructing differential light curves of high-contrast companions through direct differential spectrophotometric monitoring, using the double-grating 360°vector Apodizing Phase Plate coronagraph and the ALES integral field spectrograph. The dgvAPP360 allows high-contrast companions to be detected while also providing an image of the host star, which crucially can be used as a simultaneous photometric reference. We combine the dgvAPP360 with ALES to follow the highly successful technique of differential spectrophotometry used in exoplanet transmission spectroscopy, where light is spectrallydispersed to reduce systematic effects that otherwise dominate the variability signal we aim to measure, and then recombined into white-light flux measurements.
We demonstrated this approach using a full night of observations of substellar companion HD 1160 B. The time-series fluxes of the companion and the star in each wavelength channel were extracted simultaneously using aperture photometry. We then produced white-light measurements for both the companion and the star at each time frame by taking the median combination of the photometry in the wavelength dimension. The companion flux was then divided by that of the star to eliminate trends common to both, arising from Earth's atmosphere and other systematics, producing a differential white-light curve that only contains non-shared variations and covers a wavelength range of 3.59-3.99 µm. We find that the shape of the resulting light curve is robust against issues arising from instrumental flexure, as tested using calibration frames collected throughout the observation sequence. Using a multiple linear regression approach with eight decorrelation parameters, we modelled and removed non-shared trends from the differential whitelight curve. We find that airmass and air temperature are the most correlated parameters with the light curve. We also analyse publicly available data from the TESS mission to check for variability in the host star HD 1160 A, and confirm that it is non-varying at the 0.03% level.
We find that the detrended differential white-light curve of HD 1160 B shows sinusoidal-like variability over a short timescale. By fitting the unbinned light curve with a sinusoid, we identify that the variability has a semi-amplitude of ∼8.8% and a period of ∼3.24 hours. When binned to 18 minutes of integration time per bin, we achieve a light curve precision at the 3.7% level. After thorough investigation and rejection of systematic noise sources, we attribute this variability as likely due to heterogeneous features in the atmosphere of the companion, rotating in and out of view as it rotates. We find that if the period of this variability reflects the rotation period of HD 1160 B, physical limitations suggest that it is rotating at close to its breakup period. Alternatively, the short period variability in the light curve of HD 1160 B may arise from multi-ple periodic features in its atmosphere with different phase offsets. Furthermore, light curves in the 30 individual wavelength channels in the 3.59-3.99 µm range show tentative evidence of an increase in variability amplitude at longer wavelengths. Further observations at additional epochs will help to confirm and characterise the variability of HD 1160 B and to determine its physical explanation.
The precision that we achieve in the detrended differential white-light curve is the greatest achieved from ground-based studies of sub-arcsecond high-contrast companions to date. However, direct comparisons to other ground-based studies that instead use satellite spots to search for variability in the light curves of highcontrast companions are not straightforward due to the different magnitude and contrast of the observed systems, with HD 1160 B having a more favourable contrast Biller et al. 2021;Wang et al. 2022). The RMS of the detrended differential light curve for HD 1160 B as a function of bin size follows the same trend as the theoretical white noise expectation with no evidence of yet approaching a noise floor. This indicated that the single night of data analysed here is not yet systematic-limited, and that further observations from additional epochs could enable greater sensitivity to be reached. A deeper investigation of this type of data and its precision, including injection-recovery tests to test how effectively known variability signals can be recovered and which systematics have the greatest impact, is forthcoming (Sutlieff et al., in preparation).
While JWST will measure the variability of fainter, lower mass companions from space with unprecedented precision, its comparatively smaller aperture means it cannot outperform the largest AO-equipped ground-based telescopes at separations 0.5 at ∼3.5 µm (Girard et al. 2022), so companions such as Jupiter analogues at the closest angular separations, for all but the nearest stars, remain accessible only to ground-based monitoring techniques for the coming decade. Ground-based differential spectrophotometry with the vAPP is therefore highly complementary to space-based approaches for measuring the variability of high-contrast exoplanet and brown dwarf companions, and for searching for their transiting exomoons or binary planets. Moreover, ground-based telescopes can reach much higher spectral resolution, which then enables line profile variability studies to map atmospheric features, including storms and hurricanes like Jupiter's Great Red Spot, via Doppler imaging. These results are promising for further variability studies using vAPP coronagraphs on current and upcoming instruments and telescopes, which include ERIS on the VLT and METIS on the ELT. gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper includes data collected with the TESS mission, obtained from the Mikulski Archive for Space Telescopes (MAST) data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of NASA's Astrophysics Data System. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research made use of SAOImageDS9, a tool for data visualization supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEASARC) with support from the JWST Mission office at the Space Telescope Science Institute for 3D visualization (Joye & Mandel 2003). This work made use of the whereistheplanet 3 prediction tool (Wang et al. 2021a). This work makes use of the Python programming language 4 , in particular packages including NumPy (