# High-resolution ALMA and HST images of q1 Eri: an asymmetric debris disc with an eccentric Jupiter

J Lovell, S Marino, M Wyatt, G Kennedy, M Macgregor, K Stapelfeldt, B

Dent, J Krist, L Matrà, Q Kral, et al.

► To cite this version:

J Lovell, S Marino, M Wyatt, G Kennedy, M Macgregor, et al. High-resolution ALMA and HST images of q1 Eri: an asymmetric debris disc with an eccentric Jupiter. Monthly Notices of the Royal Astronomical Society, 2021, 506 (2), pp.1978-2001. 10.1093/mnras/stab1678. hal-03660729

HAL Id: hal-03660729

https://hal.science/hal-03660729

Submitted on 6 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# High-resolution ALMA and HST images of q<fnref n="1" /> Eri: an asymmetric debris disc with an eccentric Jupiter

J. B. Lovell<sup>1,★</sup> S. Marino<sup>1,2</sup> M. C. Wyatt<sup>1</sup> G. M. Kennedy<sup>1,3</sup> M. A. MacGregor<sup>4</sup> K. Stapelfeldt<sup>5</sup> B. Dent<sup>6</sup> J. Krist<sup>5</sup> L. Matrà<sup>7</sup> Q. Kral<sup>8</sup> O. Panić<sup>9</sup> T. D. Pearce<sup>10</sup> and D. Wilner<sup>11</sup>

Accepted 2021 June 1. Received 2021 May 27; in original form 2021 March 19

## ABSTRACT

We present Atacama Large Millimetre/sub-Millimetre Array (ALMA) 1.3 and 0.86 mm observations of the nearby (17.34 pc) F9V star q1 Eri (HD 10647, HR 506). This system, with age ~1.4 Gyr, hosts a ~2 au radial velocity planet and a debris disc with the highest fractional luminosity of the closest 300 FGK type stars. The ALMA images, with resolution ~0.′5, reveal a broad (34–134 au) belt of millimetre emission inclined by 76.7 ± 1.0 degrees with maximum brightness at 81.6 ± 0.5 au. The images reveal an asymmetry, with higher flux near the south-west ansa, which is also closer to the star. Scattered light observed with the *Hubble Space Telescope* is also asymmetric, being more radially extended to the north-east. We fit the millimetre emission with parametric models and place constraints on the disc morphology, radius, width, dust mass, and scale height. We find that the south-west ansa asymmetry is best fitted by an extended clump on the inner edge of the disc, consistent with perturbations from a planet with mass 8 $M_⊕$–11 $M_{Jup}$ at ~60 au that may have migrated outwards, similar to Neptune in our Solar system. If the measured vertical aspect ratio of $h = 0.04 ± 0.01$ is due to dynamical interactions in the disc, then this requires perturbers with sizes >1200 km. We find tentative evidence for a 0.86 mm excess within 10 au, $70 ± 22$ µJy, that may be due to an inner planetesimal belt. We find no evidence for CO gas, but set an upper bound on the CO gas mass of $4 × 10^{-6} M_⊕$ ($3σ$), consistent with cometary abundances in the Solar system.

**Key words:** techniques: interferometric–planets and satellites: dynamical evolution and stability–(stars:) circumstellar matter–stars: individual: HD 10647.

## 1 INTRODUCTION

The first detections of exoplanetary systems (systems with planets and/or planetesimal belts) were made a few decades ago (Aumann et al. 1984; Harper, Loewenstein & Davidson 1984; Wolszczan & Frail 1992; Mayor & Queloz 1995). Many hundreds of systems are now known to possess planetesimal belts and many thousands to possess planets.<fnref n="1" /> Most of these planetesimal belts are cold and reside at 10s of au, making them analogues of the Kuiper Belt in our own Solar system. They are inferred to exist from observations of emission from dust which is commonly seen to lie in belts around their

parent stars and must have been created in collisions between larger planetesimals. Such circumstellar dust and the implied planetesimals collectively form what is known as a star’s debris disc (Wyatt 2008; Hughes, Duchêne & Matthews 2018).

The number of debris discs is continually growing, and several of these are near enough and sufficiently bright to allow high-resolution imaging of the discs’ structure and substructure with mm/sub-mm instruments such as the Atacama Large Millimetre/sub-Millimetre Array (ALMA; see Matrà et al. 2018, and references therein), and at shorter optical wavelengths, for example with the *Hubble Space Telescope* (HST; see Apai et al. 2015). In many cases the observed morphologies are believed to be influenced by the presence of large perturbing bodies (e.g. planets), and many systems have been directly observed with multiple planets (e.g. HR 8799 with 4 planets, Marois et al. 2008; Su et al. 2009; Marois et al. 2010).

© 2021 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2"><sup>2</sup>Jesus College, University of Cambridge, Jesus Lane, Cambridge CB5 8BL, UK</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="4"><sup>4</sup>Department of Astrophysical and Planetary Sciences, University of Colorado, 2000 Colorado Avenue, Boulder, CO 80309, USA</footnote>

<footnote n="5"><sup>5</sup>Jet Propulsion Laboratory, M/S 321-161, 4800 Oak Grove Drive, Pasadena, CA 91109 USA</footnote>

<footnote n="6"><sup>6</sup>ESO, Alonso de Córdova 3107, Vitacura, Región Metropolitana, 763 0355, Chile</footnote>

<footnote n="7"><sup>7</sup>School of Physics, National University of Ireland Galway, University Road, Galway H91 TK33, Ireland</footnote>

<footnote n="8"><sup>8</sup>LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="9"><sup>9</sup>School of Physics and Astronomy, University of Leeds, Woodhouse, Leeds LS2 9JS, UK</footnote>

<footnote n="10"><sup>10</sup>Astrophysikalisches Institut, Friedrich-Schiller-Universität Jena, Schillergaesschen 2-3, D-07745 Jena, Germany</footnote>

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

<footnote>*E-mail: jl638@cam.ac.uk</footnote>

<footnote n="1">¹http://exoplanet.eu</footnote>
Disc-planet interactions are expected to produce a wide variety of detectable morphologies such as clumps, radial offsets from their stars, spirals and brightness asymmetries, and characterizing these can place important constraints on the architecture and evolution of the entire planetary system (Wyatt et al. 1999; Wyatt 2003; Lee & Chiang 2016; Faramaz et al. 2019).

The planets that have been observed around other stars fit into a few different classes based on their masses, radii and semimajor axes. These include objects such as short period super Earths found during transits, outer giants observed by direct imaging, and exo-Jupiters discovered by radial velocity measurements. It might be expected that the formation mechanism of different planets is reflected in the properties of their debris discs, and recent studies have explored such connections for Super Earth systems like 61 Vir (Marino et al. 2017), and for systems with known outer giants like HR 8799 and Beta Pic (see Booth et al. 2016; Matrà et al. 2017), and for populations of directly imaged giant planets (see Meshkat et al. 2017).

Exo-Jupiters have masses and semimajor axes respectively in the ranges $0.1-3 M_{Jup}$ and $1-5$ au, are observed around ~5 per cent of stars, and are commonly found to have eccentric orbits (Chiang & Laughlin 2013). The origin of their eccentric orbits is understood to arise from early stage instabilities in planetary systems which scatter planets into a broad distribution of eccentricities (Jurić & Tremaine 2008). Planetary system instabilities that excite exo-Jupiter eccentricities have been shown to deplete outer debris discs (those with radii greater than 10s of au), and in turn these outer discs have been shown to dynamically influence closer in planets (Gomes et al. 2005; Raymond et al. 2011, 2012). Whilst studies of mutual disc-planet interactions might lead us to conclude that debris belts may be more readily depleted if exo-Jupiters are present, there appears to be neither positive nor negative correlation between the presence of exo-Jupiters in systems with debris discs (Bryden et al. 2009; Moro-Martín et al. 2015; Yelverton, Kennedy & Su 2020).

$q^1$ Eri (HD 10647, HR 506) is an old (1.4 Gyr), nearby main-sequence F9V star, and is an example of a system with both a bright debris disc and an exo-Jupiter, $q^1$ Eri b (Marmier et al. 2013). The $q^1$ Eri system is at a distance of $17.344\pm0.014$ pc (Gaia Collaboration 2016, 2018), and the debris disc has the highest fractional luminosity of the closest 300 sun-like (FGK-type) stars, $f \ge 10^{-4}$ (Liseau et al. 2008; Sibthorpe et al. 2018). The planet $q^1$ Eri b has a semimajor axis $a = 2.03 \pm 0.15$ au, with a mass $M \sin(i) = 0.93 \pm 0.18 M_J$, and a low-moderate eccentricity of $e = 0.15 \pm 0.08$ (Marmier et al. 2013).

The disc of $q^1$ Eri has been resolved both in the optical (scattered light) with HST and in the far-IR (thermal emission) with Herschel. These studies showed the disc to be highly inclined ($i > 60^\circ$) to the plane of the sky (Stapelfeldt et al. 2007; Liseau et al. 2008, 2010) with emission concentrated at around ~100 au, although only weak constraints could be placed on the disc's inner edge based on the ~4 arcsec (~70 au) resolution of Herschel, and given the HST coronagraph obscuration inside ~50 au. More recent analysis with this Herschel data has demonstrated that the disc likely has a radius of $R_{disc} = 81.1_{-1.3}^{+1.8}$ au, and a broad radial width of $\Delta R_{disc} = 71.1_{-13.3}^{+1.9}$ au (Marshall et al. 2021). The HST data confirms the disc to be asymmetric, with this more extended in the NE than the SW. In addition to the belt beyond 70 au, the flux distribution indicates that an inner warm component at ~10 au is also present (Kennedy & Wyatt 2014; Schüppler et al. 2016). This component is yet to be resolved, but may be in close proximity to the exo-Jupiter $q^1$ Eri b.

In this work we present new high-resolution (subarcsec) ALMA observations of $q^1$ Eri to characterize its asymmetric outer belt and to consider the relationship of that belt to the planetary system architecture, in particular to the known exo-Jupiter. These measure-

<table><caption>**Table 1.** ALMA observational set-up, over the six different epochs. The horizontal lines separate the (top) Band 6, (middle) 2016 Band 7, and (bottom) 2018 Band 7 observation sets. All times represent time on source.</caption><thead><tr><th>D.M.Y</th><th>Time [min]</th><th>Flux calibrator</th><th>Bandpass calibrator</th><th>Phase calibrator</th></tr></thead><tbody><tr><td>26.05.16</td><td>39:58</td><td>Pallas</td><td>J2258-2758</td><td>J0210-5101</td></tr><tr><td>02.06.16</td><td>39:58</td><td>Pallas</td><td>J2258-2758</td><td>J0210-5101</td></tr><tr><td>12.07.16</td><td>48:34</td><td>Ceres</td><td>J2258-2758</td><td>J0210-5101</td></tr><tr><td>12.07.16</td><td>48:34</td><td>Ceres</td><td>J0538-4405</td><td>J0210-5101</td></tr><tr><td>27.06.18</td><td>39:58</td><td>J2258-2758</td><td>J2258-2758</td><td>J0124-5113</td></tr><tr><td>07.07.18</td><td>39:57</td><td>J0159-4546</td><td>J0159-4546</td><td>J0210-5101</td></tr></tbody></table>

ments were taken over three epochs, covering two wavelengths; one Cycle 3 (2016) observation in Band 6 (~1.25 mm) and two Band 7 observations (~856 µm) in Cycles 3 (2016) and 5 (2018). Combining these data provides a resolution of ~0.8 arcsec, five times better than that achieved by Herschel, sufficient to constrain the inner edge of the debris disc and observe the previously unresolved inner regions. In addition we present 2006 HST data, as included in Stapelfeldt et al. (2007), to complement our analysis of the mm/sub-mm observations.

We provide an overview of all our observational data sets in Section 2 and discuss our initial observational analysis of these in Section 3. We then present our methodology to model this system in Section 4, give a discussion of our findings and future work in Section 5 and summarize our key conclusions in Section 6.

# 2 OBSERVATIONS

## 2.1 Submillimetre ALMA observations

### 2.1.1 ALMA Band 6 observations

$q^1$ Eri was observed for ~80 min (on source) in two scheduling blocks with ALMA in Band 6 during Cycle 3 as shown in the top panel of Table 1, using 41 antennas with minimum and maximum baselines ranging from 15.1 to 772.8 m as part of project 2015.1.00307.S (PI: David Wilner). The correlator had three spectral windows centred on frequencies of 232.490, 244.989, and 246.989 GHz, for continuum observations with a bandwidth of 2.000 GHz and channel widths of 15.625 MHz. Also set-up was a spectral window centred on a topocentric frequency of 230.536 GHz with a bandwidth of 1.875 GHz, and channel widths of 488.281 kHz, for CO J = 2-1 spectral line observations. The visibility data set was calibrated using the CASA software version 5.1.1-5 with the standard pipeline provided by the ALMA Observatory. Additional data flagging was performed on the Band 6 data to mitigate issues introduced by several poorly performing antennas. The plotms task in CASA was used to examine the visibility amplitudes as a function of time and uv-distance. Outlying points were flagged using the flagdata task. Table 1 shows a summary of the observational set-up for ALMA data collection. Continuum imaging was conducted using the CASA tclean algorithm with natural weighting (to enhance S/N), shown in Fig. 1 (bottom) in which the disc is clearly detected. The synthesized beam size in this image is $0.67 \text{ arcsec} \times 0.54 \text{ arcsec} (\text{PA}_{\text{Beam}}=69.7^\circ)$, which at 17.34 pc corresponds to a physical size of $11.6 \times 9.4$ au.

### 2.1.2 ALMA Band 7 observations

$q^1$ Eri was observed for ~177 min (on source) by ALMA in Band 7 over four scheduling blocks in Cycles 3 and 5 as shown
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Non-primary beam corrected ALMA images of $q^1$ Eri in Band 7 (0.86 mm, top) and Band 6 (1.25 mm, bottom), with the synthesized beams shown in the lower left in white, and the stellar location marked with a red star (both in the coordinate centres). In both North is up and East is left. Unresolved emission is detected in the centre of these figures deemed to be mostly from the star. Top: A cleaned Briggs weighted (robust = 1.0) image for the combined Band 7 2016 and 2018 data sets, with contours of $+4\sigma$, $+6\sigma$, $+8\sigma$, and $+10\sigma$ significance. The beam is $0.55$ arcsec $\times$ $0.47$ arcsec, $PA_{\text{beam}} = -66.9^\circ$, and the image rms is $\sigma=13.9$ µJy beam$^{-1}$. Bottom: A cleaned natural weighted $\sim 0.4$ arcsec uv-tapered image in Band 6, with contours of $+3\sigma$, $+5\sigma$, $+7\sigma$, and $+9\sigma$ significance. The beam is $0.67$ arcsec $\times$ $0.54$ arcsec, $PA_{\text{beam}} = -67.6^\circ$, and the image rms is $\sigma=13.3$ µJy beam$^{-1}$.</figcaption></figure>

in the middle and lower panels of Table 1, with baselines ranging from 15.1 to 867.2m (2016) and 15.0 to 313.7m (2018). In 2016 (project 2015.1.01260.S, PI: Mark Wyatt) using 33 antennas, the correlator had 3 spectral windows centred on frequencies of 347.817, 335.775, and 333.817 GHz each with a bandwidth of 2.0 GHz and channel widths of 15.625 MHz (a total of 128 channels each, for continuum observations). Also set-up was a spectral window centred on a topocentric frequency of 345.817 GHz with a bandwidth of 1.875 GHz and spectral channel widths of 488.281 kHz (a total of 3840 channels), for CO J = 3–2 spectral line observations. In 2018 (project 2017.1.00167.S, PI: Mark Wyatt) using 46 antennas (in observation block 1) and 43 antennas (in observation block 2), the same correlator set-up (3 continuum, 1 spectral line) and respective channel widths were used, however the central frequencies for the continuum observations were 347.683, 335.683, and

<table><caption>**Table 2.** Stellar position of $q^1$ Eri for the six observational epochs, calculated based on Gaia DR2 data (Gaia Collaboration 2018), for which DR2 positional errors are < 60 µas.</caption><thead><tr><th>DD.MM.YY</th><th>RA</th><th>Dec</th></tr></thead><tbody><tr><td>26.05.16</td><td>01:42:29.621</td><td>-53:44:28.722</td></tr><tr><td>02.06.16</td><td>01:42:29.622</td><td>-53:44:28.723</td></tr><tr><td>12.07.16<sup>a</sup></td><td>01:42:29.624</td><td>-53:44:28.735</td></tr><tr><td>27.06.18</td><td>01:42:29.660</td><td>-53:44:28.942</td></tr><tr><td>07.07.18</td><td>01:42:29.661</td><td>-53:44:28.945</td></tr></tbody></table>

333.788 GHz, and the spectral line observation central frequency (topocentric) was 345.787 GHz (still covering the line transition frequency).

The 2016 Band 7 visibility data sets were calibrated using the CASA software version 4.7.2 with the standard pipeline provided by the ALMA Observatory, whilst the 2018 Band 7 visibility data sets were calibrated using the CASA software version 5.1.1-5 (also with the standard pipeline provided by the ALMA Observatory). Both Band 7 data sets were time averaged to 10-s widths, channel averaged into four channels per spectral window, and the CASA statwt task was run to estimate the visibility weights based on the measured variance. We note here that the Band 7 data in 2016 got a 'semi pass' rating for its quality assurance, since a smaller synthesized beam was obtained in 2016 than requested.

To improve the S/N and resolution of the Band 7 data, we combined the two Band 7 data sets. Between the 2016 and 2018 observations the proper motion of $q^1$ Eri, $\mu_{RA, Dec} = (165.83 \pm 0.10, -105.52 \pm 0.10)$ mas yr$^{-1}$, resulted in its position shifting by $\sim 0.4$ arcsec (Gaia Collaboration 2016, 2018). Since the ALMA observation phase centres were not perfectly aligned with the expected stellar position from Gaia DR2, we realigned to these coordinates using the CASA task fixvis. We then used the CASA task fixplanets to shift the coordinates of the 2018 observations to coincide with the 2016 measurement set phase centre (see Table 2 for these expected DR2 coordinates). The two resulting measurement set epochs were then combined with the CASA task concat. We note that correcting for the stellar proper motion in this way would result in background sources (if present) being smeared along the proper motion direction, however we find no evidence of any point sources which appear significantly elongated (i.e. the emission is consistent with a point source at the location of the star). Fig. 1 (top) shows the resulting Band 7 cleaned image of $q^1$ Eri with the CASA tclean algorithm with Briggs weighting (robust = 1.0), in which the beam size is $0.55$ arcsec $\times$ $0.47$ arcsec (i.e. $9.5 \times 8.1$ au). This weighting was selected here to achieve a beam size similar to the Band 6 image (naturally weighted, Fig. 1 bottom).

The star should be at the centre of the image but has a positional uncertainty of $\pm 0.05$ arcsec, due to the combination of the astrometric accuracy of Gaia DR2 ($< 60$ µas) and ALMA (which we assume to be <10 per cent of the synthesized beam based on the ALMA pipeline weblog 'timegaincal' plots; given our observation routine, choice of calibrators, and weather conditions, a reasonable target solution could be found). The peak detected at the centre of the image has a S/N=10.4 and therefore could be from a source offset from its true position by $0.065$ arcsec, i.e. due to noise and systematic errors (see equation 10.7 in Remijan et al. 2019). As such the source detected at the centre of the image is consistent with the stellar position.

<footnote><sup>a</sup>Note that all Band 7 2016 data were collected on 12.07.16, so only one position of $q^1$ Eri is reported for this epoch.</footnote>
<table><caption>**Table 3.** HST observing set-up, taken over the two consecutive orbits. Subscripts 1 and 2 in 'F606W₁' and 'F606W₂' refer to the F606W filter being used for two different exposure times within a single observation sequence. The F502N automated acquisition exposures were taken with the star behind the coronagraph's occulting spot for 0.1 s each, whilst the four F606W₂ exposures were taken with the star centred behind the 1.8 arcsec diameter occulter.</caption><thead><tr><th>Orbit</th><th>F606W<sub>1</sub> (s)</th><th>F606W<sub>2</sub> (s)</th><th>Roll (°)</th></tr></thead><tbody><tr><td>1</td><td>540 (x4) = 2160</td><td>575 (x4) = 2300</td><td>0.0</td></tr><tr><td>2</td><td>540 (x4) = 2160</td><td>575 (x4) = 2300</td><td>24.9</td></tr></tbody></table>

## 2.2 Scattered light HST observations

$q^1$ Eri was observed by the coronagraph onboard the HST's Advanced Camera for Surveys (ACS) High Resolution Camera (HRC; 0.25 arcsec pixel⁻¹), as part of HST program 10539 (PI: Karl Stapelfeldt), on 2nd September, 2006. Two observation sequences were made, as outlined in Table 3, with the F606W filter (central wavelength: 5887 Å, FWHM: 1566 Å). Each sequence took one orbit in time, and the two sequences were made in consecutive orbits. Between the two sequences, the telescope was rolled about the line-of-sight axis by ~25° to provide a mechanism to correct for systematic errors (such as instrumental artefacts) which can then be distinguished from real objects in the data. The ACS coronagraph reduces the wings of stellar diffraction patterns caused by the HST aperture and obscurations. Any remaining halo then seen around the star is caused by optical surface errors (from scattering) that the coronagraph does not suppress.

Since two images of $q^1$ Eri were collected at different roll angles, we used one as the reference point spread function (PSF) for the other. As the PSF is fixed on the detector, but the disc rotates with the roll angle, any sufficiently bright halo would appear as positive and negative signals when one image is subtracted from the other. While a face-on disc would 'self-subtract' with this method, $q^1$ Eri is sufficiently inclined for this method to be effective. An algorithm was applied to iteratively solve for the sky and PSF images (applied previously to detect debris discs by Krist 2004), producing a final sky image with residual haloes. The final science image shown in Fig. 2 is median filtered to remove cosmic rays, binned to 0.05 arcsec pix⁻¹ sampling, and shows the contour lines at the levels of +22 and +21 mag arcsec⁻². This image clearly demonstrates the extent of $q^1$ Eri's debris disc (and has been reduced identically as per Stapelfeldt et al. 2007) which we analyse further in Section 3.2.

## 2.3 Flux density distribution

Fig. 3 shows the flux density distribution of $q^1$ Eri, including fluxes derived from the ALMA data in Section 3. Values for the wavelengths and fluxes used to produce this are given in Table 4. The distribution was fitted using the methodology of Kennedy & Wyatt (2014) and Yelverton et al. (2019), finding a two-component temperature distribution (an inner warm component, and outer cool component) along with a stellar effective temperature of $T_{eff}=6100 \pm 100$ K, and a stellar luminosity and radius (with 2 per cent calibration uncertainties) of $L=1.55 L_⊙$, and $R=1.11 R_⊙$ respectively. We note that Gaia Collaboration (2018) estimated these stellar parameters (using the 'Priam' and 'FLAME' algorithms) as 6143 K, 1.59 $L_⊙$, and 1.11 $R_⊙$ respectively, all broadly in agreement.

The two-component modeRiciil finds an inner warm belt and an outer cool belt, respectively with blackbody temperatures of $T_{bb,inner}=101 \pm 4$ K and $T_{bb,outer}=41 \pm 1$ K (corresponding to uncor-

<figure><img src="image_4.png" /><figcaption>**Figure 2.** HST image of $q^1$ Eri, following F606W filtered observations (see details in Table 3), showing contours for +22 and +21 mag arcsec⁻² surface brightnesses. The inner red dashed ring indicates the scale of the coronagraph. The two outer dark-red dashed rings are at 5.7 and 6.7 arcsec from the stellar centre, showing that the emission in the NE is more extended than the SW along the major axis. To aid viewing the major axis, the black dash-dotted line is shown with a position angle of 56°. In this figure, North is up and East is left.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The flux density distribution of $q^1$ Eri, demonstrating the need to model this system with stellar emission, and both a cool (outer) and a warm (inner) component. The warm and cool emission profiles (dotted lines) are modelled as modified blackbodies. The green dash-dotted line shows the combination of both warm and cool components, and the lilac dash-dotted line the combination of the stellar emission and warm component. Note that these data include the derived flux from the ALMA measurements outlined in Section 3, and between ~5 and 36 µm the *Spitzer* IRS spectra.</figcaption></figure>

rected blackbody radii of 9.5 ± 0.7 au and 57 ± 3 au). Respectively these components have fractional luminosities of 6.3 ± 0.4 × 10⁻⁵ and 2.44 ± 0.08 × 10⁻⁴. Although due to non-photospheric emission and the possibility of additional unresolved circumstellar material such as hot dust being located near the star (and thus a strict
<table><caption>**Table 4.** Values for the fluxes used to produce the flux density distribution in Fig. 3.</caption><thead><tr><th>Source</th><th>Wavelength [µm]</th><th>Flux / magnitude</th><th>Units</th></tr></thead><tbody><tr><td>U - B<sup>a</sup></td><td>-</td><td>0.00 ± 0.12</td><td>mag</td></tr><tr><td>C<sub>1</sub><sup>b</sup></td><td>-</td><td>0.336 ± 0.012</td><td>mag</td></tr><tr><td>B<sub>T</sub><sup>b</sup></td><td>0.420</td><td>6.164 ± 0.014</td><td>mag</td></tr><tr><td>M<sub>1</sub><sup>b</sup></td><td>-</td><td>0.168 ± 0.008</td><td>mag</td></tr><tr><td>B - V<sup>a</sup></td><td>-</td><td>0.530 ± 0.026</td><td>mag</td></tr><tr><td>B - Y<sup>b</sup></td><td>-</td><td>0.354 ± 0.008</td><td>mag</td></tr><tr><td>B<sub>P</sub><sup>d</sup></td><td>0.513</td><td>5.6770 ± 0.0017</td><td>mag</td></tr><tr><td>V<sub>T</sub><sup>d</sup></td><td>0.532</td><td>5.581 ± 0.009</td><td>mag</td></tr><tr><td>H<sub>P</sub><sup>d</sup></td><td>0.542</td><td>5.638 ± 0.006</td><td>mag</td></tr><tr><td>V<sup>a</sup></td><td>0.550</td><td>5.540 ± 0.019</td><td>mag</td></tr><tr><td>G<sup>d</sup></td><td>0.642</td><td>5.3550 ± 0.0013</td><td>mag</td></tr><tr><td>R<sub>P</sub><sup>d</sup></td><td>0.780</td><td>4.9660 ± 0.0022</td><td>mag</td></tr><tr><td>J<sup>f</sup></td><td>0.791</td><td>8.757 ± 0.020</td><td>mag</td></tr><tr><td>J<sup>f</sup></td><td>1.24</td><td>4.79 ± 0.23</td><td>mag</td></tr><tr><td>H<sup>f</sup></td><td>1.65</td><td>4.40 ± 0.23</td><td>mag</td></tr><tr><td>K<sub>S</sub><sup>f</sup></td><td>2f.16</td><td>4.34 ± 0.28</td><td>mag</td></tr><tr><td>WISE W1<sup>g</sup></td><td>3.38</td><td>4.17 ± 0.38</td><td>mag</td></tr><tr><td>WISE W2<sup>g</sup></td><td>4.63</td><td>3.91 ± 0.22</td><td>mag</td></tr><tr><td>AKARI IRC<sup>i</sup></td><td>8.98</td><td>1.151 ± 0.019</td><td>Jy</td></tr><tr><td>IRAS<sup>i</sup></td><td>11.2</td><td>0.82 ± 0.06</td><td>Jy</td></tr><tr><td>WISE W3<sup>g</sup></td><td>12.3</td><td>4.22 ± 0.014</td><td>mag</td></tr><tr><td>AKARI IRC<sup>i</sup></td><td>19.2</td><td>0.312 ± 0.039</td><td>Jy</td></tr><tr><td>WISE W4<sup>g</sup></td><td>22.3</td><td>3.954 ± 0.021</td><td>mag</td></tr><tr><td>IRAS<sup>i</sup></td><td>23.3</td><td>0.34 ± 0.04</td><td>Jy</td></tr><tr><td>MIPS<sup>j</sup></td><td>23.7</td><td>196.20 ± 0.08</td><td>mJy</td></tr><tr><td>IRAS<sup>i</sup></td><td>59.4</td><td>0.85 ± 0.11</td><td>Jy</td></tr><tr><td>PACS<sup>l</sup></td><td>71.2</td><td>0.961 ± 0.010</td><td>Jy</td></tr><tr><td>MIPS<sup>j</sup></td><td>71.4</td><td>1035 ± 6</td><td>mJy</td></tr><tr><td>IRAS<sup>i</sup></td><td>100</td><td>&lt;1.08</td><td>Jy</td></tr><tr><td>PACS<sup>l</sup></td><td>101</td><td>0.925 ± 0.008</td><td>Jy</td></tr><tr><td>PACS<sup>l</sup></td><td>164</td><td>0.651 ± 0.018</td><td>Jy</td></tr><tr><td>SPIRE<sup>k</sup></td><td>249</td><td>0.312 ± 0.026</td><td>Jy</td></tr><tr><td>SPIRE<sup>k</sup></td><td>350</td><td>0.1421 ± 0.0038</td><td>Jy</td></tr><tr><td>JCMT<sup>l</sup></td><td>447</td><td>&lt;1701</td><td>mJy</td></tr><tr><td>SPIRE<sup>k</sup></td><td>504</td><td>0.0674 ± 0.0035</td><td>Jy</td></tr><tr><td>JCMT<sup>l</sup></td><td>845</td><td>20.1 ± 2.7</td><td>mJy</td></tr><tr><td>ATCA<sup>m</sup></td><td>6760</td><td>0.093 ± 0.017</td><td>mJy</td></tr><tr><td>ALMA B7<sup>n</sup></td><td>856</td><td>13.2 ± 1.3</td><td>mJy</td></tr><tr><td>ALMA B7 (inner)<sup>n</sup></td><td>856</td><td>169 ± 22</td><td>µJy</td></tr><tr><td>ALMA B6<sup>n</sup></td><td>1250</td><td>5.3 ± 0.6</td><td>mJy</td></tr><tr><td>ALMA B6 (inner)<sup>n</sup></td><td>1250</td><td>59 ± 22</td><td>µJy</td></tr></tbody></table>

Notes. All wavelengths are provided to three significant figures. The calculation of the ALMA values in the lower part of this table are derived in Section 3. $^a$Mermilliod (1987, 2006); $^b$Paunzen (2015); $^c$Høg et al. (2000); $^d$Gaia Collaboration (2018); $^e$ESA (1997); $^f$Skrutskie et al. (2006); $^g$Wright et al. (2010); $^h$Ishihara et al. (2010); $^i$Neugebauer (1984); $^j$Scott et al. (2010); $^k$Sibthorpe et al. (2018); $^l$Holland et al. (2017); $^m$Ricci et al. (2015); $^n$this work.

blackbody extrapolation may not be an accurate representation of the stellar spectrum, see White et al. (2020), we can nevertheless use these models to estimate the 856 µm and 1.25 mm emission of the star, and find these to be 99 and ~46 µJy, both with 2 per cent uncertainty respectively.

We note here that since the emission from the region consistent with the stellar location is resolved separately to the emission due to the complete $q^1$ Eri system, data points can be plotted for this inner emission (see green data points in lower-right of Fig. 3). Using modified blackbodies, these data points allowed us to constrain the inner component more tightly than would have been possible if the

images of $q^1$ Eri were unresolved. Our modified blackbodies are identical to regular blackbodies until a critical wavelength, $\lambda_0$, after which the flux density is multiplied by $(\lambda/\lambda_0)^{-\beta_0}$. The temperature (and thus blackbody radius) of this inner component is relatively well constrained (i.e. by near and mid-infrared emission near this component's peak). However, the parameters defining the blackbody modification are not well constrained (i.e. due to uncertainty on the stellar photosphere and given that the two Band 6 and Band 7 data points these model have only low signal-to-noise excesses) and thus there are a broad range of values for which $\lambda_0$ and $\beta_0$ are consistent with our data. For this reason, whilst this two-component model gives a good match to the complete flux distribution of $q^1$ Eri, and the Band 6 and 7 data points in the inner regions of the system, there still remains large uncertainty on what might be expected for the faint sub-mmm emission of the inner region at these wavelengths. We discuss the implications of this further in Section 3.1.1 and Section 5.3.

We note that there are two ~850 µm flux measurements, our ALMA data (13.2 ± 1.3mJy) and JCMT data (20.1 ± 2.7mJy). These measurements are formally consistent with each other to within 2.3 $\sigma$, but we prefer the ALMA measurement since it is consistent with the Band 6 measurement for a typical debris disc millimetre spectral index of 2.5 (Ricci et al. 2015; MacGregor et al. 2016). This consistency also suggests that the ALMA observations in the two bands with differing maximum recoverable scales do not resolve out significant levels of emission.

## 3 OBSERVATIONAL ANALYSIS

In this section we present an initial analysis of the observations outlined in Section 2 to assess the disc's morphology and any asymmetries, motivating the more detailed modelling in Section 4.

### 3.1 Continuum analysis

Prior to conducting any further image analysis we raised the S/N per beam (by ~50 per cent) of the Band 7 data by applying a uv-taper of 0.5 arcsec, effectively increasing the synthesized beam to ~0.82 arcsec × 0.75 arcsec, shown in Fig. 4. This raises the S/N in comparison to the Briggs weighted image shown in the top of Fig. 1. Throughout this work the directions and disc positions referred to as the north-east (NE), south-east (SE), north-west (NW), and south-west (SW) are shown in this diagram. We also show the major (NE-SW) and minor (NW-SE) axes shown as white dotted lines (based on the position angle of the disc, calculated in the following paragraph).

We first measure the disc position angle, PA (measured anticlockwise from North) and inclination, i (for which 90° would be edge-on), required to deproject the data. We found the position angle by plotting the total flux within 10° wedges from the stellar centre as a function of angle (iterating this procedure from 0 to 180° in 1° increments), which finds a value of PA=57.0±1.0°, consistent with previous analyses which found the position angle as ~56° and 54 ± 5° respectively (see Stapelfeldt et al. 2007; Liseau et al. 2010). We define the direction of the position angle as the disc major axis, the minor axis as the line perpendicular to this, both through the stellar position.

#### 3.1.1 Continuum brightness analysis

Fig. 5 shows the radial profiles (brightness as a function of distance from the star) in the major and minor axes (upper and lower plots respectively) for the Band 7 image. These are measured by summing
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The tapered ALMA Band 7 image (natural weighted), demonstrating the new beam size of ~0.82 arcsec × 0.75 arcsec in the lower left in white, the intensity of the disc emission with contours of +5σ, +10σ, +15σ, and +20σ significance, the directions of the major axis (NE-SW) and minor axis (NW-SE) with white dotted lines, the stellar location in the image centre as a red star, and the line dividing the disc in half by integrated flux is the white dashed line (see later Section 3.1.1).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Integrated radial profiles for the Band 7 image in Fig. 4, taken along the major axis (top) and minor axis (bottom), with $I_{\text{HM}}$, the half-maximum intensity, shown on the top plot. The shaded regions show the $1\sigma$ error. Our procedure for producing these is outlined in the first paragraph of Section 3.1.1.</figcaption></figure>

<table><caption>**Table 5.** Integrated radial profile derived radii, inner and outer edges, widths, brightnesses, and radial offsets from Fig. 5. Combined errors are found from the quadrature sum of the individual measurement values for the averages. Widths in the major-axis measurements are defined as the difference in distance between the inner and outer edges. Averages found around the ring are detailed in the lower panel of the table. Projected distances in au have assumed an inclination of $i = 77^\circ$ and the Gaia distance to $q^1$ Eri of 17.34 pc.</caption><thead><tr><th>Major axis</th><th>Direction</th><th>Value</th></tr></thead><tbody><tr><td>Peak radius</td><td>NE</td><td>4.70 ± 0.04 arcsec | 81.5 ± 0.7 au</td></tr><tr><td>Peak radius</td><td>SW</td><td>4.71 ± 0.04 arcsec | 81.7 ± 0.7 au</td></tr><tr><td>Peak brightness</td><td>NE</td><td>0.225 ± 0.013 mJy arcsec<sup>-1</sup></td></tr><tr><td>Peak brightness</td><td>SW</td><td>0.209 ± 0.013 mJy arcsec<sup>-1</sup></td></tr><tr><td>Inner edge</td><td>NE</td><td>3.62 ± 0.14 arcsec | 62.8 ± 2.4 au</td></tr><tr><td>Outer edge</td><td>NE</td><td>5.71 ± 0.18 arcsec | 99.0 ± 3.1 au</td></tr><tr><td>Inner edge</td><td>SW</td><td>3.14 ± 0.14 arcsec | 54.4 ± 2.4 au</td></tr><tr><td>Outer edge</td><td>SW</td><td>6.02 ± 0.35 arcsec | 104 ± 6 au</td></tr><tr><td>Width<sub>NE</sub></td><td>NE</td><td>2.09 ± 0.22 arcsec | 36.2 ± 3.8 au</td></tr><tr><td>Width<sub>SW</sub></td><td>SW</td><td>2.88 ± 0.38 arcsec | 50 ± 7 au</td></tr><tr><th>Minor axis</th><th>Direction</th><th>Value</th></tr><tr><td>Peak radius</td><td>NW</td><td>1.08 ± 0.10 arcsec | 77.4 ± 7.2 au</td></tr><tr><td>Peak radius</td><td>SE</td><td>1.08 ± 0.10 arcsec | 77.4 ± 7.2 au</td></tr><tr><td>Peak brightness</td><td>NW</td><td>0.247 ± 0.016 mJy arcsec<sup>-1</sup></td></tr><tr><td>Peak brightness</td><td>SE</td><td>0.258 ± 0.016 mJy arcsec<sup>-1</sup></td></tr><tr><th>Averages</th><th>Axis</th><th>Value</th></tr><tr><td>Peak-peak radius</td><td>Major</td><td>4.71 ± 0.03 arcsec | 81.6 ± 0.5 au</td></tr><tr><td>Peak-peak radius</td><td>Minor</td><td>1.16 ± 0.07 arcsec | 77.4 ± 5.0 au</td></tr><tr><td>Inner edge</td><td>Major</td><td>3.38 ± 0.10 arcsec | 58.6 ± 1.7 au</td></tr><tr><td>Outer edge</td><td>Major</td><td>5.87 ± 0.20 arcsec | 101.8 ± 3.5 au</td></tr><tr><td>Width (W<sub>HM</sub>)</td><td>Major</td><td>2.49 ± 0.22 arcsec | 43 ± 4 au</td></tr><tr><td>Peak brightness</td><td>Major</td><td>0.217 ± 0.009 mJy arcsec<sup>-1</sup></td></tr><tr><td>Peak brightness</td><td>Minor</td><td>0.253 ± 0.011 mJy arcsec<sup>-1</sup></td></tr><tr><th>Radial offset</th><th>Axis</th><th>Value</th></tr><tr><td>R<sub>inner</sub></td><td>Major</td><td>0.071 ± 0.029</td></tr><tr><td>R<sub>peak</sub></td><td>Major</td><td>-0.001 ± 0.006</td></tr><tr><td>R<sub>outer</sub></td><td>Major</td><td>-0.026 ± 0.033</td></tr></tbody></table>

the flux in pixels in columns perpendicular to the plotted axis. For the major and minor profiles, these are found in columns either ~0.3 arcsec or ~0.5 arcsec (respectively) above or below the axis being plotted, to minimize signal from off-axis emission. The peak brightnesses and associated radii found from these profiles are given in Table 5, along with the average peak to peak radial distances. This shows that the disc emission brightness peaks at an average projected distance of 81.6 ± 0.5 au from the star, that the peak brightnesses in the major axis (NE and SW) are consistent, and that the peak brightnesses in the minor axis (NW and SE) are also consistent (i.e. based on these peaks, the disc geometry is consistent with a circle centred on the star). To estimate the inclination, $i$, we used the peak to peak radial averages ($r_{\text{maj,peak-peak}} = 4.71 \pm 0.03$ arcsec and $r_{\text{min,peak-peak}} = 1.08 \pm 0.07$ arcsec), since for rings with low eccentricity, $i \approx \arccos(r_{\text{min,peak-peak}}/r_{\text{maj,peak-peak}})$, giving $i = 76.7 \pm 1.0^\circ$, also consistent with previous analyses which found this inclination to be $76^\circ$ and $>63^\circ$ respectively (see Stapelfeldt et al. 2007; Liseau et al. 2010). To attain projected radii (in au) from our radial profiles in Fig. 5 (as noted in Table 5), we assume $i = 77^\circ$.

We first define some of the features of the major axis profile (top) in Fig. 5, which is more clearly resolved and has a higher signal to noise than the minor axis profile. This profile demonstrates that there is no significant difference between the peak brightness in the NE and SW directions, nor is there a difference in the radius at which this peak is measured. Calculating the average of the two NE and
SW emission peak intensities and halving this, we also measured radii consistent with this half-maximum intensity, which we use as a measure of the disc width, and find this to be $W_{HM} = 43 \pm 4$ au (see values in Table 5, under ‘Averages’). This width measurement is affected less by noisy emission towards the outer regions of the image, and we therefore use the radial location of the half-peak emission radii to define the inner and outer edge locations in both directions, and the width of the disc ansa in the NE and SW between these radii, respectively. We note however that if the average radius in both NE and SW directions at which emission exceeds 3 σ is used to define its extent, then the disc emission can be seen to extend over a broader 34–134 au. Since the image is signal-to-noise limited, this is a lower limit, and the disc extends over at least this range.

Whereas the radial profiles beyond the emission peak (i.e. on the outer side) in the NE and SW directions are consistent within their error bars for all measured radii, this is not the case on the inner side, where there is a significant radial width over which the south-west is brighter than the north-east. We quantify the extent a given radial diagnostic *r* (i.e. the projected distance from the star) differs between the NE and SW (in $r_{NE}$ and $r_{SW}$) using what we call the radial offset, $R = (r_{NE} - r_{SW})/(r_{NE} + r_{SW})$, for which a positive value would indicate that the emission measured is offset from the stellar emission (at the coordinate centre) in the NE direction. The offset of the inner edge $R_{inner} = 0.071 \pm 0.029$ is 2.4 σ significant, but those of the location of the peak emission and of the outer edge are insignificant of $R_{peak} = -0.001 \pm 0.006$ and $R_{outer} = -0.026 \pm 0.033$. Although the significance of $R_{inner}$ does not exceed 3 σ, the SW emission exceeds the NE emission for all inner edge radii $\ge 1.2$ arcsec, and therefore there are a broad range of radii on the inner edge that such an offset measurement would yield a value larger than 0 (i.e. $R_{inner}$ can be found at a similar level over a radial extent that covers several beam widths).

We considered the minor axis emission similarly to the major axis, as plotted in the lower profile of Fig. 5, and tabulated data in Table 5. We find that the peak brightnesses do not differ significantly between the NW and SE and that the peak emission radii are strongly consistent, however in the minor axis there is no evidence of a radial offset at any radii.

Most debris discs with resolved scale heights are viewed edge-on. When considering the narrower and less inclined debris disc of HD 181327, Marino et al. (2016) showed that inclined discs ($i \sim 70°$) observed with ALMA have measurably different azimuthal brightness profiles for low and moderate scale heights. Given that we have determined $q^1$ Eri’s debris disc to be more inclined than HD 181327, we therefore investigate whether or not the dimensionless vertical aspect ratio $h = H/r$ is resolved. Herein we refer to the vertical aspect ratio as the scale height, and define $H$ as the height of dust above the mid-plane at a radius $r$.

We first note that the major axis radial profile would be largely unaffected by the vertical emission distribution, which would not be the case in the minor axis of an inclined disc. We then use the major axis radial profile to determine what the minor axis radial profile should look like if the disc was flat (i.e. $h = 0$), accounting for the inclination projection and the smoothing produced by the beam size parallel with the minor axis (i.e. by re-scaling the major axis profile to compare with the minor axis). In Fig. 6, we plot this comparison, referring to the major axis profile viewed in the minor axis as the filtered major axis (assuming this to be thin), alongside the averaged minor axis radial profile. We find that the observed minor axis profile is broader than the filtered major axis profile (if assumed to be thin) which would be the case if the disc instead has a vertical scale height.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Normalized emission intensity as a function of the distance from the stellar centre for the ‘filtered’ major axis profile (average of the SW and NE profiles, projected onto the minor axis following deconvolution and convolution to account for the change in beam extent, in solid, amber), the filtered major axis smoothed with a further Gaussian convolution (in dashed, green), and the average minor axis intensity profile (in dotted, blue).</figcaption></figure>

By convolving the filtered major axis profile with a Gaussian ($\sigma_h = 0.8$ arcsec), we broadly reproduce the minor axis profile, shown in Fig. 6 as the smoothed major axis profile. We also plot above and below the smoothed major axis profile, the result of convolving the filtered major axis instead with a Gaussian, each with $\sigma_h = 0.6$ arcsec (lower) and $\sigma_h = 1.0$ arcsec (upper). Since these also broadly fit the width of the 1 σ error region of the minor axis, this allows us to constrain the level of smoothing required, which we find to be consistent with that introduced by emission with a constant scale height $h = 0.04 \pm 0.01$. Therefore, whilst this may only be weakly constrained (and this may be biased by other substructure in the disc), such an analysis tentatively suggests that the vertical scale height of $q^1$ Eri has been resolved in Band 7, and that this is at the level of a few per cent of the disc radius.

We finally note that the peak emission coincident with the star, if assumed to be a point source, has a Band 7 flux of $F_{star, B7} = 169 \pm 22$ µJy and a Band 6 flux of $F_{star, B6} = 59 \pm 15$ µJy, with the error estimated by summing in quadrature the image rms with an assumed 10 per cent flux calibration uncertainty. For the Band 7 emission, this is higher than the ~99 µJy expected from the star by ~3.2 σ, although the ~46 µJy Band 6 emission is consistent with the stellar emission (see Table 6 and Section 2.3). This Band 7 excess may therefore comprise of emission from the stellar photosphere and an extra unresolved warm emission component, from a possible inner belt, which we discuss further in Section 5.3.

In summary, by analysing the brightness distribution of the Band 7 image, we have demonstrated that this disc is inclined, broad in extent, has a significant radial offset towards the SW (in the major axis) on the disc inner edge, is symmetric in the minor axis, and may contain detectable sub-mm emission from an inner warm component near to the star. Such a major axis offset could be the result of a larger-scale asymmetric distribution in the disc, such as the presence of a clump on the SW inner edge, due to the distribution being eccentric
<table><caption>**Table 6.** Table of flux values calculated within the defined regions from Figs 7 and 8, including the flux ratios comparing the two disc halves. Note that a 10 per cent flux calibration error has been added in quadrature to the determined $F_{tot}$ values.</caption><thead><tr><th>Flux measurements</th><th>Band 7<br>(mJy)</th><th>Band 6<br>(mJy)</th></tr></thead><tbody><tr><td>$F_{\text{star}}$</td><td>0.169 ± 0.022</td><td>0.059 ± 0.015</td></tr><tr><td>$F_{\text{tot}}$</td><td>13.2 ± 1.3</td><td>5.3 ± 0.6</td></tr><tr><td>$F_{\text{SW, 82au}}$</td><td>1.78 ± 0.06</td><td>0.74 ± 0.06</td></tr><tr><td>$F_{\text{SW, 200au}}$</td><td>7.23 ± 0.09</td><td>2.59 ± 0.09</td></tr><tr><td>$F_{\text{NE, 82au}}$</td><td>1.48 ± 0.06</td><td>0.67 ± 0.06</td></tr><tr><td>$F_{\text{NE, 200au}}$</td><td>5.90 ± 0.09</td><td>2.59 ± 0.09</td></tr><tr><td>$\delta F$</td><td>1.33 ± 0.13</td><td>0.00 ± 0.13</td></tr><tr><td>Flux ratios</td><td>Band 7</td><td>Band 6</td></tr><tr><td>$f_{82au}$</td><td>1.21 ± 0.06</td><td>1.10 ± 0.13</td></tr><tr><td>$f_{200au}$</td><td>1.23 ± 0.04</td><td>1.00 ± 0.07</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Band 7 flux profiles for the inner 200 au of $q^1$ Eri ALMA data, found by integrating the total flux within deprojected radial bins cumulatively. Top: Flux profile for the full image, with shaded regions showing the ±1σ errors. Bottom: Difference profile for the full image, for flux in the NE half minus flux in the SW half. The shaded region is the ±1σ quadrature sum of the error from the two halves, and the red-dotted line shows the flux difference at 200 au.</figcaption></figure>

(with a pericentre direction between the NW and SW), a combination of these, or simply due to noise. We explore these hypotheses further in Sections 3.1.2, 4.2, and 5.2.

### 3.1.2 Flux profile analysis

We produced flux profiles for the image by halving the Band 7 image (Fig. 4) along the minor axis, and summing the total flux within projected radial bins either side of the minor axis. The results are shown for the Band 6 and 7 images in Figs 7 and 8. Table 6 shows the values measured from these plots of the total flux, $F_{tot}$, within a projected distance of 200 au from the star. These are used to produce the flux distribution in Fig. 3, and can be used to calculate the spectral index $\alpha_{mm} = 2.34 \pm 0.29$. Given that the previous brightness profile analysis concluded that an asymmetry may exist

<figure><img src="image_10.png" /><figcaption>**Figure 8.** As per Fig. 7 but for the Band 6 data.</figcaption></figure>

on the inner edge of the disc, we sought to quantify this further with measurements of the integrated flux. By ignoring emission internal to 30 au (i.e. from the stellar photosphere and possible inner warm component), we measured the total flux on either side of the disc internal to the brightness maxima radii (i.e. internal to ~82 au), noted as $F_{SW, 82au}$ and $F_{NE, 82au}$ respectively, and internal to 200 au, noted as $F_{SW, 200au}$ and $F_{NE, 200au}$ respectively, and also the net difference between these at 200 au, $\delta F = F_{SW, 200au} - F_{NE, 200au}$. To quantify the flux enhancement on either side of the disc, we also calculated the flux ratios between the NE and SW from these same measurements, as $f_{82au} = F_{SW, 82au}/F_{NE, 82au}$ and $f_{200au} = F_{SW, 200au}/F_{NE, 200au}$ respectively. We note here that whilst the absolute values of fluxes should all include flux calibration uncertainties of 10 per cent, this uncertainty does not feature in the calculations for flux ratios, since this calibration uncertainty affects all flux in this map equally (i.e. to assess asymmetries this does not need to be included, but for individual flux measurements this does).

Considering first the Band 7 data, it can be seen in Fig. 7 that although the shapes of the integrated flux in the NE and SW are similar, the values diverge. This is more easily seen in the lower plot which quantifies the difference between these. This shows that from ~50 au (i.e. on the inner edge) these are significantly different, reaching a maximum difference at ~200 au. Although the peak flux difference value is found at this longest radius, this is consistent with no flux difference beyond ~160 au, i.e. on the outer extent of the disc (see Section 3.1.1). We quantify this difference in Table 6 as $\delta F = 1.33 \pm 0.18$ mJy, and the flux ratios along the inner edge and over the full disc extent as $f_{82au} = 1.21 \pm 0.06$ and $f_{200au} = 1.23 \pm 0.04$ respectively, all of which show strong evidence that the Band 7 disc flux is asymmetric. The Band 6 flux profile in Fig. 8 shows no significant flux asymmetry. Although the lower plot demonstrates that the SW may contain more flux, this is not significant, as quantified by the Band 6 $\delta F$ value (consistent with 0) or either the $f_{82au}$ and $f_{200au}$ values (both consistent with 1). Despite these Band 6 flux ratios and flux difference not showing evidence for
<figure><img src="image_11.png" /></figure>

**Figure 9.** A minor axis mirror-subtracted Band 7 image of $q^1$ Eri where the flux in each pixel has had subtracted that of the corresponding pixel on the opposite side of the minor axis at 0 arcsec offset (vertical in image). This also shows the continuum emission contour lines for $+6\sigma$, $+11\sigma$, and $+16\sigma$, and the significance of the emission after subtraction in the blue and red contours, respectively $\pm3$ and $\pm4\sigma'$ (where $\sigma'=\sqrt{2\sigma_{B7}}$ image). Note that the left-hand and right-hand sides of this image represent the north-eastern and south-western halves of the disc, respectively.

an asymmetry, the lower S/N of this image means that an asymmetry could still be present at broadly the same ~20 per cent level seen in Band 7. Deeper observations of $q^1$ Eri in Band 6 would be required to determine this.

### 3.1.3 Interpreting the flux and brightness profiles

If the disc was overall symmetrical, but with a clump in the SW ansa, then this would lead to a brightness asymmetry between the two ansae. To explore this, Fig. 9 shows a mirror-subtracted image in which the flux in each pixel has had subtracted that of the corresponding pixel on the opposite side of the minor axis (i.e. relative to 0 arcsec offset in the major axis). This plot shows the thin black contour lines of the original Band 7 image, and coloured contours where the mirror subtraction has resulted in significant residual emission (i.e. $>3\sigma$). This also demonstrates that the asymmetry in the SW ansa is interior to the peak brightness radius and may extend over a broad azimuthal range of the disc. We note that the two $4\sigma$ residual contours have sizes of $\approx0.7$ arcsec respectively at projected orbital radii ~60 au. No further asymmetries are evident radially beyond the emission peak, which are coincident with the disc, consistent with our analysis of the brightness. Thus, one scenario that we will investigate in Section 4 with model CL is that the sub-mm emission arises from a broad, symmetric ring of parent planetesimals with an extended clump of emission on the inner edge of the SW ansa (equivalently, this could be interpreted as the SW ansa being radially broader than the NE).

Another interpretation of the asymmetries that we will explore with model ECC is that these are not caused by a clump, but due to the disc being eccentric. If this was the case, the disc flux would be symmetrical about its pericentre direction. To explore what that pericentre direction might be, we found the angle at which the total integrated flux ratio (either side of an axis rotated around the stellar position) was equal to 1 (i.e. we found the total flux line of symmetry, accounting for the errors consistently as per our analysis of the fluxes in the two disc halves). This gave a value of $84.5 \pm 2.6^\circ$ clockwise from the minor axis (i.e. pointing between the SW and NW), shown in Fig. 4 as the thick white dashed line. We also consider a combination of models CL and ECC, i.e. an eccentric disc which also has a clump of emission on the inner edge of the SW ansa (model CL + ECC).

## 3.2 Scattered light analysis

There are two main observational conclusions that we draw from Fig. 2. First, along the major axis there is a radial offset (brightness asymmetry) in that the scattered light emission extends out to larger radii in the NE than the SW. To aid viewing this offset, the data has been annotated with dark-red dashed rings with radii of 5.7 and 6.7 arcsec (i.e. an average projected radii of ~108 au), which overlap with the $+22$ mag arcsec$^{-2}$ surface brightness contour lines in the SW and NE ansae, and a black dash-dotted line along the major axis. Similarly to the radial offset analysis in Section 3.1.1, we use these distances with an estimated error of $\pm0.1$ arcsec to find a scattered light radial offset of $R_{\text{scat}} = (r_{\text{NE}} - r_{\text{SW}})/(r_{\text{NE}} + r_{\text{SW}}) = 0.081 \pm 0.011$. This demonstrates that there is a significant offset towards the NE direction in the scattered light emission in the outer region of the disc, reminiscent of those reported for the sub-mm emission in Section 3.1.1 in the inner regions of the disc. This scattered light asymmetry is explored further in Section 5.

Secondly, the scattered light emission along the minor axis in the SE region of the disc is brighter and more azimuthally extended in comparison to the NW. This can be seen by visually comparing the emission either side of the black dash-dotted line (showing a position angle of $56^\circ$), which demonstrates there is more emission on the southern side than the northern. This type of asymmetry is usually inferred to be due to the dust emission preferentially forward scattering in this direction (Augereau et al. 1999). Therefore this may suggest that the SE of the disc is on our near-side. However, since emission within ~2.25 arcsec cannot be probed due to coronagraph obscuration, this limits our ability to estimate the significance of this minor axis asymmetry. Higher resolution scattered light imaging that can observe closer to the star than 2.25 arcsec is thus necessary to further constrain this.

## 3.3 Summary of observational analysis

Taken in conjunction, we have found four significant observational constraints that we can place on the $q^1$ Eri debris disc, the first 3 from ALMA (O1–3) and the fourth from HST (O4), which will be used to constrain our modelling in Section 4 and discussed further in Section 5: **O1 - Broad Structure:** $q^1$ Eri has a bright, broad disc, with a position angle $PA=57.0 \pm 1.0^\circ$, and an inclination $i = 76.7 \pm 1.0^\circ$. This emission peaks at a radius of $81.6 \pm 0.5$ au, and has a half-maximum width of $W_{\text{HM}} = 43 \pm 4$ au, although likely extends from at least 34–134 au. The disc vertical scale height may have also been resolved with an aspect ratio $h = 0.04 \pm 0.01$. **O2 - Major Axis Brightness Asymmetry:** the major axis emission brightness profile is asymmetric. At the same projected radii, the SW inner edge is brighter than the NE, but the outer edge is comparable on both sides. The inner edge asymmetry was quantified with the radial offset parameter $R_{\text{inner}} = 0.071 \pm 0.029$. The peak emission radius and outer edge were assessed similarly, finding values of $R_{\text{peak}} = -0.001 \pm 0.006$ and $R_{\text{outer}} = -0.026 \pm 0.033$ respectively. **O3 - Major Axis Flux Asymmetry:** the total Band 7 fluxes between 30–82 and 30–200 au are greater on the SW side of the disc than on the NE side by fractional amounts $f_{82au} = 1.21 \pm 0.06$ and $f_{200au} = 1.23 \pm 0.04$ respectively. The Band 6 images are broadly consistent with the same ~20 per cent level of asymmetry as the Band 7 images, although the lower sensitivity at the longer wavelength precludes an independent determination. **O4 - Scattered Light:** the scattered light data shows the disc brightness to be radially and azimuthally asymmetric. We showed that dust emission has a significant radial offset towards the NE of the star for the $22$ mag arcsec$^{-2}$ emission
contour lines, $R_{scat} = 0.081 \pm 0.011$, at an average projected radius of ~108 au. We also found the SE side of the disc to be brighter than the NW, suggestive that the SE edge is on the disc's near-side to us.

Three possible scenarios were proposed in Section 3.1.3 to model this system: model CL, a symmetric disc with an emission clump on the inner edge of the SW ansa; model ECC, an eccentric disc; and model CL + ECC, a combination of an eccentric disc with a clump on its SW inner edge. For comparison we will also present a symmetrical disc model SYM.

## 3.4 CO J = 3–2 spectral line analysis

Separate to the continuum analysis, for the 2018 Band 7 ALMA data, we produced a subset of the CASA measurement set with all spectral windows removed except the one containing the CO J = 3–2 spectral line ($f = 345.79599$ GHz, herein referred to as $f_{CO}$). Any circumstellar CO emission would be present across multiple data channels, as determined by the radial velocity of the star, $v_{rad}=27.82 \pm 0.15$ kms$^{-1}$ (Gaia Collaboration 2018), and the range of Keplerian velocities expected in the disc. To remove continuum emission from any CO signal, we used the CASA tool uvcontsub with a fit order of 1 external to the region defined by $f_{CO} \pm \Delta f$, where $\Delta f$ was set to a velocity width $\pm 35$ kms$^{-1}$, avoiding fitting to channels where CO could be present. We produced a data cube in the barycentric reference frame excluding all channels further than 200 kms$^{-1}$ away from the J = 3–2 line in velocity space using the CASA tclean algorithm (with zero clean iterations, since there was no significant emission per beam present in any channel).

To measure the line profiles, we applied the spectro-spatial filtering method of Matrà et al. (2015, 2017). To do so, we applied a Keplerian mask to the data cube to shift channels in both of the possible disc rotation directions (clockwise and anticlockwise) based on a stellar mass of $1.1 M_{\odot}$ (Marmier et al. 2013). We integrated all emission in rotation directions between projected radii of ~38 and 125 au from the star (i.e. within the broad region defined by the full-width half-maxima of the sub-mm emission), which are both shown as a function of radial velocity in Fig. 10. We found the rms of the two line profiles (clockwise and anticlockwise rotations) across all velocities, and measured rms=7.0 mJy in the clockwise-rotated direction, and rms=7.1 mJy in the anticlockwise rotated direction. We define our upper bound based on the larger of these two rms values, and place a $3\sigma$ limit on the CO line peak flux of $F_{CO,3\sigma}=21.3$ mJy. The velocity channel spacing of the data is $\Delta v=0.423$ kms$^{-1}$, however the effective spectral bandwidth is ~2.667 larger than this, since adjacent ALMA channels are not fully independent from each other.<fnref n="2" /> The $3\sigma$ limit on this flux is then found as $\delta F_{CO,3\sigma}=2.667 \times \Delta v \times$ rms=24.0 mJy kms$^{-1}$. Based on this, we derive a gas mass upper bound in Section 5.4.

## 4 MODELLING

Following the conclusion of the observational constraints in Section 3.3, in this section we test these against parametric models of the Band 7 data (for simplicity), first outlining our modelling methodology, then giving best-fitting model parameters and discussing their implications.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Flux as a function of radial velocity within spectrally and spatially filtered maps. The spectral shifts assume that the projected orbital motion is in the clockwise and anticlockwise directions in the top and bottom plots, respectively. The star's radial velocity is shown in grey.</figcaption></figure>

## 4.1 Model types and definitions

The models described in Section 3.1.3 are based on the same 12-parameter disc model. The 12 common parameters are the dust mass ($M_{dust}$), characteristic radius ($r_c$), scale height ($h$), inclination ($i$), position angle (PA), inner and outer power-law indices ($p_{in}$ and $p_{out}$, described further below), photospheric Band 7 flux ($F_{star, B7}$), and the phase centre offsets in RA and Dec for both Band 7 data sets ($x_{off, 16}$, $y_{off, 16}$, $x_{off, 18}$, and $y_{off, 18}$). This model can be configured to introduce (either or both) an extended 2D Gaussian source, and eccentricity. Where we include an additional extended source, this is defined by its flux ($F_S$), offset from the stellar position ($x_{offset, s}$ and $y_{offset, s}$), a position angle (PA$_S$), and major and minor axis standard deviations ($R_{maj, s}$ and $R_{min, s}$), related to the FWHM via $FWHM=2\sqrt{2\ln 2}\sigma$.

We parametrize the surface density ($\Sigma(a)$) of our model based on a two-component radial power law (Augereau et al. 1999; Kennedy et al. 2018), where on a grid of $r$ (distance from the star) and $\phi$ (azimuthal angle from the pericentre direction),

$$ \Sigma(a) \propto \left[ \left( \frac{a}{r_c} \right)^{-2p_{in}} + \left( \frac{a}{r_c} \right)^{-2p_{out}} \right]^{-1/2}, \quad (1) $$

for $a=r[1 + e\cos(\phi - \omega)]/(1 - e^2)$, where $\omega$ is the angle between the pericentre direction and the line of nodes where the disc plane crosses the sky plane (i.e. the major axis, for which $\omega = 0^{\circ}$ is the south-western direction). This parametrization of the surface density is valid for a disc containing particles with the same eccentricity and pericentre (Marino et al. 2019); note, however, that it is possible that particle eccentricities and orientations are not all identical in such a broad disc, and we discuss the implications of this in Section 4.3. The models have the same vertical Gaussian density distribution, defined by the scale height, $h = H/r$, where $H$ is the height of dust above the mid-plane at a radius $r$ (e.g. Marino et al. 2016, 2019). Where used, eccentricity $e$ is a constant. All models use the same minimum and maximum grain sizes of $a_{min}=0.9$ µm and $a_{max}=1$cm (respectively set by the grain blow out size, and by neglecting emission from larger grains at Band 7 wavelengths), a dust grain density of $2.7$ g cm$^{-3}$, a grain-size distribution with

<footnote n="2">²For a complete discussion of this, see https://safe.nrao.edu/wiki/pub/Main/ALMAWindowFunctions/Note_on_Spectral Response.pdf</footnote>
power-law exponent $\alpha = 3.5$ (Dohnanyi 1969), and a weighted mean dust opacity based on a mix of compositions with mass fractions of 70 per cent astrosilicates, 15 per cent amorphous carbon, and 15 per cent water ice (for example, identical to that used to model HD107146 by Marino et al. 2018, i.e. $1.88 \text{ cm}^2 \text{ g}^{-1}$ at $\lambda=856 \text{ µm}$). For model CL (described in detail in Section 4.2 we show the dust density and temperature profiles in Appendix B, which shows a temperature of $T = 44 \text{ K}$ at the location where the dust density peaks.

The Band 7 emission from the star was estimated to be ~99 µJy in Section 2.3. However, since stellar emission at sub-mm wavelengths is often poorly constrained, the unresolved flux at the stellar location in the Band 7 image is left as a free parameter $F_{\text{star, B7}}$.

## 4.2 Model fitting, results, and comparison

We used the RADMC-3D (Dullemond et al. 2012) package to compute the dust temperature of the disc and produce model images defined by the parameters outlined in Section 4.1, at a Band 7 wavelength of 856 µm, and compute model visibilities at the same uv-baselines as our ALMA observations using the tools developed in Marino et al. (2018). All models are produced with the same number of visibility data points as the combined 2016 and 2018 Band 7 observations, $N_{\text{dat}} = 2 \times N_{\text{vis}} = 19440912$ (9720456 for both real and imaginary components). To investigate the consistency of these models with our Band 7 data, we used the MCMC best-fitting model-parameter estimator, emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013). The posterior distribution generated by emcee is the product of the prior distribution and the likelihood function, which takes a value $\exp(-\chi^2/2)$, where for interferometric data sets

$$ \chi^2_{\text{Vis}} = \sum_{i=1}^{N_{\text{vis}}} \frac{\|V_{\text{data,i}} - V_{\text{model,i}}\|^2}{f_{\sigma}^2 \delta V_{\text{data,i}}^2}, \quad (2) $$

in which $V$ is the value of the each visibility, $\delta V$ is the intrinsic dispersion of the visibilities (as calculated by the CASA package statwt), and $f_\sigma$ a normalization parameter discussed further below.

For large interferometric data sets, individual visibility measurements can be dominated by noise (i.e. $S/N < 1$), leading to the reduced $\chi^2_{\text{Vis}}$ values ($\chi^2_{\text{red}}=\chi^2_{\text{vis}}/N_{\text{dat}}$) of even null models being close to 1 (i.e. the value expected for a perfect model). From previous ALMA modelling of debris discs, it has been shown that the CASA routine statwt (used to re-weight visibilities according to their dispersion) does not provide an accurate measure of the absolute uncertainties on visibility measurements, although it does provide a good estimate of the relative errors between visibilities. This means that the weights or variance can be off by a small factor, and subsequently affect the posterior distributions of modelled parameters (i.e. leading to incorrect parameter estimations and errors). To account for this, as per previous analyses of ALMA data (e.g. Marino et al. 2018), we apply a correction factor, $f_\sigma$ (see equation 2), such that $\chi^2_{\text{red null}}=1$ (i.e. for a model with $V_{\text{model}}=0$ everywhere). For both 2016 and 2018 data sets we find a similar $f_\sigma \sim 1.6$ (i.e. this confirms that statwt has underestimated the uncertainty in the visibilities). The need to scale uncertainties by $f_\sigma$ however means that the absolute value of $\chi^2$ cannot be used to assess the goodness of fit of an individual model (e.g. since the null model should have $\chi^2_{\text{red}} > 1$ as it should be a poor fit). However, given that $f_\sigma$ is constant, the $\chi^2_{\text{vis}}$ values can still be used to quantify the improvement between different models. We outline our assessment of the goodness of fit in the following section.

### 4.2.1 Model results

We investigated the consistency of our data with a symmetrical model SYM, defined by the 12 common parameters introduced in Section 4.1, before continuing with models ECC, CL, and CL + ECC, with 14, 18, and 20 free parameters, respectively. For these MCMC routines, between 140–200 walkers and 600–2000 steps were used to ensure parameter solutions were converged on, with burn in lengths varying between 250 and 1200. The results of our modelling procedure are shown in Table 7, where we report all best-fitting model values, along with $r_0$ (the average radius of the peak brightness found from model radial profiles) and the total Band 7 flux, $F_{\text{tot, B7}}$ (which includes $F_{\text{star, B7}}$), generated by each best-fitting model, for which we assume and quote a 10 per cent flux calibration error.

We first discuss the ability of all four models to explain the axisymmetric parameters, before their ability to interpret the asymmetric parameters. The average peak emission radius $r_0$ for all models is found to be within $<1\sigma$ from the value determined in the observational analysis. Likewise comparison of best-fitting parameters with those derived directly from the cleaned images shows the total Band 7 flux is within $1.3\sigma$, the inclination within $1.9\sigma$, the position angle within $0.4\sigma$, the stellar flux within $1.5\sigma$, and all stellar offsets consistent with 0 (no offset), within $3\sigma$ (for models SYM, CL and CL + ECC). In the case of model ECC, the offset of the 2018 data is not found at the origin (albeit nearby). These four models do however show broad agreement with these aspects. We did not derive estimates for the dust mass, or outer radial power law indices in the observational analysis, however these parameters are relatively well constrained by the four models. The inner edge power law index is however only poorly constrained, and the values presented in Table 7 show the $3\sigma$ lower limit, and the MCMC-derived best-fit value and error. Although only weakly constrained, in all models the inner edge index is shown to be sharp, and the better constrained outer edge index is shown to be shallow, which both appear consistent with the radial profiles in Fig. 5. That all models find a scale height significantly above 0 is discussed in Section 5, however this suggests that these models have at least partially resolved the vertical distribution of the disc’s emission. In Section 3.1.1 we showed that this was consistent with $h \sim 0.04$, and thus these values are all broadly in agreement with this (i.e. are at the few per cent level).

Whilst the broad disc parameters can be well modelled by each of the 4 models, we next assess which of these provides the best quantitative and qualitative fit. In Fig. 11, we show images of all four models (top) and the 2016 and 2018 epoch concatenated residual maps, re-imaged using tclean with a uvtaper of 0.5 arcsec (to match the image of Fig. 4) for each of the four models (bottom). Qualitatively we can see that whilst these models are similar, there are key differences in how these reproduce the overall disc morphology, and the relative emission brightness in different regions of the disc. In the residual maps, models SYM and ECC have multiple regions coincident with the disc where the data is fitted poorly (residuals $>3\sigma$). In model CL + ECC there is one region on the outer edge of the NE ansa where the residual map exceeds $3\sigma$ (and a few $2\sigma$ contours), whereas the residual map of model CL appears to fit the data most consistently (all residuals $<3\sigma$, with only a few $2\sigma$ contours). We note that despite model CL + ECC having 2 more free parameters than model CL, yet model CL appearing to be a qualitatively better fit, this difference is consistent with our quantitative assessment of the $\chi^2$ values (discussed later in this section).
<table><caption>**Table 7.** Model best-fitting posterior values from models SYM, CL, ECC, and CL + ECC. A dash indicates that a parameter was not a free parameter in the model this is associated with. In the case of *p*<sub>inner</sub> the 3 σ lower limit is shown along with the posterior value. The total Band 7 flux and peak emission radius, *F*<sub>tot, B7</sub> and *r*<sub>0</sub>, for each model is shown in the lowest section of this table (for the best-fitting model parameters), assuming a 10 per cent flux calibration error in *F*<sub>tot, B7</sub> and *F*<sub>star, B7</sub> (appropriate for interferometric observations with ALMA).</caption><thead><tr><th></th><th></th><th>SYM</th><th>CL</th><th>ECC</th><th>CL + ECC</th></tr></thead><tbody><tr><td>n<sub>walkers</sub></td><td></td><td>120</td><td>180</td><td>140</td><td>180</td></tr><tr><td>n<sub>steps</sub></td><td></td><td>600</td><td>1600</td><td>1000</td><td>2000</td></tr><tr><td>n<sub>burnin</sub></td><td></td><td>250</td><td>1000</td><td>400</td><td>1200</td></tr><tr><td>M<sub>dust</sub> (M<sub>☉</sub>)</td><td>(au)</td><td>0.0315 ± 0.0008</td><td>0.0275 ± 0.0010</td><td>0.0319 ± 0.0008</td><td>0.0282 ± 0.0012</td></tr><tr><td>r<sub>c</sub></td><td>(au)</td><td>75.5 ± 1.1</td><td>75.9 ± 1.1</td><td>75.1 ± 1.1</td><td>75.8 ± 1.1</td></tr><tr><td>e</td><td></td><td>–</td><td>–</td><td>0.053 ± 0.009</td><td>0.025 ± 0.012</td></tr><tr><td>ω</td><td>(°)</td><td>–</td><td>–</td><td>11 ± 15</td><td>19 ± 30</td></tr><tr><td>h</td><td></td><td>0.054 ± 0.004</td><td>0.046 ± 0.004</td><td>0.055 ± 0.004</td><td>0.048 ± 0.004</td></tr><tr><td>i</td><td></td><td>78.62 ± 0.20</td><td>78.68 ± 0.22</td><td>78.58 ± 0.20</td><td>78.63 ± 0.22</td></tr><tr><td>PA</td><td></td><td>56.67 ± 0.20</td><td>56.66 ± 0.21</td><td>56.65 ± 0.21</td><td>56.64 ± 0.22</td></tr><tr><td>p<sub>inner</sub></td><td></td><td>>5.1 (16.8 ± 3.9)</td><td>>5.6 (18.2 ± 4.2)</td><td>>5.5 (18.1 ± 4.2)</td><td>>5.9 (17.6 ± 3.9)</td></tr><tr><td>p<sub>outer</sub></td><td></td><td>-1.87 ± 0.13</td><td>-2.09 ± 0.18</td><td>-1.80 ± 0.13</td><td>-2.06 ± 0.17</td></tr><tr><td>F<sub>star, B7</sub> (mJy)</td><td>(arcsec)</td><td>0.145 ± 0.022</td><td>0.142 ± 0.022</td><td>0.141 ± 0.022</td><td>0.144 ± 0.022</td></tr><tr><td>x<sub>off, 16</sub> (arcsec)</td><td></td><td>0.057 ± 0.021</td><td>0.059 ± 0.020</td><td>0.044 ± 0.022</td><td>0.051 ± 0.022</td></tr><tr><td>y<sub>off, 16</sub> (arcsec)</td><td></td><td>-0.022 ± 0.018</td><td>-0.022 ± 0.018</td><td>-0.024 ± 0.020</td><td>-0.021 ± 0.019</td></tr><tr><td>x<sub>off, 18</sub> (arcsec)</td><td></td><td>0.004 ± 0.028</td><td>-0.024 ± 0.028</td><td>-0.11 ± 0.03</td><td>-0.03 ± 0.04</td></tr><tr><td>y<sub>off, 18</sub> (arcsec)</td><td></td><td>-0.036 ± 0.020</td><td>-0.030 ± 0.020</td><td>-0.109 ± 0.024</td><td>-0.058 ± 0.027</td></tr><tr><td>F<sub>S</sub> (mJy)</td><td></td><td>–</td><td>1.7 ± 0.6</td><td>–</td><td>1.7 ± 0.6</td></tr><tr><td>x<sub>offset, S</sub> (arcsec)</td><td></td><td>–</td><td>-4.2 ± 0.5</td><td>–</td><td>-4.3 ± 0.5</td></tr><tr><td>y<sub>offset, S</sub> (arcsec)</td><td></td><td>–</td><td>-2.6 ± 0.3</td><td>–</td><td>-2.7 ± 0.4</td></tr><tr><td>R<sub>maj, S</sub> (arcsec)</td><td></td><td>–</td><td>3.2 ± 0.5</td><td>–</td><td>3.5 ± 0.6</td></tr><tr><td>R<sub>min, S</sub> (arcsec)</td><td></td><td>–</td><td>2.3 ± 0.5</td><td>–</td><td>2.7 ± 0.6</td></tr><tr><td>PA<sub>S</sub> (°)</td><td></td><td>–</td><td>40 ± 16</td><td>–</td><td>50 ± 50</td></tr><tr><td>r<sub>0</sub> (au)</td><td>(au)</td><td>82.9 ± 1.2</td><td>83.0 ± 1.2</td><td>82.9 ± 1.2</td><td>82.9 ± 1.2</td></tr><tr><td>F<sub>tot, B7</sub> (mJy)</td><td>(mJy)</td><td>11.3 ± 1.1</td><td>12.3 ± 1.2</td><td>11.9 ± 1.2</td><td>12.5 ± 1.3</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Plots to show the best-fitting parameter models (top) and the residual maps produced when the models are subtracted from the concatenated Band 7 data sets (bottom) for model SYM (left), model CL (centre-left), model ECC (centre-right), and model CL + ECC (right). Model contours were chosen to demonstrate their inner and outer edges and how their emission varies in the two ansae. On the residual maps the black dashed ovals demonstrate the projected disc emission at the half-maximum inner and outer edges, and the contours demonstrate the ±2σ and ±3σ residuals (red positive, blue negative). On the lower left is the beam size for each of these images.</figcaption></figure>
<table><caption>**Table 8.** Model comparison table, demonstrating in the upper section of the table the $\chi^2$ values from the visibilities and images, and the difference in the visibility $\chi^2$ values (from model SYM), and then results for observational conclusions O1 (half-peak disc width, noting that all other O1 conclusions are readily generated by the models and are included in Table 7), O2 Maj (major axis radial offsets for the major axis), and O3 Maj (the flux ratio for emission compared in the major axis between 30 and 82 au and between 30 and 200 au) in the following table sections. Note that all values bold are > 3$\sigma$ from the equivalent metric determined in the Observational Analysis (or in the case of the $\chi^2_{\text{Im}}$ metric, a poor representation of the data), and those values in blue are between 2$\sigma$ and 3$\sigma$ from the same equivalent observational metric. For ease of comparison, the observational analysis parameters are shown in the rightmost column.</caption><thead><tr><th></th><th>SYM</th><th>CL</th><th>ECC</th><th>CL + ECC</th><th>Obs analysis</th></tr></thead><tbody><tr><td>$\chi^2_{\text{vis,Im}}$</td><td><strong>2.77</strong></td><td>1.13</td><td><strong>1.72</strong></td><td>1.13</td><td>-</td></tr><tr><td>$\chi^2_{\text{vis}_2}$</td><td>19 525 531.6</td><td>19 525 462.8</td><td>19 525 499.4</td><td>19 525 457.0</td><td>-</td></tr><tr><td>$\Delta\chi^2_{\text{vis}}$</td><td>0.0</td><td>-68.8</td><td>-30.2</td><td>-74.6</td><td>-</td></tr><tr><td>$\Delta\text{BIC}_{\text{vis}}$</td><td>0.0</td><td>31.9</td><td>3.4</td><td>59.7</td><td>-</td></tr><tr><td>$\Delta\text{BIC}_{\text{Im}}$</td><td>0.0</td><td>-191.8</td><td>-103.2</td><td>-181.2</td><td>-</td></tr><tr><td>$N_{\text{par}}$</td><td>12</td><td>18</td><td>14</td><td>20</td><td>-</td></tr><tr><td><strong>O1 (Width)</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>$W_{\text{HM}}$ (au)</td><td>54</td><td>55</td><td>50</td><td>54</td><td>43 ± 4</td></tr><tr><td><strong>O2 Maj</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>$R_{\text{peak}}$</td><td>0.00</td><td>0.00</td><td><strong>0.049</strong></td><td>0.00</td><td>-0.001 ± 0.006</td></tr><tr><td>$R_{\text{inner}}$</td><td>0.00</td><td>0.048</td><td>0.068</td><td>0.040</td><td>0.071 ± 0.029</td></tr><tr><td>$R_{\text{outer}}$</td><td>0.00</td><td>0.00</td><td>0.013</td><td>0.00</td><td>-0.026 ± 0.033</td></tr><tr><td><strong>O3 Maj</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>$f_{82\text{au}}$</td><td><strong>1.00</strong></td><td>1.30</td><td>1.39</td><td>1.27</td><td>1.21 ± 0.06</td></tr><tr><td>$f_{200\text{au}}$</td><td><strong>1.00</strong></td><td>1.17</td><td><strong>0.98</strong></td><td>1.14</td><td>1.23 ± 0.04</td></tr></tbody></table>

To quantify this, we report in Table 8 the parameter $\chi^2_{\text{vis,Im}}$ as the reduced $\chi^2$ of the residual maps within projected radii 38–125 au (based on the peak emission radius ± the half-peak emission width, see Table 5) such that a significant extent of the inner and outer edge emission was included. This required us to define the residual map error, which we determined from the rms of the residual maps external to the region that $\chi^2_{\text{Im}}$ is calculated within as 15.9 µJy beam$^{-1}$. Consistent with our qualitative assessment, the values shown in Table 8 demonstrate that models CL and CL + ECC are most consistent with a value of 1.0 (1.13 for both), and models ECC and SYM are both higher (1.72 and 2.77 respectively).

To assess how well these compare with respect to the mean and uncertainty in these residual maps, we start by noting that a $\chi^2_{\text{Im}}$ value of 1.0 would imply that the residual maps are consistent with noise. By producing $N = 3000$ maps of noise convolved with a beam equal to that in our residual maps (with a resulting rms equal to 15.9 µJy beam$^{-1}$), we found $\chi^2_{\text{Im}}$ exceeded 1.13 in ~4 per cent of all iterations, and in no circumstances were values found >1.72, which based on the number of runs, allowed us to place an upper bound on the frequency at which we would expect residual maps to exceed this as <0.03 per cent. We found these results were normally distributed (with a slight tail towards higher values), centred on 0.96, with a width of 0.10. We therefore assess the CL and CL + ECC residual maps as being consistent with a map comprised entirely of noise (i.e. given their $\chi^2_{\text{Im}}$ values are both 1.13), but that both the ECC and SYM residual maps are inconsistent. This implies from the $\chi^2_{\text{Im}}$ values that models CL and CL + ECC are both reasonable representations of our data, and models SYM and ECC are poor representations.

We next compare the visibilities from their entire interferometric data sets respectively (i.e. rather than between projected radii), and report the $\chi^2_{\text{vis}}$ values for each of the models in Table 8 with their respective differences, $\Delta\chi^2_{\text{vis}}$, from the value achieved with model SYM. Subtracted as such, this would mean that the most negative $\Delta\chi^2_{\text{vis}}$ would represent the best-fitting model, i.e. in order of best-to-worst fit, this would conclude that the models are ordered

CL+ECC, CL, ECC, and finally SYM. Since the difference between model CL and model CL + ECC is $\Delta\chi^2_{\text{vis}} \sim 6.8$, these two models provide a comparable fit to the data, whereas models ECC and SYM are correspondingly worse fits. Thus, we anticipate that the two additional free parameters required for model CL + ECC are not justified, and so that model CL is favoured by this analysis.

We further this assessment, by calculating the Bayesian Information Criterion (BIC), which statistically assesses the significance of additional free parameters, where

$$ \text{BIC} = \chi^2 + N_{\text{par}} \ln N_{\text{dat}}, \quad (3) $$

for a model with $N_{\text{par}}$ independent free parameters, and $N_{\text{dat}}$ independent data points. This is reported in Table 8 for each of the models relative to the baseline of model SYM, for consideration of both the visibilities and the image, i.e. the $\Delta\text{BIC}_{\text{vis}}$ and $\Delta\text{BIC}_{\text{Im}}$ values respectively. Where the $\Delta\text{BIC}$ between two models exceeds 10, this is often reported as statistically strong evidence to support the hypothesis that the model reporting the *lower* BIC is a more significant fit to the data than the model it is being compared with (Schwarz 1978; Kass & Raftery 1995).

This analysis however arrives at a different conclusion, with the model SYM being statistically favourable, with the model ECC being broadly comparable, but CL and CL+ECC being consecutively worse. This is surprising given that the residual image of model SYM shows that this is a poor fit to the observations, and that models CL and CL + ECC are far better. Thus whilst we may use the BIC parameter to distinguish between these models, it may not provide a robust metric to do so. To understand why this may have occurred, consider that the BIC parameter scales with $\ln(N_{\text{vis}})$, and so, for the number of visibility measurements in our data set, each additional model parameter would require a relative improvement in $\chi^2_{\text{vis}}$ of >16.8 to be justified. However, our data also contain a high fraction of data points that are negligibly affected by the model or the asymmetry, resulting in the BIC applying a disproportionate penalty to additional model parameters. In the image space there are
orders of magnitude fewer independent data points, therefore this issue is much less pronounced. The $\Delta$BIC<sub>lm</sub> values support this, finding that the model CL is a significantly better fit than any other models (with a value >10 lower than the next statistically favoured model CL + ECC), in line with our earlier assessment of the residual image. Nevertheless, we acknowledge that the conclusion from the BIC analysis may also reflect that our asymmetric models are also not a perfect representation of the data.

The lower three sections of Table 8 assess how well the model images reproduce the observationally derived parameters associated with points O1, O2, and O3 from Section 3.3 (when convolved with the same beam), noting that point O4 is beyond the scope of the sub-mm models considered here. As previously discussed, all models are consistent with the broad structure of the disc noted in point O1, such as the peak radius, inclination and position angle. Table 8 shows that all four models reproduced the disc width to within $3\sigma$, consistent with the Band 7 image (although all model discs were wider than the Band 7 image, measured as 43 au). We also note that the Band 7 SW-NE flux difference $\delta F = 1.33 \pm 0.18$ mJy is strongly consistent with the clump flux value, $F_S = 1.7 \pm 0.6$ mJy (i.e. within $1\sigma$ for both clump models). Only models CL and CL + ECC are able to account for the major axis offsets within $2\sigma$, however, this is not a strong conclusion given that the values for the model SYM (i.e. without any offset) would also be deemed reasonable (i.e. all within $3\sigma$), and since the clump location is not tightly constrained on the SW ansa (see Appendix C). Further, only models CL and CL + ECC are able to reproduce the major axis flux ratios $f_{82\text{au}}$ and $f_{200\text{au}}$ within $2\sigma$, with both model ECC and SYM providing poorer estimates of these parameters. This comparative analysis thus suggests that despite a few shortcomings, the CL and CL + ECC asymmetric models have been able to model the major axis asymmetries reasonably well, with the ECC model proving less adequate, though still preferable to the symmetric model.

## 4.3 Modelling summary

Of our four models, two are able to reproduce the observations to within the noise levels (see observational parameters in Tables 7 and 8), and consistency was shown with all the broad disc parameters and our measurements in Section 3 in all four cases. The asymmetric parameter analysis demonstrated that the models which introduced a clump on the SW ansa of the debris disc were favourable, and by further analysing the respective residual images and $\chi^2$ values from the model visibilities we find that the model with a single clump is the best fit to our data. Although, we note that an assessment of the BIC parameters may suggest differently to the above, with the models with fewest parameters being shown to be consecutively more favourable by this metric. Nevertheless, in conjunction we consider our model that is a symmetrical disc with a clump model (model CL) to be the most favourable scenario of those considered.

We note however that the clump modelled was a simple 2D Gaussian. A more complicated clump profile may better model this system with more success given that this distribution may be too simplistic to interpret the physical origin of the extra emission (i.e. this may be azimuthally extended around an orbit if the emission is due to a collision within the belt, which we consider in Section 5.2.2). On the other hand, whilst the disc models without a clump (either fully symmetrical or having an underlying constant eccentricity distribution) were poorer fits to the data, an eccentric model with both free and forced eccentricity terms, or one with an outwardly falling eccentricity (e.g. if due to the forced eccentricity of a planet sculpting the inner edge) may provide a better fit than the constant eccentricity modelled here (for example, similar to the models of MacGregor et al. 2017; Kennedy 2020). We suggest future work may wish to conduct detailed modelling of this scenario to explore the origin of the asymmetric emission of this debris disc, which we further discuss in Section 5.2. Thus, whilst our single clump model provides a reasonable model of our data, we leave open the possibility that more complicated clump profiles or multicomponent eccentric distributions may better and consistently interpret our observations.

# 5 DISCUSSION

In Section 3, we presented an analysis of the Band 6 and 7 images for $q^1$ Eri. From this we were prompted to investigate the disc in Band 7, given the observed asymmetries, and presented different models for the disc morphologies in Section 4. Here we interpret these results and discuss them in the context of the wider $q^1$ Eri planetary system.

## 5.1 A massive and broad debris disc

### 5.1.1 How big are the largest planetesimals, and what does this imply for the disc mass?

As outlined in Section 1, $q^1$ Eri stands out as having the brightest debris disc of the closest 300 sun-like stars, despite its age being $1.4 \pm 0.9$ Gyr (Marmier et al. 2013). This could mean that $q^1$ Eri’s disc is somehow extreme among the population (e.g. in having started with an unusually large mass in planetesimals, see Kral et al. 2015), or alternatively that there is more dust in this system than expected from steady state collisional models (e.g. due to a recent transient event). The disc in this system is broad, extending between 34 and 134 au with a peak at $\sim 81.6$ au. Despite being broad, $q^1$ Eri’s peak emission radius in the sub-mm is in accordance with the population of 26 SMA (Sub-Millimetre Array) and ALMA resolved discs in Matrà et al. (2018), which (for $L_* \sim 1.59L_\odot$) predicts a sub-mm disc radius $r_0 = 73 \pm 12$ au (i.e. within $1\sigma$ of our measurement). Further, this sub-mm resolved radius is entirely consistent with the archival far-IR Herschel radius of $R_{\text{disc}}=81.1_{-1.3}^{+1.8}$ au (Marshall et al. 2021). Whilst our measurement of the disc width is however narrower than the broad $\Delta R_{\text{disc}}=71.1_{-13.3}^{+1.9}$ au value found by Marshall et al. (2021), this is likely due to the different manner in which we assess the width (i.e. we note that the extent of the disc is much broader than the $\sim 43$ au width we define from the half-maximum emission). Further, the width measurements of Marshall et al. (2021) are also likely unresolved, which for example, may also be influenced by emission from an inner warm component that we discuss further in Section 5.3.

To constrain the disc’s total mass, we first quantify the size of planetesimals that must be feeding the disc’s collisional cascade, assuming this to be in steady state. To do so we assume the number of bodies with a given size follows a power law with slope $\alpha = -3.5$ (Dohnanyi 1969) from dust with a maximum modelled size of $D \sim 2$ cm, in which there is a mass $M_{\text{dust}}=0.028 M_\oplus$ (determined from model CL), up to a largest planetesimal size $D_{\text{max}}$, which results in a total mass of $(D_{\text{max}}/D_{\text{mm}})^{0.5} M_{\text{dust}}$. The largest planetesimal size is determined by setting the collisional lifetime $t_{\text{coll}}$ (see equation 3 of Wyatt 2008) to be equal to the age of 1.4 Gyr (although we note this age is uncertain at the $\sim 50$ per cent level). This gives a lower limit largest planetesimal size of $D_{\text{max}} \sim 1$ km and a lower-bound total disc mass estimate of $M_{\text{disc}} > 8 M_\oplus$. These are lower limits because there may be even larger planetesimals in the disc, which have not yet contributed to the collisional cascade, and so are not required to be present (e.g. Wyatt & Dent 2002). These values are found for a fixed $Q_D^* \sim 163$ J kg$^{-1}$, consistent with $\sim 1$ km Basalt planetesimals for both
Benz & Asphaug (1999) and Stewart & Leinhardt (2009), a relative velocity $v_{rel}=103.9 M_*^{0.5} h r^{-0.5}=620 \pm 30$ m s⁻¹ (using equation 10 of Matrà et al. 2019, with our model CL determined value for $h$), the stellar mass $M_* \sim 1.1 M_\odot$, at the peak brightness radius $r = 81.6 \pm 0.5$ au, and with a belt width $dr \sim 43$ au. Given the uncertainties and assumptions (e.g. constant power law slope, emission being dominated by grains with diameters smaller than 2 cm), these lower limits are only approximate, but they are consistent with the ~81.7 $M_\oplus$ disc mass determined by the detailed collisional modelling of Schüppler et al. (2016). Indeed, given that this collisional modelling included planetesimals up to 100 km in size and resulted in a disc mass of 10 higher than our ~8 $M_\oplus$ lower bound, the relationship between total disc mass and maximum planetesimal size (i.e. $M_{disc} \sim D_{max}^{0.5}$) suggests these two estimates may be strongly consistent. This shows that, while the required disc mass is certainly large, it is not unreasonably so, suggesting that the disc's unusual brightness can still be explained within the context of steady state collisional erosion.

Although we have thus far assumed that planetesimals are stirred at $t=0$ (i.e. born stirred), it is also possible that planetesimal collisions were delayed until a time $t_{delayed}$ in accordance models (e.g. self-stirring or planet-stirring Kenyon & Bromley 2004; Mustill & Wyatt 2009). Such delayed stirring may commence between 10 s of Myr to Gyr time-scales, and thus only recently in the $q^1$ Eri debris disc. fig. 5 of Wyatt (2008) shows the fractional luminosity evolution with one such delayed stirring model for a planetesimal belt extending from 30 to 150 au (i.e. similar to the extent measured here in the sub-mm for $q^1$ Eri). If this belt is stirred after just ~20 Myr, then by ~1.4 Gyr of evolution this model finds $f \sim 3 \times 10^{-4}$ (consistent with the broad planetesimal belt of $q^1$ Eri, $f = (2.44 \pm 0.08) \times 10^{-4}$, see Section 2.3) for a disc with an initial mass of one minimum mass solar nebula (1.0 MMSN). Whilst the delayed self-stirring models initiate the collisional cascade after the formation of 2000 km planetesimals within the disc, planets outside the disc can initiate the cascade in planet-stirring models. This suggests that belts as bright and old as $q^1$ Eri may only contain as much mass as the MMSN, which seems reasonable since at least this mass must have been present in the Solar system. Given the range of initial disc masses, extents and plausible range of delayed stirring time-scales that may be consistent with the observed fractional luminosity of $q^1$ Eri, the possibility that the disc was not stirred at $t=0$ adds uncertainty to any disc mass calculations. We consider the mass of planetesimals in the belt further in Section 5.1.2.

### 5.1.2 Are there further constraints from the vertical distribution?

The scale height of the disc is resolved in all of the models presented here, with an average value ~0.05, comparable with that determined for other debris discs of 0.02–0.12 (Hughes et al. 2018; Kennedy et al. 2018; Daley et al. 2019; Marino et al. 2019; Matrà et al. 2019). This agrees with both the observational and modelling conclusions that we may have resolved the vertical dust distribution. The vertical distribution in discs can be interpreted as due to dynamical stirring from massive planetesimals. Equating the relative velocity, $v_{rel}$ (as estimated in Section 5.1.1 at the peak emission radius ~81.6 au), to the escape velocity of the largest bodies (e.g. Daley et al. 2019), for an asteroid density of $\rho=2.7$ g cm⁻³, we find lower bounds on the mass and size of the bodies providing the stirring as $M_{body} > 2.3 \pm 0.5 \times 10^{21}$ kg and $D_{max} > 1200$ km, respectively ~20 per cent the mass and ~50 per cent the size of Pluto. Based on Ida & Makino (1993), we applied the same analysis as Matrà et al.

(2019) from which we use their equation (12), our model determined value of $h$, and this lower limit on $M_{body}$ to derive a surface mass density $\Sigma=0.039 M_\oplus$ au⁻². Given this value of $\Sigma$ is around 40 times larger than the MMSN at 81.6 au, this suggests that bodies larger than 1200 km indeed may be necessary in order to stir the disc within the age of the system (i.e. to reduce the surface mass density to levels that may be more plausible).

This dynamical limit on planetesimal sizes could suggest that these may exist in the belt a factor of $10^3$ larger than the ~1 km sized lower limit predicted by collisional replenishment. However, if the assumed $\alpha = 3.5$ Dohnanyi (1969) size distribution continued up to planetesimals of this size, then the disc mass would be $M_{disc} \sim 220 M_\oplus$, a factor of ~3 higher than the prediction of Schüppler et al. (2016) and over an order of magnitude larger than the lower limit derived from the collisional lifetime and age of the system. This would not violate our earlier calculation since these $\gg 1$ km planetesimals could be abundant in the disc without having collided within the age of the system. Whilst this higher total disc mass is still consistent with the dust mass measurements of protoplanetary discs (for example, see Andrews & Williams 2005; Ansdell et al. 2016; Cieza et al. 2019), high debris disc masses, i.e. those in the range of 100–1000 $M_\oplus$, become problematic since these would require a very high efficiency of primordial dust being incorporated into these larger planetesimals (an example of the so-called ‘disc mass problem’, see Krivov et al. 2018; Krivov & Wyatt 2021). Nevertheless, it might still be possible to explain the observed level of stirring by embedded bodies while retaining a lower disc mass (i.e. if the largest bodies are less frequent than the $\alpha = 3.5$ size distribution would predict). If instead the size distribution had a slope of $\alpha = 3.7$, then even with these 1200 km bodies, the total disc mass estimate would be reduced from ~220 $M_\oplus$ to ~6 $M_\oplus$, since the total number of these would be greatly reduced. On the other hand, the size distribution could be truncated, being much steeper for planetesimals larger than a few kilometres, and shallower in other regions. We estimate the effect that this can have on the derived disc mass using equation (9) of Krivov & Wyatt (2021), and find that the total mass can be reduced by an order of magnitude from ~220 $M_\oplus$ to ~17 $M_\oplus$, if based on a triple power law size distribution, with $q_{med} = 4$, $q_{big} = 3$ and $q = 3.5$. Alternatively the measured vertical scale height could be due to other dynamical interactions (e.g. stirring by a planet internal or external to the belt, or a recent stellar fly-by) or even be a remnant of the primordial disc (e.g. if this disc was born stirred, Booth & Clarke 2016). In summary, whilst $q^1$ Eri at an age of ~1.4 Gyr is an outlier in terms of its brightness, it need not be an outlier in terms of its disc mass unless many much larger planetesimals are present.

### 5.1.3 What is happening at the disc inner-edge?

Although this disc is radially broad, it also has a sharp inner edge. To within $3\sigma$, model CL has an inner edge $\Sigma(r) \propto r^{-5.6}$ (see Table 7). Kennedy & Wyatt (2010) showed that disc inner edges shaped by collisional processes have much shallower inner edges of $\Sigma(r) \propto r^{7/3}$. That our models find a steeper inner edge favours a planetary carving scenario (see, e.g. Chiang et al. 2009), though we cannot rule out that a steep edge simply reflects the primordial planetesimal distribution.

Already known in this system is the exo-Jupiter, $q^1$ Eri b, but with an orbital radius at ~2 au this planet cannot be responsible for carving the inner edge (which we place at ~57 au based on the half-maximum derived inner edge from Table 5). The lack of any gaps in the disc beyond ~57 au suggests there is no evidence for additional
~planet-mass bodies, though these would also be unlikely to explain the disc's inner edge. Rather, if an as yet unseen outer planet is sculpting this inner edge, it would likely have a semimajor axis close to ~57 au. Very massive bodies (even at this radius) can introduce long-term linear trends into RV measurements, however no such linear trend is seen by Marmier et al. (2013) for $q^1$ Eri. We therefore report an upper limit on such a planet's mass of <11 $M_{Jup}$ based on the inferred ~57 au semimajor axis, the 13 yr measurement baseline, and the radial velocity uncertainty of ~9 $m$s$^{-1}$, although we note that larger masses may be possible for unfavourable orientations.

There are two constraints we consider to set a lower-bound on the mass of a planet near the disc inner edge. The first is that the planet must be sufficiently massive to truncate the inner edge of the planetesimal belt within the ~1.4 Gyr age of the system. This results in a lower limit of $M > 2 M_⊕$ (see equation 3 of Shannon et al. 2016). The second constraint is that the planet should be more massive than the disc it is carving or else it would be forced to migrate during this interaction (Kirsh et al. 2009). This results in a stricter lower limit of $M > 8 M_⊕$.

A planet with $M \sim 8 M_⊕$ has a secular time-scale at the peak brightness radius of the disc that is shorter than the ~1.4 Gyr age of the system (i.e. the time-scale over which such a planet would impose structure on the disc, see equation 17 of Pearce & Wyatt 2014). We may therefore infer that any planet close to the inner edge is unlikely to have an eccentricity significantly greater than that of the inner edge. This implies that such a planetary orbit would be close to circular (however we note that if detailed modelling found a better fit with a more realistic eccentric profile, this conclusion would need to be revised). We can also infer that such a planet should be currently aligned with the disc mid-plane (or more accurately that secular perturbations would have caused the disc mid-plane to become aligned with the planet), and that the planet should not have started out on an orbit that was inclined to the disc by more than half the scale height (since the alignment process imposes vertical structure on the disc). These two constraints could argue against such a planet having been scattered out from a formation location closer to the star due to interactions with other planets (such as $q^1$ Eri-b), since it is unlikely that such scattering would result in a coplanar, eccentric planet. However, if such a planet was comparable in mass to the disc, its orbit could be circularized and become aligned with the disc without imprinting significant structure on the disc (e.g. Pearce & Wyatt 2015). Alternatively the planet could have formed in situ or migrated out through disc interactions while retaining a nearly circular co-planar orbit.

## 5.2 Origin of disc asymmetries

Our modelling found that the disc's asymmetry is most consistent with a broad clump on the SW ansa with a total flux $1.7 \pm 0.6$ mJy. The projected offset of the clump puts it at a radius of ~80 ± 15 au if it is in the disc mid-plane (i.e. determined from the clump position offsets taking account of the disc position angle and inclination). The clump's major and minor axis standard deviations is ~3.2 arcsec × 2.5 arcsec, and its position angle ~40° (i.e. angled towards the NW minor axis along the disc inner-edge). This is consistent with the observational analysis which inferred a clump located near the inner edge in the SW ansa with a flux of $1.33 \pm 0.29$ mJy (from the difference between the fluxes in the NE and SW ansae). We also found the inner edge to have a radial offset, and so too in the scattered light HST imaging. In this section we interpret these asymmetric measurements as either due to a disc-planet interaction, a recent collision, or due to extra-galactic background sources.

### 5.2.1 A planet-driven asymmetry?

The discussion in Section 5.1.3 concluded that the inner edge could have been sculpted by a planet. Such a planet could have formed closer in and migrated outwards via the exchange of angular momentum with planetesimals in scattering events (i.e. planetesimal-driven migration, see Fernandez & Ip 1984; Ida et al. 2000; Gomes, Morbidelli & Levison 2004; Ormel, Ida & Tanaka 2012), via interaction with the gas-disc in which it formed, or indeed by interactions with other planets in the system (e.g. with the known $0.93 M_J$ planet $q^1$ Eri b). During this migration its resonances would have swept through the planetesimal disc, and some planetesimals could have become trapped. The resonances that become populated in this process would depend on the planet's mass and migration rate and can result in the planetesimal distribution being clumpy (Wyatt 2003). The clump observed on the inner-edge could therefore be explained by planetesimal trapping in the 2:1 resonance. We previously inferred that a planet may have sharpened the inner edge, however, if an asymmetry was driven by such a planet's migration, this motion must have stalled as the planet migrated outwards. For example, this could occur if the disc ran out of angular momentum, which would imply that the planet mass is comparable with the disc mass, i.e. $M_{pl} \sim M_{disc} > 8 M_⊕$. Interesting for this scenario, is that this would be consistent with the evolution of Neptune in our own Solar system, inferred to have stalled during its migration following the resonant sweep-up of planetesimals in the Kuiper Belt (Fernandez & Ip 1984).

While an eccentric distribution of planetesimals is not needed to explain the sub-mm major axis asymmetry, smaller dust particles are more strongly affected by stellar radiation pressure. This additional force places small dust on to eccentric orbits with periastra near to the location of where they were released from their parent planetesimals (presumably in collisions). Therefore, small dust created in collisions in the clump would have apastron in the NE of the disc (see Wyatt 2006). This provides a consistent explanation for the HST image which found dust extending to greater distances in the NE (see Section 3.3). Thus a planet that arrived at the inner edge of the disc by planet migration could potentially reproduce the observed flux, size and position of the measured asymmetry, in both the sub-mm and scattered light, though further modelling of this scenario is required to investigate this.

Nevertheless, we cannot rule out a scenario in which the asymmetric belt is explained by secular perturbations from an eccentric planet. For example, if such a planet migrated towards the inner edge of the belt whilst retaining a sufficiently high eccentricity, this would force an eccentricity into the belt since the planet's secular time-scale is less than the age of the system. However, since the forced eccentricity of a belt is due to the combination of both the perturbing body's semimajor axis and eccentricity (i.e. to first order, $e_{forced} = 1.25e_{pl}a_{pl}/a_p$, where $e_{pl}$ and $a_{pl}$ are the eccentricity and semimajor axis of the planetary perturber, and $e_{forced}$ and $a_p$ are the forced eccentricity and semimajor axis of the perturbed planetesimals), there is a degeneracy between the perturbing body's semimajor axis and its eccentricity. Whilst we are therefore unable to tightly constrain these parameters, from the equations (5)–(8) of Pearce & Wyatt (2014), it can be shown that a planet orbiting at 30 au with an eccentricity of 0.1 would drive the belt's inner and outer edges to eccentricities of 0.064 and 0.036 respectively (consistent with our derived $R_{inner}$ and $R_{outer}$ parameters). Equivalent eccentricities could be produced from a planet at 60 au with an eccentricity of just 0.05, i.e. the same semimajor axis as the half-maximum derived inner edge (see Table 5). We suggest as further work that such a scenario could be
modelled in detail to better constrain required planetary parameters, and to determine if this provides a realistic model for the broad debris disc and its asymmetries.

### 5.2.2 A recent collision?

Alternatively, this clump could have formed from a recent massive collision in the belt (such as those explored in Kral et al. 2015). This would not require planetesimals to be trapped in resonance (although such a collision could have been at the resonant location). Here we consider how large a parent body would need to be for a single collision to reproduce the observed asymmetry. By scaling the total dust mass in the disc by the ratio of the 1.33 mJy excess flux in the SW ansa to the total disc flux, this implies the clump has a mm dust mass of ~0.0029 $M_⊕$ (or ~30 per cent larger if the 1.7 mJy clump flux from model CL had been used). For this amount of dust to have formed via a single collision, the parent body must have had a diameter of at least $D≈2300$ km for a density of 2.7 g cm⁻³. Since this would require the collision to fragment the parent body entirely into mm-sized grains, whereas a range of fragment sizes is more likely, the parent body would most likely be significantly larger than 2300 km.

While Section 5.1.1 concluded that no bodies larger than $D_{max}$ ≈ 1 km need to be present in the disc, if the $\alpha = 3.5$ Dohnanyi (1969) size distribution continues from $D_{max}$ to 2300 km then the collisional lifetime of 2300 km bodies can be estimated to be a factor of $\sqrt{2300 \text{ km}}/D_{max}$ longer than the age of the system, i.e. ~40 Gyr. Extrapolating from the number of 2 cm bodies present in the 0.028 $M_⊕$ disc, an $\alpha = 3.5$ size distribution implies that there would be ~10⁴ bodies of size ~2300 km present in the disc, and therefore we would expect collisions to occur every ~4 Myr. Clumps disperse after formation, and fig. 7 of Jackson et al. (2014) suggests that this process might take ~1000 orbits, i.e. ~0.5 Myr at ~60 au. Comparing this dispersal rate to the expected collision rate, this calculation suggests that although infrequent, clumps are not implausibly rare in this disc, demonstrating a ~10 per cent probability that we might observe a clump in the disc at any one time if such massive planetesimals exist in the disc and collisions are energetic enough to fragment these.

However, there are two arguments against this. Most significantly, collisions between even 2300-km-sized bodies may not be sufficient to cause the observed clump, since these would not be 100 per cent efficient at converting parent bodies into mm-sized dust, i.e. we would very likely need even larger collisional bodies. The second is that the disc mass required to get even a 10 per cent detection fraction in this scenario is ~300 $M_⊕$, which as discussed in Section 5.1.2, may be problematic. Such a high mass is inevitable if there are to be enough 2300 km bodies for collisions to occur frequently, although we note here that the presence of such large bodies might help to explain the observed scale height of the disc. If we were to instead assume the size distribution was flatter, then we could raise the clump detection probability, however this would require an even higher disc mass. For example, to raise the detection probability to order unity, there would instead need to be 10⁵–2300 km bodies; however, such a disc would then have a mass ≳3000 $M_⊕$. This scenario would then be problematic for the reason outlined in Section 5.1.2, i.e. the disc mass would then exceed the solid mass available to form planetesimals in protoplanetary discs. Therefore, whilst a collisional origin for this asymmetry could explain the observed asymmetry in both the sub-mm and scattered light, it is perhaps unlikely.

### 5.2.3 Clumpiness from extra-galactic emission?

An extra-galactic sub-mm galaxy (SMG) in the SW ansa is one possible explanation for the asymmetry seen in the ALMA images. Perhaps the strongest argument against this is that such an SMG would not explain the extended scattered light emission in the NE, and the simplest explanation for the asymmetries measured in the sub-mm and in scattered light is that they are caused by a single phenomenon. Nevertheless, a more complicated scenario involving two phenomena cannot be ruled out. Thus, we also consider the probability of detecting an SMG coincident with the disc. From the SMG counts of Simpson et al. (2015), we find a detection probability (within a projected 200 au from the image centre with $\delta F \sim 1.33$ mJy) of $P(\text{SMG}) \sim 0.2$ per cent. Whilst not implausibly rare, this casts doubt on this interpretation. We also found this clump to be resolved or azimuthally broadened, which would require this SMG to either be an extended source, or be due to multiple sources nearby on the sky. However, the frequency of SMGs appearing nearby around this flux level, and the number of large SMGs (i.e. >0.7 arcsec, see Section 3.1.3) are both rare, additionally arguing against an SMG interpretation.

The spectral index has also been used to discriminate between possible origins of sub-mm emission (Su et al. 2017; Booth et al. 2019), since this differs between optically thin dust from planetesimal collisions, $\alpha_{mm} \approx 2-2.5$, and extra-galactic emission, $\alpha_{mm} \approx 3-4$. In Section 3.1.2 we derived the spectral index for the full disc of $\alpha_{mm} = 2.34 \pm 0.29$, which is consistent with measurements of other debris discs Ricci et al. (2015). By instead considering the spectral index due to the flux difference between the SW and NE ansae, $\delta F$, in Bands 6 and 7 (see Table 6), we find only a weak constraint on the spectral index of the clump $\alpha_{mm} \sim 9 \pm 4$ (i.e. consistent with the spectral index for thin dust or an SMG). We therefore cannot claim the spectral index of $q^1$ Eri as inconsistent with either optically thin dust or extra-galactic emission, either from its total sub-mm fluxes, or those of the clump. Whilst we cannot exclude the possibility that the clump is due to extra-galactic emission, this interpretation seems less likely than a planet-driven asymmetry (see Section 5.2.1) since this cannot simultaneously explain the scattered light asymmetry and such an SMG would have a very low occurrence rate. Further ALMA observations could ascertain if this inner edge asymmetry is co-moving with the belt in the future and definitively rule this out as due to background SMG emission based on the signal-to-noise of the clump. For example, if this can be measured with an SNR>3 with the same beam in 2018, given the proper motion of the system is ~200 mas yr⁻¹, and the precision to which a Gaussian can be centred on this emission is 0.5 × beamsize/SNR, measuring a >3 σ change in this location will require waiting a further ~2 yr from the 2018 observations (i.e. such a change may now be measurable).

### 5.2.4 Interactions with the ISM?

Although interactions with the interstellar medium (ISM) are not expected to affect the orbits of the sub-mm grains, the ISM can significantly influence the orbits of smaller micron-sized dust observed in the scattered light. Indeed, this has been shown to be the case with HD 15115, which shows a swept-back asymmetric disc in scattered light as observed with HST (for which the east-west emission shows a radial offset of 2, see Kalas, Fitzgerald & Graham 2007), however when observed with ALMA, there are no observed radial offset asymmetries (see MacGregor et al. 2019). This suggests that there does not to be a single mechanism to interpret sub-mm and scattered light asymmetries simultaneously. Whilst a single mechanism may
still however be preferable for probabilistic reasons, we cannot rule
out the possibility that the asymmetries observed in the scattered light
and the sub-mm are independent, i.e. with both ISM interactions and
planet-disc interactions affecting these respectively.

5.3 What is the inner warm component?

q¹ Eri has previously been found to be consistent with having a
multicomponent debris disc, comprising an outer cool belt with r<sub>bb</sub>
~60 au, and an inner warm belt with r<sub>bb</sub> ~ 10 au (Kennedy & Wyatt
2014). Schüppler et al. (2016) likewise modelled this system, and
found two disc components, an inner warm belt between 3 and 10 au,
close to the known exo-Jupiter at ~2 au, and an outer cool belt
between 75 and 125 au. Although our ALMA images clearly resolved
the disc's outer belt, the only other significant Band 7 emission
detected was unresolved and coincident with the star with a flux
of 169±22 µJy. Given that our flux distribution model predicts the
856 µm stellar photospheric flux to be ~99 µJy (with a 2 per cent
uncertainty), this suggests that at the ~3 σ level, we may have
detected additional emission coincident with the location of the star
with a flux F<sub>inner</sub>=70±22 µJy. This is consistent with what might be
expected for the emission of an inner warm component (see Fig. 3),
however we note that the uncertainties associated with this modelled
inner component are large since the slope of the inner modified
blackbody are determined entirely by the two ALMA Band 6 and
Band 7 data points (i.e. there are no resolved measurements of
the inner component in the near, mid- or far-infrared). Since the
angular size of emission at ~10 au from the star at a distance of
17.34 pc would be ~0.6 arcsec, i.e. comparable with our beam size,
emission internal to this radius would be largely unresolvable from
the stellar emission in our images. This emission could arise from a
planetary belt inside 10 au, be due to additional stellar emission,
or image noise coincident with the stellar position. We discuss each
of these in turn.

We first consider whether this emission is evidence of a warm
planetesimal belt internal to 10 au, which if present, may resemble
the Solar system's asteroid belt, though with higher mass and external
to the orbit of the known exo-Jupiter, q¹ Eri b. Belts consistent with
this have been previously inferred towards q¹ Eri by Kennedy &
Wyatt (2014) with a blackbody radius r<sub>bb</sub> =10 au (and consistent
with the flux distribution presented here in Fig. 3,) and at 3 au by
the flux distribution modelling of Schüppler et al. (2016) from the
mid-IR data. Even though this blackbody radius is relatively well
constrained, such an assessment does not come without uncertainty.
The value determined by Kennedy & Wyatt (2014) is dominated
by emission from the mid-IR flux measurements (i.e. not the
longer wavelength ALMA data), and such radius estimations can
significantly underestimate the true radii of discs (see equation
8 of Pawellek & Krivov 2015). For example, for a 50 per cent
astrosilicate grain-50 per cent ice dust composition around a star with
a luminosity of L=1.59 L<sub>☉</sub> (consistent with q¹ Eri, see Section 2.3),
the true belt radius could be larger than the blackbody radius by up
to a factor of ~4.6. If so, the warm emission analysed by Kennedy &
Wyatt (2014) could have a contribution from dust as far out as ~50 au,
i.e. dust on the inner edge of the disc, and so be potentially associated
with the clump rather than an inner planetesimal belt.

However, for two reasons we believe this to be highly unlikely.
First, we used the IRS data to consider the possibility of the warm
emission being associated with the inner edge clump at ~60 au (3.5
arcsec). This clump would lie within the IRS slit for the 'long'
modules (λ >14 μm), and could feasibly shift the centroid of IRS
emission as it changes from star to disc-dominated (as seen in the

flux distribution in Fig. 3). However, the spatial profiles and PSF-
subtracted residual images in the CASSIS data base (Lebouteiller
et al. 2011) appear consistent with a single unresolved source across
all wavelengths, so it is unlikely that any clump has contributed
flux to the IRS spectrum. Secondly, we note that applying the same
correction factor to the main belt would lead us to predict its belt to
be at ~276 au, whereas it is observed at 81.6 au, a factor of just ~1.4
larger than r<sub>bb</sub>. Although the inner and outer belt correction factors
may be independent (e.g. if they are compositionally different), this
suggests that such a high correction factor may still be unrealistic.
If applied accordingly, such a correction factor would predict the
inner planetesimal belt to have a radius of ~14 au, and (at the flux
observed) be difficult to resolve from the star with this ALMA data.
Therefore, since this emission seems unlikely to have been confused
with the outer planetesimal belt, we cannot rule out the presence of
a high mass planetesimal belt internal to 10 au, which like in Marino
et al. (2018), we have shown may be detectable in the sub-mm.

If such a belt was present, and approximated a narrow Gaussian
ring, equation (5) of Matrà et al. (2020) predicts that in visibility space
this would produce a Bessel function with a first null at ~130 kλ. This
might suggest that such a function may be visible near such baselines
(for example, in Fig. A1). However on inspecting this visibility data,
given the SNR of emission at these longest baselines (i.e. shortest
angular scales) no such function is visible with significant emission
at the ~70 µJy level, even when the bin sizes are increased to raise
the SNR per bin. Given this warm emission is faint, this analysis
therefore cannot exclude the existence of a narrow inner belt inside
10 au.

One alternative to explain this emission could be that it has origi-
nated in sub-mm stellar flaring events. Although common around M-
type stars (see MacGregor, Osten & Hughes 2020), second-minute
time-scale sub-mm variability can be detected towards solar-type
main-sequence stars at the level of tens-hundreds of µJy, consistent
with the warm excess measured here (Burton et al. in preparation).
Although between the two Band 7 epochs (separated by two years) we
found consistent flux measurements (within 10 µJy) at the location
of the star, we did not explore stellar variability in detail in this
analysis. Therefore, we equally cannot rule out the possibility that
this warm excess is due to stellar variability from sub-mm flares.

Finally, we note that we cannot rule out this emission as due to
noise coincident with the stellar location, although at the 3 σ level, or
contamination with the disc flux. Although this appears less plausible
than the other two interpretations, if the inner edge of emission on
the minor axis is only marginally resolved from the star, then this
could have contributed to the measured excess. Therefore, to better
understand the nature of this inner component, we note that further
high-resolution scattered light, mid-IR and sub-mm measurements
(e.g. with SPHERE, JWST and higher resolution and deeper ALMA
imaging) are necessary.

5.4 CO Gas mass

In Section 3.4 we demonstrated that CO spectral signatures for
the J = 3–2 transition line are not present in our ALMA data.
This non-detection of CO resulted in an upper limit on the CO
flux of 24.0 mJy kms⁻¹. This can be used to derive an upper limit
on the gas mass from the excitation conditions of the gas set by
the radiation environment, electron density and kinetic temperature
(Matrà et al. 2015). We compute the level populations of the CO
J = 3 rotational level using an NLTE (Non-Local Thermodynamic
Equilibirum) code including fluorescence (Matrà et al. 2015, 2018).
For three temperatures 10 K, 50 K and 250 K (covering a range of
<figure><img src="image_14.png" /><figcaption>**Figure 12.** CO gas mass as a function of the electron density, for the three temperatures of 10, 50, and 250 K.</figcaption></figure>

values appropriate for debris discs), together with the stellar flux (see Fig. 3) and the peak emission radius of ~81.6 au, Fig. 12 shows the estimated CO gas mass as a function of the electron density in the disc. We plot a range of collision-partner densities, from low densities where molecular excitation is dominated by radiation (e.g. fluorescence, from electron excitation by starlight, followed by decay through higher J rotational levels), through to larger densities where line excitation is instead dominated by collisions, and as such the line populations are in local thermodynamic equilibrium where the line fluxes depend only on temperature. From this we can then set an NLTE upper bound CO gas mass of $M_{CO} < 4 \times 10^{-6} M_{\oplus}$, the constraint for low collider densities.

Based on the assumption that icy planetesimals create second-generation CO gas through planetesimal collisions, Kral et al. (2017) predicted a CO gas mass in $q^1$ Eri's disc to be $1.1 \times 10^{-7} M_{\oplus}$, based on the predicted collision rate of bodies from mid- and far-IR flux measurements, a stellar luminosity and temperature consistent with values in this work, and a disc radius of 105 au. Such a predicted gas mass is consistent with the upper limit derived here, being two orders of magnitude lower than the upper bound mass given above. From equation (2) of Matrà et al. (2017) we can use this upper limit CO mass alongside the stellar mass, luminosity, disc fractional luminosity, disc peak emission radius and belt width, with a CO photodissociation time of 120 yr (see Visser, van Dishoeck & Black 2009) to estimate an upper bound on the fraction of planetesimals composed of CO + CO₂, which is $f_{CO+CO_2} < 95$ per cent. This upper bound is consistent with the Solar system ($f_{CO+CO_2} \sim 10$ per cent). Nevertheless, higher sensitivity ALMA data could place tighter constraints on this upper limit CO gas mass, and if measured with sufficient depth, detect any CO, if present. Furthermore, the presence of different chemical species can be constrained by ALMA (for example, CI, OH, HCN, $N_2H+$). With measurements of such other gaseous species, in conjunction these may provide us with a far better understanding of the composition of this debris disc.

## 6 CONCLUSIONS

We have presented Band 6 and 7 ALMA observations, and HST scattered light observations of the $q^1$ Eri debris disc. Until now the debris disc architecture has been interpreted from low-resolution thermal imaging (~4 arcsec) and scattered light HST imaging with

a poor inner working angle (~2.25 arcsec). Here we have explored the full extent of the $q^1$ Eri debris disc at subarcsecond resolution, placing bounds on disc features such as the inner and outer edges, peak emission radius, extent, inclination, position angle, the vertical dust distribution, and placed an upper limit on the CO gas mass and on the ice mass fraction of CO + CO₂ present.

From the 1.4 Gyr age of the system, we placed lower limits on the size of the largest planetesimals in this disc from their collisional lifetimes, finding $D_{max} > 1$ km, leading to an estimate of the mass of this disc of $M_{\text{disc}} > 8 M_{\oplus}$, suggesting that whilst this disc is an outlier in brightness, it need not be an outlier in terms of its mass. From observations and modelling that tentatively determined this disc to be resolved in its vertical direction, we investigated the size and mass of bodies that could be responsible for stirring the disc. This interpretation required bodies with a size at least 1200 km significantly above the lower limit on the largest planetesimal size, although the disc could alternatively have been born stirred (or planet-stirred), in which case such large planetesimals would not be necessary.

We have demonstrated through image analysis and modelling that a series of asymmetric features exist in the sub-mm emission: with a total flux asymmetry between the SW and NE sides, and a radial offset towards the SW ansa in the major axis. Our modelling found this disc is most consistent with an axisymmetric disc with an emission clump on its SW inner edge.

At just ~2 au, the known planet $q^1$ Eri b is too close in to affect the observed outer main belt, and we discuss how the broad disc may be connected to other aspects of this planetary system, such as other potential planets and planetesimal collisions. One interpretation of this system is that the belt of planetesimals had its inner edge carved by a planet which formed a clump on the inner edge whilst migrating, similar to the inferred evolution of Neptune. This scenario may simultaneously explain the inner edge sharpness, the SW ansa sub-mm clump and the scattered light asymmetry on the outer edge in the NE.

By assessing the emission coincident with the star, and by comparing this with previous modelling and our flux distribution, we show that there is tentative evidence for an inner warm component. This may be due to a belt of planetesimals closer to the planet $q^1$ Eri b, between radii 3 and 10 au, and thus a high mass Asteroid Belt analogue.

Throughout this work we have suggested further analysis that could be undertaken to better understand this disc. New discoveries and confirmations of interpretations laid out here could be made with more detailed modelling, and future observations with instruments/observatories such as ALMA, HST, SPHERE, and JWST.

## ACKNOWLEDGEMENTS

We thank the anonymous reviewer for their comments that improved the quality of this work. JBL is supported by an STFC postgraduate studentship. SM is supported by a research fellowship from Jesus College, Cambridge. GMK is supported by the Royal Society as a Royal Society University Research Fellow. The research of OP is funded through the Royal Society Dorothy Hodgkin Fellowship. TDP is supported by DFG grants Kr 2164/14-2 and Kr 2164/15-2.

## DATA AVAILABILITY STATEMENT

This work makes use of the following ALMA data: ADS/JAO.ALMA 2017.1.00167.S, 2015.1.01260.S, and 2015.1.00307.S. ALMA is a partnership of ESO (representing its member states), NSF (USA)
and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA), program 10539. This work has made use of data from the European Space Agency (ESA) mission Gaia,<fnref n="3" /> processed by the Gaia Data Processing and Analysis Consortium (DPAC<fnref n="4" />). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The Combined Atlas of Sources with Spitzer/IRS Spectra (CASSIS) is a product of the Infrared Science Center at Cornell University, supported by NASA and JPL. This work has made use of IRS data with a programme ID: 20463. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

## REFERENCES

Andrews S. M., Williams J. P., 2005, ApJ, 631, 1134

Ansdell M. et al., 2016, ApJ, 828, 46

Apai D., Schneider G., Grady C. A., Wyatt M. C., Lagrange A.-M., Kuchner M. J., Stark C. J., Lubow S. H., 2015, ApJ, 800, 136

Augereau J. C., Lagrange A. M., Mouillet D., Papaloizou J. C. B., Grorod P. A., 1999, A&A, 348, 557

Aumann H. H. et al., 1984, ApJ, 278, L23

Benz W., Asphaug E., 1999, Icarus, 142, 5

Booth M. et al., 2016, MNRAS, 460, L10

Booth M. et al., 2019, MNRAS, 482, 3443

Booth R. A., Clarke C. J., 2016, MNRAS, 458, 2676

Bryden G. et al., 2009, ApJ, 705, 1226

Chiang E., Kite E., Kalas P., Graham J. R., Clampin M., 2009, ApJ, 693, 734

Chiang E., Laughlin G., 2013, MNRAS, 431, 3444

Cieza L. A. et al., 2019, MNRAS, 482, 698

Daley C. et al., 2019, ApJ, 875, 87

Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531

Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, Astrophysics Source Code Library, record ascl:1202.015

ESA 1997, The HIPPARCOS and TYCHO Catalogues. Astrometric and Photometric Star Catalogues Derived from the ESA HIPPARCOS Space Astrometry Mission, ESA SP-1200. ESA, Noordwijk

Faramaz V. et al., 2019, AJ, 158, 162

Fernandez J. A., Ip W. H., 1984, Icarus, 58, 109

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

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

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

Gomes R. S., Morbidelli A., Levison H. F., 2004, Icarus, 170, 492

Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Nature, 435, 466

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Harper D. A., Loewenstein R. F., Davidson J. A., 1984, ApJ, 285, 808

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

Holland W. S. et al., 2017, MNRAS, 470, 3606

Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541

Ida S., Bryden G., Lin D. N. C., Tanaka H., 2000, ApJ, 534, 428

Ida S., Makino J., 1993, Icarus, 106, 210

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

Jackson A. P., Wyatt M. C., Bonsor A., Veras D., 2014, MNRAS, 440, 3757

Jurić M., Tremaine S., 2008, ApJ, 686, 603

Kalas P., Fitzgerald M. P., Graham J. R., 2007, ApJ, 661, L85

Kass R. E., Raftery A. E., 1995, J. Am. Stat. Assoc., 90, 773

Kennedy G. M., 2020, R. Soc. Open Sci., 7, 200063

Kennedy G. M., Marino S., Matrà L., Panić O., Wilner D., Wyatt M. C., Yelverton B., 2018, MNRAS, 475, 4924

Kennedy G. M., Wyatt M. C., 2010, MNRAS, 405, 1253

Kennedy G. M., Wyatt M. C., 2014, MNRAS, 444, 3164

Kenyon S. J., Bromley B. C., 2004, AJ, 127, 513

Kirsh D. R., Duncan M., Brasser R., Levison H. F., 2009, Icarus, 199, 197

Kral Q., Matrà L., Wyatt M. C., Kennedy G. M., 2017, MNRAS, 469, 521

Kral Q., Thébault P., Augereau J. C., Boccaletti A., Charnoz S., 2015, A&A, 573, A39

Krist J. E., 2004, Proc. SPIE, 5487, 1284

Krivov A. V., Ide A., Löhne T., Johansen A., Blum J., 2018, MNRAS, 474, 2564

Krivov A. V., Wyatt M. C., 2021, MNRAS, 500, 718

Lebouteiller V., Barry D. J., Spoon H. W. W., Bernard-Salas J., Sloan G. C., Houck J. R., Weedman D. W., 2011, ApJS, 196, 8

Lee E. J., Chiang E., 2016, ApJ, 827, 125

Liseau R. et al., 2008, A&A, 480, L47

Liseau R. et al., 2010, A&A, 518, L132

MacGregor A. M., Osten R. A., Hughes A. M., 2020, ApJ, 891, 80

MacGregor M. A. et al., 2016, ApJ, 823, 79

MacGregor M. A. et al., 2017, ApJ, 842, 8

MacGregor M. A. et al., 2019, ApJ, 877, L32

Marino S. et al., 2016, MNRAS, 460, 2933

Marino S. et al., 2017, MNRAS, 465, 2595

Marino S. et al., 2018, MNRAS, 479, 5423

Marino S., Yelverton B., Booth M., Faramaz V., Kennedy G. M., Matrà L., Wyatt M. C., 2019, MNRAS, 484, 1257

Marmier M. et al., 2013, A&A, 551, A90

Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J., Lafrenière D., Doyon R., 2008, Science, 322, 1348

Marois C., Zuckerman B., Konopacky Q. M., Macintosh B., Barman T., 2010, Nature, 468, 1080

Marshall J. P., Wang L., Kennedy G. M., Zeegers S. T., Scicluna P., 2021, MNRAS, 501, 6168

Matrà L. et al., 2017, MNRAS, 464, 1415

Matrà L. et al., 2020, ApJ, 898, 146

Matrà L., Marino S., Kennedy G. M., Wyatt M. C., Öberg K. I., Wilner D. J., 2018, ApJ, 859, 72

Matrà L., Panić O., Wyatt M. C., Dent W. R. F., 2015, MNRAS, 447, 3936

Matrà L., Wyatt M. C., Wilner D. J., Dent W. R. F., Marino S., Kennedy G. M., Milli J., 2019, AJ, 157, 135

Mayor M., Queloz D., 1995, Nature, 378, 355

Mermilliod J.-C., 1987, A&AS, 71, 413

Mermilliod J. C., 2006, VizieR Online Data Catalog, p. II/168

Meshkat T. et al., 2017, AJ, 154, 245

Moro-Martín A. et al., 2015, ApJ, 801, 143

Mustill A. J., Wyatt M. C., 2009, MNRAS, 399, 1403

Neugebauer G., et al., 1984, ApJL, 278, L1

Ormel C. W., Ida S., Tanaka H., 2012, ApJ, 758, 80

Paunzen E., 2015, A&A, 580, A23

Pawellek N., Krivov A. V., 2015, MNRAS, 454, 3207

Pearce T. D., Wyatt M. C., 2014, MNRAS, 443, 2541

<footnote n="3">³https://www.cosmos.esa.int/gaia</footnote>

<footnote n="4">⁴https://www.cosmos.esa.int/web/gaia/dpac/consortium</footnote>
Pearce T. D., Wyatt M. C., 2015, MNRAS, 453, 3329

Raymond S. N. et al., 2011, A&A, 530, A62

Raymond S. N. et al., 2012, A&A, 541, A11

Remijan A. et al., 2019, ALMA Doc. 7.3, ver. 1.0, ISBN 978-3-923524-66-2

Ricci L., Maddison S. T., Wilner D., MacGregor M. A., Ubach C., Carpenter
J. M., Testi L., 2015, ApJ, 813, 138

Schüppler C., Krivov A. V., Löhne T., Booth M., Kirchschlager F., Wolf S.,
2016, MNRAS, 461, 2146

Schwarz G., 1978, Anna. Stat., 6, 461

Scott K. S. et al., 2010, ApJS, 191, 212

Shannon A., Bonsor A., Kral Q., Matthews E., 2016, MNRAS, 462, L116

Sibthorpe B., Kennedy G. M., Wyatt M. C., Lestrade J.-F., Greaves J. S.,
Matthews B. C., Duchêne G., 2018, MNRAS, 475, 3046

Simpson J. M. et al., 2015, ApJ, 807, 128

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

Stapelfeldt K., Krist J., Bryden G., Chen C., 2007, In the Spirit of Bernard
Lyot: The Direct Detection of Planets and Circumstellar Disks in the 21st
Century. Univ. California, Berkeley, CA, p. 47

Stewart S. T., Leinhardt Z. M., 2009, ApJ, 691, L133

Su K. Y. L. et al., 2009, ApJ, 705, 314

Su K. Y. L. et al., 2017, AJ, 154, 225

Visser R., van Dishoeck E. F., Black J. H., 2009, A&A, 503, 323

White J. A. et al., 2020, ApJ, 894, 76

Wolszczan A., Frail D. A., 1992, Nature, 355, 145

Wright E. L. et al., 2010, AJ, 140, 1868

Wyatt M. C., 2003, ApJ, 598, 1321

Wyatt M. C., 2006, ApJ, 639, 1153

Wyatt M. C., 2008, ARA&A, 46, 339

Wyatt M. C., Dent W. R. F., 2002, MNRAS, 334, 589

Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., Holmes
E. K., Piña R. K., 1999, ApJ, 527, 918

Yelverton B., Kennedy G. M., Su K. Y. L., 2020, MNRAS, 495, 1943

Yelverton B., Kennedy G. M., Su K. Y. L., Wyatt M. C., 2019, MNRAS, 488, 3588

APPENDIX A: VISIBILITY DATA

For completeness, we include here the visibility data for all three ALMA epochs. Using the derived values for *i* and PA, Fig. A1 shows the interferometric visibilities as a function of the deprojected baselines. Although these profiles differ (their absolute values of their maxima and minima are inconsistent), they show consistently located nulls at ~10 kλ, ~35 kλ, and ~60 kλ, and consistently located peaks at ~20 kλ and ~45 kλ. This suggests that although these data sets may have different profiles, their radial emission has brightness peaks and minima at similar radial locations. The imaginary data (seen in the lower plots, which probe the azimuthal structure) show a number of significant features (particularly at the largest angular scales with uv-distances below 10 kλ for the 2018 Band 7 data). Although the majority of imaginary visibility data points are consistent with zero (i.e. axisymmetry), for the very shortest baselines (i.e. those associated with larger scale structure) there are departures from this, symptomatic of the type of departure from symmetry as measured in Section 3.1.1.

<figure><img src="image_15.png" /><figcaption>**Figure A1.** Composite plots of the binned visibility data for the three ALMA epochs. Error bars here are calculated from the standard deviation divided by the number of independent data points within each bin, where we have included 120 bins per measurement set.</figcaption></figure>
## APPENDIX B: DUST DENSITY AND TEMPERATURE DISTRIBUTION

Here we provide the dust and temperature distributions as calculated by RADMC – 3D for the best-fitting model CL in Fig. B1.

## APPENDIX C: POSTERIOR DISTRIBUTION OUTPUTS

Here we provide a visual of the MCMC posterior distribution corner plots for the model with a clump in the SW ansa (model CL), in two plots showing the disc parameters (in Fig. C1) and for the SW clump (in Fig. C2).

<figure><img src="image_16.png" /><figcaption>**Figure B1.** Dust temperature and dust density as a function of radial distance from the star, as computed by RADMC – 3D.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure C1.** MCMC posterior distribution of the 18 parameter 'Model CL': the symmetric belt with clump in the SW ansa. Posterior distributions show the range of values determined for each parameter (see each histogram), and how each parameter varies as a function of each of the other parameters. This allows us to interpret whether the fits are well constrained, and if there are any degeneracies between variables. Here, we show the 12 parameters that describe the underlying disc. It can be seen that the distribution of the parameter *p*<sub>inner</sub> is flat beyond a value of 18, meaning that this is poorly constrained.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure C2.** MCMC posterior distribution, as per Fig. C1. Here we show the six parameters that describe the SW clump.</figcaption></figure>

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