# Scattering and sublimation: a multiscale view of µm-sized dust in the inclined disc of HD 145718

Claire L. Davies, Evan A. Rich, Tim J. Harries, John D. Monnier, Anna S. E. Laws, Sean M. Andrews, Jaehan Bae, David J. Wilner, Narsireddy Anugu, Jacob Ennis, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Claire L. Davies, Evan A. Rich, Tim J. Harries, John D. Monnier, Anna S. E. Laws, et al., Scattering and sublimation: a multiscale view of µm-sized dust in the inclined disc of HD 145718. Monthly Notices of the Royal Astronomical Society, 2022, 511, pp.2434-2452. 10.1093/mnras/stac149. insu-03705355

HAL Id: insu-03705355

https://insu.hal.science/insu-03705355

Submitted on 24 Mar 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Scattering and sublimation: a multiscale view of µm-sized dust in the inclined disc of HD 145718

Claire L. Davies<sup>1,★</sup> Evan A. Rich<sup>1,†</sup>,<sup>2</sup> Tim J. Harries,<sup>1</sup> John D. Monnier<sup>1,‡</sup>,<sup>2</sup> Anna S. E. Laws,<sup>1</sup>

Sean M. Andrews,<sup>3</sup> Jaehan Bae<sup>1,4,5</sup> David J. Wilner,<sup>3</sup> Narsireddy Anugu,<sup>6,7</sup> Jacob Ennis,<sup>2</sup>

Tyler Gardner,<sup>2</sup> Stefan Kraus,<sup>1</sup> Aaron Labdon,<sup>1,8</sup> Jean-Baptiste le Bouquin,<sup>2,9</sup> Cyprien Lanthermann,<sup>7</sup>

Gail H. Schaefer,<sup>7</sup> Benjamin R. Setterholm,<sup>2</sup> Theo ten Brummelaar<sup>7</sup> and the G-LIGHTS collaboration

Accepted 2022 January 17. Received 2022 January 17; in original form 2021 August 27

## ABSTRACT

We present multi-instrument observations of the disc around the Herbig Ae star, HD 145718, employing geometric and Monte Carlo radiative transfer models to explore the disc orientation, the vertical and radial extent of the near-infrared (NIR) scattering surface, and the properties of the dust in the disc surface and sublimation rim. The disc appears inclined at 67–71°, with position angle, PA = −1.0 to 0.6°, consistent with previous estimates. The NIR scattering surface extends out to ~75 au and we infer an aspect ratio, $h_{scat}(r)/r \sim 0.24$ in J band; ~0.22 in H band. Our Gemini Planet Imager images and VLTI + CHARA NIR interferometry suggest that the disc surface layers are populated by grains $\gtrsim\lambda/2\pi$ in size, indicating these grains are aerodynamically supported against settling and/or the density of smaller grains is relatively low. We demonstrate that our geometric analysis provides a reasonable assessment of the height of the NIR scattering surface at the outer edge of the disc and, if the inclination can be independently constrained, has the potential to probe the flaring exponent of the scattering surface in similarly inclined ($i \gtrsim 70°$) discs. In re-evaluating HD 145718’s stellar properties, we found that the object’s dimming events – previously characterized as UX Or and dipper variability – are consistent with dust occultation by grains larger, on average, than found in the ISM. This occulting dust likely originates close to the inferred dust sublimation radius at 0.17 au.

**Key words:** accretion, accretion discs – radiative transfer – techniques: high angular resolution – circumstellar matter – stars: formation – stars: individual: HD 145718.

## 1 INTRODUCTION

The planet formation process requires sub-µm sized particles, typical of the interstellar medium (ISM), to grow by >12 orders of magnitude to produce planetesimals and planets. Moreover, given the relatively short lifetimes of protoplanetary discs (~3–11 Myr; Ribas, Bouy & Merín 2015), such growth has to be highly efficient. The process is complicated – with grain evolution involving coagulation, vertical settling, radial drift, and fragmentation – and dependent on the local structure of the disc and the properties of the dust therein (Testi et al. 2014).

Highly inclined ($i \gtrsim 70°$), dust rich protoplanetary discs provide unique opportunities to study vertical and radial disc structure. In

particular, the optically thick nature of protoplanetary discs across the optical/near-infrared (NIR) allows observations at these wavelengths to directly trace the disc surface layers. The dust in these regions is expected to be dominated by sub-µm sized grains and be well-coupled to the gas. Meanwhile, large grains are expected to preferentially settle towards the disc mid-plane (Dubrulle, Morfill & Sterzik 1995).

Observational evidence for such vertical stratification is seen in the wavelength dependence of the vertical extent of near-edge-on discs (e.g. Glauser et al. 2008; Duchêne et al. 2010; Villenave et al. 2019; Wolff et al. 2021). However in apparent contrast, the presence of ‘large’ aggregates (radius, $a \gtrsim \lambda/2\pi$, where $\lambda$ is the observed wavelength) in the surface layers of protoplanetary discs has been inferred from (i) asymmetric brightness patterns in polarized differential imaging data (e.g. Mulders et al. 2013; Stolker et al. 2016a, b; Avenhaus et al. 2018; Garufi et al. 2020)

<footnote n="1"><sup>1</sup>Astrophysics Group, Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="2"><sup>2</sup>Astronomy Department, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="3"><sup>3</sup>Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="4"><sup>4</sup>Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="6"><sup>6</sup>Steward Observatory, Department of Astronomy, University of Arizona, 933 North Cherry Avenue, Tucson, AZ, 85721, USA</footnote>

<footnote n="7"><sup>7</sup>The CHARA Array of Georgia State University, Mount Wilson Observatory, Mount Wilson, CA 91023, USA</footnote>

<footnote n="8"><sup>8</sup>European Southern Observatory, Casilla 19001, Santiago 19, Chile</footnote>

<footnote n="9"><sup>9</sup>Institut de Planetologie et d'Astrophysique de Grenoble, F-38058 Grenoble, France</footnote>

<footnote>* E-mail: c.davies3@exeter.ac.uk</footnote>
and (ii) the colours and polarization of stars exhibiting photometric variability attributable to occultation by circumstellar dust (namely 'dippers', Bouvier et al. 1999; Stauffer et al. 2015; Bredall et al. 2020 and UX Ors, Huang et al. 2019). Numerical simulations and lab experiments show that such large, porous dust grains may be key to overcoming the bouncing barrier (Wada et al. 2011; Kothe et al. 2013; Brisset et al. 2017) and the radial drift barrier (Okuzumi et al. 2012; Kataoka et al. 2013) to dust grain growth.

Bright, relatively close-by ($d \lesssim 300$ pc) young stellar objects (YSOs) permit the detailed study of disc structure using 8-m class telescopes and infrared (IR) and millimetre (mm) interferometers. Here, we focus on HD 145718 (common aliases include PDS 80 and V718 Sco), an intermediate mass YSO ($\sim 1.5-3 M_\odot$) - i.e. a Herbig Ae star - in the Upper Scorpius (USco) association (Rizzuto, Ireland & Robertson 2011; Pecaut, Mamajek & Bubar 2012; Luhman et al. 2018). In particular, we combine new and archival observations of the dusty circumstellar environment around HD 145718 probing sub- to tens of au scales. We combine these observations to constrain the nature of the dust grains in the surface layers of the innermost and outermost disc regions and assess whether circumstellar dust obscuration is likely responsible for the photometric variability observed in this object.

Our paper is organized as follows. Section 2 provides an overview of previous studies involving HD 145718 and our knowledge of its circumstellar environment to-date. Section 3.1 describes our Gemini Planet Imager (GPI) observations of HD 145718, conducted as part of the Gemini Large Imaging with GPI Herbig/T-Tauri Survey (G-LIGHTS; Monnier et al. 2017; 2019; Laws et al. 2020; Rich et al., in preparation). Sections 3.2 and 3.3 describe our complementary CHARA/MIRC-X NIR interferometric observations, VLTI/PIONIER and VLTI/GRAVITY archival NIR interferometric data sets, and archival multiband, multi-epoch photometry, and IR spectroscopy. Our combination of NIR interferometric and polarized scattered light imaging probes sub-µm- to µm-sized dust grains on multiple angular scales. We first employ a simple geometric model to explore the disc orientation and extent. The methodology and results of this part of our investigation are presented in Section 4. We further build on this in Section 5 using full Monte Carlo radiative transfer models to simultaneously model the GPI images, NIR interferometry, and spectral energy distribution (SED). This includes a re-evaluation of HD 145718's stellar luminosity, radius, and visual extinction in Section 5.1. We discuss our results in the context of HD 145718's photometric variability in Section 6.1, assess the robustness of our geometric modelling in Section 6.2, and summarize our findings in Section 7.

## 2 HD 145718

Previous studies of HD 145718 have reported the existence of an inclined ($i \gtrsim 70^\circ$; Guimarães et al. 2006; Garufi et al. 2018; Gravity Collaboration 2019; Ansdell et al. 2020), gas-rich (Dent, Greaves & Coulson 2005; Ansdell et al. 2020) and dust-rich (e.g. Gregorio-Hetem et al. 1992; Oudmaijer et al. 1992; Friedemann, Guertler & Loewe 1996) Meeus et al. (2001) Group II disc (Keller et al. 2008). Dent et al. (2005) obtained a marginal $J = 3 - 2^{12}\text{CO}$ detection towards HD 145718 and (accounting for the different stellar distances adopted between their study and ours – see Section 5.1) estimated the gaseous disc extends out to $70 \pm 35$ au.

The star itself has a spectral type of A5 (Carmona et al. 2010). Its identification as photometrically variable saw it classified as an eclipsing binary throughout the 1900s. However, by comparing their radial velocities to earlier measurements by Carmona et al. (2010), Ripepi et al. (2015) found no evidence of multiplicity

in the system. Adaptive optics imaging has also ruled out the presence of companions at 20–200 milliarcsecond (mas) separations down to $\Delta L' = 2.6-4.8$ mag (Ansdell et al. 2020). More recently, HD 145718's photometric variability has been re-attributed to inherent stellar variability ($\delta$ Scuti-type pulsations; Ripepi et al. 2015) and obscuration by circumstellar dust (dipper and UX Ori variability; Poxon 2015; Ansdell et al. 2018; Cody & Hillenbrand 2018; Rebull et al. 2018). Temporal variations like those seen in the blueshifted portion of HD 145718's $H\alpha$ line profile (compare, for example, the line profiles presented in Vieira et al. 2003; Carmona et al. 2010; Ripepi et al. 2015; Wichittanakom et al. 2020) and higher order Balmer series lines (Guimarães et al. 2006) can also be attributed to accretion and outflow processes in YSOs (e.g. Muzerolle et al. 2004) and are likely also linked to the dipper/UX Ori variability.

The evolutionary status of HD 145718 has been debated in the literature and isochronal age estimates for the object vary from 5.7 to 20 Myr (Alecian et al. 2013; Fairlamb et al. 2015; Vioque et al. 2018; Arun et al. 2019; Wichittanakom et al. 2020). However, estimating the age of individual disc-hosting, photometrically variable young stars using isochrone fitting is fraught with difficulty (e.g. Davies, Gregory & Greaves 2014). Strong evidence for HD 145718 being pre-main-sequence rather than an evolved star is found in the p-mode frequencies of the object's $\delta$ Scuti-type pulsations: the highest p-mode frequency observed – which scales linearly with stellar age (Zwintz et al. 2014) – is consistent with those of disc-less, non-accreting USco members with isochronal ages of $\sim 10$ Myr (Ripepi et al. 2015).

## 3 OBSERVATIONS AND SUPPLEMENTARY ARCHIVAL DATA

### 3.1 GPI data

J- and H-band polarimetry mode observations of HD 145718 with GPI (Macintosh et al. 2014), situated at the Gemini South telescope, were obtained on 2018 June 07 and 08, respectively (program ID GS-2018A-LP-12). The 32 frames, each with two co-adds, were observed with exposure time of 29 s (total integration time = 1862 s per wave-band). Between each frame, the half-wave plate was rotated between $0^\circ, 22.5^\circ, 45^\circ$, and $67.5^\circ$, creating eight independent sequences with four equally spaced half-wave plate angles. Additionally, J- and H-band coronagraphs were used with focal plane diameters of 0.184 and 0.246 arcsec, respectively.

To reduce the data, we used the GPI Data Reduction Pipeline (DRP) version 1.5 (Maire et al. 2010; Perrin et al. 2014), with modifications to the flux calibration of the polarized images and the removal of stellar and instrumental polarization (see Monnier et al. 2019; Laws et al. 2020; Rich et al., in preparation). In summary, the GPI DRP was used to subtract dark background, extract the polarization spots, correct for bad pixels, remove microphonics noise, flat-field the frames using a low-frequency flat, and measure the star position using a radon transformation of the satellite spots.<fnref n="1" /> This resulted in 32 polarized images: four polarized images within each of the eight polarization sets. Each of the four polarized images were combined together using the double-differencing technique, creating eight Stokes cubes containing *I*, *Q*, and *U* images. The polarization sets were rotated so that the top of the image pointed north and the stellar and instrumental polarization were removed (Appendix B).

<footnote n="1">¹These so-called satellite spots are created by diffractive elements in the pupil plane of the GPI instrument. Their radial structure point to the star's location behind the coronographic mask.</footnote>
We then projected the eight polarization sets of Q and U (oriented with North up; East left) to local $Q_\phi$ and $U_\phi$, based on the stellar position determined above. Specifically, for a pixel grid of (X, Y) coordinates with centre pixel ($X_0, Y_0$) and coordinate position angle, $\gamma$:

$$Q_{\phi} = -U \sin(2\phi) - Q \cos(2\phi) \quad (1)$$

$$U_{\phi} = Q \sin(2\phi) - U \cos(2\phi), \quad (2)$$

where $\phi$ is the polar angle

$$\phi = \tan^{-1} \left( \frac{Y - Y_0}{X - X_0} \right) + \gamma \quad (3)$$

(see appendix A of Monnier et al. 2019). The eight polarization sets were then averaged together to create the combined $I$, $Q_\phi$, and $U_\phi$ images.

The satellite spots in the polarized images were used in the flux calibration. The 32 polarized images were averaged together to increase the signal-to-noise ratio of the satellite spots, as described in Laws et al. (2020). We used only the second-order and only the first-order satellite spots for the J- and H-band observations, respectively. Using 2MASS photometry (Skrutskie et al. 2006; Table A1), we measured J- and H-band flux scaling factors of $3.387 \pm 0.75$ mJy arcsec$^{-1}$/(ADU/sec/coadd) and $2.492 \pm 0.37$ mJy arcsec$^{-1}$/(ADU/sec/coadd), respectively.

### 3.1.1 Visual inspection of flux-calibrated images

Our flux-calibrated Stokes $I$, $U_\phi$, and $Q_\phi$ images of HD 145718 are shown in the top, middle, and bottom rows of Fig. 1, respectively. Fig. 2 highlights the main features in the $Q_\phi$ images. The black circle at the centre of each image indicates the size of the inner working angle (IWA) of the coronographic mask (radii of ~9 and ~11 pixels for the J- and H-band images, respectively; Rich et al., in preparation).

Both the J- and H-band $Q_\phi$ images feature a central ellipse (marked $E_J$ and $E_H$, respectively, in Fig. 2), elongated along a north-south direction. The $U_\phi$ images feature positive flux to the north-west and south-east and negative flux in the north-east and south-west. The brightest regions (negative and positive flux) predominantly extend along a north-south direction with minimal extension to the east and west of the IWA. Taken together, these $Q_\phi$ and $U_\phi$ features are consistent with the presence of an inclined circumstellar disc around HD 145718, with a major axis position angle, $\mathrm{PA} \approx 0^\circ$. This is consistent with previous results from K band and mm continuum interferometry: $\mathrm{PA} = 2 \pm 2^\circ$ (Gravity Collaboration 2019) and $\mathrm{PA} = 1 \pm 1^\circ$ (Ansdell et al. 2020), respectively.

The major axes of $E_J$ and $E_H$ appear offset to the west of the image centre. While we cannot directly measure how well the coronograph is centred on the star in HD 145718, we can use other targets in our G-LIGHTS sample with companion detections within the field of view of GPI to comment on the significance of this apparent offset. For example, in the eight cycles we observed for G-LIGHTS target HD 50138, we estimate a standard deviation for the centroid of its companion to be ~0.22 pixels (~3.1 mas), much smaller than the offset we observe. Thus, we do not expect the coronograph centring to be a large contributing factor to the offset we observe for HD 145718. Instead, the offset of the ellipse centre from the image centre is likely symptomatic of a flared disc structure where the ellipse traces the open face of the disc towards the observer.

The H-band $Q_\phi$ image features a bright arc to the east of $E_H$ (marked $A_H$ in Fig. 1). This likely traces scattering events close

to the outer edge of the surface of the disc facing away from the observer. The drop in surface brightness between $E_H$ and $A_H$ is likely a result of the opaque disc mid-plane. The J-band $Q_\phi$ image lacks a similar arc feature. Instead, two dark features ($Q_\phi$ flux significantly below the background level) are observed immediately to the east of the coronographic mask (marked $D_{J1}$ and $D_{J2}$ in Fig. 1). These may also be attributable to the opaque disc mid-plane or may be artefacts of imperfect stellar polarization subtraction (Appendix B).

## 3.2 Complementary near-infrared interferometry

Fully reduced and calibrated VLTI/PIONIER data were retrieved from the Optical Interferometry Database (OIDB). These probe H-band emission from HD 145718 on mas scales. Details of the reduction and calibration procedure are provided in Lazareff et al. (2017). K-band NIR interferometric data, obtained using VLTI/GRAVITY and originally published in Gravity Collaboration (2019), were retrieved from the European Southern Observatory archive. The data were reduced and calibrated using the GRAVITY pipeline (version 1.1.2) using default settings. We restrict our analysis to the low spectral dispersion ($R = \Delta\lambda/\lambda \sim 30$) GRAVITY fringe tracker (FT) data. The standard star HD 145809 (uniform disc diameter, UDD = $0.402 \pm 0.002$ mas; Bourges et al. 2017) was used to estimate the transfer function and calibrate the visibilities and closure phases. The bluest spectral channel of the GRAVITY FT was not used as this is known to be corrupted by the metrology laser operating at 1.908 µm (Lippa et al. 2016).

To probe smaller scale circumstellar emission, we also obtained a single snapshot observation of HD 145718 using the MIRC-X instrument (Anugu et al. 2018, 2020; Kraus et al. 2018) of the Center for High Angular Resolution Astronomy (CHARA) Array on UT date 2021 May 12. The (u, v)-plane coverage of our observations are compared to those of the GRAVITY and PIONIER interferometric data sets in Fig. 3. The CHARA Array comprises six 1-m class telescopes arranged in a Y-shaped array. Its operational baselines between 34 and 330 m (ten Brummelaar et al. 2005) provide ~0.5 mas resolution<fnref n="2" /> across the H band. MIRC-X is capable of combining light from all six CHARA telescopes. However, fibre injection issues on one beam associated with a telescope focus problem caused by the mounting mechanism of the primary mirror limited our observations to the five-telescope configuration E1-W2-W1-S2-E2. We used the PRISM 50 spectral mode that provides five spectral channels across the H band ($\lambda \sim 1.4 - 1.7$ µm). The data were reduced using pipeline version 1.3.5,<fnref n="3" /> described in Anugu et al. (2020). We applied the bispectrum bias correction, set the number of coherent coadd to 10, and adopted a flux threshold of 10. Otherwise, we adopted default reduction settings. Standard stars HD 145965 (UDD = $0.209 \pm 0.005$ mas; Bourges et al. 2017) and HD 139487 (UDD = $0.305 \pm 0.009$ mas; Bourges et al. 2017) were observed either side of HD 145718 in a CAL-SCI-CAL concatenation. These data were inspected for signatures of binarity but none were found. They were used to estimate the transfer function to calibrate the visibilities and closure phases.

### 3.2.1 Inspection of the NIR interferometry

The CHARA/MIRC-X interferometry was inspected for (i) consistency with the VLTI/PIONIER data and (ii) signatures of binarity.

<footnote n="2">²$\approx \lambda/2B$, with B denoting the baseline length and $\lambda$ the operational wavelength.</footnote>

<footnote n="3">³https://gitlab.chara.gsu.edu/lebouquj/mircx_pipeine.git</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** GPI images: Total intensity, *I* (top row), radial Stokes *U*<sub>φ</sub> (middle row), and *Q*<sub>φ</sub> (bottom row). The left-hand (right-hand) column shows the *J*-band (*H*-band) images (north is up; east is left). The physical scale is shown in the bottom right corner of all the images. The IWA of the coronographic mask (radii of ~9 and ~11 pixels for the *J*- and *H*-band, respectively) is indicated by a black ring in the centre of each image.</figcaption></figure>

The shorter baseline MIRC-X visibilities showed good consistency with those obtained by PIONIER. However, the different angular scales and position angles probed by the two data sets – together with the variations in brightness that HD 145718 is known to exhibit – make direct comparison of the data difficult. The MIRC-X closure

phases are consistent with zero, indicating that the underlying bright-
ness distribution is centro-symmetric. The VLTI data were known
to show a similar lack of deviation from centro-symmetry (Lazareff
et al. 2017; Gravity Collaboration 2019). We thus find no indication
of binarity and restrict our analysis in Section 5 to the visibilities.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Annotated J-band (top) and H-band (bottom) Stokes Q<sub>φ</sub> images highlighting the main features E<sub>J</sub>, D<sub>J1</sub>, D<sub>J2</sub>, E<sub>H</sub>, and A<sub>H</sub> (see Section 3.1.1 for details). The IWA of the coronographic mask is indicated by the black ring in the centre of the image.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** (u, v)-plane coverage of the interferometry. Left: H-band VLTI/PIONIER (dark green data points) and CHARA/MIRC-X (lime green data points). Right: K-band VLTI/GRAVITY (brown data points). North is up; east is left.</figcaption></figure>

## 3.3 Archival multiwavelength photometry and infrared spectroscopy

Complementary multiwavelength archival photometry and flux measurements were retrieved using the Spectral Energy Distribution Builder for Young Stars (SEDBYS; Davies 2021). The WISE W4-band magnitude was flagged and removed due to its discrepant low flux and we added the 1.3 mm ALMA flux reported in Garufi et al. (2018)

to the collated data set. The full list of flux-calibrated photometry, together with their references, are provided in Appendix A. Flux-calibrated Spitzer Infrared Spectrograph (IRS; Houck et al. 2004) Short-Low and Long-High module spectra were also retrieved from the Infrared Science Archive (IRSA). Details regarding the reduction of these data are provided in Keller et al. (2008).

## 4 GEOMETRIC MODELLING OF THE GPI IMAGES

We fit elliptical ring models to isophotes of surface brightness, S<sub>v</sub>, tracing the E<sub>J</sub> and E<sub>H</sub> features in the Q<sub>φ</sub> images. These allow us to assess the radial and vertical extents, as well as the orientation, of the disc around HD 145718 prior to the more computationally expensive and time-consuming radiative transfer modelling (Section 5). Our elliptical ring model was prescribed as a circular ring of radius, r, inclined by angle, i, rotated through position angle, PA,<sup>4</sup> and translated in right ascension, RA, and declination, Dec., by coordinates (δRA, δDec.) from the image centre. Assuming our observations trace light scattered by dust close to the disc surface, we expect the centres of E<sub>J</sub> and E<sub>H</sub> to be offset from the image centre along a vector which, when projected on to the sky, lies perpendicular to the disc PA. Thus, we can relate PA, δRA, and δDec. to the height, h<sub>scat</sub>(r), of the scattering surface above the disc mid-plane at radius r:

$$PA = \tan^{-1} \left( \frac{\delta \text{RA}}{\delta \text{Dec}} \right) + \frac{\pi}{2}, \quad (4)$$

and

$$h_{\text{scat}}(r) = d((\delta \text{RA})^2 + (\delta \text{Dec})^2)^{1/2}, \quad (5)$$

where d is the stellar distance. We follow Vioque et al. (2018) and adopt d = 152.5<sup>+3.2</sup><sub>-3.0</sub> pc, corresponding to the inverse of the Gaia data release (DR) 2 parallax<sup>5</sup> (Gaia Collaboration 2016, 2018). Our elliptical ring model is then fully prescribed using four parameters: r, i, PA, and h<sub>scat</sub>(r).

Before extracting the S<sub>v</sub> isophotes, we first masked the Q<sub>φ</sub> images to exclude the central pixels within the IWA of the coronograph. The pixel coordinates of the S<sub>v</sub> isophotes tracing the E<sub>J</sub> and E<sub>H</sub> features were then isolated from the full list returned by the contour function of MATPLOTLIB.PYPLOT (Hunter 2007). The western side of each ellipse-tracing S<sub>v</sub> isophote deviated from an elliptical shape, likely due to the combined effects of (i) the relatively narrow vertical extent of the scattering surface compared to the east-west extent of the coronographic mask; and (ii) the scattering phase function of the dust grains in the disc resulting in a lower back-scattered than forward-scattered flux (e.g. Stolker et al. 2016b; Tazaki et al. 2019, see Section 5.3). Meanwhile, the eastern side of each S<sub>v</sub> isophote was shaped by features D<sub>J1</sub>, D<sub>J2</sub>, and A<sub>H</sub>. To isolate the ellipse-tracing portion of each S<sub>v</sub> isophote, the collated coordinate arrays were inspected by-eye and cuts were applied to the horizontal and vertical pixel coordinates. Additionally, isophotes with S<sub>v</sub> > 1.25 mJy arcsec<sup>-2</sup> (J band) or S<sub>v</sub> > 1.85 mJy arcsec<sup>-2</sup> (H band) and S<sub>v</sub> < 0.70 mJy arcsec<sup>-2</sup> (J band) or S<sub>v</sub> < 0.80 mJy arcsec<sup>-2</sup> (H band) were not used as they did not sufficiently trace the apexes of the

<footnote n="4"><sup>4</sup>Position angles are those of the disc major axis, measured East of North.</footnote>

<footnote n="5"><sup>5</sup>This is within the range of the Bailer-Jones et al. (2018) estimate of 151.9 ± 1.9 pc (which accounts for the non-linear nature of the parallax-distance transformation), based on the Gaia DR2 parallax, and the Gaia early DR3 inverse parallax estimate of 154.7 ± 0.5 pc (Gaia Collaboration 2021; Lindegren et al. 2021).</footnote>
<table><caption>**Table 1.** Results of our geometric modelling of the $S_v$ isophotes of the J- and H-band $Q_\phi$ images. Column 1: isophote surface brightness; columns 2, 3, and 4: elliptical ring radius, inclination, and position angle; column 5: height of the scattering surface above the disc mid-plane at radius $r$.</caption><thead><tr><th>S<sub>v</sub><br>(mJy/arcsec<sup>2</sup>)</th><th>r<br>(au)</th><th>i<br>(°)</th><th>PA<br>(°)</th><th>h<sub>scat</sub>(r)<br>(au)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center; font-weight:bold;">J band</td></tr><tr><td>0.70</td><td>71.4<sup>+0.6</sup><sub>-0.6</sub></td><td>68.1<sup>+0.6</sup><sub>-0.7</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.4</sub></td><td>9.5<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>0.75</td><td>71.0<sup>+0.6</sup><sub>-0.6</sub></td><td>68.3<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.5<sup>+0.3</sup><sub>-0.4</sub></td><td>9.6<sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>0.80</td><td>70.7<sup>+0.6</sup><sub>-0.6</sub></td><td>68.2<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.4</sub></td><td>9.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>0.85</td><td>69.5<sup>+0.6</sup><sub>-0.6</sub></td><td>67.1<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.1<sup>+0.3</sup><sub>-0.4</sub></td><td>10.3<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>0.90</td><td>69.1<sup>+0.7</sup><sub>-0.6</sub></td><td>66.9<sup>+0.8</sup><sub>-0.9</sub></td><td>-0.6<sup>+0.4</sup><sub>-0.5</sub></td><td>10.8<sup>+0.6</sup><sub>-0.5</sub></td></tr><tr><td>0.95</td><td>66.9<sup>+0.6</sup><sub>-0.6</sub></td><td>67.7<sup>+0.8</sup><sub>-0.9</sub></td><td>-0.1<sup>+0.3</sup><sub>-0.4</sub></td><td>9.2<sup>+0.6</sup><sub>-0.5</sub></td></tr><tr><td>1.00</td><td>66.6<sup>+0.6</sup><sub>-0.6</sub></td><td>68.4<sup>+0.8</sup><sub>-0.9</sub></td><td>0.3<sup>+0.3</sup><sub>-0.3</sub></td><td>9.0<sup>+0.6</sup><sub>-0.6</sub></td></tr><tr><td>1.05</td><td>66.0<sup>+0.6</sup><sub>-0.6</sub></td><td>69.0<sup>+0.8</sup><sub>-0.8</sub></td><td>0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>8.5<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.10</td><td>65.4<sup>+0.6</sup><sub>-0.6</sub></td><td>68.6<sup>+0.8</sup><sub>-0.9</sub></td><td>0.2<sup>+0.3</sup><sub>-0.4</sub></td><td>8.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.15</td><td>64.7<sup>+0.6</sup><sub>-0.6</sub></td><td>68.6<sup>+0.8</sup><sub>-0.9</sub></td><td>0.2<sup>+0.3</sup><sub>-0.3</sub></td><td>9.1<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.20</td><td>63.6<sup>+0.5</sup><sub>-0.5</sub></td><td>67.8<sup>+0.8</sup><sub>-0.8</sub></td><td>0.1<sup>+0.3</sup><sub>-0.3</sub></td><td>9.4<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.25</td><td>62.9<sup>+0.5</sup><sub>-0.6</sub></td><td>67.8<sup>+0.9</sup><sub>-1.0</sub></td><td>0.2<sup>+0.3</sup><sub>-0.4</sub></td><td>9.4<sup>+0.7</sup><sub>-0.6</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">H band</td></tr><tr><td>0.80</td><td>71.6<sup>+0.5</sup><sub>-0.5</sub></td><td>67.8<sup>+0.6</sup><sub>-0.7</sub></td><td>-0.1<sup>+0.3</sup><sub>-0.3</sub></td><td>8.0<sup>+0.5</sup><sub>-0.4</sub></td></tr><tr><td>0.85</td><td>71.3<sup>+0.4</sup><sub>-0.4</sub></td><td>68.6<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.8<sup>+0.2</sup><sub>-0.3</sub></td><td>8.3<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>0.90</td><td>70.4<sup>+0.4</sup><sub>-0.4</sub></td><td>69.0<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.7<sup>+0.2</sup><sub>-0.3</sub></td><td>7.7<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>0.95</td><td>70.0<sup>+0.5</sup><sub>-0.5</sub></td><td>68.6<sup>+0.6</sup><sub>-0.6</sub></td><td>-1.0<sup>+0.3</sup><sub>-0.3</sub></td><td>8.0<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.00</td><td>69.5<sup>+0.5</sup><sub>-0.5</sub></td><td>68.5<sup>+0.6</sup><sub>-0.7</sub></td><td>-1.0<sup>+0.3</sup><sub>-0.4</sub></td><td>8.3<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.05</td><td>69.4<sup>+0.4</sup><sub>-0.4</sub></td><td>69.6<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.6<sup>+0.2</sup><sub>-0.3</sub></td><td>7.8<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.10</td><td>68.3<sup>+0.4</sup><sub>-0.4</sub></td><td>69.4<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.3<sup>+0.2</sup><sub>-0.3</sub></td><td>7.7<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.15</td><td>68.2<sup>+0.4</sup><sub>-0.4</sub></td><td>69.4<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.2<sup>+0.2</sup><sub>-0.3</sub></td><td>7.7<sup>+0.4</sup><sub>-0.3</sub></td></tr><tr><td>1.20</td><td>67.5<sup>+0.4</sup><sub>-0.4</sub></td><td>69.4<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.3<sup>+0.2</sup><sub>-0.2</sub></td><td>7.7<sup>+0.4</sup><sub>-0.3</sub></td></tr><tr><td>1.25</td><td>67.5<sup>+0.4</sup><sub>-0.4</sub></td><td>68.8<sup>+0.5</sup><sub>-0.6</sub></td><td>-0.5<sup>+0.2</sup><sub>-0.3</sub></td><td>8.3<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.30</td><td>66.7<sup>+0.4</sup><sub>-0.4</sub></td><td>68.3<sup>+0.6</sup><sub>-0.7</sub></td><td>-0.6<sup>+0.2</sup><sub>-0.3</sub></td><td>8.5<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.35</td><td>65.6<sup>+0.4</sup><sub>-0.4</sub></td><td>69.9<sup>+0.6</sup><sub>-0.6</sub></td><td>-0.2<sup>+0.2</sup><sub>-0.2</sub></td><td>6.9<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.40</td><td>65.4<sup>+0.4</sup><sub>-0.4</sub></td><td>70.2<sup>+0.5</sup><sub>-0.6</sub></td><td>-0.1<sup>+0.2</sup><sub>-0.2</sub></td><td>7.0<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.45</td><td>65.1<sup>+0.4</sup><sub>-0.4</sub></td><td>70.7<sup>+0.5</sup><sub>-0.5</sub></td><td>-0.2<sup>+0.2</sup><sub>-0.2</sub></td><td>6.8<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.50</td><td>64.6<sup>+0.4</sup><sub>-0.4</sub></td><td>70.5<sup>+0.6</sup><sub>-0.6</sub></td><td>-0.4<sup>+0.2</sup><sub>-0.3</sub></td><td>7.0<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.55</td><td>64.5<sup>+0.4</sup><sub>-0.4</sub></td><td>70.4<sup>+0.6</sup><sub>-0.6</sub></td><td>-0.5<sup>+0.2</sup><sub>-0.3</sub></td><td>7.1<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>1.60</td><td>63.8<sup>+0.4</sup><sub>-0.4</sub></td><td>69.6<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>7.7<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.65</td><td>62.8<sup>+0.4</sup><sub>-0.4</sub></td><td>69.0<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>7.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>1.70</td><td>62.3<sup>+0.4</sup><sub>-0.4</sub></td><td>70.0<sup>+0.7</sup><sub>-0.8</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>7.0<sup>+0.6</sup><sub>-0.5</sub></td></tr><tr><td>1.75</td><td>62.1<sup>+0.4</sup><sub>-0.4</sub></td><td>70.0<sup>+0.7</sup><sub>-0.9</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>7.0<sup>+0.7</sup><sub>-0.6</sub></td></tr><tr><td>1.80</td><td>61.9<sup>+0.4</sup><sub>-0.4</sub></td><td>70.0<sup>+0.8</sup><sub>-0.9</sub></td><td>-0.6<sup>+0.3</sup><sub>-0.3</sub></td><td>7.0<sup>+0.7</sup><sub>-0.6</sub></td></tr><tr><td>1.85</td><td>61.6<sup>+0.5</sup><sub>-0.4</sub></td><td>69.8<sup>+0.8</sup><sub>-0.9</sub></td><td>-0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>7.1<sup>+0.7</sup><sub>-0.7</sub></td></tr></tbody></table>

ellipse. Limiting the range of $S_v$ isophotes used in this way results in a limited range of disc radii (and therefore $h_{scat}(r)$) being explored.

## 4.1 Inferred disc geometry and potential flaring of the disc scattering surface

The minimize function of LMFIT (version 1.0.1; Newville et al. 2020) was used to fit our elliptical ring model to the data. Specifically, we minimized the radial distance between the cylindrical coordinates of the isophote and the elliptical ring at the same azimuthal angle. Our results are summarized in Table 1 and we show example fits

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example isophotes extracted from the J- (left) and H-band (right) $Q_\phi$ images at a surface brightness of 0.7 and 1.3 mJy arcsec$^{-2}$, respectively, and their corresponding best-fitting elliptical ring model.</figcaption></figure>

in Fig. 4. To estimate the uncertainties on our fits, we used the MINIMIZER.EMCEE package of LMFIT (Foreman-Mackey et al. 2013; Newville et al. 2020) to explore the posterior probability distribution of each of the parameters in our model. The values and uncertainties quoted in Table 1 correspond to the median and $1\sigma$ quantiles of these probability distributions. We do not incorporate the uncertainty on the distance (Table 2) in our $r$ and $h_{scat}(r)$ uncertainties.

Our fits to the J- and H-band $S_v$ isophotes provide broadly consistent results for the disc geometry: $67^\circ \le i \le 69^\circ$ and $-0.7^\circ \le PA_{major} \le 0.6^\circ$ in J band; $68^\circ \le i \le 71^\circ$ and $-1.0^\circ \le PA_{major} \le -0.1^\circ$ in H band. As such, the disc appears more inclined and oriented closer to 0° PA in our GPI images than has previously been inferred from geometric modelling and image reconstruction of PIONIER interferometry (Lazareff et al. 2017; Kluska et al. 2020). Our results are more consistent with findings from Fourier-plane analysis of mm continuum and K-band interferometry: $i = 70.4 \pm 1.2^\circ$ and PA = $1 \pm 1^\circ$ inferred from 1.3 µm ALMA observations (Ansdell et al. 2020); $68^\circ < i < 72^\circ$ and PA = $2 \pm 2^\circ$ inferred from GRAVITY observations (Gravity Collaboration 2019). We discuss these results further, particularly in relation to HD 145718’s photometric variability, in Section 6.1.

The values in Table 1 suggest a scattering surface aspect ratio ($h_{scat}(r)/r$) in the range ~0.13–0.16 (J band) and ~0.10–0.13 (H band). This is at the lower end of the range found by e.g. Avenhaus et al. (2018) and Ginski et al. (2016) who used NIR scattered light images of six protoplanetary discs with concentric ring features to estimate scattering surface aspect ratios, finding $h_{scat}(r)/r = 0.09–0.25$. The results from this simple isophote fitting procedure also suggest that the J-band scattering surface of HD 145718 may be more vertically extended than the H-band scattering surface. This is not wholly unexpected as longer wavelengths of light should penetrate deeper into the disc due to the grain size dependent vertical stratification of the disc as a result of vertical settling (e.g. Pinte et al. 2007; Duchêne et al. 2010). We examine this in more detail using Monte Carlo radiative transfer models in Section 6.2.

Looking closely at our best-fitting $h_{scat}(r)$ and $r$ values, the height of the scattering surface also appears to increase with radius suggesting we may be probing the degree of flaring of the disc scattering surface in the outer disc. We explore this further, and assess the robustness of the relatively simplistic fitting procedure we have employed in Section 6.2.
<table><caption>**Table 2.** Adopted stellar parameters. The effective temperature (column 2) and surface gravity (column 3) are from Fairlamb et al. (2015). The distance (column 4) is the inverse of the Gaia DR2 parallax (Gaia Collaboration 2016, 2018). The visual extinction (column 5), radius (column 6), and luminosity (column 7) were re-evaluated herein using photometry from Lazareff et al. (2017, see Appendix C).</caption><thead><tr><th></th><th><i>T</i><sub>eff</sub><br>(K)</th><th>log (g)</th><th><i>d</i><br>(pc)</th><th><i>A</i><sub>V</sub><br>(mag)</th><th><i>R</i><sub>*</sub><br>(<i>R</i><sub>☉</sub>)</th><th><i>L</i><sub>*</sub><br>(<i>L</i><sub>☉</sub>)</th></tr><tr><th></th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th></tr></thead><tbody><tr><td>HD 145718</td><td>8000 ± 250</td><td>4.37 ± 0.15</td><td>152.5<sup>+3.2</sup><sub>-3.0</sub></td><td>0.89<sup>+0.34</sup><sub>-0.08</sub></td><td>1.97<sup>+0.12</sup><sub>-0.11</sub></td><td>14.3<sup>+3.9</sup><sub>-3.1</sub></td></tr></tbody></table>

## 5 RADIATIVE TRANSFER MODELLING OF THE SED, GPI IMAGES, AND INTERFEROMETRIC VISIBILITIES

We build on our analysis above using Monte Carlo radiative transfer modelling with TORUS (Harries 2000; Harries et al. 2019). This allows us to probe the vertical and radial structure of the disc, and its scattering surface, in a more physically motivated and self-consistent manner. To help constrain our models, we complement our GPI images with new and extant NIR interferometry (Section 3.2) and multiwavelength spectrophotometry (Section 3.3), allowing us to probe the surface layers of the disc over its full radial extent.

The circumstellar environment of HD 145718 was modelled as a passive disc, illuminated by a central star (see Section 5.1), and is built on a two-dimensional, cylindrical adaptive mesh refinement (AMR) grid. The density structure of the gaseous portion of the disc is prescribed following Shakura & Sunyaev (1973):

$$ \rho_{\text{gas}}(r, z) = \frac{\Sigma_{\text{gas}}(r)}{h_{\text{gas}}(r)\sqrt{2\pi}} \exp\left\{-\frac{1}{2}\left[\frac{z}{h_{\text{gas}}(r)}\right]^2\right\}, \quad (6) $$

where $r$ and $z$ are the radial distance from the star into the disc and the vertical distance from the disc mid-plane, respectively. The pressure scale height, $h_{\text{gas}}(r)$, and surface density, $\Sigma_{\text{gas}}(r)$, of the gas are prescribed to follow simple radial power laws:

$$ h_{\text{gas}}(r) = h_{0,\text{gas}} \left( \frac{r}{r_0} \right)^{\beta} \quad (7) $$

and

$$ \Sigma_{\text{gas}}(r) = \Sigma_{0,\text{gas}} \left( \frac{r}{r_0} \right)^{-p}. \quad (8) $$

Here, $h_{0,\text{gas}}$ and $\Sigma_{0,\text{gas}}$, are the pressure scale height and surface density of the gas, respectively, evaluated at canonical radius, $r_0 = 100$ au. We keep the surface density power law exponent, $p = 1.0$ fixed in all models.

The Lucy (1999) algorithm is used to iteratively solve for radiative equilibrium and dust sublimation. The disc is populated with two populations of dust, comprising ‘surface’ and ‘settled’ grains (see Section 5.2). Dust is added to grid cells whose temperature is cooler than the dust sublimation temperature after the fourth Lucy iteration. Convergence is typically achieved after seven iterations.

We further used TORUS to generate model SEDs, 4 × 4 arcsec model Stokes I, Q<sub>φ</sub>, and U<sub>φ</sub> images at $\lambda = 1.25$ µm (J band) and 1.65 µm (H band), and 24 × 24 mas model total intensity images at $\lambda = 1.65$ µm and 2.13 µm (K band). We followed the procedure outlined in Davies et al. (2018) to extract visibilities from our mass-scale images at the (u, v) coordinates of our interferometric data set (Fig. 3).

### 5.1 Stellar parameters

The stellar parameters we adopt as input parameters for our TORUS models are listed in Table 2. We follow Vioque et al. (2018) and adopt the effective temperature, $T_{\text{eff}}$, and surface gravity, log (g), estimates from Fairlamb et al. (2015). However, we choose to re-evaluate the stellar luminosity, $L_*$, and visual extinction, $A_V$, rather than adopt the values in Vioque et al. (2018). Our reasons for this are two-fold:

(i) the Vieira et al. (2003) BVRI photometry used by Vioque et al. (2018) to estimate $L_*$ trace a fainter epoch than the Høg et al. (2000), Lazareff et al. (2017), and Gaia Collaboration (2018) photometry (see Table A1), suggesting the star may be inherently brighter than the Vieira et al. (2003) photometry suggests;

(ii) HD 145718's B-V colour is bluer during fainter epochs, consistent with increased scattering during obscuration by circumstellar dust. Indeed, HD 145718 has previously been identified as displaying UX Ori-type (Poxon 2015) and dipper variability (Ansdell et al. 2018; Cody & Hillenbrand 2018; Rebull et al. 2018). If the dust grains responsible for the occultations are larger, on average, than those in the ISM, the total-to-selective extinction, $R_V$, may be closer to 5.0 (Hernández et al. 2004) than the value of 3.1 adopted by Vioque et al. (2018).

Using the brighter epoch BVRI photometry from Lazareff et al. (2017), we follow the methodology outlined in Fairlamb et al. (2015) to re-estimate $A_V$, $R_*$, and $L_*$ (Table 2). In doing so, we find that consistent values of $R_*$, and $L_*$ can be used to fit the bright and faint epochs of photometry if $R_V$ changes from 3.1 to 5.0 during dimming events (see Appendix C). This suggests that dust grains larger, on average, than those found in the ISM are present in the surface layers of the disc. We discuss this further in Section 6.2.

### 5.2 Disc mass and dust prescription

We prescribe populations of ‘surface’ and ‘settled’ grains in our models, both of which are prescribed as comprising solely of Draine (2003) astronomical silicates.<fnref n="6" /> Our settled grains are larger in size and dominate the disc in terms of its mass. Meanwhile, the grains in the disc surface dominate our GPI and interferometric data.

Our surface grain population is assumed to be well coupled to the gas and therefore follow the vertical and radial density prescriptions in equations (7) and (8). They are prescribed to sublimate when they exceed a temperature,

$$ T_{\text{sub},1} = G \rho_{\text{gas}}^{\gamma}(r, z), \quad (9) $$

<footnote n="6">⁶PAH emission is evident in the IR spectrum of HD 145718 (Keller et al. 2008) but our assumption is a reasonable approximation as the prominent 10 µm spectral feature indicates silicate grains are readily abundant in the surface layers of the disc.</footnote>
where $G = 2000$ and $\gamma = 1.95 \times 10^{-2}$ (Pollack et al. 1994). As the density of disc material is most concentrated around the disc mid-plane and tapers off at larger scale heights, the dependence of $T_{\text{sub,1}}$ on $\rho_{\text{gas}}(r, z)$ results in a curved sublimation rim (Isella & Natta 2005). How far from the star a grain sublimates also depends on how efficiently it can cool and larger grains cool more efficiently than smaller grains. The location and radial extent of the sublimation rim depends on the size of the largest grains in the mixture (Tannirkulam, Harries & Monnier 2007) as these will shield smaller grains from incident stellar radiation. We allow our surface grains to range in size between a fixed minimum value, $a_{\text{min}} = 0.01$ µm, and a maximum that is varied between models: $0.14 \le a_{\text{max}} \le 1.30$ µm. These values of $a_{\text{max}}$ reflect the range over which an increase in grain size produces an increase in cooling efficiency and associated decrease in sublimation radius (Isella & Natta 2005; Davies et al. 2020; Davies & Harries, in preparation) and provide a range of scattering phase functions (Stolker et al. 2016b; Tazaki et al. 2019).

Larger, mm-sized grains are expected to have settled closer to the disc mid-plane and are therefore absent from the disc surface layers. We restrict the vertical extent of these settled grains to a fraction, $f$, of $h_{\text{gas}}(r)$, and vary $f$ between models. This increases the density of material in the disc mid-plane, further affecting the shape and radial extent of the sublimation rim (Tannirkulam et al. 2007). How the settling of dust grains larger than a few microns in size influences the location, shape, and extent of the sublimation rim has not been well explored and is beyond the scope of this paper. Instead, after exploring a range of sublimation temperatures for the mm-sized grains, we set the sublimation temperature of the settled grains to a density-independent value of $T_{\text{sub,2}} = 1200$ K. This ensured that the settled grains were contained within the sublimation rim structure forged by our surface grain population.

Our mm-sized settled grains also dominate the flux at mm-wavelengths. We used the 1.3 mm flux, $F_{\nu}$, reported in Garufi et al. (2018; see Table A1), to estimate a total disc mass (gas + dust),

$$M_{\text{disc}} = \frac{F_{\nu}d^2}{\kappa_{\nu}B_{\nu}(T_{\text{dust}})} = 0.0097 M_{\odot}. \quad (10)$$

Here, $B_{\nu}(T_{\text{dust}})$ is the blackbody radiation at frequency, $\nu$, for dust at temperature, $T_{\text{dust}}$. We assumed $T_{\text{dust}} = 20$ K, and an opacity, $\kappa_{\nu} = 0.1(\nu/10^{12} \text{ Hz})^{\beta_{\kappa}} \text{ cm}^2 \text{ g}^{-1}$ with $\beta_{\kappa} = 1.0$ (Beckwith et al. 1990), which accounts for the adopted 100:1 gas-to-dust ratio. Assuming the density of dust grains follows $n(a) \propto a^{-3.5}$ (where $a$ represents the grain size), the larger, settled grains will contribute 96.7 per cent of the dust mass budget (with the smaller, surface grains contributing the remaining 3.3 per cent).

## 5.3 TORUS parameter grid exploration

We computed a grid of TORUS models in which we varied the gas scale height, $h_{0,\text{gas}}$; the settling height of the 1 mm-sized grains, $f$; the flaring exponent of the gas pressure scale height, $\beta$; the disc outer radius, $R_{\text{out}}$; the disc orientation ($i$ and PA); and the maximum size of the dust grains in the disc surface layers, $a_{\text{max}}$ (see Table 3). For each model, we visually inspected the fit to the SED, the GPI $Q_{\phi}$ and $U_{\phi}$ images, and the NIR visibilities.

In assessing the goodness of fit of each model to the SED, the flux across optical to IR wavelengths produced by the model was allowed to range between the bounds of the dereddened and non-dereddened data (black and grey data points in the top panel of Fig. 9, respectively). In this way, we assumed that the dereddened observed photometry (using $A_{\nu} = 0.89$ and $R_{\nu} = 3.1$) provided

<table><caption>**Table 3.** Best-fitting disc parameters found from our radiative transfer modelling with TORUS (see Section 5.3 for the meanings of each of the symbols). Where applicable, physical sizes in au were converted to angular scales using $d = 152.5$ pc (see Section 4).</caption><thead><tr><th>Parameter</th><th>Values explored</th><th>Best model</th></tr></thead><tbody><tr><td>$h_{0,\text{gas}}$ (au)</td><td>5, 6, 7, 8, 9, 10, 11</td><td>10</td></tr><tr><td>$f(h_{0,\text{gas}})$</td><td>0.05, 0.1, 0.2, 0.3, 0.4, 0.5</td><td>0.1</td></tr><tr><td>$\beta$</td><td>1.07, 1.08, 1.09, 1.11, 1.13, 1.15, 1.17</td><td>1.15</td></tr><tr><td>$a_{\text{max}}$ (µm)</td><td>0.14, 0.30, 0.40, 0.50, 0.60, 0.70, 1.30</td><td>0.50</td></tr><tr><td>$R_{\text{out}}$ (au)</td><td>70, 75, 80</td><td>75</td></tr><tr><td>$i$ (°)</td><td>48, 65, 68, 70, 72, 74, 76, 78, 80</td><td>72</td></tr><tr><td>PA (°)</td><td>-2, 0, +2</td><td>0</td></tr></tbody></table>

an upper limit to the allowed model flux and the non-dereddened observed photometry (i.e. $A_{\nu} = 0$) provided an extreme lower limit.

In the sections that follow, we briefly discuss the impact of these parameters on the shape of the disc and on the synthesized observables before exploring the model providing the best overall fit in greater detail in the next subsection.

### 5.3.1 Gas scale height

Increasing $h_{0,\text{gas}}$ inflates the vertical extent of the gaseous disc across all radii. This results in increased excess emission in the SED across NIR to mm wavelengths and, in general, brighter $Q_{\phi}$ and $U_{\phi}$ images. The height of the scattering surface is also increased, resulting in a broader gap between the elliptical and arc features in the $Q_{\phi}$ images. The NIR visibilities are sensitive to the contrast between the stellar and circumstellar flux components as well as the shape of the circumstellar NIR emitting region. An increase in $h_{0,\text{gas}}$ results in a greater fraction of circumstellar flux, decreasing the visibilities at spatial frequencies which at least partially resolve the inner disc rim. Above a certain value of $h_{0,\text{gas}}$ (which is dependent on the values of $\beta$, $a_{\text{max}}$, and $i$), the inner disc regions are tall enough to result in direct line-of-sight occultation of the central star (as seen previously for RY Tau; Davies et al. 2020). This decreases the SED flux across optical wavelengths and increases the fraction of total NIR flux emanating from the circumstellar regions in the synthesized image, further lowering the visibilities.

### 5.3.2 Flaring exponent on the gas pressure scale height

As $\beta$ is increased between models, the spectral index between ~20 and 100 µm flattens out while the flux across the NIR to mid-IR (MIR) decreases. The visibility level on any given baseline increases due to the reduction in circumstellar NIR flux while the $Q_{\phi}$ and $U_{\phi}$ images brighten as more material is distributed to larger radii.

### 5.3.3 Settled height of mm grains

The $Q_{\phi}$ and $U_{\phi}$ images are unaffected by changes in $f$. We used a ray-tracing algorithm in TORUS to trace the $\tau_{\nu} = 1.0$ scattering surface at $J$ and $H$ band and found that the height of the scattering surface remained unchanged beyond $\approx 1.5 - 2$ au (Fig. 5) when increasing the value of $f$ from 0.1 $h_{0,\text{gas}}$ to 0.4 $h_{0,\text{gas}}$. Differences between these models are observed in the SED (shortward of ~20 µm) and the visibilities. As the value of $f$ is increased, the total NIR to MIR flux in the SED drops. For $f \gtrsim 0.3$, the inner disc rim broadly resembles the curved rim in the models of e.g. Isella & Natta (2005), Tannirkulam et al. (2007), and Kama, Min & Dominik (2009; see the surface traced
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: Comparison between the *J*- and *H*-band scattering surfaces resulting from TORUS models with different values of *f*. All other model parameters were identical in the two models. Bottom: zoom-in on the inner ~2 au of the *H*-band scattering surface traced by our ray-tracing algorithm. The horizontal spread in data points highlights the partially optically thick nature of the upper layers of the dust sublimation rim.</figcaption></figure>

by the magenta data points in Fig. 5). For smaller values of *f*, the rim has a 'stepped' feature (see the surface traced by blue data points in Fig. 5), more akin to the disc rim model of McClure et al. (2013). Here, a more tightly curved sublimation rim forms with inner edge at $r_{in}$ and a more loosely curved surface emerges above this at $r > r_{in}$. The emergence of the small grains out of the settled rim extends the mas-scale NIR brightness distribution to larger scales. This improves the fit to the shorter baseline visibilities that are overresolved in the models using larger *f* values.

### 5.3.4 Maximum size of surface-layer dust grains

The value of $a_{max}$ affects the radius at which the small grains emerge above the settled disc rim. Isella & Natta (2005) previously showed that larger grains can survive at higher temperatures as they are more efficient at cooling. This same process is responsible for the effect we see here. Increasing $a_{max}$ therefore affects the NIR flux level and the shape of the visibilities in much the same way as seen for single grain size models (Isella & Natta 2005; Davies et al. 2018, 2020). Changing $a_{max}$ also affects the scattering phase function: at the relatively high inclinations we explore, we see an asymmetric $Q_{\phi}$ brightness distribution and increasing $a_{max}$ results in a decrease in back-scattering efficiency, relative to forward scattering (Fig. 6). The extent of this difference is consistently more marked in the *J*-band image than the *H*-band. This is associated with the dependence

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Impact of changing $a_{max}$ on *J*- and *H*-band model $Q_{\phi}$ images (left-hand and right-hand panels, respectively). From top to bottom, the maximum grain size of non-settled grains is increased from 0.3 to 0.6 µm. Each image is 2 × 2 arcsec.</figcaption></figure>

of the scattering phase function on the grain size (Stolker et al. 2016b; Tazaki et al. 2019). One can also see from Fig. 6 that the arc features on the eastern side of the $Q_{\phi}$ images are dimmer for models with larger $a_{max}$. The $U_{\phi}$ images also decrease in brightness with increasing $a_{max}$.

### 5.3.5 Best-fitting TORUS model

The TORUS model providing the best overall fit was found to have $h_{0,gas} = 10$ au, $\beta = 1.15$, $R_{out} = 75$ au, $a_{max} = 0.50$ µm, $i = 72^{\circ}$, and PA = 0$^{\circ}$. The corresponding model $Q_{\phi}$ and $U_{\phi}$ images are compared to our GPI observations in Figs 7 and 8, respectively. The model SED and mas-scale images are compared to the respective observational data in Fig. 9.

The model $Q_{\phi}$ image is able to broadly replicate the $S_{v}$ level, location and extent of the main $E_{H}$ and $A_{H}$ features in the observed *H*-band image (Fig. 7). Reducing *i* or increasing the height of the scattering surface by increasing $h_{0,gas}$ and/or $\beta$ causes the separation between the elliptical and the arc feature in the model image to increase. Increasing $a_{max}$ decreases the overall $S_{v}$ in the $Q_{\phi}$ and $U_{\phi}$ images and reduces the back-scattering efficiency, resulting in a larger dark portion on the west side of the elliptical feature in the *J*- and *H*-band model $Q_{\phi}$ images, like those seen in Fig. 6. The weaker back-scattering we observe in the disc of HD 145718 has also been observed in scattered light imaging of other inclined discs
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: J- (top) and H-band (bottom) $Q_\phi$ images for our best TORUS model ($h_{0.gas} = 10$ au, $\beta = 1.15$, $R_{\text{out}} = 75$ au, $a_{\text{max}} = 0.50$ µm, $i = 72^\circ$, and PA = 0$^\circ$). Middle: observed $Q_\phi$ images overlaid with 0.8 and 4.0 mJy arcsec$^{-2}$ surface brightness contours extracted from these model images (white solid and dashed lines, respectively). Right: zoomed-in view comparing the 0.8 and 4.0 mJy arcsec$^{-2}$ contours extracted from the model (black solid and dashed lines, respectively) and the 0.8 and 1.8 mJy arcsec$^{-2}$ contours extracted from the observed $Q_\phi$ images (cyan solid and dashed lines, respectively).</figcaption></figure>

(e.g. DoAr 25, Garufi et al. 2020; IM Lup, Avenhaus et al. 2018), suggesting that the surface layers of discs may be routinely populated by grains of size, $a \gtrsim \lambda/2\pi$. These larger grains possibly have an aggregate structure that provides them with aerodynamic support against settling. Alternatively, this may indicate a relative dearth of smaller grains ($a \ll \lambda/2\pi$). For instance, Wolff et al. (2021) found smaller dust grains were confined to a more diffuse region above the disc surface in their modelling of SSTC2D J163131.2–242627. We do not see evidence of such a diffuse region around HD 145718.

None of the models we explored were able to reproduce the four quadrants of emission seen in our $U_\phi$ images. Instead, we note that the eastern quadrants of emission in the $U_\phi$ image in Fig. 1 appear to overlap with the dark lane between the front and rear sides of the disc. As we discuss in Appendix B, the shape and $S_v$ levels in our $U_\phi$ are sensitive to our method of stellar and instrument polarization correction. Without higher angular resolution observations from instruments with improved instrument polarization characterization, we are unable to assess whether this discrepancy is due to some underlying astrophysical process or an artefact of imperfect calibration.

The model SED is displayed as the black dashed line in the top panel of Fig. 9 while the photometric data are shown as black filled circles (de-reddening applied assuming $A_V = 0.89$ mag with $R_V = 3.1$) or grey open circles (no de-reddening applied). Our model SED reproduces the shape of the full SED well. We also show the same model SED computed at a lower inclination of 48° (solid black line). The slightly lower flux across optical wavelengths between the low and intermediate inclination model SEDs indicates that the circumstellar disc slightly occults the star along our line of sight in this model. Multicolour photometric monitoring (e.g. Petrov et al. 2019) and/or contemporaneous photometric and polarimetric/interferometric monitoring is required to confirm whether dust in the surface layers of disc obscures the star even during bright epochs.

The visibility profiles in Fig. 9 (middle and bottom panels) are displayed as a function of deprojected baseline

$$B_{\text{de-proj}} = B \left[ \sin^2(\phi) + \cos^2(i) \cos^2(\phi) \right]^{1/2}. \qquad (11)$$

Here, B is the baseline length, i is the disc inclination, and φ is the difference between the baseline position angle and the disc minor axis
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Model $U_\phi$ images for our best TORUS model ($h_{0,gas}$ = 10 au, $\beta$ = 1.15, $R_{\text{out}}$ = 75 au, $a_{\text{max}}$ = 0.50 µm, $i$ = 72°, and PA = 0°). Top: J-band; bottom: H-band.</figcaption></figure>

position angle. By displaying the visibility as a function of $B_{\text{de-proj}}$, we account for the foreshortening of the brightness distribution along baseline position angles that trace the disc minor axis. Any vertical spread in visibility still present in the plot should then reflect a wavelength-dependence in the data, azimuthal variations in the flux contrast, imperfect calibration of the data, or noise. As our mas-scale model images are computed at a single wavelength (1.65 µm for H band; 2.13 µm for K band), we implicitly assume a greybody approximation to the visibilities in each waveband. To account for the azimuthal variations in image brightness seen in our mas-scale model images, we extracted visibility curves along baseline position angles that trace the disc major and minor axes. These are shown as grey dashed and dot-dashed lines in the middle and bottom panels of Fig. 9, respectively.

None of the models we explored could simultaneously reproduce both the sharp drop in visibility for $B_{\text{de-proj}} \lesssim 40$ M$\lambda$ and the relatively flat visibility profile beyond. However, models with $f < 0.3 h_{0,gas}$ are able to reproduce the general shape of the visibilities at longer baselines and recover the visibility level on the longest H-band baselines well. The value of $f$ we infer from our best-fitting model (0.1 $h_{0,gas}$) is consistent with recent results from Villenave et al. (2020) where the vertical extents of mm grains in a sample of edge-on discs were found to be on the order of a few au at 100 au.

Our inability to recover the sharp drop in visibility at short baselines suggests the presence of more extended circumstellar NIR emission than we are able to produce with our current models, as previously indicated by Lazareff et al. (2017) and Kluska et al. (2020). This may indicate the presence of a photoevaporative or magnetohydrodynamically driven disc wind, like that inferred for SU Aur based on similar analyses (Labdon et al. 2019). However, unlike for SU Aur (Ginski et al. 2021), we see no evidence of an outflow on larger scales in our GPI images but we note that the image regions along the disc minor axis are most affected by the

<figure><img src="image_11.png" /><figcaption>**Figure 9.** SED (top) and visibilities (middle and bottom) of our best TORUS model ($h_{0,gas}$ = 10 au, $\beta$ = 1.15, $R_{\text{out}}$ = 75 au, $a_{\text{max}}$ = 0.50 µm, $i$ = 72°, and PA = 0°) compared to the observed spectrophotometry and NIR interferometry. Top: model SEDs at $i$ = 72° (dashed black line) and $i$ = 48° (solid black line) are compared to the *Spitzer* spectrum (blue line) and photometry (open grey circles have zero de-reddening applied; black filled circles have been de-reddened using $A_V$ = 0.89 and $R_V$ = 3.1 – see Table 2). Middle and bottom: deprojected visibility profiles extracted from the 1.65 µm (middle) and 2.13 µm (bottom) model total intensity images along baseline position angles tracing the major (grey dashed line) and minor (grey dot-dashed line) disc axes. Visibilities on the shortest baselines are overestimated by the model, indicating an additional extended emission component, unaccounted for in our models, is present. Colours represent the different beam combiners: see Fig. 3. The 24 × 24 mas-scale model total intensity images are shown inset at the lower left corner of each plot.</figcaption></figure>
choice of stellar polarization subtraction (Appendix B). Alternatively, it may indicate that our non-settled, 0.01 µm ≥ *a* ≥ 0.50 µm-sized dust grain mixture does not fully describe the disc surface, at least in the innermost disc regions. If we were to further segregate our dust prescription by grain size into three or more populations, we would likely see a smaller grain population emerge above this surface at larger radii, thus extending the NIR emitting region. Further investigation of this necessitates detailed theoretical work to simulate the combined effects of coagulation, settling, radial drift, collisional fragmentation, and sublimation to predict where differently sized grains exist in the sublimation rim.

# 6 DISCUSSION

## 6.1 Disc orientation and dust obscuration

Our GPI images of HD 145718 reveal an inclined disc with its major axis oriented along a north-south direction. In Section 4, we fit elliptical ring models to coordinates tracing isophotes of surface brightness in the *J*- and *H*-band *Q*<sub>Φ</sub> images. Specifically, the full list of coordinates was trimmed to avoid features within (or close to the edge of) the inner working angle of the coronograph and the western portion of the *E*<sub>J</sub> and *E*<sub>H</sub> features, which deviated from an elliptical shape (Fig. 4). From our radiative transfer analysis in Section 5, we saw that this deviation is attributable to the scattering phase function of the dust grains in the surface layers of the disc (Fig. 6). The strong forward scattering and weak back-scattering we observe is typical of grain mixtures dominated by grains of size *a* ≳ λ/2π. These large grains are likely porous, providing them with aerodynamic support against settling (see e.g. Mulders et al. 2013).

We infer a disc inclination in the range 67–71°, with major axis position angle between 1.0° west of north and 0.6° east of north. These results are consistent with previous assessments based on mm continuum and *K*-band interferometry (Gravity Collaboration 2019; Ansdell et al. 2020). Why Lazareff et al. (2017) and Kluska et al. (2020) infer lower inclinations from *H*-band VLTI/PIONIER data is unclear: the NIR continuum emitting regions at *H* and *K* band are expected to be roughly coincident and therefore strongly aligned. If the lower inclination of 48° measured by Kluska et al. (2020) is used to de-project the baseline (equation 11), the PIONIER and MIRC-X are not observed to follow the approximately Gaussian profile we see in the middle panel of Fig. 9 when we use $i = 72°$, regardless of the PA. We note that Kluska et al. (2020) estimate the disc geometry from reconstructed images while Gravity Collaboration (2019) and Ansdell et al. (2020) perform their analysis in the Fourier plane. To reliably recover the inclination of highly inclined discs when using image reconstructions, one must ensure that the minor axis is well resolved. Otherwise, as in this case, the emission along the minor axis is smoothed out by the interferometric beam and the disc will appear less inclined.

At the inclination we infer, the surface layers of the inner disc of HD 145718 partially obscure the star along the observer's line of sight. This is most clearly seen in the top panel of Fig. 9 where the low inclination ($i = 48°$) model SED (black solid line) is slightly higher than the $i = 72°$ model (black dashed line) across optical wavelengths. This difference is within the allowed range of the de-reddened multi-epoch photometry (filled data points). In estimating the extinction, stellar luminosity, and radius, we had found that differences between the bright and faint epoch photometry could be explained by a difference in the value of the total-to-selective extinction (Appendix C). Similarly, HD 145718 is known to exhibit UX Ori and dipper variability, typically attributed to aperiodic stellar

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Apparent height of the scattering surface, $h_{scat}(r)$, as a function of the elliptical ring radius, *r* (data are from Table 1: magenta and green data points are from fits to the *J*-band and *H*-band isophotes, respectively). The dashed lines show the best-fitting power law profiles (equation 7) to the *H*-band (green line), the *J*-band (magenta line), and the combined *J*- and *H*-band data (black line). The amplitude, $h_{65}$, and exponent, $\beta$, of the power law fits are provided in the legend.</figcaption></figure>

occultation by circumstellar dust close to the star (e.g. Dullemond et al. 2003; Tambovtsева & Grinin 2008). What the data considered herein are not able to definitively assess is whether the stellar surface is always at least partially obscured by circumstellar material. Multicolour photometric monitoring of HD 145718 such as that undertaken for RY Tau and SU Aur (Petrov et al. 2019) would be useful to establish this.

## 6.2 Assessing the robustness of our isophote fitting procedure

In addition to assessing the geometry of the disc, our isophote fitting procedure (Section 4) allowed us to explore the vertical extent of the *J*- and *H*-band scattering surfaces in a much more time and computationally efficient manner than our Monte Carlo radiative transfer modelling allows. Our results from these simple models suggested that (i) the *J*-band scattering surface is more extended than the *H*-band scattering surface and (ii) there appears to be a slight increase in $h_{scat}(r)$ with *r* in both wavebands. We show both of these effects more clearly in Fig. 10 by plotting our best-fitting $h_{scat}(r)$ values against the associated best-fitting *r* values for each $S_v$ in Table 1: magenta and green data points represent the *J*- and *H*-band data, respectively. In this section, we provide some background on why probing the vertical structure of the disc is important, assess the degree of flaring inferred from our isophote fits, and examine the robustness of our isophote fitting by repeating the procedure on the synthetic observations computed from our best-fitting TORUS model.

The vertical extent of the gaseous component of circumstellar discs is expected to follow a flared profile (equation 7). For instance, based on theoretical predictions for centrally irradiated, steady-state accretion discs, the power-law exponent on the gas scale height, $\beta_{gas}$, is expected to lie in the range 1.125–1.3 (e.g. Kenyon & Hartmann 1987; Chiang & Goldreich 1997). Indeed, we find that the value of $\beta$ in our best-fitting TORUS model (=1.15) lies in this range. However, even in scenarios where dust grains are well coupled to the gas, these values of $\beta$ are not expected to also describe the flaring exponent of the scattering surface (see e.g. Avenhaus et al. 2018). This is because the gas pressure scale height depends only on the gas temperature, while the height of the scattering surface depends on the dust properties (e.g. opacity, scattering phase function etc.). Thus, probing the vertical disc structure allows us to probe the properties of
the dust content of the disc. In recent years, observational constraints on the degree of flaring in the surface layers of dusty and gaseous components of protoplanetary discs have begun to be made (Ginski et al. 2016; Avenhaus et al. 2018; Pinte et al. 2018; Villenave et al. 2020; Rich et al. 2021). Elliptical models, similar to those we employ in Section 4, have been used to determine the height of the scattering surface in concentric ring features of discs exhibiting substructure, with power-law profiles (equation 7) then used to constrain the flaring exponent of the scattering surface, $\beta_{scat}$. At the resolution of current observations, HD 145718 does not show indications of disc substructure and so these same methods are not applicable.

Assuming, first of all, that our surface brightness isophotes do trace concentric disc annuli, we use the results of our isophote fitting in Table 1 to derive an initial estimate for the flaring of the J- and H-band scattering surfaces close to the apparent outer edge of the disc. Specifically, we employ the power-law parametrization for the vertical height of the scattering surface (e.g. equation 7). Based on the range of r values in Table 1, we adopt a canonical radius, $r_0 = 65$ au such that $h_0$ is the height at 65 au (which we denote $h_{65}$). To estimate $\beta_{scat}$ and $h_{65}$, we performed least-squares fitting to the linear relation

$$ \log_{10}(h_{\text{scat}}(r)) = \log_{10}(h_{65}) + \beta_{\text{scat}} \log_{10}(r/65 \text{ au}). \quad (12) $$

Specifically, we drew 10 000 realizations of $h_{scat}(r)$ and r from split-normal distributions, based on their lower and upper bounded errors, and repeated the fit each time. The best-fitting values and uncertainties on $\beta_{scat}$ and $h_{65}$ are then the median and $1\sigma$ upper and lower quartiles from these fits. Fitting the data in this manner allowed us to account for the asymmetric errors on $h_{scat}(r)$ and r. The resulting profiles are shown by the dashed lines in Fig. 10 while the shaded regions illustrate the uncertainties on the best-fitting values: $h_{65} = 9.2_{-1.0}^{+1.0}$ with $\beta_{scat} = 0.7_{-0.4}^{+0.4}$ for the J-band data (magenta dot-dashed line and shaded region); $h_{65} = 7.4_{-1.0}^{+1.0}$ with $\beta_{scat} = 1.0_{-0.3}^{+0.3}$ for the H-band data (green dot-dashed line and shaded region); and $h_{65} = 7.4_{-1.0}^{+1.0}$ with $\beta_{scat} = 1.0_{-0.3}^{+0.3}$ for the combined J- and H-band data (black dashed line and grey shaded region). These values of $\beta_{scat}$ are within the range of previous observational results determined by fitting power laws to the heights of concentric ring features in discs with substructure: Avenhaus et al. (2018) found $\beta_{scat} = 1.605 \pm 0.132$ for V4046 Sgr, $\beta_{scat} = 1.116 \pm 0.095$ for RX J 1615.3-3255, and $\beta_{scat} = 1.271 \pm 0.197$ (IM Lup) while Ginski et al. (2016) found $\beta_{scat} = 1.73$ for CU Cha.

We used our synthetic $Q_{\phi}$ TORUS images to examine whether our isophote fitting procedure is indeed tracing the height of the scattering surface over multiple disc radii (and therefore the flaring of the disc scattering surface). Using our synthetic $Q_{\phi}$ images (top row of Fig. 7), we repeat the procedure outlined in Section 4 to extract the coordinates of $S_v$ isophotes. We trimmed the data using the same constraints as before so as to reduce the effect of the low back-scattering efficiency. We kept i and PA as free parameters and plot the resultant r and $h_{scat}(r)$ values in Fig. 11 (we use magenta and green crosses to signify J- and H-band data, respectively). Again, we use equation (7), to estimate $h_{65}$ and $\beta_{scat}$ from the results of the isophote fits. We find that the inferred J- and H-band scattering surfaces are coincident and the apparent flaring is much more pronounced, with $\beta \sim 2.4-2.5$, casting doubt on the applicability of our isophote fitting procedure.

To further inspect our isophote modelling procedure, we used a ray-tracing algorithm to compute the $\tau_{scat} = 1.0$ scattering surface at 1.25 and 1.65 µm. We compare this to the apparent surface traced by our isophote fits to the best-fitting synthetic $Q_{\phi}$ images in Fig. 11. The $\tau_{scat} = 1.0$ scattering surfaces at 1.25 and 1.65 µm are indicated by the dashed magenta and green lines, respectively. Two things are

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Application of our isophote fitting procedure (Section 6.1) to the synthetic J- and H-band $Q_{\phi}$ images computed from our best-fitting TORUS model. Radii and heights extracted from the J- and H-band images are plotted in magenta and green, respectively. The power-law fits to these data are shown by the dot-dashed lines, with the shaded areas representing the $1\sigma$ uncertainties on $h_{65}$ and $\beta_{scat}$ (see plot legend). The $\tau_{scat} = 1$ scattering surfaces at 1.25 µm (J-band) and 1.65 µm (H-band), as measured by TORUS, are shown by the dashed lines. The cross-hair in the lower right corner represents the pixel scale in our GPI and TORUS images.</figcaption></figure>

immediately clear: (i) our isophote fitting procedure does not recover the degree of flaring in the $\tau_{scat} = 1.0$ scattering surface; and (ii) the J-band scattering surface does extend to larger scale heights than the H-band scattering surface. The J-band surface has an aspect ratio of ~0.24 while that of the H-band is ~0.22. Our best-fitting TORUS model has a gas scale height of 10 au at 100 au and, with $\beta_{gas} = 1.15$, this corresponds to a height of ~6.1 au at 65 au, indicating that the NIR scattering surface lies at a height of ≈2.5 gas pressure scale heights.

Close to the outer edge of the disc, our isophote fitting method does provide a reasonable estimate of the vertical extent of the scattering surface, given the pixel resolution (indicated by the cross-hair in the lower right corner of Fig. 11) and the centring accuracy of the GPI coronograph (equivalent to ~0.47 au at a distance of 152.5 pc; see Section 3.1).

To understand why our isophote fitting procedure does not recover the scattering surface traced by the ray-tracing algorithm, we determined the elliptical ring parameters of the ellipse drawn out by the $\tau_{scat} = 1.0$ scattering surface at r = 55, 65, and 75 au and compared these to the isophotes extracted from the model image at the same disc radii. These are shown in Fig. 12, overlaid with the $S_v$ isophote coordinates and corresponding best-fitting ellipse.

We find that for the smallest values of $S_v$ (left-most panel in Fig. 12), the coordinates tracing the eastern side of the ellipse appear offset further to the east compared to the $\tau_{scat} = 1.0$ scattering surface at r = 75 au. This is right at the outer edge of the disc ($R_{out} = 75$ au; see Table 3), indicating that these $S_v$ isophotes trace scattering events at the radial edge of the disc, rather than in the disc surface layers. This results in an overestimation of the width of the ellipse along its minor axis, therefore affecting the inferred values of $h_{scat}(r)$ and i. Indeed, fits using larger values of $S_v$ tend to larger values of i (see Fig. 13). The second panel from the left in Fig. 12 corresponds to the ellipse fit where the inferred inclination matches the prescribed value of 72°. Even here, the eastern portion of the ellipse is well traced but the isophote does not sufficiently trace the western portion of the apex of the ellipse, resulting in an underestimation of $h_{scat}(r)$. The apparent truncation of the ellipse
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Coordinates of $S_v$ isophotes extracted from synthetic TORUS model $Q_\phi$ images (red crosses) and the resultant elliptical ring fits (solid black ellipse) are compared to ellipses traced by the wavelength-dependent $\tau_{scat} = 1.0$ scattering surface at 55 au (dashed blue ellipse), 65 au (dashed orange ellipse), and 75 au (dashed green ellipse). Top row: *J*-band; bottom row: *H*-band. The values of $S_v$ are provided above each subplot window.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Trend observed between the surface brightness of the isophote, $S_v$, and the inferred inclination, *i*, when fitting elliptical ring models to the isophotes extracted from the synthetic $Q_\phi$ images computed using TORUS.</figcaption></figure>

apex to the west is worse for larger $S_v$ isophotes, giving rise to the apparent $S_v-i$ relation we observe, and is attributable to the low back-scattering efficiency for the grain mixture used in this model (see Fig. 6). The isophote fitting technique used in Section 4 may therefore be improved by further requiring the ellipse apex to pass through the coordinates with the largest radial extent from the image centre and/or by independently constraining the disc inclination (e.g. using ALMA).

# 7 SUMMARY

We present a multi-instrument study of the intermediate mass YSO, HD 145718. Specifically, we complemented GPI *J*- and *H*-band polarized differential imaging from G-LIGHTS (Monnier et al. 2017, 2019; Laws et al. 2020) with NIR interferometry (archival VLTI/GRAVITY and VLTI/PIONIER data plus new CHARA/MIRC-X observations), IR spectroscopy, and multiband, multi-epoch photometry.
Our GPI images reveal an inclined disc with major axis position angle close to 0°. The offset between the centroid of emission in the open face of the disc and the image centre indicates the disc scattering surface is elevated above the disc mid-plane. Further, the strong forward scattering and weak back-scattering evident in the $Q_\phi$ images indicate that large grains (size, $a \gtrsim \lambda/2\pi$) are abundant in the surface layers of the disc.

Inspection of the multi-epoch BVRI photometry indicated that previously published estimates of HD 145718's stellar luminosity and radius have been affected by direct line-of-sight occultation by circumstellar dust. Using the brightest epoch photometry retrieved by SEDBYS (Davies 2021) and a total-to-selective extinction, $R_V = 3.1$, we re-estimated the visual extinction ($A_V = 0.89_{-0.08}^{+0.34}$ mag), stellar radius ($R_* = 1.97_{-0.11}^{+0.12} R_\odot$), and stellar luminosity ($L_* = 14.3_{-3.1}^{+3.9} L_\odot$) for HD 145718. Furthermore, we find consistent estimates of $R_*$ and $L_*$ when applying the same method to faint epoch photometry from Vieira et al. (2003) and using $R_V = 5.0$. This further supports the idea that the UX Ori and dipper photometric variability displayed by HD 145718 is associated with obscuration by dust grains and that, during periods of occultation, the dust grains are larger, on average, than those in the ISM (Hernández et al. 2004). Further multicolour photometric monitoring, similar to that undertaken for RY Tau and SU Aur by Petrov et al. (2019), is required to determine whether smaller circumstellar grains still contribute to the extinction during the brightest epochs.

We used an off-centre elliptical ring model to fit isophotes of surface brightness in the $Q_\phi$ images, finding:

(i) the optically thick disc extends out to a radius of ~75 au, assuming a distance to HD 145718 of 152.5 pc (Gaia Collaboration 2018);

(ii) the disc is oriented with major axis position angle between -1.0° and 0.6° east of north and inclined at 67° to 71°, consistent with previous measurements based on mm continuum and K-band interferometry (Gravity Collaboration 2019; Ansdell et al. 2020).

We used detailed radiative transfer modelling to self-consistently investigate the radial and vertical disc structure and to assess the extent to which our isophote fitting could be used to probe the disc scattering surface. From our radiative transfer modelling, we found that a model comprising a centrally illuminated passive disc with gas pressure scale height, $h_{0,gas} = 10$ au, flaring exponent, $\beta = 1.15$, outer disc radius, $R_{\text{out}} = 75$ au, maximum size of non-settled grains, $a_{\text{max}} = 0.50$ µm, large-grain settling factor, $f = 0.1 h_0$, inclination, $i = 72°$, and major axis position angle, $\text{PA} = 0°$ provides a good fit to:

(i) the surface brightness, location, and extent of the ellipse and arc features in the J- and H-band $Q_\phi$ images;

(ii) the shape of the full SED from optical to millimetre wavelengths;

(iii) the general shape and stellar-to-circumstellar flux contrast level traced by the H-band visibilities on the largest baselines probed;

(iv) and the general shape of the K-band visibilities.

However, the model could not account for the immediate drop in visibility on the shortest baselines in both H and K bands. As previously suggested by e.g. Kluska et al. (2020), this likely indicates the presence of more extended NIR emission, potentially in the form of an outflow. We find no robust evidence of this outflow on the larger scales probed by our GPI images. Further assessment of this requires better assessment of the instrument polarization which would allow us to improve our stellar polarization subtraction.

By comparing synthetic images to the $Q_\phi$ images we obtained with GPI, we find that fitting ellipses to isophotes of $Q_\phi$ surface brightness recovers the general elliptical shape of the emission but, due to the azimuthally dependent scattering efficiency, cannot reliably recover both the disc inclination and scattering surface height. We propose that our simple isophote fitting method could be improved by independently constraining the disc inclination using e.g. ALMA continuum observations to counter the preference we observe for lower (higher) surface brightness isophotes to appear less (more) inclined.

## ACKNOWLEDGEMENTS

The authors thank the anonymous referee for their comments that improved the content of this paper. We thank René Oudmaijer for useful discussions regarding the estimation of visual extinction and Fred C. Adams for a detailed reading of the manuscript. CLD acknowledges financial support from the College of Engineering, Mathematics and Physical Sciences at the University of Exeter. EAR and JDM acknowledge funding from a National Science Foundation grant (reference NSF-AST1830728). ASEL acknowledges financial support from a Science Technology and Facilities Council (STFC) studentship (reference 1918673). S.K. acknowledges support from an ERC Starting Grant (Grant Agreement No. 639889). A.L. received support from a Science Technology and Facilities Council (STFC) studentship (No. 630008203). JB acknowledges support by NASA through the NASA Hubble Fellowship grant \#HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This study is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério Ciência, Tecnologia e Inovação (Brazil) and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). This work is based upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1636624 and AST-1715788. Institutional support has been provided from the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. MIRC-X received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant No. 639889). The calculations for this paper were performed on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter. This work uses data obtained from the ESO Science Archive Facility. This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology; the Jean-Marie Mariotti Center OiDB service;<fnref n="7" /> the SIMBAD database, operated at CDS, Strasbourg, France; the VizieR catalogue access tool, CDS, Strasbourg, France; NASA's Astrophysics Data System Bibliographic Services.

<footnote n="7">⁷Available at http://oidb.jmmc.fr.</footnote>
## DATA AVAILABILITY

The GPI data are available from the Gemini Observatory Archive at https://archive.gemini.edu/searchform and can be accessed using proposal number GS-2018A-LP-12. The MIRC-X data will be made available through the OiDB (http://oidb.jmmc.fr) following publication. VLTI/GRAVITY and VLTI/PIONIER data are available in the ESO archive (http://archive.eso.org/cms.html) and the OiDB and can be accessed using the target name 'HD 145718'. The photometry are accessible through SEDBYS, available at https://gitlab.com/clairedavies/sedbys. The IR spectra are available in the IRSA (https://irsa.ipac.caltech.edu).

## REFERENCES

* Alecian E. et al., 2013, MNRAS, 429, 1001

* Ansdell M. et al., 2018, MNRAS, 473, 1231

* Ansdell M. et al., 2020, MNRAS, 492, 572

* Anugu N. et al., 2018, in Creech-Eakman M. J., Tuthill P. G., Mérand A., eds, Proc. SPIE Conf. Ser. Vol. 10701, Optical and Infrared Interferometry and Imaging VI. SPIE, Bellingham, p. 1070124

* Anugu N. et al., 2020, AJ, 160, 158

* Arun R., Mathew B., Manoj P., Ujjwal K., Kartha S. S., Viswanath G., Narang M., Paul K. T., 2019, AJ, 157, 159

* Avenhaus H. et al., 2018, ApJ, 863, 44

* Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58

* Beckwith S. V. W., Sargent A. I., Chini R. S., Guesten R., 1990, AJ, 99, 924

* Bourges L., Mella G., Lafrasse S., Duvert G., Chelli A., Le Bouquin J. B., Delfosse X., Chesneau O., 2017, in Manset N., Forshay P., eds, ASP Conf. Ser. Vol. 485. Astron. Soc. Pac., San Francisco, p. 223

* Bouvier J. et al., 1999, A&A, 349, 619

* Bredall J. W. et al., 2020, MNRAS, 496, 3257

* Brisset J., Heißelmann D., Kothe S., Weidling R., Blum J., 2017, A&A, 603, A66

* Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

* Carmona A., van den Ancker M. E., Audard M., Henning T., Setiawan J., Rodmann J., 2010, A&A, 517, A67

* Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds, Proc. IAU Symp. 210, Modelling of Stellar Atmospheres. Kluwer, Dordrecht, p. A20

* Chiang E. I., Goldreich P., 1997, ApJ, 490, 368

* Cody A. M., Hillenbrand L. A., 2018, AJ, 156, 71

* Cutri R. M. et al., 2012, Explanatory Supplement to the WISE All-Sky Data Release Products

* Davies C. L., 2021, SoftwareX, 14, 100687

* Davies C. L. et al., 2018, ApJ, 866, 23

* Davies C. L., Gregory S. G., Greaves J. S., 2014, MNRAS, 444, 1157

* Davies C. L. et al., 2020, ApJ, 897, 31

* Dent W. R. F., Greaves J. S., Coulson I. M., 2005, MNRAS, 359, 663

* Draine B. T., 2003, ApJ, 598, 1026

* Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237

* Duchêne G. et al., 2010, ApJ, 712, 112

* Dullemond C. P., van den Ancker M. E., Acke B., van Boekel R., 2003, ApJ, 594, L47

* Fairlamb J. R., Oudmaijer R. D., Mendigutía I., Ilee J. D., van den Ancker M. E., 2015, MNRAS, 453, 976

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Friedemann C., Guertler J., Loewe M., 1996, A&AS, 117, 205

* Gaia Collaboration, 2016, A&A, 595, A1

* Gaia Collaboration, 2018, A&A, 616, A1

* Gaia Collaboration, 2021, A&A, 649, A1

* Garufi A. et al., 2018, A&A, 620, A94

* Garufi A. et al., 2020, A&A, 633, A82

* Ginski C. et al., 2016, A&A, 595, A112

* Ginski C. et al., 2021, ApJ, 908, L25

* Glauser A. M., Ménard F., Pinte C., Duchêne G., Güdel M., Monin J. L., Padgett D. L., 2008, A&A, 485, 531

* Gravity Collaboration, 2019, A&A, 632, A53

* Gregorio-Hetem J., Lepine J. R. D., Quast G. R., Torres C. A. O., de La Reza R., 1992, AJ, 103, 549

* Guimarães M. M., Alencar S. H. P., Corradi W. J. B., Vieira S. L. A., 2006, A&A, 457, 581

* Harries T. J., 2000, MNRAS, 315, 722

* Harries T. J., Haworth T. J., Acreman D., Ali A., Douglas T., 2019, Astron. Comput., 27, 63

* Hernández J., Calvet N., Briceño C., Hartmann L., Berlind P., 2004, AJ, 127, 1682

* Høg E. et al., 2000, A&A, 355, L27

* Houck J. R. et al., 2004, ApJS, 154, 18

* Huang P. C. et al., 2019, ApJ, 871, 183

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Isella A., Natta A., 2005, A&A, 438, 899

* Ishihara D. et al., 2010, A&A, 514, A1

* Kama M., Min M., Dominik C., 2009, A&A, 506, 1199

* Kataoka A., Tanaka H., Okuzumi S., Wada K., 2013, A&A, 557, L4

* Keller L. D. et al., 2008, ApJ, 684, 411

* Kenyon S. J., Hartmann L., 1987, ApJ, 323, 714

* Kluska J. et al., 2020, A&A, 636, A116

* Kothe S., Blum J., Weidling R., Güttler C., 2013, Icarus, 225, 75

* Kraus S. et al., 2018, in Creech-Eakman M. J., Tuthill P. G., Mérand A., eds, Proc. SPIE Conf. Ser. Vol. 10701, Optical and Infrared Interferometry and Imaging VI. SPIE, Bellingham, p. 1070123

* Labdon A. et al., 2019, A&A, 627, A36

* Laws A. S. E. et al., 2020, ApJ, 888, 7

* Lazareff B. et al., 2017, A&A, 599, A85

* Lindegren L. et al., 2021, A&A, 649, A4

* Lippa M. et al., 2016, in Malbet F., Creech-Eakman M. J., Tuthill P. G., eds, Proc. SPIE Conf. Ser. Vol. 9907, Optical and Infrared Interferometry and Imaging V. SPIE, Bellingham, p. 990722

* Lucy L. B., 1999, A&A, 344, 282

* Luhman K. L., Herrmann K. A., Mamajek E. E., Esplin T. L., Pecaut M. J., 2018, AJ, 156, 76

* Macintosh B. et al., 2014, Proc. Natl. Acad. Sci., 111, 12661

* Maire J. et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. SPIE, Bellingham, p. 773531

* Mamajek E. E., 2012, ApJ, 754, L20

* McClure M. K. et al., 2013, ApJ, 775, 114

* Meeus G., Waters L. B. F. M., Bouwman J., van den Ancker M. E., Waelkens C., Malfait K., 2001, A&A, 365, 476

* Monnier J. D. et al., 2017, ApJ, 838, 20

* Monnier J. D. et al., 2019, ApJ, 872, 122

* Mulders G. D., Min M., Dominik C., Debes J. H., Schneider G., 2013, A&A, 549, A112

* Murakami H. et al., 2007, PASJ, 59, 369

* Muzerolle J., D'Alessio P., Calvet N., Hartmann L., 2004, ApJ, 617, 406

* Newville M. et al., 2020, lmfit/lmfit-py 1.0.1

* Okuzumi S., Tanaka H., Kobayashi H., Wada K., 2012, ApJ, 752, 106

* Oudmaijer R. D., van der Veen W. E. C. J., Waters L. B. F. M., Trams N. R., Waelkens C., Engelsman E., 1992, A&AS, 96, 625

* Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ, 746, 154

* Perrin M. D. et al., 2014, in Ramsay S. K., McLean I. S., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V. SPIE, Bellingham, p. 91473J

* Perrin M. D. et al., 2015, ApJ, 799, 182

* Petrov P. P. et al., 2019, MNRAS, 483, 132

* Pinte C., Fouchet L., Ménard F., Gonzalez J. F., Duchêne G., 2007, A&A, 469, 963

* Pinte C. et al., 2018, A&A, 609, A47

* Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., Fong W., 1994, ApJ, 421, 615

* Poxon M., 2015, J. Am. Assoc. Var. Obs., 43, 35
Rebull L. M., Stauffer J. R., Cody A. M., Hillenbrand L. A., David T. J.,
Pinsonneault M., 2018, AJ, 155, 196

Ribas Á., Bouy H., Merín B., 2015, A&A, 576, A52

Rich E. A. et al., 2021, ApJ, 913, 138

Ripepi V., Balona L., Catanzaro G., Marconi M., Palla F., Giarrusso M., 2015, MNRAS, 454, 2606

Rizzuto A. C., Ireland M. J., Robertson J. G., 2011, MNRAS, 416, 3108

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Stauffer J. et al., 2015, AJ, 149, 130

Stolker T. et al., 2016a, A&A, 595, A113

Stolker T., Dominik C., Min M., Garufi A., Mulders G. D., Avenhaus H., 2016b, A&A, 596, A70

Tambovtseva L. V., Grinin V. P., 2008, Astron. Lett., 34, 231

Tannirkulam A., Harries T. J., Monnier J. D., 2007, ApJ, 661, 374

Tazaki R., Tanaka H., Muto T., Kataoka A., Okuzumi S., 2019, MNRAS, 485, 4951

ten Brummelaar T. A. et al., 2005, ApJ, 628, 453

Testi L. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. Univ. Arizona Press, Tucson, p. 339

van der Veen W. E. C. J., Waters L. B. F. M., Trams N. R., Matthews H. E., 1994, A&A, 285, 551

Vieira S. L. A., Corradi W. J. B., Alencar S. H. P., Mendes L. T. S., Torres C. A. O., Quast G. R., Guimarães M. M., da Silva L., 2003, AJ, 126, 2971

Villenave M. et al., 2019, A&A, 624, A7

Villenave M. et al., 2020, A&A, 642, A164

Vioque M., Oudmaijer R. D., Baines D., Mendigutía I., Pérez-Martínez R., 2018, A&A, 620, A128

Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2011, ApJ, 737, 36

Wichittanakom C., Oudmaijer R. D., Fairlamb J. R., Mendigutía I., Vioque M., Ababakr K. M., 2020, MNRAS, 493, 234

Wolff S. et al., 2021, AJ, 161, 238

Zwintz K. et al., 2014, Science, 345, 550

## APPENDIX A: PHOTOMETRY

The full list of photometry returned by SEDBYS is provided in Table A1.

<table><caption>**Table A1.** Collated photometry. Individual references are provided in column 3.</caption><thead><tr><th>λ<br>(μm)<br>(1)</th><th>λF<sub>λ</sub><br>(10<sup>-13</sup> W m<sup>-2</sup>)<br>(2)</th><th>Reference<br>(3)</th></tr></thead><tbody><tr><td>0.36</td><td>14</td><td>Vieira et al. (2003)</td></tr><tr><td>0.36</td><td>22</td><td>Lazareff et al. (2017)</td></tr><tr><td>0.43</td><td>51 ± 1</td><td>Høg et al. (2000)</td></tr><tr><td>0.44</td><td>44</td><td>Vieira et al. (2003)</td></tr><tr><td>0.44</td><td>55</td><td>Lazareff et al. (2017)</td></tr><tr><td>0.51</td><td>51.3 ± 0.7</td><td>Gaia Collaboration (2018)</td></tr><tr><td>0.53</td><td>58.8 ± 0.9</td><td>Høg et al. (2000)</td></tr><tr><td>0.54</td><td>61</td><td>Lazareff et al. (2017)</td></tr><tr><td>0.54</td><td>46</td><td>Vieira et al. (2003)</td></tr><tr><td>0.64</td><td>49.5 ± 0.3</td><td>Gaia Collaboration (2018)</td></tr><tr><td>0.65</td><td>56</td><td>Lazareff et al. (2017)</td></tr><tr><td>0.65</td><td>43</td><td>Vieira et al. (2003)</td></tr><tr><td>0.78</td><td>47.4 ± 0.5</td><td>Gaia Collaboration (2018)</td></tr><tr><td>0.79</td><td>39</td><td>Vieira et al. (2003)</td></tr><tr><td>0.79</td><td>48</td><td>Lazareff et al. (2017)</td></tr><tr><td>1.25</td><td>29</td><td>Lazareff et al. (2017)</td></tr><tr><td>1.25</td><td>32.1 ± 0.7</td><td>Skrutskie et al. (2006)</td></tr><tr><td>1.60</td><td>22.5</td><td>Lazareff et al. (2017)</td></tr><tr><td>1.65</td><td>23.1 ± 0.6</td><td>Skrutskie et al. (2006)</td></tr></tbody></table>

<table><caption>**Table A1** – continued</caption><thead><tr><th>λ<br>(μm)<br>(1)</th><th>λF<sub>λ</sub><br>(10<sup>-13</sup> W m<sup>-2</sup>)<br>(2)</th><th>Reference<br>(3)</th></tr></thead><tbody><tr><td>2.15</td><td>19.6 ± 0.4</td><td>Skrutskie et al. (2006)</td></tr><tr><td>2.18</td><td>17.3</td><td>Lazareff et al. (2017)</td></tr><tr><td>3.35</td><td>13.6 ± 0.7</td><td>Cutri et al. (2012)</td></tr><tr><td>4.60</td><td>12.9 ± 0.5</td><td>Cutri et al. (2012)</td></tr><tr><td>9.00</td><td>13.88 ± 0.04</td><td>Ishihara et al. (2010)</td></tr><tr><td>11.60</td><td>11.99 ± 0.09</td><td>Cutri et al. (2012)</td></tr><tr><td>12.00</td><td>14.2</td><td>Oudmaijer et al. (1992)</td></tr><tr><td>18.00</td><td>6.35 ± 0.08</td><td>Ishihara et al. (2010)</td></tr><tr><td>25.00</td><td>7.1</td><td>Oudmaijer et al. (1992)</td></tr><tr><td>60.00</td><td>2.5</td><td>Oudmaijer et al. (1992)</td></tr><tr><td>65.00</td><td>1.7 ± 0.1</td><td>Murakami et al. (2007)</td></tr><tr><td>90.00</td><td>1.51 ± 0.09</td><td>Murakami et al. (2007)</td></tr><tr><td>100.00</td><td>0.9</td><td>Oudmaijer et al. (1992)</td></tr><tr><td>450.00</td><td>&lt;0.06</td><td>van der Veen et al. (1994)</td></tr><tr><td>800.00</td><td>0.0016 ± 0.0006</td><td>van der Veen et al. (1994)</td></tr><tr><td>1100.00</td><td>0.0011 ± 0.0004</td><td>van der Veen et al. (1994)</td></tr><tr><td>1300.00</td><td>0.0011 ± 0.0001</td><td>Garufi et al. (2018)</td></tr></tbody></table>

## APPENDIX B: STELLAR POLARIZATION SUBTRACTION

Improper or incomplete removal of polarization from the instrument, ISM, or unresolved central star can result in $Q_{\phi}$ and $U_{\phi}$ artefacts that can be misinterpreted as disc structure. To remove the stellar polarization, we measure the fractional $Q/I$ ($f_Q$) and fractional $U/I$ ($f_U$). We then multiply this fractional polarization with the intensity cube and subtract the new image from the $Q$ and $U$ frames, respectively (see e.g. Laws et al. 2020). The instrumental polarization will vary with the parallactic angle. Thus, stellar and instrumental polarization are removed for each individual Stokes frame where the change in parallactic angle is minimal (Perrin et al. 2015). Since the stellar and instrumental polarization are constant values of $Q$ and $U$ for a given observation cycle, the Stokes images are rotated into the $Q_{\phi}$ and $U_{\phi}$ frames. Any residual stellar or instrumental polarization will form a quadrupole pattern and will be seen in the image (first column of Fig. B1).

Measuring $f_Q$ and $f_U$ can be complicated as we need to disentangle the stellar and instrumental polarization signal from the unknown disc structure. We explored four different methods to measure the $f_Q$ and $f_U$ with the resulting mean-combined $Q_{\phi}$ and $U_{\phi}$ images shown in Fig. B1. We demonstrate these methods with the J-band data only. The H-band data showed similar results.

### B1 Method 1

This method is used in the GPI Data Reduction Pipeline (DRP) v1.5. All of the $Q$, $U$, and $I$ counts within the coronagraphic spot (< 6 pix) are used and the calculated $f_Q$ and $f_U$ values are based on the summed $Q$, $U$, and $I$ values. The resultant image (first column of Fig. B1) has a strong quadrupole signal, highlighting that this method is insufficient to remove the stellar and instrumental value. Laws et al. (2020) drew similar conclusions: using the light within the coronagraphic spot does not result in the best subtraction of the stellar polarization.

### B2 Method 2

We employed an algorithm that utilizes the full field of view of the image while masking out individual pixels of the image where the
<figure><img src="image_16.png" /></figure>

**Figure B1.** J-band mean-combined $Q_{\phi}$ (top row) and $U_{\phi}$ images (bottom row) for the four different stellar and instrumental polarization removal techniques (see text for details). All images are plotted with the same symlog spread and smoothed with a Gaussian kernel to highlight the quadrupole pattern in the background seen in the first three columns.

disc is bright ($f_Q$ or $f_U$ > 0.05) or where the signal-to-noise ratio in the pixel is low. All the non-masked pixels are summed, giving summed values of Q, U, and I from which to calculate $f_Q$ and $f_U$. This is the main method employed by the G-LIGHTS team (Rich et al., in preparation). The resulting $Q_{\phi}$ and $U_{\phi}$ images are shown in the second column of Fig. B1. The quadrupole pattern is not as strong as for Method 1 but it is still present in the $U_{\phi}$ image. Additionally, the $Q_{\phi}$ image appears to indicate the presence of scattering along the apparent disc minor axis.

## B3 Method 3

Laws et al. (2020) summed the Q, U, and I intensities in a $70 < r < 80$ pixel radius around the target star. The resultant mean-averaged $Q_{\phi}$ and $U_{\phi}$ images are shown in the third column of Fig. B1. These show a further reduction in quadrupole structure, though some is still visible in the $U_{\phi}$ image. However, the negative $Q_{\phi}$ flux values along the apparent disc minor axis produced using this method are not expected to be real.

## B4 Method 4

This involves by-eye selection of the best values of $f_Q$ and $f_U$ to simultaneously minimize the quadrupole structure in the outer portions of the image and remove the positive/negative $Q_{\phi}$ flux along the apparent disc minor axis in the eight individual Stokes $Q_{\phi}$ and $U_{\phi}$ frames. The rationale here was to see whether the apparent minor axis structure found in the $Q_{\phi}$ image when using Method 2 (column 2 of Fig. B1) could solely be a result of residual stellar and instrumental polarization. Ultimately, we can find a combination of $f_Q$ and $f_U$ for each of the eight Stokes frames that fully remove the residual quadrupole pattern and the minor axis feature with this method, as shown in the final column of Fig. B1. As this method results in no extraneous structures in the image, we use this by-eye minimization of $Q_{\phi}$ and $U_{\phi}$ to remove stellar and instrumental polarization herein.

In Fig. B2, we plot the stellar and instrumental polarization angle and the percent polarization for each of the eight Stokes cycles using Method 3 and Method 4. We measure an average percent polarization of 1.0 per cent and an average polarization angle of 102°. We note

<figure><img src="image_17.png" /></figure>

**Figure B2.** Percent polarization (top panel) and polarization angle (lower panel) of stellar and instrumental polarization removed for each of the eight cycles of observations (method three is shown in blue; method four in orange). The average values for each of the methods are shown as dashed lines and the average values are written in the upper left of the panel.

that both the polarization angle and percent polarization have strong deviations from the average with Method 4 resulting in standard deviations of 0.12 per cent and 5.1°, respectively. We also see that the average per cent polarization and polarization angle for both methods do not differ wildly (0.1 per cent and 1°) while some individual Stokes frames have very different polarization values (e.g. Stokes cycle 3 and 6: Fig. B2). Deviations from the average are expected as the magnitude and direction of the instrument polarization changes as a function of parallactic angle. However, for continuous observations such as these, the stellar and instrumental polarization should follow a sinusoidal function while the values shown in Fig. B2 do not.

One potential explanation for the change in variation is related to HD 145718 photometric variability. While this effect has not been studied in great detail to the best of our knowledge, dust from
circumstellar occultation events should have a polarization signal from surface scattering which would vary with orbital phase. Indeed, Perrin et al. (2015) concluded from their instrumental polarization measurements of β Pic that GPI should have a polarization accuracy of 0.1 per cent, slightly below the 0.12 per cent standard deviation we measure for HD 145718. However, since there has been no systematic modelling of the instrument polarization for GPI and our observations do not trace a sufficiently long time-scale to allow us to undertake a similar analysis for HD 145718 as Perrin et al. (2015) did for β Pic, we cannot definitively conclude whether we are, in fact, observing polarization variation from dust obscuration events.

## APPENDIX C: RE-ESTIMATING STELLAR LUMINOSITY AND VISUAL EXTINCTION USING BRIGHT EPOCH PHOTOMETRY

We followed Fairlamb et al. (2015) and computed a grid of reddened Castelli & Kurucz (2003) model atmospheres with $T_{\text{eff}} = 7750, 8000$, and $8250$ K; $\log(g) = 4.37$; and $A_V$ ranging from 0.50 to 2.50 in steps of 0.01 mag. The value of $\log(g)$ makes little difference to the fit so its value is not changed. Each reddened model was fit to the Lazareff et al. (2017) BVRI photometry and the best-fitting model was identified from the minimum of all the $\chi^2$ values. This procedure was run twice: once using the Cardelli, Clayton & Mathis (1989) reddening law with total-to-selective extinction, $R_V = 3.1$ and once using $R_V = 5.0$. We also re-ran the fitting to the fainter epoch photometry from Vieira et al. (2003) to highlight the differences we observe. In each case, we follow Wichittanakom et al. (2020) and use the models providing $\chi^2$ values twice the minimum-$\chi^2$ to estimate the uncertainty in $A_V$.

The reddened model atmosphere providing the best fit to the BVRI photometry was then scaled to the V-band photometry. This scaling factor – which accounts for fitting models of surface flux to photometry – corresponds to $(d/R_*)^2$. We use this and the Gaia DR2 stellar distance to estimate the stellar radius, $R_*$, The quoted uncertainties on $R_*$ take into account the uncertainties on $A_V$ and on d. Finally, we use the Stefan–Boltzmann law to calculate $L_*$ from $T_{\text{eff}}$ and $R_*$, assuming an effective temperature for the Sun of 5771.8 K (Mamajek 2012). Our results are presented in Table C1 and displayed in Fig. C1.

The models reddened using $R_V = 3.1$ provide an improved fit to the Lazareff et al. (2017) photometry over those reddened using $R_V = 5.0$. Conversely, a better fit to the Vieira et al. (2003) photometry is provided by the models reddened using $R_V = 5.0$.

<table><caption>**Table C1.** Results from fitting Castelli & Kurucz (2003) stellar atmospheres with to BVRI photometry. Column 1: the source of the photometry; column 2: the adopted total-to-selective extinction; columns 3 and 4: the best-fitting $A_V$ and $R_*$, assuming $d = 152.5$ pc; column 5: inferred $L_*$.</caption><thead><tr><th>Photometry source<br>(1)</th><th><i>R</i><sub>V</sub><br>(mag)<br>(2)</th><th><i>A</i><sub>V</sub><br>(mag)<br>(3)</th><th><i>R</i><sub>*</sub><br>(<i>R</i><sub>☉</sub>)<br>(4)</th><th><i>L</i><sub>*</sub><br>(<i>L</i><sub>☉</sub>)<br>(5)</th></tr></thead><tbody><tr><td>Vieira et al. (2003)</td><td>3.1</td><td>0.89<sup>+0.31</sup><sub>-0.05</sub></td><td>1.70<sup>+0.08</sup><sub>-0.07</sub></td><td>10.7<sup>+2.6</sup><sub>-2.0</sub></td></tr><tr><td>Vieira et al. (2003)</td><td>5.0</td><td>1.22<sup>+0.04</sup><sub>-0.03</sub></td><td>1.98<sup>+0.07</sup><sub>-0.08</sub></td><td>14.5<sup>+3.1</sup><sub>-2.7</sub></td></tr><tr><td>Lazareff et al. (2017)</td><td>3.1</td><td>0.89<sup>+0.34</sup><sub>-0.08</sub></td><td>1.97<sup>+0.12</sup><sub>-0.11</sub></td><td>14.3<sup>+3.9</sup><sub>-3.1</sub></td></tr><tr><td>Lazareff et al. (2017)</td><td>5.0</td><td>1.21<sup>+0.52</sup><sub>-0.21</sub></td><td>2.28<sup>+0.35</sup><sub>-0.25</sub></td><td>19.2<sup>+9.7</sup><sub>-5.8</sub></td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure C1.** Comparison of Castelli & Kurucz (2003) model atmospheres with two epochs of previously published BVRI photometry for HD 145718. The solid lines show our best fits to the Vieira et al. (2003) photometry (filled circles) while dashed lines show our best fit to the Lazareff et al. (2017) photometry (filled diamonds). In the top and bottom plots, we show the fits adopting $R_V = 3.1$ and $R_V = 5.0$, respectively (see Table C1).</figcaption></figure>

This paper has been typeset from a TeX/LATeX file prepared by the author.