© The Authors 2025

1 Introduction

Direct imaging of exoplanets is a challenging type of observation in modern astronomy due to the high star-planet contrast and small angular separation. Although only ~1% of currently known exoplanets have been directly imaged, these observations permit detailed characterization of their atmospheres. Present and upcoming improvements in high-angular-resolution capabilities, including 30 m-class telescopes (Quanz et al. 2015; Bowens et al. 2021) and interferometry (ground-based: Defrère et al. 2018; Wagner et al. 2021a,b; Ertel et al. 2022; or spacebased: Kammerer et al. 2022), open up the Habitable Zone (HZ) of nearby main-sequence stars for direct imaging. The HZ commonly refers to the circumstellar region where planetary surfaces may sustain liquid water given an adequate atmosphere (Kopparapu et al. 2013). This region is relevant for the search of potential exoplanets located at the Earth-equivalent insolation distance (EEID) - the orbital distance where objects receive the same irradiance as Earth - and the assessment of their habitability. However, this region can also be rich in debris material, similar to the zodiacal cloud of the Solar System, known as exo-zodiacal dust. Like other debris disks found at larger separations, exozodiacal disks are expected to appear at the end of planetary formation, during the disk-clearing phase (between 10 Myr and a few Gyr age, Hinz et al. 2009; Ercolano & Pascucci 2017). Collisions and evaporation of planetesimals serve as sources of exozodiacal dust (Kennedy & Piette 2015; Marboeuf et al. 2016; Faramaz et al. 2017; Rigley & Wyatt 2022) while other phenomena tend to deplete it (e.g., Poynting-Robertson (P-R) drag: Wyatt & Whipple 1950; Burns et al. 1979; planetary sweeping: Bonsor et al. 2018; or radiation pressure: Burns et al. 1979. The study of exozodiacal disks around nearby main-sequence stars serves several main purposes. It enables investigation of the properties and dynamics of the dust (Mennesson et al. 2014; Defrère et al. 2015; Lebreton et al. 2016), provides insights into the architecture and formation mechanisms of planetary systems (Stark & Kuchner 2008), and allows evaluation of cometary bombardment and its impact on the habitability of inner planetary systems (Chyba et al. 1990; O’Brien et al. 2014; Rotelli et al. 2016; Kral et al. 2018; Wyatt et al. 2020; Ritson et al. 2020). Studying these disks is also important for estimating their impact on rocky planet imaging, both from the shot noise they generate in the infrared (Defrère et al. 2010; Stark et al. 2014, 2019), and from their scattering of stellar light in the visible (Beichman et al. 2006; Defrère et al. 2012; Roberge et al. 2012).

A survey of exozodiacal disks was carried out by the Hunt for Observable Signatures of Terrestrial Systems (HOSTS: Ertel et al. 2018, 2020a) which quantified the presence of warm dust (~300K, i.e., located in the HZ) and its surface density around nearby (<30 pc) main-sequence stars. The survey was performed with the nulling mode of the Large Binocular Telescope Interferometer (LBTI: Defrère et al. 2015; Hinz et al. 2016; Ertel et al. 2020b). As part of the survey, the F7V type star θ Boo was observed. This star, located at 14.5 pc, showed signs of HZ dust emission estimated at 148.2 ± 27.7 zodi (1 zodi being the dust surface density at the EEID of the Solar System, i.e., 1 au; see definition in Weinberger et al. 2015; Kennedy et al. 2015). Table 1 summarizes the stellar parameters of θ Boo used in this study. Among these, the stellar luminosity L_★ is used to calculate the local EEID following the method presented in Weinberger et al. (2015), where the radius r_EEID is proportional to $\sqrt{{\text{L}}_{\star }}$. For θ Boo, we find r_EEID ≃ 2 au, which corresponds to a separation of ~ 138 mas. The star was previously observed by the space telescopes Spitzer (2004) and Herschel (2010) at 24, 70, and 100 μm wavelength to search for warm (~600 K) and cold dust (~120 K) emission. These surveys did not show significant emission excess within their detection limits (Bryden et al. 2006; Trilling et al. 2008; Gáspár et al. 2013; Montesinos et al. 2016). Figure 1 shows the observed spectral energy distribution (SED) of θ Boo, and the upper limits for the dust emission obtained from space surveys¹. The upper limits are estimated following the method used in Yelverton et al. (2019), in which the observed data are fitted by a model with both a stellar and a dust component. The stellar component is a PHOENIX BT-Settl model (Allard et al. 2012), and the dust component is a modified blackbody. The error bars of the observed data are then propagated to the results of the dust model.

This raises the question of the origin of the exozodiacal dust in the system. The presence of such warm dust is strongly correlated with the presence of a cold debris disk (Ertel et al. 2020a). Statistical results from the HOSTS survey support this correlation, finding a detection rate of HZ dust of ${78}_{-18}^{+8}$% in systems with a known cold debris disk. Only 3 of 28 stars have known HZ dust without cold debris, with θ Boo being one of them. To investigate the nature of the exozodiacal disk of θ Boo, a new observation in nulling mode was performed by the LBTI in 2023 that covered a wider range of parallactic angles than those obtained during the HOSTS survey. This work is part of an ongoing effort to study individual targets of the HOSTS survey (βLeo: Defrère et al. 2021; αLyr: Faramaz et al., in prep.; 110 Her: Rousseau et al., in prep.; εEri: Weinberger et al. in prep.; βUma: Bryden et al. in prep.), and to investigate the nature of these disks.

In this paper, we present a descriptive model for the spatial distribution of dust around θ Boo. We also constrain the presence of companions using observations from the LBTI/L- and M-band InfraRed Camera (LMIRCam) and the System for coro-nagraphy with High-order Adaptive optics from R to K band (SHARK)-Near InfraRed (NIR) at the Large Binocular Telescope (LBT). Section 2 describes the different observations of θ Boo used in this study. Section 3 explains the data reduction processes for these observations and their results. The constraints on the presence of possible companions in the system are presented in Sect. 4. The models of the dust distribution and its variability at the EEID are then shown in Sect. 5. Finally, we discuss the constraints on the dust distribution, its possible origins, and the source of its variability in Sect. 6.

Fig. 1 Observed spectral energy distribution (SED) of θ Boo. The black points indicate the observed photometric data. The gray line represents the best fit stellar model, and the green dashed line shows the dust model. The green points indicate the upper limits of the dust emission.

2 Observations of θ Boo

2.1 Instrument description

The nulling observations of θ Boo were conducted at the LBT on UT 2017 April 11 and UT 2018 May 23 with the LBTI/Nulling Optimized Mid-Infrared Camera (NOMIC) in the N′-band (9.81-12.41 μm) as part of the HOSTS survey. An additional nulling observation with NOMIC was obtained on UT 2023 May 25. High-contrast adaptive optics (AO) imaging with LMIRCam in the L′-band (3.41-3.99 μm) and with SHARK-NIR in the H-band (1.38-1.82 μm) was performed on UT 2024 February 24. Table 2 gives an overview of these observations.

Hill et al. (2012) describe the LBT, which consists of two 8.4m diameter apertures, each corrected with AO. The AO system achieves Strehl ratios of 80%, 95%, and 99% at 1.6 μm, 3.8 μm, and 10 μm, respectively, using Pyramid Wave Front Sensing (WFS) with the visible portion of the light (Bailey et al. 2010; Bailey et al. 2014; Pinna et al. 2016). The infrared portion of the light is directed toward the Nulling and Imaging Camera (NIC), where it is split and focused onto different cameras as explained in Hinz et al. (2008).

The LBTI is the instrument that combines the light from the two LBT apertures to perform interferometric observations (Hinz et al. 2016). LMIRCam and NOMIC can be used in this configuration by imaging both beams separately or by overlapping them. Figure 2 provides an overview of the LBTI architecture. A K-band fast readout camera, called Phasecam, monitors the phase difference between the two apertures (Defrère et al. 2014). It measures the differential piston and corrects it using two piston correctors located upstream on the tip-tilt mirrors.

The SHARK-NIR instrument is located near the Gregorian focus of the left aperture and performs observations in the near-infrared (Farinato et al. 2022; Marafatto et al. 2022). The LMIRCam instrument covers wavelengths from the J- to M-bands, while NOMIC operates in the N-band.

The only relevant change made at the LBT between the initial observations of 2017 and 2018 and those in 2024 was the AO upgrade (Pinna et al. 2016). This upgrade is not expected to have a major impact on the nulling observations, as the previous AO system already achieved a similar Strehl ratio in the N-band on the HOSTS targets. Moreover, any change in AO performance would have resulted in changes in instrumental null or phase jitter, which are considered in the data reduction pipeline and the error estimations.

Fig. 2 Block diagram of the LBTI architecture. Light from the two apertures enters from the top left and right and passes through the beam combiner subsystem (purple) and the NIC cryostat (blue). Both the beam combiner and the NIC are maintained at cryo-temperatures. At the entrance of the beam combiner subsystem, visible light is reflected by dichroics and directed to the wavefront sensing (WFS) units. Infrared light is transmitted to the beam combiner where two mirrors adjust the tip-tilt and the optical path difference (OPD) of the beams for nulling interferometry. After recombination of the two LBTI pupils by the roof mirror, the light is split in the NIC. The L-band (3-5 μm) is directed to LMIRCam for high-contrast imaging of exoplanets, whereas the K-band (1.5-2.5 μm) is directed to Phasecam to monitor the differential tip-tilt and phase between the two apertures. Finally, the N-band (8-13 μm) light is directed to NOMIC which can be used in either imaging interferometry or nulling interferometry modes.

2.2 N′-band imaging observations

In the N′-band, the NOMIC camera was used in a configuration designed for nulling observations. In this mode, the phase difference between the two apertures is locked at π. When combined in the pupil plane, light from a point source at the center of the field of view (FOV) destructively interferes and is not observed when reimaged at the detector plane. However, off-axis sources add an extra phase difference between the apertures, depending on their on-sky position. This results in emission from off-axis sources and any sufficiently extended stellar components (termed “stellar leakage” or “instrumental null depth”) being imaged on the NOMIC camera.

The throughput of off-axis sources as a function of their on-sky position is described by the “transmission map” T. For the LBTI, which has two apertures separated by a baseline B = 14.4 m, T = sin² (πφB/λ), where φ is the on-sky position projected onto the baseline vector and λ is the wavelength (Kennedy et al. 2015). The transmission of an object on the camera therefore depends on its on-sky position and the parallactic angle rotation during the night (i.e., the baseline rotation). The coherent FOV θc - the on-sky region where light sources can interfere - is defined as $${\theta }_c\equiv \frac{{\lambda }^2}{B\times \mathrm{\Delta }\lambda },$$(1)

where Δλ is the spectral bandwidth (Thompson et al. 1986). For nulling in the N′-band (9.81-12.41 μm), the coherent FOV is θ_c ≃ 670 mas.

The amount of instrumental null depth can be estimated using a set of calibrator stars (labeled “CAL”) observed before and after the science target θ Boo. These calibrators were selected following Mennesson et al. (2014), using the two catalogs from Bordé et al. (2002) and Mérand et al. (2005), with additional stars from the JSDC catalog and the SearchCal tool from Chelli et al. (2016). Table 3 lists the calibrators used in this study and their properties.

Observations consisted of several observation blocks (OBs), each comprising 2000 × 42.7 ms frames. Between nulling OBs, the telescope was offset by ±2.3 ″ at the detector in a vertical nod movement. The instrument uses an Aquarius detector, with a FOV of 18″ × 18″ and a pixel scale of 17.9 mas. This allows for measurement of the mid-infrared thermal background noise at pixels that may have received a signal in previous observations. Each observation includes multiple nulling OBs (with overlapping beams), one photometric OB (with separated beams), and one background OB (with beams nodded off the detector). The N′-band observations on 2017 April 11, 2018 May 23, and 2023 May 25 targeted exozodiacal dust at a separation range of 40-700 mas (i.e., 0.6-10.4 au for θ Boo). Table 2 summarizes the three observing nights. The observation sequences for each night were as follows: CAL1-θBoo-CAL2 (2017 April 11); CAH-θ Boo-CAL3 (2018 May 23); and CAL1-θBoo-CAL4θ Boo-CAL5 (2023 May 25). Each sequence began with CAL1, enabling comparison of results across the three nights.

2.3 L′ -band imaging observations

We obtained L′-band observations with LMIRCam on the LBT right arm (with one 8.4m aperture) using Director’s Discretionary Time (DDT). These observations, conducted on 2024 February 24, targeted a potential companion at a separation range of 0.3-4″ (i.e., 4-58 au for θ Boo, see Sect. 3.2). They consisted of28 sequences of exposures, each with 2000 images with 55 ms integration time, totaling 3080 s (~51 minutes). We did not use a coronograph and centering was achieved using a rotationalbased algorithm following Morzinski et al. (2015). The Teledyne H2RG detector provided a FOV of 20″ × 20″ and a pixel scale of 10.9 mas. We acquired observations in pupil-stabilized mode that allowed the rotation of the FOV to perform angular differential imaging (ADI: Marois et al. 2006). During observations, the parallactic angle rotated by 104°, and the seeing fluctuated between 1.05″ and 2″. To monitor the thermal sky background and detector drifts, the stellar image was nodded up and down on the detector by an offset angle of 4.5″ between every observing block.

2.4 NIR imaging observations

In parallel with LMIRCam, we conducted simultaneous observations with SHARK-NIR operating on the LBT left arm. Using a broadband H filter (central wavelength 1.6 μm, bandwidth 0.218 μm) and a Gaussian Lyot coronagraph with an inner working angle (IWA) of 150 mas, we acquired data with a Teledyne H2RG detector (FOV: 18″ × 18″; pixel scale: 14.5 mas). To avoid detector saturation near the star and minimize exposure time, we selected an observing mode where only a central region of dimensions 2048 × 200 pixels was read out. This limits the usable region for companion detection, and for the definition of contrast, to a separation of 1.3″ (see Sect. 3.3). Each frame had a 0.84 s exposure time, and we acquired 6855 frames. This gave a total exposure time of 5758.2 s (~96 minutes) on the target. The observations were also taken to allow ADI during postprocessing. To correctly estimate the obtained contrast, images of the stellar point spread function (PSF) without the coronagraph were taken at the start and end of the observations. These frames were taken using an appropriate neutral density filter (ND3) to avoid detector saturation.

Fig. 3 Null depths obtained after the data reduction pipeline described in Sect. 3.1. These points are computed using the reduction pipeline with an aperture radius of 24 pixels (i.e., 430 mas). The OBs are plotted as a function of their parallactic angle (left) and UT of observation (right), for each of the three nights: 2017 April 11 (top), 2018 May 23 (middle), and 2023 May 25 (bottom). In each panel, light blue points and light red crosses show the measured null depth of the OBs for the calibrators and for θ Boo, respectively. Dark blue points and dark red crosses indicate the null depth values for each pointing, corresponding to the unweighted mean values of the OBs included in the pointing. The instrumental null depths measured from the calibrators are averaged to calculate the null floor, which represents the reference level below which any excess cannot be reliably detected. The null floor and its error bar are shown by the solid and dotted black lines, respectively. The error bar corresponds to the 1σ standard deviation.

3 Data reduction

3.1 N′-band nulling data

We used the LBTI nulling pipeline from Defrère et al. (2016) for the N′ -band nulling data reduction and calibration. The pipeline first corrects the raw images by removing bad pixels and subtracting the mid-infrared thermal background. Because the background cannot be measured at the star position when the star is present, it is estimated from the series of OBs with an opposite nod position. The background is then subtracted using either the mean value of the series or a principle component analysis (PCA) approach (Rousseau et al. 2024). In this work, we used the mean background subtraction. The pipeline then computes the flux contained within an aperture radius using the background-subtracted images. The PSF of the telescope at 11 μm (λ/D), where D is the diameter of an individual LBT aperture, has a full width at half maximum (FWHM) of ~286mas, which is oversampled by 16 pixels (each 17.9 mas). We therefore selected a standard photometric aperture radius of ~0.5λ/D (i.e., 8 pixels, or 143 mas), which maximizes the signal-to-noise ratio (S/N) for a point-like source’s measured flux.

To probe the dust spatial distribution, we performed photometric measurements using several aperture radii: 8, 16, 24, 32, and 40 pixels (i.e., 143, 286, 430, 573, and 716 mas). After computing the flux for each radius, the pipeline calculates and calibrates the null for each OB using the nulling self-calibration (NSC) approach, originally developed for the Palomar Fiber Nuller (Hanot et al. 2011; Mennesson et al. 2011) and adapted for the LBTI (Defrère et al. 2016; Mennesson et al. 2016). This approach removes errors in the nulling setpoint between the science star and its calibrators. The computed null is also corrected for the instrumental null depth generated by the science target and calibrators (see Sect. 2.2). The limb-darkened angular diameters θLD and their 1σ uncertainties are listed in Table 3. The null depths obtained for the calibrators and θ Boo are shown in Fig. 3. The instrumental null depths from the calibrators define a “null floor”, corresponding to their mean value. The difference between the θ Boo nulls and the null floor is attributed to dust emission passing through the transmission pattern T (see Sect. 2.2). This difference is referred to as the “calibrated null depth.” Small sample statistics (as described in Mawet et al. 2014) are not used in the NSC for small apertures. The statistical noise is estimated from the time series of images and the flux measurement in each image. With approximately 2000 images per OB and several OBs per observation, the noise estimated by the NSC is not affected by small sample statistics.

After calculating the calibrated null depth for each OB of the science target with a given aperture radius, we take their mean value as the calibrated null depth for that radius. The calibrated null depths for all apertures jointly probe the excess flux from the dust and its radial distribution. The left panel of Fig. 4 shows the data reduction results for the three nulling observation nights. The null depths of the OBs associated with each point are listed in Appendix A. The difference between each aperture radius is 8 pixels (143 mas), corresponding to 0.5λ/D. The different calibrated nulls are correlated due to the PSF size, so a model is required to determine the dust location in the system. Modeling results for these observations are presented in Sect. 5.

To assess whether these results depend on the parallactic angle, we analyzed the two pointings from the 2023 May 25 observation individually. Figure 3 shows the parallactic angle ranges for the two pointings, and the right panel of Fig. 4 plots the calibrated null excess for the pointings individually and for the combined data. The results show no significant differences in calibrated nulls compared to the combined results, indicating that the excess can be assumed to be independent of the on-sky parallactic angle. This suggests that the exozodiacal disk is likely close to face-on if the disk is vertically thin, which is useful for modeling the dust distribution in Sect. 5. Additionally, this allows us to compare results from the three nights of observation despite their different parallactic angle ranges.

We used the calibrated nulls from Fig. 4 to estimate the total infrared excess from the dust measured by the LBTI at 11 μm. Estimating this excess without prior assumptions on the dust angular distribution is challenging due to the transmission map T. Since T follows a sin² function with the on-sky position parallel to the baseline vector, we assume that approximately 50% of the excess generated by the dust in the FOV is reimaged on the camera (Rigley & Wyatt 2020). To estimate the total excess for each observing night, we take the aperture radius with the highest calibrated null, and double it. This yields total emission excesses at 11 μm of 1.83±0.49%, 2.49±0.66%, and 2.92±0.51% of the stellar emission for the 2017, 2018, and 2023 observing nights, respectively. The theoretical blackbody emission of θ Boo is ~2315 mJy, which gives an excess of 42± 11, 58± 15, and 68±12 mJy from the dust for the 2017, 2018, and 2023 observing nights, respectively. These estimates are consistent with the upper limit shown in Fig. 1. An increase in the total excess by ~2σ is observed between 2017 and 2023. The mean excess value is 56±7.4 mJy.

Fig. 4 Left: calibrated null depths for the different aperture radii and for the three nights of observation with NOMIC. The aperture radii are 8, 16, 24, 32, and 40 pixels (i.e., 143, 286, 430, 573, and 716 mas, respectively). Right: calibrated null depths for the 2023 May 25 observation as a function of aperture radius. Results for the two individual pointings are shown as dotted and dashed lines, and the combined results are shown as crosses. The first pointing corresponds to the OBs taken from 04:00 to 06:00 UT and the second pointing is from 06:30 to 07:30 UT.

3.2 L′-band direct imaging data

We processed the LMIRCam data using the LBTI Exozodi Exoplanet Common Hunt (LEECH) survey pipeline (Stone et al. 2018). The pipeline first replaces bad pixels with the median of the eight closest pixels. As with the N′-band observations, the thermal background is subtracted using the median value of the opposite nod images from the closest OB. The images are then coarsely corrected for distortion using the dewarp coefficients of Maire et al. (2015) which remain sufficiently precise, as astrometric precision is not required. The telescope PSF has a FWHM of 95 mas at 3.7 μm. The detector pixels have an on-sky size of 10.7 mas, so the PSF is oversampled by a factor of about 4. Each image is therefore binned 2 × 2 to remove cosmic rays or any remaining bad pixels.

The pipeline then fits and removes the diffracted light from the star using PCA (Amara & Quanz 2012; Soummer et al. 2012; Gonzalez et al. 2017) before derotating and stacking the images. To optimize high contrast capability, the PCA method was performed annulus by annulus with each annulus having a width of 9 pixels (~2λ/D) and a subtraction region annulus of 1 pixel width. Artificial planets were injected at different angular radii, and the number of PCA components was chosen to maximize their S/N. The best contrast was achieved after several iterations. The planet injection and contrast estimation procedures are detailed in Wagner et al. (2019). The two nod positions were reduced separately, but their images were recombined using a weighted mean. Weights were set for each annulus to optimize the detection of the artificial planets. Regions of poor sensitivity due to dust diffraction within LMIRCam were down-weighted accordingly.

The left panel of Fig. 5 shows the reduced image of θ Boo. No feature resembling a companion can be observed between 0.3″ and 4″ . The detection limits for this observation are presented in Sect. 4. To account for the reduced number of resolution elements near the star - and its impact on the confidence level of the detection thresholds - the detection limits were calculated including small sample statistics (Mawet et al. 2014).

Fig. 5 Left: ADI image of θ Boo obtained with LMIRCam after processing with the LEECH-survey pipeline. The image has a FOV of 4″. A central mask of 5λ/D radius is used to improve the visualization. No feature indicating the presence of a companion is detected between 0.3″ and 4". Right: ADI image of θ Boo from SHARK-NIR after processing with the SHARP pipeline using five PCA components. The image has a FOV of 1.3″. For both images, the outer circle indicates the coherent FOV of NOMIC in the N′-band (~670 mas; see Eq. (1)), and the inner dashed circle represents the EEID (~ 138 mas). No feature indicating the presence of a companion is detected between 0.15″ and 1.3″.

3.3 NIR direct imaging data

The SHARK-NIR data were reduced using the Python pipeline called SHARP, currently under development. The pipeline first subtracts the dark and divides the image by the flat. The internal deformable mirror generates four satellite spots that are images of the star with lower brightness and symmetric with respect to the position of the target, arranged in a cross configuration that can be used to determine the position of the star behind the coronagraph. Fitting two lines through two opposite spots, their intersection accurately defines the central position of the stellar PSF. Frames of poor quality, for example those not adequately masked by the coronagraph due to bad weather conditions, were removed. After this procedure, 6422 frames remained. To improve the reduction time, the frames were grouped in sets of four and averaged, forming a data cube of 1605 frames. Speckle subtraction was then performed on this cube, using the PCA algorithm, followed by image derotation and stacking. The PCA was run with 1, 5, 10, 15, 20, and 25 components. A reduced image was produced for each case; however, no significant differences were observed between the different images.

The right panel of Fig. 5 shows the resulting SHARK-NIR image using five PCA components. No feature resembling a companion is observed between 0.15″ and 1.3″. We then calculate the brightness contrast by measuring the standard deviation in one-pixel-wide rings and dividing by the normalization factor obtained using the non-coronagraphic images described in Sect. 2.4. The detection limits corresponding to this observation are presented in Sect. 4. As with the LMIRCam data, the detection limits were calculated accounting for small sample statistics.

4 Detection limits in the L′-band and NIR

The images obtained from LMIRCam and SHARK-NIR after post-processing (Fig. 5) show no apparent point-like feature, but they can still be used to compute companion detection limits. Figure 6 presents these detection limits (corresponding to S/N=5) as contrast curves relative to angular separation from the star. Constraints on planetary masses are derived using the COND evolutionary model described in Baraffe et al. (2003). The model assumes a “hot start” formation scenario, in which gravitational instabilities in the protoplanetary disk cause dust and gas to collapse and form a planet (Spiegel & Burrows 2012). The collapsing gas conserves its entropy, resulting in a high initial entropy for the planet. The COND model neglects dust opacity by assuming that atmospheric dust immediately falls below the photosphere after formation. Given the contrast curves and stellar parameters, the model provides planetary masses corresponding to various contrast levels. The best detection limit in the NIR is achieved by SHARK-NIR with a limit of 11 M_Jup around 1.3″. The LMIRCam is unable to probe the planetary mass regime, with a limit of 18 M_Jup for separations >1.5″.

The mass constraints for these observations are calculated assuming an age of 500 Myr for θ Boo from Gáspár et al. (2013). However, the age of the system is uncertain, with different methods yielding a range of estimates. Table 4 summarizes the current estimates for each method, which span from 400 Myr to several Gyrs. Rachford & Foight (2009) reported that θ Boo is approximately twice as luminous as a zero-age main-sequence star of the same color. This implies that the star is either at the end of its main-sequence phase or is a binary with equal components. Using the Stefan-Boltzmann law with the stellar parameters listed in Table 1 yields a stellar radius of 1.77 R_⊙. Interferometric observations at the Center for High Angular Resolution Astronomy (CHARA) measured a stellar radius of 1.73 ± 0.01 R_⊙ (Boyajian et al. 2013). The theoretical and measured stellar radii are about 30% larger than that of a zero-age F7V star (i.e., 1.324R, Pecaut & Mamajek 2013). This indicates that θ Boo is likely at the end of its main-sequence phase and not a binary with equal components. Its age is therefore likely toward the higher end of the range in Table 4, between 3 and 4.7 Gyr. For ages beyond 500 Myr, the internal heat resulting from planet formation has already dissipated, and there should be no significant difference in the evolutionary model. The mass limits are therefore still relevant even for this age range. It should also be noted that θ Boo has a binary companion detected at a separation of 70″ (~1000 au, Eggen 1956). However, this separation is too large to impact the presence of a debris disk in the system (Yelverton et al. 2019).

Fig. 6 Detection limits for companions around θ Boo in the L′-band with LMIRCam (left) and in the H-band with SHARK-NIR (right). The blue line shows the contrast as a function of angular separation. For four contrast levels, the corresponding detection limit is given in terms of mass (in units of Jupiter mass). The best mass constraint is obtained in the H-band around 1″, with a limit of 11 MJup. The system is assumed to be 500 Myr for these mass constraints.

5 Dust modeling

In this section, we constrain the spatial distribution of the exo-zodiacal dust. We first describe the dust emission models used to represent exozodi emission. These models are then combined with an MCMC fitting procedure, which allows us to recover information about the shape of the exozodiacal dust disk around θ Boo at different epochs, and thus explore its temporal evolution.

5.1 Model description

Here we describe how we generate exozodiacal dust emission models for comparison with the LBTI observations shown in Fig. 4. Our models are based on the framework presented in Kennedy et al. (2015) - which serves as the baseline for the HOSTS survey at the LBTI.

In this framework, exozodiacal dust disks are assumed to be both optically and vertically thin. The face-on surface density of the dust, Σ, at a distance r from the star is assumed to follow a power-law distribution: $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0\times r^{\alpha },$$(2)

between inner and outer radii R_in and R_out, respectively, with r expressed in au. The disk surface density of cross-sectional area is represented by Σ, and hence is expressed in au²/au², and is analogous to optical depth. The surface density at 1 au is denoted by Σ₀, and α is the surface density slope of the disk. The dust surface brightness S at wavelength λ is then given by $$S(r)=2.35\times {10}^{-11}\mathrm{\Sigma }(r)B_{\nu }(\lambda ,T_{\text{d}}(r)),$$(3)

with B_ν being the dust grain emission at temperature T_d. The numerical factor of 2.35 × 10⁻¹¹ converts the surface brightness to units of Jy.arcsec⁻². The temperature T_d is that of a blackbody and is defined by $$T_{\text{d}}=278.3L_{\star }^{0.25}{r}^{-0.5}K.$$(4)

We summarize the parameters describing the dust spatial distribution in Table 5. The stellar parameters used to compute the dust emission can be found in Table 1.

Geometry of the system. The face-on dust distributions in the model are limited to centrosymmetric geometries with the following parameters: the inner radius R_in, the outer radius R_out, the dust density at 1 au, Σ₀, and the dust density slope, α. We define the mean radius of the belt as R = (R_in + R_out)/2. Three different geometries are considered when modeling the observations. The first is a “wide ring” model with free R_in and R_out and fixed α = 0. The second is a “thin ring” model, which is a sub-model of the wide ring. It has a free R_in, a fixed (R_out - R_in)/R = 0.1², and fixed α = 0. This model is useful to constrain the region where most of the dust emission is located. The third geometry is a “disk” model where R_in is fixed at the sublimation radius of silicates (corresponding to T_BB = 1500 K), R_out fixed at an arbitrarily large distance (typically for T_BB = 88 K³), and α is free. For θ Boo, R_in and R_out are calculated to be 0.07 and 20 au respectively (i.e., 4.7 and 1380 mas).

The parameter Σ₀ is systematically included as a free parameter for each geometry and is constrained in the range [0, 10⁻²]. R_out is constrained in the range [0, 20″], and R_in in the range [0, R_out[. The range of α is [−10, 10]. For the “wide ring” and “thin ring” models, we set α = 0 to simulate a dust belt with a uniform optical depth radial distribution. This assumption is used to study the collisional evolution of cold debris disks, as in Wyatt et al. (2011) (see the P-R drag model description in Sect. 6.2).

Orientation of the system. The orientation of the disk is defined by two parameters: the disk inclination, i, and the disk position angle, PA. As mentioned in Sect. 3.1, the 2023 null depths show no visible dependence on the parallactic angles - and hence on the transmission map orientation - which indicates that the disk inclination is expected to be close to face-on (i = 0°), under the assumption of a vertically thin disk. Several preliminary test runs established that the MCMC procedure did not converge within reasonable computational times when fitting the geometry (R_in, R_out, Σ₀, α) and orientation parameters (i, PA) simultaneously. We therefore assume the disk to be face-on, and perform the MCMC procedure with i systematically set to 0° for the different datasets. The value for PA is then irrelevant.

Fig. 7 Results of the 2023 May 25 night modeling, showing the different steps of the procedure. The geometry used is a disk with Σ0 = 3 × 105 and α = 0, the orientation is face-on. From (a) to (d) show: the angular distribution of Fdust/F★, the impact of the LBTI transmission map, the convolution with the PSF of the LBT, and the comparison between the observed excess and the modeled excess. Panel (c) indicates the photometric apertures used in the reduction with the dashed circles. The coherent FOV is indicated with the solid orange circle. The maps have a pixel size of 10 mas.

5.2 Modeling strategy

Our goal is to model the dust emission as measured with NOMIC. The next step is therefore to compute the dust raw radiance, F_dust, by integrating Eq. (3) over the passband of NOMIC (N′-band: 9.81-12.41 μm), taking into account the transmission of the NOMIC filter. The value of F_dust is then normalized by the stellar radiance, F_★, calculated from its blackbody emission. Finally, the model accounts for the impact of the PSF and the transmission map on the detected flux: the raw radiance map is first multiplied by the LBTI transmission pattern (Sect. 2.2) and then convolved with the modeled PSF of the LBT. The measured excess can then be estimated by integrating the pixels over the different apertures defined in Sect. 3.1, and compared to the results from Fig. 4. This comparison uses an MCMC approach: we generate 30 walkers with flat priors across all parameters, each allowed to perform 500 sampling steps (after a burn-in phase of 1000 models). For each model, we calculate its likelihood function with the LBTI observation and its corresponding χ². The MCMC then explores the parameter space to minimize the χ² value of the models. The model with the minimum χ² is finally returned as the best-fit model. Figure 7 shows this modeling procedure as applied to the 2023 May 25 observation. The different steps are detailed, from the computation of raw radiance to the retrieval of the excess signature on the detector. The various photometric apertures used to constrain the dust distribution are also shown, with the coherent FOV calculated from Eq. (1).

5.3 Results

First, assuming that the system is in a steady state, we use the three geometries considered to model the dust spatial distribution (thin ring, wide ring, and disk) to fit the combined 2017, 2018, and 2023 datasets. Figure 8 shows the best-fit models, and Table 6 summarizes the median values of the model parameters with their standard deviation, obtained from their MCMC posterior distributions. All MCMC posteriors are presented in Appendix B. The median values of χ² and the reduced χ² (${\chi }^2({\chi }_{\text{red}}^2={\chi }^2/\nu$), where ν is the degree of freedom) are also indicated.

The different geometries are then used to fit the individual 2017, 2018, and 2023 datasets. Figure 9 shows the best-fit models. The best-fit values of ${\chi }_{\text{red}}^2$ are significantly worse when fitting the combined nights compared to the individual nights. The steady-state hypothesis for the system is therefore not favored, but cannot be fully rejected. It is worth noting that, although the thin ring is a submodel included within the wide ring model, their results are not systematically compatible. Adding R_out as an additional free parameter allows exploration of more geometries and can yield models with better χ² - although not significantly better. For the individual fits, the χ² values obtained for each geometry are consistent across the three observing nights. Therefore, our data do not distinguish between the three possible dust distributions, assuming a face-on, optically thin exozodi.

Figure 10 presents the 2D maps of dust radiance for the individual models of 2017, 2018, and 2023, using the median values of the model parameters from Table 6.

5.4 Analysis of the dust models

Dust temperature. Equation (3) shows a degeneracy between the dust temperature and the dust density Σ(r). For the dust brightness detected by NOMIC, this degeneracy means that the dust density cannot be precisely estimated without constraining the dust temperature - and, by extension, the dust composition and grain size distribution. The values of Σ₀ in Table 6 therefore depend on the assumption that the dust grains behave as ideal blackbodies and should be interpreted with caution. Since the nulling observations performed with the LBTI are photometric, they cannot constrain the dust temperature using spectral information. Nonetheless, using the same temperatures T_d(r) between the 2017, 2018, and 2023 datasets allows for a comparison of their Σ values and an estimate of dust variability. Estimates of the dust size distribution, and hence its equilibrium temperature profile, are discussed in Sect. 6.

Dust variability at the EEID. As a probe for the amount of exo-zodiacal dust predicted by the models in the HZ of θ Boo, we consider the zodi level, z. This zodi level is defined as the dust surface density, Σ, at the EEID (~2au or 138 mas for θ Boo) divided by the dust surface density at 1 au in the Solar System (Weinberger et al. 2015). For a disk model, Σ can be calculated self-consistently for each combination of parameters Σ₀ and α using Eq. (2). We use the MCMC chain samples to build posterior distributions of z for the 2017, 2018, and 2023 disk models — z₂₀₁₇, z₂₀₁₈, and z₂₀₂₃, respectively. The resulting median values for the 2017, 2018, and 2023 datasets, together with their 1σ confidence intervals, are shown in Table 6. It should be noted that the zodi level obtained for the 2017 dataset is consistent with that calculated by the HOSTS survey (148.2±27.7 zodis).

To further quantify the increase in HZ dust brightness found in our models between 2017 and 2023, we calculate the posterior distribution of the ratio rz = (z₂₀₂₃ - z₂₀₁₇)/z₂₀₁₇. From the two MCMC runs for 2017 and 2023, we sample 800 models from each run and randomly pair MCMC models from the two epochs 640 000 times to construct the posterior histogram of the ratio. The posterior can then be integrated to determine posterior probabilities for the zodi level increase between epochs. For example, the integral over r_z > 2 denotes the probability that the zodi level at least doubled between 2023 and 2017. To perform this integration efficiently and analytically, we approximate the posterior using a Gaussian kernel density estimate (KDE: Rosenblatt 1956; Parzen 1962), which is shown together with the posterior histogram in Fig. 11. The resulting posterior probability that there is an increase in zodi level is 93.7% (~1.86σ).

The posterior probability thus suggests a tentative increase in the HZ dust brightness between 2017 and 2023. We note that the probability calculations still make implicit assumptions about the correctness and quality of both the data acquisition (including calibration and error model) and the dust modeling, and depend slightly on the details of the adopted kernel density estimate. These results should therefore be interpreted with caution. The zodi levels calculated for 2017, 2018, and 2023 are also useful for providing an order-of-magnitude estimate of the absolute dust brightness in the HZ of θ Boo. However, these values depend on the assumed dust temperature, which is not well constrained. Estimates of the zodi level for realistic dust grains with different sizes are given in Sect. 6.

Fig. 8 Best-fit models for the combined 2017, 2018, and 2023 datasets. Solid dark lines represent the best-fit models using a disk geometry, as defined in Sect. 5.1. Dashed gray lines represent the best-fit models using a wide ring geometry and dash-dotted silver lines represent the best-fit models using a thin ring geometry. Parameters of the models are listed in Table 6.

Fig. 9 Best fit models for the 2017 (left), 2018 (middle), and 2023 (right) datasets. Solid dark lines represent the best-fit models using a disk geometry, as defined in Sect. 5.1. Dashed gray lines represent the best-fit models using a wide ring, and dash-dotted silver lines represent the bestfit models using a thin ring geometry. The parameters for the models are listed in Table 6.

Fig. 10 2D maps of the dust surface brightness (Fdust) for the models of 2017 (left), 2018 (middle), and 2023 (right). Maps correspond to the disk (top), wide ring (middle), and thin ring (bottom) geometries. Fdust is normalized by the stellar radiance (F★). Both Fdust and F★ are integrated over the NOMIC passband. The dashed circle indicates the EEID (~138 mas). The maps have a pixel size of 10 mas.

Fig. 11 MCMC posterior distribution sample histogram of rz, representing the relative difference in zodi level between the 2017 and 2023 models. A Gaussian kernel density estimate (KDE) is used to fit the sample and calculate the posterior probabilities of the dust variability.

6 Discussion

6.1 Constraints from the models

SED of the dust models. The LBTI observations around 11 μm cannot discriminate between the different dust distribution models obtained in Sect. 5 or the properties of the dust grains. One way to further constrain their properties is to simulate the SED from each model and compare it with the upper limits from Fig. 1. Since the emitting dust may be considerably smaller than the wavelengths, we must account for the reduced emissivity Q_abs of realistic grains compared to ideal blackbodies. We assume compact spherical grains composed of amorphous silicates and refractory organics with volume fractions of 1/3 and 2/3, respectively (Li & Greenberg 1998). We then compute Q_abs and the grains’ resulting equilibrium temperatures using the Mie-type optical properties code of Sommer et al. (2025), which utilizes the optical constants for the constituent materials from Li & Greenberg (1997). Figure 12 shows the resulting Q_abs. We consider grains with diameters of 3, 5, and 10 μm. This choice is motivated by the assumption that the geometric optical depth is dominated by the smallest grains bound to the system, with an estimated blowout size of approximately 2.5 μm (calculated using the same grain model).

We retrieve the radial profile of the dust surface brightness S from the models shown in Fig. 10. For the three grain sizes considered, the radial profile of the dust’s optical depth, Σ, is recomputed to match its surface brightness with S, using the relation $$\mathrm{\Sigma }(r)=\frac{S(r)}{Q_{\text{abs}}(\lambda )\times B_{\nu }(\lambda ,T_{\text{d}}(r))}.$$(5)

Since the surface brightness S is constant between the models with ideal blackbodies (Sect. 5) and the models with realistic grains, all models should fit the data in a similar way. Finally, we choose to truncate the optical depth profile of the dust at 10 au (i.e., ~700 mas) for the SED calculation, since the dust located outside this distance is beyond our coherent FOV. This putative outer dust also has an infrared emission that does not comply with the upper limits, particularly for models with a slope α > 1. We argue that if such quantities of outer dust were present at these distances in a static way, they would have been detected by previous far-infrared observations.

Figure 13 shows the SEDs corresponding to the recomputed optical depth profiles. Only the SEDs for the disk and wide ring models are presented for comparison, as the SEDs for the thin ring models yield similar results. The estimate of the mean infrared excess at 11 μm from the current LBTI/NOMIC observations, calculated in Sect. 3.1 is also included. This data point serves as a consistency check, confirming that all models are compatible with our observations at 11 μm. No increase is observed at 11 μm between the 2017 and 2023 models. This outcome arises because the apparent increase in dust brightness is localized to the inner region of the system, while the dust emission integrated over the entire FOV remains constant for the three observations. This is evident in Fig. 4, where the data points corresponding to the largest apertures are consistent.

The SED results indicate that for a grain size of 3 μm, all models are consistent with the upper limits. For a grain size of 5 μm, some models show a discrepancy with the upper limits near a wavelength of 20 μm. For a grain size of 10 μm, all models exhibit a significant infrared excess between 20 μm and 50 μm. Therefore, the simulated SEDs do not clearly discriminate between the different dust distribution models obtained, but favor a dust grain size distribution between 3 and 5 μm in order to satisfy the various upper limits. Table 7 provides the updated estimates of zodi level from the recomputed optical depth profiles. The zodi levels change by a factor of 2.8 for dust grains between 3 and 10 μm.

Further constraints on the dust SED could be obtained using observations from the James Webb Space Telescope (JWST: Rigby et al. 2023; Gardner et al. 2023), particularly with the Medium Resolution Spectrometer (MRS: Wells et al. 2015; Argyriou et al. 2023) of the Mid-InfraRed Instrument (MIRI: Wright et al. 2023). The measured far-infrared emission could then be compared with our modeled SEDs to identify not only the dust distribution but also the dust properties, such as size, temperature, and composition.

Orientation of the system. The results presented in Sect. 3.1 support the hypothesis that the planetary system of θ Boo has a low inclination, possibly even a face-on orientation. This hypothesis is central for the modeling strategy described in Sect. 5; however, no other independent observations of the planetary system are currently available to further constrain this orientation. It is possible to constrain the stellar inclination from spectroscopic measurements of θ Boo. The inclination of debris disks is often found to match the stellar inclination (Watson et al. 2011; Guilloteau et al. 2011). To determine the stellar inclination, we follow the method used in Faramaz et al. (2014) for ζ² Reticuli. This method uses the color index (B - V), the stellar radius, R_★, and the activity indicator, R′_HK, to estimate the rotational period at the equator of the star using the activity-rotation diagram from Fig. 6b in Noyes et al. (1984). This period is then compared with measurements of ν sin(i), from which a value of i can be deduced. The stellar properties of θ Boo are summarized in Table 1. We estimate the rotational period to be P = 3.22 ± 0.83 days and the equatorial rotational velocity to be ν = 27.2 ± 7.1 km/s. This is consistent with observations giving ν sin(i) = 31.8 ± 2 km/s, which suggests that the star is likely viewed edge-on (i.e., with i ~ 90°). If the exozodiacal disk shares the same inclination as the star, this is incompatible with our LBTI observations for a vertically thin disk. To be compatible, the disk would need to be sufficiently thick so that no significant changes in null depths would have been observed across the different parallactic angles in Fig. 3. However, this hypothesis cannot currently be tested as no models are available for the LBTI to simulate observations of vertically thick disks.

The disk may also be misaligned with the stellar rotation axis, forming a polar circumstellar disk. This configuration has proven to be stable over the stellar lifetime and can arise from interactions in a young stellar cluster (Kennedy et al. 2012). The wide companion of θ Boo could potentially have caused such a misalignment. The disk may also be tilted or diffused by giant planets orbiting with a planet-disk misalignment (Pearce & Wyatt 2014; Brady et al. 2023). Additional observations are therefore required to further constrain the orientation of the planetary system around θ Boo.

Detection of giant planets. If the orbits of potential giant planets in the system share an inclination close to face-on with the exozodi, the constraints from radial velocity and transit observations become less stringent. Alternatively, the proper-motion anomaly measured by Gaia Early Data Release 3 (EDR3) could be analyzed to model what types of planets might generate such an anomaly (Kervella et al. 2022). Observations of θ Boo with Gaia are showing a renormalized unit weight error (RUWE) of 3.38, which is sufficient to indicate the possible presence of a massive body in the system, but this anomaly might also be explained by the wide companion of θ Boo.

Fig. 12 Simulated emissivities for dust grains of size 3, 5, and 10 μm. The volume fraction of the grains is 1/3 amorphous silicates and 2/3 refractory organics.

6.2 Origin of the dust and its variability

Poynting–Robertson (P-R) drag hypothesis. One of the main hypotheses to explain the presence of warm exozodiacal dust is the P-R drag mechanism. In this scenario, dust forms through planetesimal collisions in an outer belt located further out in the system and is dragged into the inner region by interaction with stellar radiation, eventually evaporating upon reaching sublimation temperature. To simulate this phenomenon, we applied the model presented in Rigley & Wyatt (2020) to the θ Boo system. In this analytical model, the stellar parameters and the radius of the outer belt serve as input parameters to predict the size distribution of particles in the disk. The size distribution is first simulated in the outer belt using the method described in Wyatt et al. (2011). The dust particles are then assumed to evolve independently inward until the system reaches a steady state where the optical depth distribution changes by less than 1% in a logarithmic time step. The smallest grains that dominate the optical depth typically reach that state after ~10Myr. The radial profile for each particle is determined using the approach of Wyatt (2005). The model then calculates the 2D optical depth distribution as a function of particle size and distance from the star. Finally, the model adjusts the width and mass of the outer belt to fit the surface brightness model from the observations. To determine the emission from the grains in the modeled disk, we again assume the optical properties used to simulate their SEDs in Sect. 6.1.

However, the simulations from this model are only applicable to a system in a steady state, on timescales of years. The time required for dust particles of blow-out size to drift from 8 au to 2au is estimated to be 35 000 yr (Wyatt & Whipple 1950; Burns et al. 1979). We therefore compared the P-R drag model with the model for the combined observing nights of the LBTI (Fig. 8) using a disk geometry. The wide ring and thin ring geometries are not compatible with P-R drag alone and require the presence of a massive object (e.g., a giant exoplanet) to carve their edges. Figure 14 shows the result of the P-R drag model, including the best-fit surface brightness profile and its corresponding SED. The result is obtained for an outer belt extending from 6 to 9 au, with a mass of 7 x 10⁻⁶ M_⊕, considering all dust particles with sizes ≤ 1 cm. We find that P-R drag can explain our model, with a discrepancy appearing only for <1 au separations. However, the outer belt required to explain this distribution produces an infrared excess that exceeds the upper limits in Fig. 1 around 25 μm. This suggests that such an outer belt should have been detected by Spitzer.

Since our models also do not favor the steady-state hypothesis (Sect. 5.3), we conclude that other phenomena are likely to influence the dust distribution and generate yearly variability.

Extreme debris disk hypothesis. The models presented in Sect. 5 indicate possible variability of the exozodiacal disk around θ Boo, both in its spatial distribution and brightness. Photometric variability of warm dust on monthly to yearly timescales has previously been reported for extreme debris disks (EDDs, Balog et al. 2009). These disks are detected around FGK-type main-sequence stars with ages between ~100 to ~200Myr (Melis et al. 2012; Meng et al. 2012, 2015; Rieke et al. 2021; Moór et al. 2022), and even > 1 Gyr (Thompson et al. 2019). This variability is commonly attributed to major collisions of planetesimals during the terrestrial planet formation phase (Su et al. 2019). Balog et al. (2009) also define EDDs as systems with total emission at least four times greater than the stellar photospheric output at 24 μm. For θ Boo, the age of the system is estimated to be 3-4.7 Gyr (Sect. 4), which is well after the terrestrial planet formation phase. Upper limits also indicate that no dust emission greater than one-tenth of the stellar photosphere has been detected at 24 μm. These findings cast doubt on the hypothesis of an EDD for θ Boo, with massive planetesimal collisions as the source of the dust and its variations.

Catastrophic event hypothesis. An alternative hypothesis to explain this variability is the occurrence of a catastrophic event such as a massive collision, fragmentation of a massive comet, or the passing of a massive object with highly eccentric orbit-that would have occurred in the inner region. After several years, such an event could generate new dust or disrupt the dust distribution, resulting in an observed brightness variability. For example, debris from a planetary collision occurring at 2 au would have a circular orbit of 2.5 years around θ Boo. After five years, the debris would begin to form a belt in the HZ that could be detected by the LBTI. Direct imaging observations in the H- and L-bands could not detect the presence of giant planets with masses lower than 11 M_Jup. An eccentric planet could, however, disrupt a potential outer belt and periodically scatter exocomets inward (Faramaz et al. 2017). From the increase in optical depth in our models, we estimate the additional mass of material needed to explain this variability. This material could originate from a massive collision or comet fragmentation. For a grain size of 3 μm, the required mass ranges from 10⁻³ M_Vesta (4 × 10⁻⁸ M_⊕) for our disk models, to 3 × 10⁻³ M_Vesta (1.2 × 10⁻⁷ M_⊕) for our wide ring models. For a grain size of 10 μm, the estimated mass is tripled.

Fig. 13 Simulated dust SEDs for realistic dust grains with sizes of 3 μm (left), 5 μm (middle), and 10 μm (right). The SEDs correspond to the disk (top) and wide ring (bottom) models for 2017, 2018, and 2023.

Fig. 14 Left: radial distribution of the dust surface brightness models for the combined 2017, 2018 and 2023 dataset. Models are calculated at a wavelength of 11 μm, with 1σ error bars derived from the MCMC posteriors. The solid line represents the predicted surface brightness from the P-R drag model for an outer belt extending from 6 au to 9 au, with a mass of 7 × 10−6 M⊕. The dust composition is considered to be 1/3 silicate and 2/3 carbon. Right: SED of the dust from the P-R drag model.

7 Summary

The HOSTS survey (Ertel et al. 2020a) demonstrated a correlation between the presence of cold outer dust and exozodiacal dust. θ Boo is one of the three exceptions where exozodiacal dust has been detected by the LBTI without prior observations of significant mid- to far-infrared excess emission. In this study, we analyzed multi-epoch observations of the exozodiacal disk around θ Boo using nulling interferometry with the LBTI in the N′-band. Three observations were obtained on UT 2017 April 11, UT 2018 May 23, and UT 2023 May 25, respectively. The data were reduced using the LBTI nulling pipeline and analyzed with a model dedicated to exozodiacal disk observations with the LBTI (Kennedy et al. 2015).

A prior assumption for our modeling procedure is that the planetary system of θ Boo is close to face-on (i.e., i ~ 0°). This is supported by our extended nulling observation in 2023, which shows no significant relation between the detected emission and the parallactic angle. However, calculations of the stellar inclination indicate that θ Boo is likely to be edge-on (i.e., i ~ 90°).

For the exozodiacal disk to also be edge-on, it would need to be sufficiently thick or diffuse. Because our modeling tool can only simulate a vertically thin disk, the face-on assumption is retained for this work. The results and conclusions presented here are therefore highly dependent on this assumption, and future observations are necessary to constrain the system’s inclination with a higher confidence.

From the models obtained for the individual and combined 2017, 2018, and 2023 datasets, we propose several spatial distributions of the dust that can explain our observations. The models for the combined nights fit the data less well than those for the individual nights, supporting the hypothesis of an unsteady state system. The different dust distribution geometries fit each night similarly; therefore, we are unable to discriminate between them with the LBTI observations. Using the upper limits on dust excess emission in the mid- to far-infrared from the Spitzer (2004) and Herschel (2010) observations, we further constrain the dust properties. Although our simulated SEDs do not clearly discriminate between the different geometries, they favor a representative dust grain size distribution of 3-5 μm. Using our models for the dust spatial distribution of 2017 and 2023, we also estimate the zodi level and its variability at the EEID. Considering our disk models, which assume a broad dust distribution in the orbital plane and that dust grains behave as ideal blackbodies, we estimate the dust brightness to be ${158}_{-67}^{+62}$ and ${443}_{-170}^{+98}$ zodis for 2017 and 2023, respectively. Considering the probability distribution of the two values, we find a tentative increase in dust emission with 1.86σ significance at the EEID. This result suggests possible brightness variability of exozodiacal disks in the HZ of main-sequence stars. Additional observations are needed to confirm this trend for θ Boo. However, the obtained zodi levels depend on the chosen dust temperature, and consequently on the size distribution of the dust grains. For realistic grains with a size of 3 μm, for example, the obtained zodi levels are ${63}_{-26}^{+24}$ and ${173}_{-67}^{+39}$ zodis for 2017 and 2023, respectively. Further constraints on the dust properties are therefore needed to estimate the absolute zodi level of the system.

Several hypotheses are considered for the origin of the exo-zodiacal dust and its variability: (1) P-R drag, (2) the EDD, and (3) a catastrophic event. The first hypothesis, P-R drag, is one of the main mechanisms explaining the presence of warm dust, continuously generated by a collisionally active asteroid belt analogue. However, nondetections of far-infrared excess limit the significance of such a belt. Using our steady-state model, with the combined 2017, 2018, and 2023 nights, we retrieve the dust distribution from the P-R drag for an outer parent belt extending from 6 au to 9 au with a mass of 7 × 10⁻⁶ M_⊕. However, since our steady-state model is not favored, another mechanism is likely to influence the dust distribution besides P-R drag, and generate a yearly variability.

The second hypothesis, considering an EDD, could explain the variability, but is unlikely, given the age of θ Boo and the level of excess infrared emission from the dust. The third hypothesis involves a catastrophic event, such as a massive collision, fragmentation of a massive comet, or passage of a massive object with a highly eccentric orbit. We estimate the amount of additional material injected into the system to explain this variability. Considering dust particles ranging from 3 μm to 10 μm, and for our different model geometries, we find an additional mass of 10⁻³−10⁻² M_Vesta (M_vesta ≃ 4 × 10⁻⁵ M_⊕). The presence of a substellar companion that could disrupt the dust distribution has also been constrained using high-contrast AO observations in the H- and L′-bands with the SHARK-NIR and LMIRCam instruments, respectively. These observations were performed on UT 2024 February 24, but did not show any features indicating the presence of a point-like companion source in the system. Nonetheless, the contrast limits of the images allowed us to constrain the presence of giant planets down to 11 M_Jup at 1.3″ angular separation, based on a “hot start” evolutionary model.

To overcome the current challenges in studying the exozo-diacal disk of θ Boo, additional observations will be required. In particular, detections at different wavelengths in the mid-infrared (e.g., using the JWST/MIRI instrument) would further constrain the dust properties and its temperature profile. More sensitive far-infrared observations would also search for emissions from a cold outer debris disk that could feed the warm dust, providing clues to the system’s orientation. Additional LBTI observations with reduced noise would also enable a more precise modeling of the distribution and orientation of the warm dust, particularly by considering the dependence of the null measurements on the parallactic angle. The implementation of the PCA method for background subtraction (Rousseau et al. 2024), along with the upcoming new detector for NOMIC, should reduce instrumental noise and improve future observations.

Acknowledgements

G.G., D.D., and M.A.M. have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG - 866070). S.E. is supported by the National Aeronautics and Space Administration through the Exoplanet Research Program (Grant No. 80NSSC21K0394) and the Astrophysics Decadal Survey Precursor Science program (Grant No. 80NSSC23K1473). V.F. acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grants No. 80NSSC21K0394 and 80NSSC23K1473 (PI: S. Ertel), and Grant No 80NSSC23K0288 (PI: V. Faramaz). T.D.P. acknowledges support of the Research Foundation - Flanders (FWO) under grant 11P6I24N (Aspirant Fellowship). The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation Members are: The University of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, The Leibniz Institute for Astrophysics Potsdam, and Heidelberg University; The Ohio State University, representing OSU, University of Notre Dame, University of Minnesota and University of Virginia. Part of the observing time for this program was granted as Director’s Discretionary Time by the LBT director. We acknowledge the use of the Large Binocular Telescope Interferometer (LBTI) and the support from the LBTI team, specifically from Amali Vaz, Jordan Stone, Jenny Power, Phil Willems, Grant West, Jared Carlson, Andrew Cardwell, and Alexander Becker. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Multilateral Agreement.

Appendix A Null depths

Fig. A.1 Nulls with respect to UT (left) and parallactic angle (right) for 8 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.

Fig. A.2 Nulls with respect to UT (left) and parallactic angle (right) for 16 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.

Fig. A.3 Nulls with respect to UT (left) and parallactic angle (right) for 24 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.

Fig. A.4 Nulls with respect to UT (left) and parallactic angle (right) for 32 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.

Fig. A.5 Nulls with respect to UT (left) and parallactic angle (right) for 40 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.

Appendix B Corner plots

Fig. B.1 Corner plots showing the MCMC posteriors for the model fitting of the 2017 (left) and 2018 (right) datasets. The posteriors are given for a “disk” (top), “wide ring” (middle), and “thin ring” (bottom) geometry.

Fig. B.2 Corner plots showing the MCMC posteriors for the model fitting of the 2023 (left) and the combined 2017/2018/2023 (right) datasets. The posteriors are given for a “disk” (top), “wide ring” (middle), and “thin ring” (bottom) geometry.

References

All Tables

All Figures

	Fig. 1 Observed spectral energy distribution (SED) of θ Boo. The black points indicate the observed photometric data. The gray line represents the best fit stellar model, and the green dashed line shows the dust model. The green points indicate the upper limits of the dust emission.
In the text

Fig. 2 Block diagram of the LBTI architecture. Light from the two apertures enters from the top left and right and passes through the beam combiner subsystem (purple) and the NIC cryostat (blue). Both the beam combiner and the NIC are maintained at cryo-temperatures. At the entrance of the beam combiner subsystem, visible light is reflected by dichroics and directed to the wavefront sensing (WFS) units. Infrared light is transmitted to the beam combiner where two mirrors adjust the tip-tilt and the optical path difference (OPD) of the beams for nulling interferometry. After recombination of the two LBTI pupils by the roof mirror, the light is split in the NIC. The L-band (3-5 μm) is directed to LMIRCam for high-contrast imaging of exoplanets, whereas the K-band (1.5-2.5 μm) is directed to Phasecam to monitor the differential tip-tilt and phase between the two apertures. Finally, the N-band (8-13 μm) light is directed to NOMIC which can be used in either imaging interferometry or nulling interferometry modes.

In the text

Fig. 3 Null depths obtained after the data reduction pipeline described in Sect. 3.1. These points are computed using the reduction pipeline with an aperture radius of 24 pixels (i.e., 430 mas). The OBs are plotted as a function of their parallactic angle (left) and UT of observation (right), for each of the three nights: 2017 April 11 (top), 2018 May 23 (middle), and 2023 May 25 (bottom). In each panel, light blue points and light red crosses show the measured null depth of the OBs for the calibrators and for θ Boo, respectively. Dark blue points and dark red crosses indicate the null depth values for each pointing, corresponding to the unweighted mean values of the OBs included in the pointing. The instrumental null depths measured from the calibrators are averaged to calculate the null floor, which represents the reference level below which any excess cannot be reliably detected. The null floor and its error bar are shown by the solid and dotted black lines, respectively. The error bar corresponds to the 1σ standard deviation.

In the text

Fig. 4 Left: calibrated null depths for the different aperture radii and for the three nights of observation with NOMIC. The aperture radii are 8, 16, 24, 32, and 40 pixels (i.e., 143, 286, 430, 573, and 716 mas, respectively). Right: calibrated null depths for the 2023 May 25 observation as a function of aperture radius. Results for the two individual pointings are shown as dotted and dashed lines, and the combined results are shown as crosses. The first pointing corresponds to the OBs taken from 04:00 to 06:00 UT and the second pointing is from 06:30 to 07:30 UT.

In the text

Fig. 5 Left: ADI image of θ Boo obtained with LMIRCam after processing with the LEECH-survey pipeline. The image has a FOV of 4″. A central mask of 5λ/D radius is used to improve the visualization. No feature indicating the presence of a companion is detected between 0.3″ and 4". Right: ADI image of θ Boo from SHARK-NIR after processing with the SHARP pipeline using five PCA components. The image has a FOV of 1.3″. For both images, the outer circle indicates the coherent FOV of NOMIC in the N′-band (~670 mas; see Eq. (1)), and the inner dashed circle represents the EEID (~ 138 mas). No feature indicating the presence of a companion is detected between 0.15″ and 1.3″.

In the text

Fig. 6 Detection limits for companions around θ Boo in the L′-band with LMIRCam (left) and in the H-band with SHARK-NIR (right). The blue line shows the contrast as a function of angular separation. For four contrast levels, the corresponding detection limit is given in terms of mass (in units of Jupiter mass). The best mass constraint is obtained in the H-band around 1″, with a limit of 11 MJup. The system is assumed to be 500 Myr for these mass constraints.

In the text

Fig. 7 Results of the 2023 May 25 night modeling, showing the different steps of the procedure. The geometry used is a disk with Σ0 = 3 × 105 and α = 0, the orientation is face-on. From (a) to (d) show: the angular distribution of Fdust/F★, the impact of the LBTI transmission map, the convolution with the PSF of the LBT, and the comparison between the observed excess and the modeled excess. Panel (c) indicates the photometric apertures used in the reduction with the dashed circles. The coherent FOV is indicated with the solid orange circle. The maps have a pixel size of 10 mas.

In the text

	Fig. 8 Best-fit models for the combined 2017, 2018, and 2023 datasets. Solid dark lines represent the best-fit models using a disk geometry, as defined in Sect. 5.1. Dashed gray lines represent the best-fit models using a wide ring geometry and dash-dotted silver lines represent the best-fit models using a thin ring geometry. Parameters of the models are listed in Table 6.
In the text

	Fig. 9 Best fit models for the 2017 (left), 2018 (middle), and 2023 (right) datasets. Solid dark lines represent the best-fit models using a disk geometry, as defined in Sect. 5.1. Dashed gray lines represent the best-fit models using a wide ring, and dash-dotted silver lines represent the bestfit models using a thin ring geometry. The parameters for the models are listed in Table 6.
In the text

Fig. 10 2D maps of the dust surface brightness (Fdust) for the models of 2017 (left), 2018 (middle), and 2023 (right). Maps correspond to the disk (top), wide ring (middle), and thin ring (bottom) geometries. Fdust is normalized by the stellar radiance (F★). Both Fdust and F★ are integrated over the NOMIC passband. The dashed circle indicates the EEID (~138 mas). The maps have a pixel size of 10 mas.

In the text

	Fig. 11 MCMC posterior distribution sample histogram of rz, representing the relative difference in zodi level between the 2017 and 2023 models. A Gaussian kernel density estimate (KDE) is used to fit the sample and calculate the posterior probabilities of the dust variability.
In the text

	Fig. 12 Simulated emissivities for dust grains of size 3, 5, and 10 μm. The volume fraction of the grains is 1/3 amorphous silicates and 2/3 refractory organics.
In the text

	Fig. 13 Simulated dust SEDs for realistic dust grains with sizes of 3 μm (left), 5 μm (middle), and 10 μm (right). The SEDs correspond to the disk (top) and wide ring (bottom) models for 2017, 2018, and 2023.
In the text

Fig. 14 Left: radial distribution of the dust surface brightness models for the combined 2017, 2018 and 2023 dataset. Models are calculated at a wavelength of 11 μm, with 1σ error bars derived from the MCMC posteriors. The solid line represents the predicted surface brightness from the P-R drag model for an outer belt extending from 6 au to 9 au, with a mass of 7 × 10−6 M⊕. The dust composition is considered to be 1/3 silicate and 2/3 carbon. Right: SED of the dust from the P-R drag model.

In the text

	Fig. A.1 Nulls with respect to UT (left) and parallactic angle (right) for 8 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.
In the text

	Fig. A.2 Nulls with respect to UT (left) and parallactic angle (right) for 16 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.
In the text

	Fig. A.3 Nulls with respect to UT (left) and parallactic angle (right) for 24 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.
In the text

	Fig. A.4 Nulls with respect to UT (left) and parallactic angle (right) for 32 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.
In the text

	Fig. A.5 Nulls with respect to UT (left) and parallactic angle (right) for 40 pixels aperture radius. (top) 2017 April 11, (middle) 2018 May 23, (bottom) 2023 May 25.
In the text

	Fig. B.1 Corner plots showing the MCMC posteriors for the model fitting of the 2017 (left) and 2018 (right) datasets. The posteriors are given for a “disk” (top), “wide ring” (middle), and “thin ring” (bottom) geometry.
In the text

	Fig. B.2 Corner plots showing the MCMC posteriors for the model fitting of the 2023 (left) and the combined 2017/2018/2023 (right) datasets. The posteriors are given for a “disk” (top), “wide ring” (middle), and “thin ring” (bottom) geometry.
In the text