# Probing the eccentricity in protostellar discs: Modelling kinematics and morphologies

Enrico Ragusa, Elliot Lynch, Guillaume Laibe, Cristiano Longarini, Simone Ceppi

► To cite this version:

Enrico Ragusa, Elliot Lynch, Guillaume Laibe, Cristiano Longarini, Simone Ceppi. Probing the eccentricity in protostellar discs: Modelling kinematics and morphologies. Astronomy & Astrophysics - A&A, 2024, 686, 10.1051/0004-6361/202449583. insu-04838866

HAL Id: insu-04838866

https://insu.hal.science/insu-04838866v1

Submitted on 15 Dec 2024

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
# Probing the eccentricity in protostellar discs:
Modelling kinematics and morphologies

Enrico Ragusa<fnref n="1" />,<fnref n="2" />,<fnref n="3" />✉, Elliot Lynch<fnref n="2" />, Guillaume Laibe<fnref n="2" />, Cristiano Longarini<fnref n="3" />,<fnref n="4" />, and Simone Ceppi<fnref n="3" />✉

Received 13 February 2024 / Accepted 2 April 2024

## ABSTRACT

**Context.** Protostellar discs are mostly modelled as circular structures of gas and dust orbiting a protostar. However, a number of physical mechanisms, for example, the presence of a (sub)stellar companion or initial axial asymmetry, can cause the gas and dust orbital motion to become eccentric. Theoretical studies have revealed that, when present, disc eccentricity is expected to occur with predictable profiles that can be long-lasting and potentially observable in protostellar systems.

**Aims.** We construct an analytical model predicting the typical features of the kinematics and morphology of eccentric protostellar discs, with the final goal of characterising the observational appearance of eccentricity in discs.

**Methods.** We validate the model using a numerical simulation of a circumbinary disc (where the binary makes the disc eccentric). We finally post-process the simulation with Monte Carlo radiative transfer to study how eccentric features would appear through the 'eyes' of ALMA.

**Results.** Besides the motion of the material on eccentric Keplerian orbits in the disc orbital plane, the most characteristic eccentric feature emerging from the analytical model is strong vertical motion with a typical anti-symmetric pattern (with respect to the disc line of pericentres). A circumbinary disc with a ≈ 40 au eccentric cavity ($e_{cav} = 0.2$), carved by an $a_{bin} = 15$ au binary, placed at a distance $d = 130$ pc, is expected to host in its upper emission surface vertical oscillations up to $v_z \sim 400$ m s⁻<fnref n="1" /> close to the cavity edge, that is to say, well within ALMA spectral and spatial resolution capabilities. A residual spiral pattern in the vertical velocity $\Delta v_z \sim 150$ m s⁻<fnref n="1" /> of the simulation cannot be captured by the theoretical model, we speculate it to be possibly linked to the presence of a companion in the system.

**Key words.** protoplanetary disks – planet-disk interactions – binaries: general – stars: formation – stars: pre-main sequence

## 1. Introduction

In celestial mechanics, particles orbiting spherically symmetric objects follow Keplerian ellipses. This classical physics result, at the heart of solar and extrasolar planet mechanics for centuries, comes from Newton's $1/r$ potential. These orbits naturally extend to continuous distributions of matter orbiting point-like masses, known as accretion discs, where the horizontal fluid motion follows a set of nested confocal, Keplerian ellipses.

The mathematical complexity of describing the evolution of fluid flows in eccentric discs has meant that circular discs are almost always preferred as 'spherical cows' to derive models that can be compared with observations (e.g. Shakura & Sunyaev 1973; Pringle 1981), filtering out effects such as variable scale heights and non-axisymmetric orbital compression. In that sense, the widely used terminology 'Keplerian discs', abusively used for 'circular Keplerian discs', is misleading, since it propagates the idea that discs have no eccentricity by default.

The dynamics of eccentric discs has been the subject of a number of theoretical studies (e.g. Ogilvie 2001; Goodchild & Ogilvie 2006; Ogilvie & Barker 2014; Teyssandier & Ogilvie 2016; Ogilvie & Lynch 2019; Lee et al. 2019a,b), which have shown that the eccentricity propagates in the disc as a form of density waves. In these formalisms the (generally non-linear)

radial and azimuthal velocity perturbations from circular motion can be expressed in terms of orbital eccentricity and the orientation of the apsidal angle. Disc eccentricity can be excited by a companion in the disc (binary star, planet, e.g. Lubow 1991a; Kley & Dirksen 2006). It has also been shown to be damped by bulk viscosity, but it is generally excited by shear viscosity (e.g. Syer & Clarke 1992; Kley et al. 1993).

However, it is not clear that viscous models of turbulence are even applicable to eccentric discs, making the evolution of eccentric discs a promising approach to constraining the behaviour of real disc turbulence – see for example Ogilvie (2001, 2003), Ogilvie & Proctor (2003), and Lynch & Ogilvie (2021) for theoretical works, and Papaloizou (2005), Pierens et al. (2020), Dewberry et al. (2020), Oyang et al. (2021), and Chan et al. (2022) for numerical works concerning turbulence in eccentric discs.

Disc eccentricity also places a chronological constraint on scenarios of planet formation, with recent work showing that discs are likely born eccentric (Hennebelle et al. 2020; Lebreuilly et al. 2021; Lovascio et al. 2024), consistent with the young non-axisymmetric discs observed in the eDisk survey (Ohashi et al. 2023; Thieme et al. 2023; van't Hoff et al. 2023). This contrasts with the large fraction of more evolved discs that are observed to be nearly circular (Long et al. 2018;

<footnote n="1">¹ Dipartimento di Matematica, Università degli Studi di Milano, Via Saldini 50, 20133 Milano, Italy
 e-mail: enrico.ragusa@unimi.it</footnote>

<footnote n="2">² ENS de Lyon, CRAL UMR5574, Université Claude Bernard Lyon 1, CNRS, Lyon, 69007, France</footnote>

<footnote n="3">³ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, MI, Italy</footnote>

<footnote n="4">⁴ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK</footnote>
Andrews et al. 2018; Öberg et al. 2021), with some exceptions represented by some discs with central cavities (often referred to as transition discs) where kinematic or morphological measurements of protostellar disc eccentricity were possible (Dong et al. 2018; Kuo et al. 2022; Garg et al. 2021; Yang et al. 2023; Kozdon et al. 2023); other systems showed non-axisymmetric residuals when a circular Keplerian velocity map was subtracted from their kinematics (Wölfer et al. 2023), possibly hinting at their eccentric nature. The future availability of new high spatial and spectral resolution datasets (e.g. the coming exoALMA survey) creates the need for new detailed models of disc kinematics and also enables the possibility of testing them directly with observations.

Hence, in this study, we aim to characterise kinematic signatures of eccentricity in protostellar accretion discs. We show how the kinematics and morphology of eccentric discs is fully described by fixing three functions: the eccentricity profile $e(a)$, the pericentre longitude profile $\omega(a)$, the disc mass distribution $M_a(a)$, acting as a sort of density profile that is constant on orbits, and an assumption about the disc thermodynamics, here done prescribing a locally isothermal sound speed profile $\langle c_s \rangle(a)$. In the definition of the profiles we just discussed, $a$ is the disc semi-major axis coordinate from celestial mechanics.

We finally note, for completeness, that eccentric disc dynamics has been proven to be relevant in other astrophysical contexts than protostellar discs: for example, disc eccentricity has been directly measured and discussed in debris discs in the late stages of star formation (MacGregor et al. 2013; Olofsson et al. 2019; Faramaz et al. 2019; Booth et al. 2021; MacGregor et al. 2022; Lynch & Lovell 2022; Lovell & Lynch 2023); similarly, debris discs formed after the tidal disruption of asteroids around white dwarfs feature changes in the shape of spectral lines that have been interpreted as the precession of an eccentric disc due to pressure and/or to general relativistic effects (Wilson et al. 2015; Manser et al. 2016a,b; Cauley et al. 2018; Miranda & Rafikov 2018; Dennihy et al. 2018, 2020); and finally, in tidal disruption events of stars around black holes where the disc formed after the stellar disruption is expected to be eccentric (Liu et al. 2017; Cao et al. 2018; Hung et al. 2020; Wevers et al. 2022).

This paper is divided as follows: in Sect. 2, we provide a simplified analytical model that, taking as input $e(a)$, $\omega(a)$, $M_a(a)$, and $\langle c_s \rangle(a)$, predicts the following: the radial motion $v_R$, the orbital motion $v_\theta$, the vertical motion $v_z$, and disc morphology, in the form of $H/R$ and $\Sigma(a, \phi)$ – that is, predicting over-dense regions due to the eccentric nature of the disc. We compare the eccentric disc theoretical model developed to a numerical simulation of an eccentric circumbinary disc in Sects. 3 and 4. Finally, in Sect. 5, we examine the observability of these features in protostellar discs using ALMA. To do so, we rescaled the system and performed Monte Carlo radiative transfer simulations to quantify the amplitude of the velocity perturbations. We draw our conclusions in Sect. 6.

## 2. Analytic representation of the disc kinematics and structure

Eccentric discs can be described as a continuous set of nested, confocal, ellipses, orbiting the baricentre located at the common focus, with an eccentricity $e(a)$ and pericentre phase $\omega(a)$ profiles (see Fig. 1), where $a$ is the semi-major axis of the ellipse. Below we summarise some key aspects concerning the formalism used for describing their structure and dynamics. Details concerning eccentric disc evolution (i.e. how the eccentricity

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of the morphology of a disc composed by a set of nested ellipses with both eccentricity $e(a)$ and pericentre phase $\omega(a)$ varying for different values of the semi-major axis.</figcaption></figure>

profile evolves) are beyond the scope of the paper, however we refer the reader to Appendix A and references therein for a quick overview of the theory.

To do so, as is common in eccentric disc literature (e.g. Ogilvie 2001; Ogilvie & Barker 2014), we introduce in Sect. 2.1 a new coordinate system specialised to deal with eccentric orbits, then putting together celestial mechanics results, hydrodynamic considerations we provide a description of eccentric discs morphology and kinematics. The model presented summarises results from the existing literature, here self-consistently re-derived to unify the notation, for an $(a, \phi, z)$ coordinate system in Eulerian form in Sects. 2.1–2.4. We also introduce a new effective prescription for pressure corrections of the azimuthal velocity to simplify the equations in 2.5.

### 2.1. Eccentric discs coordinates

In this section we introduce a convenient coordinate system for describing the kinematics and structure of eccentric discs. We first note that using cylindrical coordinates $(R, \theta, z)$ when dealing with eccentric disc geometry is very inconvenient: a single value of $R$ spans through multiple elliptical orbits for different values of $\theta$. As mentioned above, since the orbital structure of an eccentric disc is defined by $e(a)$ and $\omega(a)$, a much better choice is an elliptical coordinate system using $(a, \phi)$, where $a$ is the semi-major axis of the ellipse, while the azimuthal coordinate $\phi$ is related to the true anomaly by $f = \phi - \omega(a)$. Since we assume that the focus of the ellipses is fixed at the origin, $\phi$ is coincident with the polar azimuthal coordinate $\theta$.

Assuming that the focus of all the confocal ellipses is at the origin of the cartesian coordinate system $(x, y, z)$, the transformation rules between $(a, \phi, z) \to (x, y, z)$ coordinates are

$$x(a, \phi) = R(a, \phi) \cos[\phi] \quad (1)$$

$$y(a, \phi) = R(a, \phi) \sin[\phi], \quad (2)$$

$$z = z, \quad (3)$$

where $R(a, \phi)$ is

$$R(a, \phi) = \frac{a[1 - e(a)^2]}{1 + e(a) \cos[\phi - \omega(a)]}. \quad (4)$$
The Jacobian matrix from this transformation is

$$J = \frac{\partial(x, y, z)}{\partial(a, \phi, z)} \qquad (5)$$

Its determinant, $J(a, \phi) = \det(J)$, will be useful for a number of purposes below, it has the dimension of a length and and it is given by

$$J(a, \phi) = R(a, \phi)^3 \frac{1 - e(a)[e(a) + 2ae_a(a)]}{a^2[1 - e(a)^2]^2} \cdot \{1 - q \cos[\phi - \varpi(a) - \alpha]\}, \qquad (6)$$

where the angle $\alpha$ and variable $q$ are introduced to simplify the form of the equations, by exploiting the relation for $\cos(a - b)$, and are related to $e$, $ae_a$, and $a\varpi_a$ by

$$q \cos(\alpha) = \frac{[1 + e(a)^2]ae_a(a) - [1 - e(a)^2]e(a)}{1 - e(a)[e(a) + 2ae_a(a)]}, \qquad (7)$$

$$q \sin(\alpha) = \frac{e(a)[1 - e(a)^2]a\varpi_a(a)}{1 - e(a)[e(a) + 2ae_a(a)]}. \qquad (8)$$

The quantities $e_a(a)$ and $\varpi_a(a)$ indicate the derivative with respect to $a$ of $e(a)$ and $\varpi(a)$, and are often referred to as ‘eccentricity gradient’ and ‘twist’, respectively. Equations (7) and (8) also imply that the angle $\alpha$ in Eq. (6) satisfies

$$\tan(\alpha) = \frac{ae(a)[1 - e(a)^2]\varpi_a(a)}{[1 + e(a)^2]ae_a(a) - [1 - e(a)^2]e(a)}. \qquad (9)$$

We note that since the coordinate $z$ does not depend on $a$ and $\phi$, the expression of $J$ is the same both when we deal with a 2D disc using the $(a, \phi)$ coordinates only, and for a full 3D disc $(a, \phi, z)$.

In order to provide a physical interpretation to $J$, we note that it constitutes the multiplicative factor to preserve element volume $dV = dxdydz$ or area $dA = dxdy$ while transforming to eccentric coordinates as follows

$$dV = Jd\alpha d\phi dz, \qquad (10)$$

$$dA = Jd\alpha d\phi. \qquad (11)$$

In the light of this, we note that a vanishing value of $J$ indicates a coordinate singularity (one value of $a$ characterises multiple ellipses at that location), the nested ellipses overlap causing an orbit intersection (e.g. Statler 2001). From Eq. (6), it is easy to verify that this instance occurs when $q \ge 1$. Similarly $\alpha$ represents the true anomaly, that is the angle from the pericentre longitude $\varpi$, where the maximum orbital compression occurs.

## 2.2. Orbital velocity field: Keplerian motion

In this section we define the velocities in the ‘eccentric’ $(a, \phi)$ coordinate system.

From celestial mechanics we know that the radial and azimuthal orbital velocities, in cylindrical coordinates $(R, \theta)$, of a test mass in Keplerian motion around a central body are (Murray & Dermott 1999)

$$v_R(a, \phi) = \dot{R} = a\Omega_0 \frac{e(a)}{\sqrt{1 - e^2(a)}} \sin[\phi - \varpi(a)], \qquad (12)$$

$$v_{\theta,K}(a, \phi) = R\dot{\theta} = a\Omega_0 \frac{1 + e(a) \cos[\phi - \varpi(a)]}{\sqrt{1 - e^2(a)}}, \qquad (13)$$

where the dot notation indicates total time derivative, $d/dt$, and $\Omega_0$ is the mean motion

$$\Omega_0 = \sqrt{\frac{GM}{a^3}}. \qquad (14)$$

From Eq. (13) we define the angular orbital velocity along an eccentric orbit as

$$\Omega(a, \phi) \equiv \dot{\theta} = \Omega_0 \frac{\{1 + e(a) \cos[\phi - \varpi(a)]\}^2}{[1 - e^2(a)]^{3/2}}. \qquad (15)$$

The standard analytical treatment of differential operators in the eccentric coordinate system requires the definition of contravariant velocities as $v^a \equiv \dot{a}$ and $v^\phi \equiv \dot{\phi}$. With these definitions in mind we can define $v^a$ and $v^\phi$ consistently with the orbital velocities in Eqs. (12) and (13). Along an eccentric orbit the semi-major axis does not change, so that

$$v^a = 0, \qquad (16)$$

while, given the coincidence between $\phi$ and $\theta$, $v^\phi$ is

$$v^\phi = \Omega. \qquad (17)$$

We remark that $v^\phi$ is a contravariant variable, with the dimension of an angular velocity. While $v_\theta$ is the azimuthal velocity, with the dimension of a standard spatial velocity, not the covariant velocity component that might be expected from the notation. In this context, $v_\theta$ and $v^\phi$ have different dimensions and should not be confused.

## 2.3. Vertical motion

During an eccentric orbit the material undergoes vertical oscillations driven by the periodic variation of the vertical gravitational potential along the orbit (the distance from the central source of gravity varies through the orbit Ogilvie & Barker 2014). Since we use a new coordinate system, we write the continuity equation in contravariant form,

$$\frac{\partial \rho}{\partial t} = -\frac{1}{J} \frac{\partial}{\partial a} (J v^a \rho) - \frac{1}{J} \frac{\partial}{\partial \phi} (J v^\phi \rho) - \frac{\partial}{\partial z} (v^z \rho), \qquad (18)$$

obtained from the expression of the divergence of a vector field $u$ in a generic coordinate system $\{x^i\}$,

$$\nabla_i (u^i) = \frac{1}{J} \frac{\partial}{\partial x_i} J u^i. \qquad (19)$$

For our purpose we also need the equation for momentum conservation in the vertical direction,

$$\frac{\partial v^z}{\partial t} + v^a \frac{\partial v^z}{\partial a} + v^\phi \frac{\partial v^z}{\partial \phi} + v^z \frac{\partial v^z}{\partial z} = - \frac{GM}{R^2} \frac{z}{R} - \frac{1}{\rho} \frac{\partial p}{\partial z}. \qquad (20)$$

We now look for steady solutions and substitute $v^a = 0$ and $v^\phi = \Omega(a, \phi)$ – as we prescribed in Sect. 2.2. After these operations, Eqs. (18) and (20) reduce to

$$\frac{\partial}{\partial z}(v^z \rho) = -\frac{1}{J} \frac{\partial}{\partial \phi}(J \Omega \rho), \qquad (21)$$

$$\Omega \frac{\partial v^z}{\partial \phi} + v^z \frac{\partial v^z}{\partial z} = - \frac{GM}{R^2} \frac{z}{R} - \frac{1}{\rho} \frac{\partial p}{\partial z}. \qquad (22)$$
We assume the disc to be locally isothermal so that $\rho$ and $p$ have Gaussian form (e.g. Lodato 2008 for a review of classical disc dynamics)

$$ \rho = \frac{\rho_0}{\sqrt{2\pi}} \exp\left(-\frac{z^2}{2H^2}\right), \qquad (23) $$

$$ p = \frac{p_0}{\rho_0} \rho, \qquad (24) $$

where $p_0$ and $\rho_0$ represent the density value in the midplane. We note that $p_0/\rho_0 = c_s^2$ is the squared local isothermal sound speed at the disc midplane; $H(a, \phi)$, is the standard deviation of the vertical density distribution, here representing the disc local vertical scale-height<fnref n="1" />. The form of the final equations does not change for different assumptions on $\rho$ and $p$, with some caveats discussed in Appendix C.

In order to study the evolution of $H$ along the orbit, we need a reasonable ansatz for $v^z$. We impose the quantity $z/H$ to be a Lagrangian variable – that is, $z/H$ does not change along the orbit – that follows the disc expansion

$$ \frac{d}{dt} \frac{z}{H} = \frac{\partial}{\partial t} \frac{z}{H} + \Omega \frac{\partial}{\partial \phi} \frac{z}{H} + v^z \frac{\partial}{\partial z} \frac{z}{H} = 0. \qquad (25) $$

This ansatz is physically the homogeneous vertical expansion of the gas column. When the disc is stationary ($\partial/\partial t = 0$) such a condition is satisfied if $v^z$ has the form<fnref n="2" />

$$ v^z = \Omega \frac{\partial H}{\partial \phi} \frac{z}{H}. \qquad (27) $$

We use Eq. (27) as an ansatz for the solution of $v^z$. By substituting this into Eqs. (21) and (22), they can be re-written as

$$ \Omega \frac{\partial \rho_0}{\partial \phi} = -\rho_0 \left[ \frac{1}{J} \frac{\partial}{\partial \phi} (J\Omega) + \frac{\Omega}{H} \frac{\partial H}{\partial \phi} \right], \qquad (28) $$

$$ \Omega^2 \frac{\partial^2 H}{\partial \phi^2} = -\Omega \frac{\partial H}{\partial \phi} \frac{\partial \Omega}{\partial \phi} - \frac{GM}{R^3} H + \frac{c_s^2}{H}. \qquad (29) $$

A similar set of equations can be obtained from the Lagrangian viewpoint (e.g. Ogilvie & Barker 2014; Lynch & Ogilvie 2021).

We first solve Eq. (29), as it does not depend on Eq. (28). In particular, Fig. 2 shows the solutions to Eq. (29) for the azimuthal dependence of $H$ for different values of eccentricity and fixed $\sigma = 0$, and the resulting $v^z$ from Eq. (27). The value of $\rho_0$ solving Eq. (28) can then be easily obtained by noting that Eq. (28) can be rewritten as

$$ \frac{1}{J} \frac{\partial}{\partial \phi} (J \Omega H \rho_0) = 0, \qquad (30) $$

implying that $J\Omega H\rho_0 \equiv \mathcal{F}(a)$ depends on the semi-major axis only and is a conserved mass flux along the orbit. Further considerations on $\mathcal{F}(a)$ will be used in the next section to discuss the disc morphology.

$$ v^z \equiv c_s \left[ \omega(\phi) \frac{\partial h}{\partial \phi} \right] \frac{z}{H}. \qquad (26) $$

where $c_s = H_0\Omega_0$, highlighting that $c_s$ sets the natural scale for the velocity $v^z$.

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

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Vertical height $H$ and vertical velocity $v_z$ as a function of $\phi$. Top panel: $H/H_0$ as a function of $\phi$, where $H_0 = c_s/\Omega_0$, obtained solving Eq. (29). Bottom panel: $v_z/c_s$ as a function of $\phi$, obtained from Eq. (27) after solving Eq. (29), calculated at an elevation $z = H$ above the midplane – $c_s$ is the most natural normalisation<fnref n="2" /> for $v_z$. Solution to Eq. (29) was obtained solving for periodic solutions across one orbit, i.e. $H(a, 0) = H(a, 2\pi)$, using a shooting method for solving the boundary value problem.</figcaption></figure>

To conclude this section, we comment on the validity of the assumption that the disc is locally isothermal. This assumption implies that the gas temperature does not experience changes in its temperature due to compression or expansion, even though for large eccentricities the disc might undergo a very strong vertical expansion (compression) at the apocentre (pericentre). Relaxing the assumption of the disc being locally isothermal implies treating also the internal energy or assuming an adiabatic equation of state. This would cause pressure forces to increase more steeply during the compression at the pericentre, reducing the ratio between maximum and minimum disc scale-height at apocentre and pericentre $H_{\text{apo}}/H_{\text{peri}}$. A full derivation with an adiabatic equation of state is beyond the scope of the paper, but can be found in Ogilvie & Barker (2014) or Lynch & Ogilvie (2021).

Beyond values of eccentricity of $e \approx 0.5-0.6$ the value of the $c_s/H$ term in Eq. (29) can become very large close to pericentre, causing integration issues.

<footnote n="1">¹ In the general framework (e.g. when the disc is not locally isothermal), $H$ is defined by the second vertical moment of the density distribution as $H^2 = \int_z^2 \rho dz / \int_z^2 \rho dz$.</footnote>

<footnote n="2">² By making explicit the azimuthal dependence of $\Omega = \Omega_0\omega(\phi)$ and $H = H_0h(\phi)$ it is possible to rewrite</footnote>
## 2.4. Morphology

Eccentric discs might present steady over-densities, caused by orbital compression (i.e. $e_a \neq 0$ and $\bar{\omega}_a \neq 0$), that result in a characteristic non-axisymmetric morphology. In this section we put together the previous geometrical considerations and conservation of flux along the orbit to derive the expression of $\Sigma(a, \phi)$.

We obtained in Eq. (30) conservation of mass flux along elliptical orbits. We note that $\Sigma = \int_z \rho dz = H\rho_0$, thus we can rewrite $\mathcal{F}(a)$ as<fnref n="3" />

$$ \mathcal{F}(a) = J(a, \phi)\Omega(a, \phi)\Sigma(a, \phi). \quad (31) $$

We note that integrating the r.h.s. of Eq. (31) over an orbital time $t_{\text{orb}}$ implies

$$ \int_{0}^{t_{\text{orb}}} J(a, \phi)\Sigma(a, \phi)\Omega(a, \phi)dt = \int_{0}^{2\pi} J(a, \phi)\Sigma(a, \phi)d\phi = M_{a}(a) \quad (32) $$

where $M_a(a)$ is the derivative of the total disc mass enclosed within the semi-major axis $a$, $M(a)$, with respect to the semi-major axis – that is, $M_a da$ is the mass of an infinitesimal eccentric annulus of width $da$.

Since $\mathcal{F}(a)$ does not depend on $\phi$, we have $t_{\text{orb}}\mathcal{F}(a) = M_a(a)$, which can be rewritten as $\mathcal{F}(a) = M_a(a)\Omega_0/2\pi$, where we recall $\Omega_0$ is the mean motion (Eq. (14)). In the light of these considerations, Eq. (31) can be rewritten as (Ogilvie 2001; Statler 2001)

$$ \Sigma(a, \phi) = \frac{M_a(a)\Omega_0}{2\pi J(a, \phi)\Omega(a, \phi)}. \quad (33) $$

If $e_a(a) = 0$ and $\bar{\omega}_a(a) = 0$ – that is, the disc has vanishing eccentricity gradient and twist – the expression of $\Sigma(a, \phi)$ simplifies to an expression that does not depend on $\phi$.

This implies that an eccentric disc can, in principle, show no morphological features at all for arbitrary large eccentricities, provided its eccentricity and pericentre phase profiles are constant. Conversely, a negative (positive) eccentricity gradient $e_a(a) < 0$ ($e_a(a) > 0$) will produce an azimuthal overdensity located at the apocentre (pericentre); if the disc also presents some orbit twisting ($\bar{\omega}_a(a) \neq 0$), the overdensity is located at angle $\phi$ satisfying $\phi - \bar{\omega}(a) = \alpha$, where $\alpha$ is presented in Eq. (9).

It is important to note that any $e(a) = \text{const}$, $\bar{\omega}(a) = \text{const}$ profile is not steady. Such a disc will evolve in time producing a strong eccentricity gradient and twist, due to the transport throughout the disc of angular momentum and energy in the form of eccentric waves. If some dissipation is present, the disc will finally settle in a steady eccentricity configuration. More considerations concerning the evolution of eccentricity profiles can be found in Appendix A.1.

## 2.5. Approximated pressure corrections

In a pressure supported disc the centrifugal balance is reached with a sub/super-Keplerian orbital velocity. In analogy with classic disc accretion theory, we apply a generalisation of the standard pressure correction $\delta v_{\theta,P}$ required for a circular accretion

disc, that is, (e.g. Lodato 2008 for a review of classical disc dynamics):

$$ v_{\theta} = v_{\theta,K} \sqrt{1 + \eta}, \quad (34) $$

$$ \eta = (\Omega_0 a)^{-2} \left( \frac{a}{\bar{\Omega}_0} \frac{dP}{da} \right) \quad (35) $$

$$ P = \langle c_s \rangle^2 \Sigma_0, \quad (36) $$

where $v_{\theta,K}$ is the Keplerian velocity (Eq. (13)), $\Sigma_0$ constitutes an averaged value along the orbit defined as

$$ \Sigma_0 = \frac{M_a}{2\pi a}, \quad (37) $$

and $\langle c_s \rangle$ is an averaged value of $c_s$ along an orbit with semi-major axis $a$.

Equation (34) constitutes a simplification with respect to solving the fluid dynamics equations accounting for the pressure term. The analysis in Appendix A of Ogilvie & Lynch (2019) allows for a self-consistent derivation of the pressure correction within the Hamiltonian hydrodynamics framework, however for the clarity of presentation we prefer to stick to the simplified approach provided in this section, that provides a good match with numerical simulations anyway, as will be discussed in Sect. 3.

## 2.6. Disc eccentricity evolution: growth and damping

We summarise here for completeness the main physical mechanisms responsible for disc eccentricity growth and damping, despite remarking that the model we developed is independent of how the system arrived in its eccentricity configuration. More considerations concerning the existence and evolution of secularly evolving eccentric eigenmodes is discussed in Appendix A.1.

Concerning growth, two main mechanisms can be found: eccentricity pumping by a companion mass and primordial growth of the disc eccentricity.

In the first, a second mass in the system<fnref n="4" />, above a certain mass-ratio threshold ($q \ge 10^{-3}$, Kley & Dirksen 2006, even though larger values of $q$ appear to be required in 3D, Li et al. 2023) injects angular momentum and energy at resonant locations that cause the disc eccentricity to grow (Lubow 1991a; Goldreich & Sari 2003). This has been confirmed by numerical studies explicitly investigating the evolution of the disc eccentricity assuming both low mass-ratio companions (Papaloizou et al. 2001; Kley & Dirksen 2006; D’Angelo et al. 2006; Dunhill et al. 2013; Ragusa et al. 2018; Teyssandier & Lai 2019; Dempsey et al. 2021; Tanaka et al. 2022) and high mass-ratio companions, both in circumbinary discs (MacFadyen & Milosavljević 2008; Marzari et al. 2009; Shi et al. 2012; Dunhill et al. 2015; Miranda et al. 2017; Ragusa et al. 2020; Muñoz & Lithwick 2020; Pierens et al. 2020; Dittmann & Ryan 2022; Siwek et al. 2023) and in the discs surrounding the individual masses (Lubow 1991a,b; Whitehurst 1994; Murray 1996; Kley et al. 2008; Regály et al. 2011). Numerical simulations of flybys – that is, the two stars are not gravitationally bound – in star formation regions have also measured a substantial growth of the disc eccentricity after the close encounter between the two stars (Cuello et al. 2019). Furthermore numerous numerical works on circumbinary discs,

<footnote n="3">³ Equation (31) can be directly obtained in this form by vertically integrating the continuity equation Eq. (18) using $\Sigma = \int_z \rho dz$ and looking for stationary solutions.</footnote>

<footnote n="4">⁴ Binary stars or planets for protostellar systems – but also moons, for planets with rings like Saturn, or compact objects, for X-ray binaries, cataclysmic variables and stellar/supermassive binary black hole.</footnote>
despite not measuring directly the disc eccentricity, highlighted the formation of visibly eccentric cavities (e.g. D'Orazio et al. 2013; Farris et al. 2014; Ragusa et al. 2016, 2017; Price et al. 2018b; Calcino et al. 2019, 2020; Heath & Nixon 2020; Tiede et al. 2020; Franchini et al. 2023) or eccentric gaps for planets and low mass-ratio companions (e.g. Ataiee et al. 2013; Zhu & Stone 2014; Scardoni et al. 2022).

A second relevant scenario for eccentric discs concerns the primordial deviations from circular motion of the material in the discs that were excited and impressed during its formation – that is, the disc was born eccentric.

In the process of star formation, it is reasonable to expect that the phase of cloud collapse is not spherically symmetric, seeding an initial eccentricity in the protostellar disc. Numerical simulations studying this phase appear to produce, as a natural outcome, eccentric discs, even in presence of symmetric initial conditions (Hennebelle et al. 2020; Lebreuilly et al. 2021; Lovascio et al. 2024). Similarly, tidal disruption events (TDE) of eccentric (or parabolic) bodies around compact objects are expected to form eccentric discs surrounding the central body – be it a white dwarf (Trevascus et al. 2021), neutron star (Kurban et al. 2023) or a black hole (Shiokawa et al. 2015; Bonnerot et al. 2016; Zanazzi & Ogilvie 2020; Lynch & Ogilvie 2021; Cufari et al. 2022). In this scenario, the debris are expected to inherit the eccentricity of the disrupted body – in the event of parabolic orbits as in the context of TDEs of stars on supermassive black holes, half of the disrupted body remains bound to the central object with an eccentricity $e \le 1$.

Concerning disc eccentricity damping, we can identify a resonant and a viscous/dissipative mechanisms. On the one hand, while some resonances can pump the disc eccentricity, others can damp it (Goldreich & Sari 2003). Both eccentricity damping and pumping resonances can be present in a disc. Whether the disc eccentricity is growing or decreasing depends on the balance of the contributions of individual resonances, akin to what determines the direction of satellite migration: when some resonances push the satellite inwards and others outwards.

On the other hand, viscous/dissipative processes decrease the disc eccentricity. Viscosity is often included in theoretical models to describe turbulence, that, among its effects, causes the material to accrete on to the central mass. Pure shear, as for the Shakura & Sunyaev (1973) α-prescription, is not expected to produce eccentricity damping. In contrast, in some instances it can even trigger a viscous-overstability (e.g. Syer & Clarke 1992; Kley et al. 1993; Lyubarskij et al. 1994; Latter & Ogilvie 2006) that results in a growth of the disc eccentricity. Bulk viscosity instead always produce eccentricity damping (Ogilvie 2001; Goodchild & Ogilvie 2006; Lynch 2022).

Numerical dissipation, intrinsically present in any numerical scheme, produces a spurious form of bulk viscosity that causes eccentricity damping. In general, this implies that any numerical simulation is affected by a spurious decrease of the disc eccentricity due to unavoidable numerical dissipation. This produces a damping rate that is most likely to be higher than what would be reasonable to expect if the fluid was inviscid and turbulent. To date, no reliable estimate about eccentricity damping rates is available. A study about the interplay between disc eccentricity and turbulence has been attempted by Wienkers & Ogilvie (2018).

## 2.7. How the disc eccentric nature affects the dynamics

In the previous sections we provided equations that describe the kinematics and structure of eccentric discs. As obvious from

celestial mechanics considerations, compared with the circular case where $v_{\theta,circ} = \text{const}$ and $v_R = 0$ along an orbit, an eccentric disc has $v_{\theta}$ and $v_R$ varying around the orbit. The azimuthal velocity $v_{\theta}$ (Eqs. (13), (34) for the pressure corrected) reaches its maximum at the phase of pericentre $\varpi$ and minimum at the phase of the orbit apocentre. Similarly, $v_R$ (Eq. (12)) oscillates around $v_R = 0$, reaching its absolute maxima ($|v_R|$) at true anomaly $f = \pi/2$ and $f = 3\pi/2$.

However, eccentric discs are also characterised by motion out of the disc plane, in the vertical direction. This motion is the result of the change of the vertical gravitational field along the orbit that produces a change in the vertical scale height $H$ of the disc (Eq. (29)). This effect causes $H$ to reach its minimum at pericentre and maximum at apocentre. The ratio between $H_{apo}$ at apocentre and $H_{peri}$ at pericentre depends on the orbital eccentricity and can easily reach $H_{apo}/H_{peri} \ge 6$ for $e \ge 0.4$ (see Fig. 2). It is useful to note that the ratio $z/H$ remains constant along the orbit (Eq. (25)), implying that the material moves along the orbit following eccentric streamlines with $z \propto H$.

This change of altitude along the orbit implies that the material develops a vertical velocity $v_z$ (Eq. (27)) that oscillates between $v_z = 0$ at the pericentre and apocentre and reaches its maximum absolute value $|v_z|$ at $f \sim \pi/2$ and $f \sim 3\pi/2$.

The amplitude of these vertical $v_z$ oscillations scales with $z/H$ (Eq. (27)). This implies that the oscillations increase in amplitude the higher a fluid element is in the disc atmosphere, and are anti-symmetric with respect to the midplane – that is, top and bottom layers of the disc have opposite sign of $v_z$. In particular, at $z = H$, the amplitude of the vertical oscillations has as natural scale factor the sound speed $c_s$ (see footnote 2). Given the linear scaling with $z$, $v_z$ can easily become supersonic in higher disc layers with respect to the midplane ($z > H$). This type of oscillations are often referred to as a ‘breathing mode’.

Changes in the disc geometry due to eccentricity and pericentre phase variations across the disc result in overdensities and depletions of the disc surface density, $\Sigma$ (Eq. (33)). However, one should always keep in mind that a change in the disc thickness might affect the volume density, $\rho$, at the midplane without affecting the surface density: for example, a disc with $e = \text{const}$ and $\varpi = \text{const}$ across its entire domain, produces no variations in the surface density morphology, that is, $\Sigma(\theta) = \text{const}$. But, since the volume density at the midplane is $\rho = \Sigma/H$, $\rho$ changes as the inverse of $H$ – that is, $\rho$ is maximum at pericentre and minimum at apocentre even when $\Sigma$ is constant along the orbit.

## 3. Numerical simulation

In this section, we present a 3D hydrodynamical simulation of a circumbinary disc that throughout its evolution becomes eccentric (as expected due to binary-disc interaction) and perform a comparison with the analytical model predictions. We note that for the remainder of the paper we do not use contravariant notation – from now on we refer to cylindrical (non-covariant) velocities as $v_R$, $v_{\theta}$, $v_z$ and to cartesian velocities as $v_x$, $v_y$, $v_z$.

We approached the problem as follows (Sect. 3.3 for more details): we first extracted from the simulations the profiles of $e(a)$, $\varpi(a)$ and $M(a)$, on which the analytical model depends. From these profiles, using the formalism discussed in Sect. 2, we generated a 3D analytical model attributing to a grid of pixels (x, y, z), representing the coordinates of the disc surface, the corresponding ($v_x$, $v_y$, $v_z$). We finally compared the predictions of the theoretical model with the simulations.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Rendering of the surface density of the dump used as a reference numerical simulation for testing the analytical models. Dump from simulation 4A (mass-ratio $M_2/M_1 = 0.1$, $H/R = 0.05$) in Ragusa et al. (2020) after $t = 500t_{\text{bin}}$. The dump features a cavity with an eccentricity of $e_{\text{cav}} \approx 0.2$ and a cavity with a semi-major axis $a_{\text{cav}} \approx 2.7a_{\text{bin}}$.</figcaption></figure>

We note that quantities generated using the theoretical model are indicated with a subscript ‘th’ (e.g. $x_{\text{th}}$) while those calculated directly from the simulation with subscript ‘sim’ (e.g. $x_{\text{sim}}$).

### 3.1. SPH numerical setup

As a benchmark simulation to test the analytical model presented in Sect. 2 we use one of the simulations from the numerical set presented in Ragusa et al. (2020), that is, a 3D hydrodynamical disc surrounding a binary (here used as a source for the disc eccentricity).

Simulations in Ragusa et al. (2020) were meant to study the evolution of the disc eccentricity exploring various choices of binary mass-ratios and disc parameters. All simulations with $M_2/M_1 > 0.05$ presented in Ragusa et al. (2020) show a rapid growth of the disc eccentricity after $t \approx 400–500t_{\text{bin}}$ due to binary-disc interaction, and constitute the perfect testbed for the analytical model developed in this paper.

In general, we note that the dynamics in the simulation is a priori richer than the one described by the theoretical model: viscosity, waves, quadrupolar potential, accretion and disc spreading are not at all considered in the derivation we provided above.

From that set, we select simulation 4A (mass-ratio $M_2/M_1 = 0.1$). In particular, its $t = 500t_{\text{bin}}$ dump (shown in Fig. 3), where $t_{\text{bin}}$ is the binary orbital timescale represents a good candidate for our direct comparison: i) because at that time the disc eccentricity has reached a saturation value of $e_{\text{cav}} \approx 0.2$ at the cavity edge, which is qualitatively consistent with the few observationally measured ones in transition discs (e.g. Garg et al. 2021; Yang et al. 2023); ii) because the simulations with binary mass-ratios $M_2/M_1 > 0.1$ feature prominent density spirals and an orbiting overdense lump of material that alter the kinematics. The

motion of the material in those simulations is in fact eccentric, but it is also characterised by additional perturbations that are not captured by the model in Sect. 2.

Simulation 4A of Ragusa et al. (2020) was performed using the code PHANTOM (SPH, Price et al. 2018a), it used $N_p = 10^6$ particles, and consisted of a circular, live binary with mass-ratio $M_2/M_1 = 0.1$ surrounded by a circumbinary disc with initial mass $M_d = 5 \cdot 10^{-3} M_{\text{bin}}$ ($M_d = 4.8 \cdot 10^{-3} M_{\text{bin}}$ after $t = 500t_{\text{bin}}$), where $M_{\text{bin}} = M_1 + M_2$. The equation of state was chosen to be locally isothermal – with radial temperature profile producing $c_s \propto R^{-0.25}$ and $H/R = 0.05$ at $R_{\text{in}} = 2a_{\text{bin}}$. A circular cavity with $R = R_{\text{in}}$ was initially excised at $t = 0$. The mechanisms responsible for angular momentum transport were modelled using SPH artificial viscosity, producing an effective Shakura & Sunyaev (1973) $\alpha_{\text{ss}} \approx 5 \times 10^{-3}$ (obtained using $\alpha_{\text{AV}} = 0.2$, $\beta = 2$ Price et al. 2018a)$^5$ at the beginning of the simulation.

### 3.2. Obtaining maps of $v_R$, $v_{\theta}$, $v_z$, $H$, $\Sigma$ from the simulation

For the comparison with the analytical model, we start generating projected $x-y$ maps of $v_{x,\text{sim}}$, $v_{y,\text{sim}}$, $v_{z,\text{sim}}$, $H_{\text{sim}}$, and $\Sigma_{\text{sim}}$ using vertical particle integration and density weighted integration – that is, an SPH equivalent of a density weighted average along the vertical direction. The calculation of these quantities is performed using SPLASH as detailed below, the implementation of the integration procedures is described in Price (2007).

Velocities in the disc orbital plane are defined as $v_{x,\text{sim}} = \langle v_x \rangle_z$ and $v_{y,\text{sim}} = \langle v_y \rangle_z$, where $\langle \cdot \rangle_z$ indicates density weighted integration along the vertical direction and $v_x$ and $v_y$ the velocities of individual SPH particles in the density weighted integration. Then, $v_{R,\text{sim}}$ and $v_{\theta,\text{sim}}$ are obtained projecting $v_{2D} = \{v_{x,\text{sim}}, v_{y,\text{sim}}\}$ on radial and azimuthal unit vectors, respectively.

The map for $H_{\text{sim}}$ is obtained through a density weighted integration along the vertical direction as $H_{\text{sim}} = \langle z^2 \rangle_z^{1/2}$ – that is, based on the definition of scale-height $H_{\text{sim}}$ as the second moment of the particles vertical density distribution. Again, here $z^2$ indicates that the value of $z^2$ of individual SPH particles has been used in the density weighted integration. Since the disc is locally isothermal with a sound speed $c_s(R)$, the vertical density profile has the form of a Gaussian with standard-deviation $\sigma = H_{\text{sim}}$.

Similarly, the map for $v_{z,\text{sim}}$ is obtained performing a density weighted integration of the quantity $v_{z,\text{sim}} = \sqrt{\pi/2}\langle v_z \text{sign}(z) \rangle_z$ along the vertical direction. This choice regarding the functional form and normalising pre-factor of the quantity to be vertically averaged has the following motivations: i) concerning the functional form, since $v_z$ is anti-symmetric with respect to the midplane, the quantity $v_z \text{sign}(z)$ has the nice property of providing a $v_{z,\text{sim}}$ maintaining the sign of the velocities of the $z > 0$ half of the disc during the average across the midplane. In contrast, $\langle |v_z| \rangle_z$ or $\langle v_z^2 \rangle_z$ would produce a map that is everywhere positive, while the vertical average of $\langle v_z \rangle_z$, would produce $v_{z,\text{sim}} = 0$ across the whole disc. ii) The pre-factor $\sqrt{\pi/2}$ is introduced because, since the vertical density distribution is Gaussian, $\langle z \text{sign}(z) \rangle_z = H \sqrt{2/\pi}$. Under the assumption that $v_z \propto z/H$, as derived in Eq. (27), the introduction of the pre-factor makes sure that $v_{z,\text{sim}}$ traces the vertical velocity of particles at $z \sim H_{\text{sim}}$ above the midplane. We anticipate that the results presented in the following sections confirm the validity of this assumption.

<footnote>$^5$ Such an implementation implies an unavoidable amount of bulk viscosity (Lodato & Price 2010), in contrast with the pure shear nature of the Shakura & Sunyaev (1973) prescription.</footnote>
Finally, the map of Σ is obtained using the standard vertical integration of the particles' masses performed by SPLASH.

### 3.3. Constructing a 3D analytical model from the hydrodynamical simulation

The entire set of equations constituting our analytical model to predict the kinematics and the morphology of eccentric discs rely exclusively on the profiles $e(a)$, $\bar{\omega}(a)$, their derivatives, $e_a$ and $\bar{\omega}_a$, $M_a$ and an assumption on the sound speed $c_s$, that here we choose to have a radial dependence $c_s(R)$ for consistency with what is prescribed in the simulation. In this section, starting from the aforementioned profiles, we generate a 3D analytical model attributing a velocity vector ($v_x, v_y, v_z$) to each ($x, y, z$), where $z$ represents the reference altitude of the disc surface – this will be further clarified below.

From the reference snapshot discussed in Sect. 3.1 we extract the profiles of $e(a)$, $\bar{\omega}(a)$, and $M_a(a)$ (details about how semi-major axis profiles are calculated can be found in Teyssandier & Ogilvie 2017 and Ragusa et al. 2020), we filter each profile using the Savitzky–Golay filter (Savitzky & Golay 1964) to reduce the noise of the datasets and we interpolate them using cubic spline interpolation so that we can compute the derivatives $e_a(a)$, $\bar{\omega}_a(a)$ with respect to the semi-major axis. We set the filter window to be $w_{sg} = 50$ profile points (corresponding to a $\Delta a \sim 1.6a_{\text{bin}}$) and polynomial order $p_{sg} = 2$. We note that filtering is a very important step, since noisy datasets of $e(a)$, $\bar{\omega}(a)$ would produce unusable outputs when taking derivatives with respect to $a$.

In Fig. 4, we show the profiles extracted from our reference snapshot and the interpolated profiles we use for generating the analytical model. No appreciable differences can be found changing the parameter choice of the filter: as shown in Fig. 4, the filtering process has the only effect of reducing the noise of the dataset but it does not alter appreciably its qualitative form.

We use Eqs. (12), (34) to calculate $v_{R,\text{th}}$ and $v_{\theta,\text{th}}$. We then solve Eq. (29) to obtain the vertical scale height $H_{\text{th}}$ and use it for calculating $v_{z,\text{th}}$ from Eq. (27). We finally use Eq. (33) to calculate the surface density morphology $\Sigma_{\text{th}}$. Results are plotted as 2D maps where the region with $a < a_{\text{cav}}$ has been excised – where $a_{\text{cav}} = 2.7a_{\text{bin}}$ is the semi-major axis where $M_a$ is maximum. This is done in order to exclude from our analysis the cavity area in the simulation, where the dynamics of the gas is strongly affected by the binary.

## 4. Numerical simulation versus the analytical model

### 4.1. Eccentric model vs. circular

We show a comparison between the simulation and the 3D analytical model discussed in the previous section in Figs. 5–8. Velocities $v_R$, $v_\theta$, and $v_z$ are presented in the form of residuals obtained subtracting the theoretical model from the corresponding maps computed from the simulation. For comparison we plot also the residuals obtained subtracting a circular Keplerian velocity field in order to highlight how large would be the systematic error of not accounting for the disc being eccentric, and the typical residual patterns obtained. We note that subtraction of kinematic models different from circular, non-tilted discs is not yet a common practice when analysing line emission observations.

More specifically, Fig. 5 shows a comparison between the normalised residuals for $v_R$ and $v_\theta$. In the same figure, the residuals obtained subtracting $v_{R,\text{circ}} = 0$ and $v_{\theta,\text{circ}} = \sqrt{GM_{\text{bin}}/R}$. We

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

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

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Profiles of $e(a)$ (top panel), $\bar{\omega}(a)$ (middle panel), and $M_a$ (bottom panel) from our reference simulation (blue curves) and corresponding filtered-interpolated profiles (orange curves). For clarity, we remark that the orange curves are obtained by filtering and interpolating with cubic splines the blue curve. The orange curve is used as an input to generate the theoretical model based on the numerical simulation.</figcaption></figure>

normalise the residuals for $v_\theta$ with $v_{\theta,\text{circ}}$, that is, the circular Keplerian velocity at the $x-y$ pixel radius; we normalise $v_R$ using $e_{\text{cav}}v_{\theta,\text{circ}}$ where $e_{\text{cav}} = 0.2$, which represent the maximum $v_R$ achieved by the material at the cavity edge.

Figure 6 shows a comparison between the disc aspect ratio $H/R$ in the simulation and the one calculated from Eq. (29).
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Maps of $v_R$ and $v_\theta$. Top panels: comparison of the velocity fields $v_R$ and $v_\theta$ in the form of residuals between the simulation and the theoretical predictions in Eqs. (12) and (34): $(v_{R,\text{sim}} - v_{R,\text{th}})/v_{\theta,\text{circ}}$ (top left panel) and $(v_{\theta,\text{sim}} - v_{\theta,\text{th}})/v_{\theta,\text{circ}}$ (top right panel), respectively. Bottom panels: similar to top panels but subtracting a circular Keplerian velocity profile ($v_{\theta,\text{circ}}$) from the $v_{R,\text{sim}}$ and $v_{\theta,\text{sim}}$ in order to highlight the impact on the residuals of not accounting for the eccentric nature of the disc and of pressure corrections in the radial/azimuthal motion; $(v_{R,\text{sim}} - v_{R,\text{circ}})/v_{\theta,\text{circ}}$ (bottom left panel), where obviously $v_{R,\text{circ}} = 0$, and $(v_{\theta,\text{sim}} - v_{\theta,\text{circ}})/v_{\theta,\text{circ}}$ (bottom right panel). In the bottom right panel it can be clearly seen that the simulation azimuthal velocity $v_{\theta,\text{sim}}$ at large radii is sub-Keplerian, as $v_{\theta,\text{circ}}$ does not account for the pressure support term.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Maps of disc aspect-ratio $H/R$. Left panel: disc $H_{\text{sim}}/R$ from simulation, calculated as $H_{\text{sim}}/R = \langle z^2(R) \rangle^{1/2}/R$; we note that $\langle z^2(R) \rangle = H_{\text{sim}}$ by definition, with it being the vertical second moment of the density distribution. Right panel: disc $H_{\text{th}}/R$ from solving Eq. (29) using $e(a)$, $\omega(a)$, and $c_5(R)$ from the simulation as input.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Maps of disc vertical velocity $v_z$. Top panels: $v_{z,\text{sim}}/c_s$ velocity field from the simulation (top left panel) compared to the theoretical model $v_{z,\text{th}}/c_s$ obtained solving Eq. (29) (top right panel). Bottom panels: residuals of the simulation against the theoretical model $(v_{z,\text{sim}} - v_{z,\text{th}})/c_s$ (bottom left panel), i.e. the top left panel minus the top right panel of this figure, and with respect to a circular disc $(v_{z,\text{sim}} - v_{z,\text{circ}})/c_s$ (bottom right panel), i.e. $v_{z,\text{circ}} = 0$; this last plot is in fact the same as top left panel but with rescaled colourbar for a direct comparison with the other residuals plot in the bottom left panel. We note that in the simulation the *m* = 1 spiral shaped vertical velocity feature standing out even after the subtraction of the eccentric theoretical model suggests that some additional physical processes are affecting the disc dynamics.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Σ density field from the simulation (left panel) and relative analytical model (right panel) from Eq. (33).</figcaption></figure>
Figure 7 shows a comparison between the vertical velocity field in the simulation $v_{z,sim}$ against the one obtained from the theoretical model $v_{z,th}$. In Fig. 7, we also show the residuals from the subtraction of the two, and one where $v_{z,circ} = 0$ is subtracted. This allows us to perform a direct comparison of the order of magnitude of the residuals when vertical oscillations are neglected, e.g. when considering a circular disc model. We normalise the residuals for $v_z$ with $c_s$: as mentioned in footnote 2, that is a natural scale for eccentric vertical velocity perturbations at $z = H$.

Finally, in Fig. 8, we show a comparison between the surface density $\Sigma$ for both the simulation and the theoretical model, obtained from Eq. (33).

## 4.2. Agreement between the simulation and the model

We here discuss more quantitatively the agreement between the numerical simulation and the theoretical model results presented in the previous sections.

The residuals shown in Fig. 5 highlight that the azimuthal and radial motion in the simulation are well described by the theoretical model.

Residuals of the $v_R$ field show a maximum residual value of $|\Delta v_R|_{max}/(e_{cav}v_{\theta,circ}) \sim 15\%$ located at the cavity pericentre, where a tidal stream of material is launched from the tidal interaction with the binary (which is not captured by the theoretical model). The average residual value is $|\Delta v_R|_{avg} \le 5\%$, the pattern of these residuals suggest that the disc wobbles around the centre of the numerical domain. This is expected, since the disc has finite mass, both the centre of mass of the disc and that of the binary orbit around the centre of mass of the binary + disc system (which by construction is located in the centre of the domain in $(x, y) = (0, 0)$). Not considering the corrections due to the disc eccentric nature for $v_R$, that is, in a circular disc, produces $|\Delta v_{R,circ}|_{max}/(e_{cav}v_{\theta,circ}) \sim 100\%$. The comparison with the circular case, shows as expected an increase of the radial velocity reaching its maximum in absolute value $|v_R|$ at true anomaly $f = \pi/2$ and $f = 3\pi/2$.

The $v_\theta$ field similarly has a maximum residual value of $\Delta v_\theta/v_{\theta,circ} \sim 2\%$, again located at the cavity pericentre where the tidal interaction of the binary is strongest, resulting in a transfer of angular momentum and energy to the disc. In the comparison with the circular disc case, the residuals grow to $|\Delta v_{\theta,circ}|/v_{\theta,circ} \sim 15\%$. Also in this case, the comparison with the circular case, shows that the azimuthal velocity reaches its maximum at pericentre and minimum at apocentre, as expected.

Figure 6 shows a comparison between the disc vertical scale height in the simulation and the theoretical prediction in the form of the disc aspect ratio $H/R$. Both the numerical simulation and theoretical model highlight the presence of a 'breathing' mode that causes the disc to compress towards the midplane at pericentre and an expansion at apocentre. The difference between the disc thickness at pericentre and apocentre at the edge of the cavity reaches $H_{apo}/H_{peri} \sim 3$, as visible in Fig. 2 for eccentricity $e \approx e_{cav} = 0.2$.

The vertical velocity field $v_z$ associated with this expansion and compression at apocentre and pericentre, respectively, is plotted in Fig. 7. We can recognise a good overall agreement between the simulation and the theoretical model. When compared with the residuals for the circular case, where $v_{z,circ} = 0$, the theoretical eccentric disc model appears to perform much better in reproducing the numerical simulation, producing overall a reasonable agreement.

However, some differences are evident and should not be ignored. In particular, an $m = 1$ spiral feature in the numerical simulation appears not to be captured by the theoretical model. Furthermore, the $v_z$ field in the numerical simulation shows a small level of asymmetry with respect to the pericentre longitude: that is, the positive $v_z$ lobe looks slightly different in shape with respect to the negative one. Finally, the simulation $v_z$ map appears to show a slight shift with respect to the pericentre and apocentre of the location where $v_z$ vanishes, in contrast with what predicted by the theoretical model.

We identify two possible origins of these effects, that are not mutually exclusive: on the one hand, the binary nature of the central source of gravity in the simulation is responsible for the additional features; on the other hand, the additional features in the simulation might be captured by the eccentric disc theoretical model if we considered also the effect of viscosity when calculating the disc vertical dynamics (Eqs. (27) and (29)).

We explored this second possibility by including the bulk viscosity stress term in Eq. (20) – the derivation of the final equation including bulk viscosity can be found in Appendix B. A bulk viscosity component is expected to be spuriously present in any numerical scheme due to numerical dissipation. As explained in the appendix, the effect of this additional viscous term produce some asymmetry in the vertical velocity field and a shift of the location where $v_z$ vanishes, consistently with the one observed in the simulation. However, the new theoretical prescription for $v_z$ fails in any case to reproduce the $m = 1$ spiral pattern visible in the residuals of $v_z$. For this reason we decided to keep the model in a simpler form, and not to include the comparison with the viscous theoretical model for $v_z$, which is beyond the scope of the current paper.

In conclusion, we cannot identify the precise origin of the discrepancies in the predicted $v_z$ and the simulated one. Such a spiral feature in the $v_z$ field of the simulation appears to be related to the presence of the binary or to spurious dissipation resulting in an effective bulk viscosity. We postpone a thorough investigation about the origin of this discrepancy to a future work.

Spiral features in the vertical $v_z$ have been previously described in planet-disc simulations by Bae et al. (2021), and have been referred to as ‘buoyancy spirals’. Even though the spiral features in our simulations possibly share a similar dynamical ‘trigger’, that is, the presence of a secondary mass in the system, buoyancy effects cannot be excited in our vertically isothermal numerical simulations, implying that vertical thermal stratification is not a fundamental requirement for these features to arise, at least in the $q \gtrsim 0.1$ mass-ratio regime. To our knowledge, such spiral features in the vertical $v_z$ have not been documented before in the literature of locally isothermal circumbinary disc simulations.

More generally, a new analysis of the numerical set in Ragusa et al. (2020) reveals that other simulations show similar spiral features in the disc vertical velocity field. We missed this at first because our standard approach to produce 2D maps of 3D simulations requires vertically averaging the fluid properties. Given the anti-symmetric nature of the perturbation with respect to the midplane, vertical integration of the $v_z$ 3D field produces $v_z = 0$ maps in the entire disc $x-y$ domain.

We remark that this perspective is particularly intriguing, in the light of the detection of a spiral kinematic feature possibly associated with vertical motion in the $^{12}$CO emission of the system HD142527 (Garg et al. 2021), that features a well known circumbinary eccentric disc.
# 5. Synthetic observations

In this section we investigate the physical magnitude of the eccentric kinematic features discussed in the previous sections trying to assess their observability. To do so we rescale the hydrodynamical simulation and perform Monte Carlo radiative transfer on it, in order to explicitly check whether the velocity perturbations are within the current observational capabilities.

## 5.1. Monte Carlo radiative transfer

We perform Monte Carlo radiative transfer simulations using the code MCFOST (Pinte et al. 2006, 2009) to produce channel maps of the line emission from the 3–2 transition of the $^{12}$CO. For this purpose, we rescale the numerical simulations for the binary to have a separation of 15 au, by multiplying the spatial quantities by $x_{sc} = 15$. This results in a cavity size in the disc of ≈ 40 au. The two stars are assumed to have masses $M_1 = 0.91 M_\odot$ and $M_2 = 0.09 M_\odot$ and to radiate as two black bodies with temperature $T_1 = 4200$ K and $T_2 = 2900$ K consistently with an age of 3 Myr (Siess et al. 2000), respectively. Following this assumption about the binary mass, we rescale the velocities by $v_{sc} = 29.7/\sqrt{15}$ to have the velocities in km s$^{-1}$. The total gas mass is $5 M_1$, while the dust component is included reflecting the gas density distribution rescaled by a factor 100. We remark that our purpose is to produce mock synthetic observations of the eccentric kinematic features discussed in the previous sections. As a consequence, a self-consistent evolution of both gas and dust is beyond the scope of the paper. We assume a chemical abundance ratio$^6 \rho_{12\text{CO}}/\rho_{\text{H}} = 5 \times 10^{-5}$.

We here perform the radiative transfer analysis with the sole purpose of investigating the qualitative appearance of the gas kinematics, which is instead not affected by this assumption. A self-consistent treatment of the dynamics of the gas + dust mixture is beyond the scope of the paper.

Dust opacities are calculated using Mie theory assuming astrosilicate grains (Weingartner & Draine 2001), with sizes ranging between $s = 0.03–1000$ µm and a distribution $dn(s)/ds \propto s^{7/2}$.

We use $N_\gamma = 10^7$ photons for calculating the temperature structure, which is then used for determining the line emission properties; the final channel maps are then calculated through ray-tracing with $N_\gamma = 10^7$ photons.

Angular scales and fluxes are calculated placing the source at a distance $d = 130$ pc, consistently with the distance of the Taurus region. We study the appearance of the system for two inclinations: $i = 0°$ – that is, face-on – and $i = 30°$. These values are chosen in order to study the detectability of the vertical oscillations and of the other eccentric features predicted in the previous sections for both face-on and for moderately inclined discs.

The result of the radiative transfer is a datacube with ‘infinitive’ spatial resolution (set by the hydrodynamical simulation), while spectral resolution bins in the cube are $\Delta v = 0.004$ km s$^{-1}$ for the face-on case and $\Delta v = 0.01$ km s$^{-1}$ for $i = 30°$.

## 5.2. Observability of the eccentric features

With the purpose of extracting a projected velocity map of the disc surface from the radiative transfer data cube, we produce moment 1 maps (M1, i.e. channels are collapsed attributing to

each pixel the brightness-weighted average of the velocity from brightness-velocity profile) and moment 9 (M9, i.e. maps where the channels are collapsed attributing the value of the ‘brightest’ velocity in the brightness-velocity profile of each pixel) maps, shown in Fig. 9. A selection of channel maps can be found for completeness in Fig. D.1. In order to qualitatively reproduce a real observation, we convolve the image spatially with a circular, Gaussian beam $\Delta s = 0.05$ arcsec and the velocity with $\Delta v = 0.05$ km s$^{-1}$ (i.e. within ALMA observational capabilities).

It can be easily noted that the collapse of the channel maps into the moment 1 map significantly underestimates the vertical velocity component in the $i = 0°$ case compared to the moment 9. This can be understood by looking at the brightness-velocity plot in the bottom panels of Fig. 9. The emission in each spatial pixel of the data cube has a broad spectrum that depends on the projected velocity and temperature in the different disc layers encountered along the line of sight – the temperature stratification typical of protostellar discs (i.e. hot on the disc surface, cold close to the disc midplane) explains also the secondary peaks that are visible in both $i = 0°$ and $i = 30°$ pixel spectra (Pinte et al. 2018). Producing an M1 map implies taking the brightness averaged velocity as the representative velocity of each pixel: with brightness profiles such as those shown in the right panels of Fig. 9, M1 maps strongly underestimate the velocity of the upper emission surface in the pixel. M9 maps better capture the fact that the velocity in the upper layers is $v_z \neq 0$, but remain far from a precise estimate.

To strengthen this statement, in Fig. 10 we show the brightness map of the $0.72$ km s$^{-1}$ channel from the radiative transfer data cube of the $i = 0°$ simulation. We can clearly see that, even though the M9 map shows a maximum vertical velocity of $v_{z,max} \approx 0.25$ km s$^{-1}$, higher velocity channels can be quite bright and detectable using ALMA. We also note that in the top layer, the two lobes have opposite velocities, while in Fig. 10 they appear almost equally bright with positive velocities: that is, each pixel contains emission from both the top and bottom layers. See also the following section for a direct comparison of the moment maps with the upper emission surface kinematics.

Despite the poor performance of moment maps at portraying the upper emission surface kinematics, M9 maps show recognizable patterns of the disc’s eccentric nature both for face-on discs and for inclined ones: in the form of vertical motion and asymmetry between the projected velocity at apocentre and pericentre, respectively.

In conclusion, both M9 and M1 maps are not representative of the velocities in the upper emission surface, so we cannot provide a direct comparison of the radiative transfer with the theoretical model. However, the typical eccentric features, that is, double anti-symmetric vertical velocity pattern and asymmetry between apocentre and pericentre are recognisable. More sophisticated tools are required to properly reconstruct the geometry of the upper emission surface from channel maps (e.g. DISCMINER, Izquierdo et al. 2021, 2023), in order to carry out a detailed comparison with the model.

In the light of these considerations, we do not attempt a direct comparison of the radiative transfer results with the theoretical model. Instead, we limit our analysis to inclining the simulation data and compare them with an inclined theoretical model generated for this purpose in the next section.

## 5.3. Modelling the disc upper emission surface

In the previous section we saw we cannot directly compare the results from the radiative transfer moment maps with the

<footnote n="6">⁶ Our choice of $\rho_{12\text{CO}}/\rho_{\text{H}}$ instead of the typical $\rho_{12\text{CO}}/\rho_{\text{H}} = 10^{-4}$ constitutes a conservative estimate for the vertical velocity. An optically thinner disc receives less contribution in the line emission from the highest disc layers, that have larger vertical velocities.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Moment maps from the RT simulations. Top panels: Moment 1 maps from Monte Carlo radiative transfer simulations *i* = 0° (left) and *i* = 30° (right). The grey circle in the bottom left corner represents the corresponding beam size. Middle panels: Moment 9 maps, again obtained from the Monte Carlo radiative transfer of the simulation dump, *i* = 0° (left) and *i* = 30° (right). Bottom panels: Brightness-velocity maps, *i* = 0° (left) and *i* = 30° (right), for selected pixels for highlighting the large span of different velocity contributions in individual pixels. The maps are obtained after convolving spatially the image with a Gaussian beam of Δ*s* = 0.05 arcsec (marked as a grey circle in the bottom left corner of each image) and Δ*v* = 0.05 km s⁻¹. We note that the convention about velocities in observations is that blue-shift of lines is associated with negative velocities, while red-shift is associated with positive velocities – for example, this implies that *v*<sub>z</sub> in this plot is opposite in sign with respect to that in Fig. 7. The case *i* = 0° shows well how M1 maps tend to underestimate the vertical velocity contribution compared to M9 maps.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Brightness map of the 0.72 km s⁻¹ channel from the data cube of the Monte Carlo radiative transfer, face-one *i* = 0° inclination, performed on the simulation dump. The image clearly shows that, even though the M9 map in Fig. 9 shows maximum velocity of ~ 0.25 km s⁻¹ for vertical velocities, a bright contribution in higher velocities channels is also expected to be observable.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Map of the ratio between $H_{\tau=1}$, obtained from radiative transfer, tracing the upper emission surface and $H_{sim}$, representing the simulation scale height. The ratio $H_{\tau=1}/H_{sim} \approx 2-3$ across the disc. We consider $H_{\tau=1} \approx 2.5H_{sim}$ a conservative estimate of the location of the upper emission surface – ‘conservative’ since the higher in the disc atmosphere, the larger the vertical velocity.</figcaption></figure>

theoretical model, as it is not possible to define precisely the velocities of the upper emission surface. Let us assume that we observe a circumbinary disc and that we have the tools to properly reconstruct from the channel maps the velocity projected on the observer’s line of sight in the upper disc layers – this is motivated by the current development of new observational tools for extracting the kinematics of the upper emission surface (DISCMINER, Izquierdo et al. 2021). In this section, we investigate how the residuals from the subtraction of an eccentric disc model compare to those from a circular model, and how strongly the eccentric nature of the disc stands out compared to the circular model.

To answer these questions, we create a 3D model of the disc for both the theoretical models (eccentric and circular) and the numerical simulation, and we perform a direct comparison. The 3D model consists of spatial coordinates $S = \{x, y, z\}$ defining the 3D disc upper emission surface, and velocities $v = \{v_x, v_y, v_z\}$ as the velocity vector field of the material on it.

We first define the altitude of the disc upper surface using the altitude of the $\tau = 1$ surface above the midplane $H_{\tau=1}$ output from the Monte Carlo radiative transfer simulation. The value of $\tau$ is calculated integrating along the line of sight the optical depth using the local maximum opacity $\max_v[k_v(x, y, z)]$ (i.e. the centre of the local line), where $k_v(x, y, z)$ is the local opacity, that is the surface $z(x, y)$ satisfying

$$ \tau = \int_{+\infty}^{z(x,y)} \rho(x, y, z) \max_v[k_v(x, y, z)] dz = 1. \quad (38) $$

Such a surface qualitatively represents the highest elevation in the disc atmosphere at which the material contributes to the line emission, assuming the disc is face-on. We plot the altitude of the upper emission surface as $H_{\tau=1}/H_{sim}$ for a comparison with the vertical scale height of the simulation in Fig. 11. The average value $\langle H_{\tau=1}/H_{sim} \rangle \approx 2.5$, and it appears to be relatively constant throughout the entire disc. For this reason, we define our disc surface to be $z = 2.5H_{sim}$ for the simulation and $z = 2.5H_{th}$ for the circular and eccentric theoretical models.

Since Eq. (38) is based on the maximum local opacity, it does not imply necessarily that the material at $z < H_{\tau=1}$ is not

visible, as that would be the case only if $v_z$ = const throughout the entire gas column. Indeed, $\tau$ defined in Eq. (38) should not be confused with $\tau_v$, which instead is a function of the frequency, and is obtained by integrating the opacity using the local $k_v$ (and not its local maximum as done for $\tau$). Since the opacity $k_v$ of the material changes with the material velocity, $\tau_v = 1$ qualitatively defines a set of iso-velocity surfaces (one for each value of $v$, that is, one for each velocity channel), probing different altitudes above the mid-plane. In fact, different layers of the disc atmosphere along the line of sight contribute to the line emission at different frequencies, despite lying at altitudes $z < H_{\tau=1}$.

The disc has a top and a bottom surface, thus for our 3D eccentric and circular model we define them to be $S_{top} = \{x, y, 2.5H_{th}\}$ for the top surface and $S_{bottom} = \{x, y, -2.5H_{th}\}$.

We define the velocity field on the disc surfaces as $v_{top} = \{v_x, v_y, v_z\}$ for the top layer and $v_{bottom} = \{v_x, v_y, -v_z\}$ for the bottom layer, where $v_x$ and $v_y$ are obtained projecting $v_R$ and $v_\theta$ from Eqs. (12) and (34) on to the $x$ and $y$ unit vectors, while $v_z$ is calculated from Eq. (27), calculated at $z = 2.5H_{th}$. We produce a model with the same characteristics also for a circular disc, by projecting $v_{circ}$ on to Cartesian unit vectors.

Similarly, after rescaling the numerical simulation as described in Sect. 5.1 ($a_{bin} = 15$ au, $M_{bin} = 1 M_\odot$, $d = 130$ pc, that is, $x_{sc} = 15$, $v_{sc} = 29.7/\sqrt{15}$), we define the upper surface in the simulation to be $S_{sim} = x_{sc}\{x, y, 2.5H_{sim}\}$ and $v = v_{sc}\{v_{x,sim}, v_{y,sim}, 2.5v_{z,sim}\}$, where, we note, we applied an additional scaling of a factor 2.5 to the original $z_{sim}$ and $v_{z,sim}$, to account for the higher elevation of the upper emission surface compared to the vertical scale-height set by pressure – this operation exploits the fact that $v_z \propto z/H_{th}$.

## 5.4. Morphology of projected velocity residuals in circumbinary discs

Consider an observer located at $z = +\infty$, and define the projected velocity $v_{proj} = -v_z$ (following the convention that redshifted velocities are positive). At this point, the observer, looking down on the $x-y$ plane observing $v_{proj}$ will see the kinematics and geometry of a face-on disc. We do this for both the eccentric and
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Residuals of the subtraction of the eccentric theoretical model (left panels) and circular model (right panels) from the projected velocity in the simulation for *i* = 0° (top panels) and *i* = 30° (bottom panels). These maps were obtained assuming the upper emission surface to be *H*<sub>*τ*</sub>=1 = 2.5*H* (as discussed in Sect. 5.3), and rescaling the simulations to have a binary with *a*<sub>bin</sub> = 15 au, *M*<sub>bin</sub> = 1 *M*<sub>☉</sub>, placed at a distance of *d* = 130 pc. Similarly to Fig. 9, we note that the convention about velocities in observations is that blue-shift of lines is associated with negative velocities, while red-shift is associated with positive velocities – for example, this implies that *v*<sub>*z*</sub> in this plot is opposite in sign with respect to that in Fig. 7. We note that the top right panel residual is obtained by subtracting *v*<sub>*z*,circ</sub> = 0 (as the circular case does not have vertical motion) from *v*<sub>*z*</sub>, thus this residuals map portrays the vertical velocity map on the *τ* = 1 surface for the case *i* = 0°. This highlights that the amplitude of the vertical oscillations reaches *v*<sub>*z*</sub> ~ 0.4 km s<sup>-1</sup> in the upper emission layer.</figcaption></figure>

circular models, which constitute our face-on models. For the
eccentric model the line of apses inherits its orientation from the
simulation.

Inclined models are generated by rotating the face-on discs (described by **S** and **v**), keeping the observer (and associated coordinate system) fixed. First we perform a rotation around the z-axis, which sets the longitude of pericentre (this step is omitted in this work, meaning the inclined discs inherit the longitude of pericentre of the face-on disc). We then rotate around the x-axis by *i* = 30°, resulting in an inclined disc with its line of nodes aligned with the x-axis. The longitude of pericentre is unchanged by this rotation, however the projected line of apses will be shifted towards the line of nodes as a result of foreshortening. Finally, without loss of generality, we set the position angle, PA = −90° (calculated from the north anti-clockwise), as the position angle has no effect on the imaged disc.

In general, the kinematics of inclined eccentric disc is
affected by two angles: the inclination, which sets the relative
contribution of the vertical and horizontal velocity; and the argu-
ment of periapsis, which sets the relative contribution of the

on-apse and off-apse velocities to the horizontal velocities. For
symmetric beams, considered here, the position angle has no
effect on the image.

We now have three different models of the geometry and
projected velocity kinematics of the upper emission surface of
a disc inclined by *i* = 30°. One for the circumbinary disc simula-
tion, the other two for the eccentric and circular disc theoretical
models.

For the analysis described in this section, we consider only
*S*<sub>top</sub>, under the assumption we made above that we can isolate the
kinematics coming exclusively from the top emission layer. We
subtract both the ‘eccentric theoretical model’ and the ‘circular
theoretical model’ from the simulation, the residuals of these two
operations are shown in Fig. 12.

A large, trailing, *m* = 1 spiral with a velocity residual of 0.15 km s⁻¹ appears in the *v*<sub>proj</sub> kinematics (i.e. detectable with ALMA) when subtracting the eccentric theoretical model and the circular model. The eccentric theoretical model however appears to perform well in capturing all the other eccentric kinematic features that are not reproduced by a purely circular model.
Indeed, residuals from the subtraction of the circular model from the simulation show instead a broad single-lobed azimuthal kinematic feature at the edge of the disc cavity with $\Delta v > 0.3 \text{km s}^{-1}$. Interestingly, such a feature does not have the appearance one might naively expect when subtracting a circular model from an eccentric disc: a double lobed feature (two lumps in the residuals with opposite signs) due to the faster and slower velocity of the material at pericentre and apocentre, respectively. However, the structure of the residuals will change when orienting the system in a way that results in the disc having a different argument of periapsis (e.g. it is possible to get a two lobed feature with a different disc orientation).

In the light of the results presented here, it is important to consider the possible eccentric nature of the disc for explaining single and double-lobed structures in the projected velocity residuals. We also remark that we cannot be conclusive about the origin of the spiral-shaped velocity feature, not captured by the eccentric disc model. However, when observed, it might hint at the presence of a hidden (sub)stellar companion (if undetected).

## 6. Conclusion

In this paper we discussed the morphology and kinematics of eccentric discs, with the final goal of providing observational signatures, aiding their identification in observational campaigns. To do so, we presented a theoretical model that takes as input information only three functions: namely, the disc eccentricity profile $e(a)$, the longitude of pericentre profile $\omega(a)$, the derivative of the disc mass distribution $M_a(a)$, and an assumption about the locally isothermal disc sound speed ($c_s(a)$). Eccentric discs show a number of characteristic features:

(i) The azimuthal and radial velocities $v_\theta v_R$, along with the surface density $\Sigma$ vary around the orbit (Eqs. (12), (34) and (33)).

(ii) The changing vertical gravity around the eccentric orbit excites a vertical oscillation, corresponding to a variation of the disc scale height $H$ around the orbit (Eq. (29)), and results in a vertical velocity $v_z \propto z/H$ (Eq. (27)).

(iii) The theoretical model explains many features of a numerical simulation of an eccentric circumbinary disc. However, a residual, $m = 1$, spiral pattern in $v_z$ is not captured by the theoretical model. We speculate this could be produced by the binary potential or the effect of viscosity, neither included in our analytical treatment. The origin of this feature will be explored in future studies. We estimate this residual velocities in the spirals to be of the order of $\Delta v_z \approx 0.15 \text{ km s}^{-1}$.

(iv) We perform radiative transfer Monte Carlo on our numerical simulation snapshot to obtain information about $^{12}$CO line emission for inclination $i = 0^\circ$ (face-on) and $i = 30^\circ$. We find that moment 1 maps are not suitable to highlight the vertical motion in the $i = 0^\circ$ case, since $v_z$ is anti-symmetric with respect to the mid-plane. Better results can be obtained with moment 9 maps.

(v) Both $i = 0^\circ$ and $i = 30^\circ$ radiative transfer M9 maps show features of eccentric disc evolution that are within ALMA observational capabilities: the first showing a vertical velocity contribution that is not present in circular face-on discs; the second showing an asymmetric butterfly diagram of the isovelocity contours, indicating different velocity at pericentre and apocentre.

(vi) In the face-on case ($i = 0^\circ$), the M9 map reveals the vertical velocity reaches a maximum of $v_{z,max} \sim 0.25 \text{ km s}^{-1}$. However, channel maps show significant signal at higher velocities from the upper layers of the disc, making $v_z$ eccentric

signatures potentially detectable in face-on discs even for eccentricities $e \ll 0.2$ (e.g. $v_{z,max} \sim 0.4 \text{ km s}^{-1}$ on the $\tau = 1$ surface in Fig. 12).

(vii) The subtraction of a circular Keplerian disc model from an eccentric disc's kinematics can, counterintuitively, produce a single lobed pattern in the projected velocity map. Eccentric disc kinematics should be considered when patterns such as those visible in Fig. 12 are observed. However, other patterns are possible for different disc orientations.

In future works we will further investigate the origin of the spiral velocity pattern discussed in (iii), determining whether such a feature can be considered a smoking gun signature of the presence of a binary (sub)stellar companion or a feature of any eccentric disc. In any case, the results presented in this paper represent an intriguing starting point for the interpretation of similar spiral patterns in the vertical velocity field that are starting to be observed (e.g. Garg et al. 2021).

**Acknowledgements.** We thank the anonymous referee for their comments to the manuscript. E.R. acknowledges Feng Long for valuable feedback and suggestions to the manuscript and Giovanni Rosotti, Stefano Facchini, Giuseppe Lodato and Claudia Toci for fruitful discussion. E.R. acknowledges financial support from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101102964 (ORBIT-D). E.R., E.L. and G.L. acknowledge financial support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 864965, PODCAST). This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 823823 (DUSTBUSTERS). C.L. acknowledges financial support from the UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1. The simulation performed for this paper used the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). Figure 3 was created using SPLASH (Price 2007). All the other figures were created using MATPLOTLIB python library (Hunter 2007).

## References

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Ataiee, S., Pinilla, P., Zsom, A., et al. 2013, A&A, 553, L3

Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

Bonnerot, C., Rossi, E. M., Lodato, G., & Price, D. J. 2016, MNRAS, 455, 2253

Booth, M., Schulz, M., Krivov, A. V., et al. 2021, MNRAS, 500, 1604

Calcino, J., Price, D. J., Pinte, C., et al. 2019, MNRAS, 490, 2579

Calcino, J., Christiaens, V., Price, D. J., et al. 2020, MNRAS, 498, 639

Cao, R., Liu, F. K., Zhou, Z. Q., Komossa, S., & Ho, L. C. 2018, MNRAS, 480, 2929

Cauley, P. W., Farihi, J., Redfield, S., et al. 2018, ApJ, 852, L22

Chan, C.-H., Piran, T., & Krolik, J. H. 2022, ApJ, 933, 81

Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

Cufari, M., Coughlin, E. R., & Nixon, C. J. 2022, ApJ, 924, 34

D'Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698

D'Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, MNRAS, 436, 2997

Dempsey, A. M., Muñoz, D. J., & Lithwick, Y. 2021, ApJ, 918, L36

Dennihy, E., Clemens, J. C., Dunlap, B. H., et al. 2018, ApJ, 854, 40

Dennihy, E., Xu, S., Lai, S., et al. 2020, ApJ, 905, 5

Dewberry, J. W., Latter, H. N., Ogilvie, G. I., & Fromang, S. 2020, MNRAS, 497, 451

Dittmann, A. J., & Ryan, G. 2022, MNRAS, 513, 6158

Dong, R., Liu, S.Y., Eisner, J., et al. 2018, ApJ, 860, 124

Dunhill, A. C., Alexander, R. D., & Armitage, P. J. 2013, MNRAS, 428, 3072

Dunhill, A. C., Cuadra, J., & Dougados, C. 2015, MNRAS, 448, 3545

Faramaz, V., Krist, J., Stapelfeldt, K. R., et al. 2019, AJ, 158, 162

Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014, ApJ, 783, 134

Franchini, A., Lupi, A., Sesana, A., & Haiman, Z. 2023, MNRAS, 522, 1569

Garg, H., Pinte, C., Christiaens, V., et al. 2021, MNRAS, 504, 782

Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024

Goodchild, S., & Ogilvie, G. 2006, MNRAS, 368, 1123

Heath, R. M., & Nixon, C. J. 2020, A&A, 641, A64

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67
Hung, T., Foley, R. J., Ramirez-Ruiz, E., et al. 2020, ApJ, 903, 31

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

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021, A&A, 650, A179

Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, A&A, 674, A113

Kley, W., & Dirksen, G., 2006, A&A, 447, 369

Kley, W., Papaloizou, J. C. B., & Lin, D. N. C. 1993, ApJ, 409, 739

Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671

Kozdon, J., Brittain, S. D., Fung, J., et al. 2023, AJ, 166, 119

Kuo, I. H. G., Yen, H.-W., Gu, P.-G., & Chang, T.-E. 2022, ApJ, 938, 50

Kurban, A., Zhou, X., Wang, N., et al. 2023, MNRAS, 522, 4265

Latter, H. N., & Ogilvie, G. I. 2006, MNRAS, 372, 1829

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lee, W.-K., Dempsey, A. M., & Lithwick, Y. 2019a, ApJ, 872, 184

Lee, W.-K., Dempsey, A. M., & Lithwick, Y. 2019b, ApJ, 882, L11

Li, Y.-P., Chen, Y.-X., & Lin, D. N. C. 2023, MNRAS, 526, 5346

Liu, F. K., Zhou, Z. Q., Cao, R., Ho, L. C., & Komossa, S. 2017, MNRAS, 472, L99

Lodato, G. 2008, New A Rev., 52, 21

Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Lovascio, F., Commerçon, B., Lynch, E., & Ragusa E. 2024, A&A, submitted

Lovell, J. B., & Lynch, E. M. 2023, MNRAS, 525, L36

Lubow, S. H. 1991a, ApJ, 381, 259

Lubow, S. H. 1991b, ApJ, 381, 268

Lynch, E. M. 2022, MNRAS, 510, 3460

Lynch, E. M., & Lovell, J. B. 2022, MNRAS, 510, 2538

Lynch, E. M., & Ogilvie, G. I. 2021, MNRAS, 501, 5500

Lyubarskij, Y. E., Postnov, K. A., & Prokhorov, M. E. 1994, MNRAS, 266, 583

MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83

MacGregor, M. A., Wilner, D. J., Rosenfeld, K. A., et al. 2013, ApJ, 762, L21

MacGregor, M. A., Hurt, S. A., Stark, C. C., et al. 2022, ApJ, 933, L1

Manser, C. J., Gänsicke, B. T., Marsh, T. R., et al. 2016a, MNRAS, 455, 4467

Manser, C. J., Gänsicke, B. T., Koester, D., Marsh, T. R., & Southworth, J. 2016b, MNRAS, 462, 1461

Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493

Miranda, R., & Rafikov, R. R. 2018, ApJ, 857, 135

Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170

Muñoz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106

Murray, J. R. 1996, MNRAS, 279, 402

Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (UK: Cambridge University Press)

Öberg, K. I., Guzmán, V. V., Walsh, C., 2021, ApJS, 257, 1

Ogilvie, G. I. 2001, MNRAS, 325, 231

Ogilvie, G. I. 2003, MNRAS, 340, 969

Ogilvie, G. I. 2008, MNRAS, 388, 1372

Ogilvie, G. I., & Barker, A. J. 2014, MNRAS, 445, 2621

Ogilvie, G. I., & Lynch, E. M. 2019, MNRAS, 483, 4453

Ogilvie, G. I., & Proctor, M. R. E. 2003, J. Fluid Mech., 476, 389

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Olofsson, J., Milli, J., Thébault, P., et al. 2019, A&A, 630, A142

Oyang, B., Jiang, Y.-F., & Blaes, O. 2021, MNRAS, 505, 1

Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263

Pierens, A., McNally, C. P., & Nelson, R. P. 2020, MNRAS, 496, 2849

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47

Price, D. J. 2007, PASA, 24, 159

Price, D. J., Wurster, J., Tricco, T. S., et al. 2018a, PASA, 35, e031

Price, D. J., Cuello, N., Pinte, C., et al. 2018b, MNRAS, 477, 1270

Pringle, J. E. 1981, ARA&A, 19, 137

Ragusa, E., Lodato, G., & Price, D. J. 2016, MNRAS, 460, 1243

Ragusa, E., Dipierro, G., Lodato, G., Laibe, G., & Price, D. J. 2017, MNRAS, 464, 1449

Ragusa, E., Rosotti, G., Teyssandier, J., et al. 2018, MNRAS, 474, 4460

Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362

Regály, Z., Sándor, Z., Dullemond, C. P., & Kiss, L. L. 2011, A&A, 528, A93

Savitzky, A., & Golay, M. J. E. 1964, Analyt. Chem., 36, 1627

Scardoni, C. E., Clarke, C. J., Rosotti, G. P., et al. 2022, MNRAS, 514, 5478

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

Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118

Shiokawa, H., Krolik, J. H., Cheng, R. M., Piran, T., & Noble, S. C. 2015, ApJ, 804, 85

Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593

Siwek, M., Weinberger, R., Muñoz, D. J., & Hernquist, L. 2023, MNRAS, 518, 5059

Statler, T. S. 2001, AJ, 122, 2257

Syer, D., & Clarke, C. J. 1992, MNRAS, 255, 92

Tanaka, Y. A., Kanagawa, K. D., Tanaka, H., & Tanigawa, T. 2022, ApJ, 925, 95

Teyssandier, J., & Lai, D. 2019, MNRAS, 490, 4353

Teyssandier, J., & Ogilvie, G. I. 2016, MNRAS, 458, 3221

Teyssandier, J., & Ogilvie, G. I. 2017, MNRAS, 467, 4577

Thieme, T. J., Lai, S.-P., Ohashi, N., et al. 2023, ApJ, 958, 60

Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102

Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2020, ApJ, 900, 43

Trevascus, D., Price, D. J., Nealon, R., et al. 2021, MNRAS, 505, L21

van't Hoff, M. L. R., Tobin, J. J., Li, Z.-Y., et al. 2023, ApJ, 951, 10

Weingartner, J. C., & Draine, B. T., 2001, ApJ, 548, 296

Wevers, T., Nicholl, M., Guolo, M., et al. 2022, A&A, 666, A6

Whitehurst, R. 1994, MNRAS, 266, 35

Wienkers, A. F., & Ogilvie, G. I. 2018, MNRAS, 477, 4838

Wilson, D. J., Gänsicke, B. T., Koester, D., et al. 2015, MNRAS, 451, 3237

Wölfer, L., Facchini, S., van der Marel, N., et al. 2023, A&A, 670, A154

Yang, H., Fernández-López, M., Li, Z.-Y., et al. 2023, ApJ, 948, L2

Zanazzi, J. J., & Ogilvie, G. I. 2020, MNRAS, 499, 5562

Zhu, Z., & Stone, J. M. 2014, ApJ, 795, 53
# Appendix A: Other useful pieces of information about eccentric discs

## A.1. Eccentric disc eigenmodes

In this section, we discuss the main aspects relating to the evolution of the eccentricity $e(a, t)$ and pericentre longitude profiles $\varpi(a, t)$ with time.

Under the assumption that no perturbations to the central, point-like, gravitational potential are present and under the assumption that the material in the disc does not self-interact, neither gravitationally nor through hydrodynamical effects (e.g. pressure/viscosity), $e(a)$ and $\varpi(a)$ do not depend on time.

However, for gaseous astrophysical discs the effects of pressure and the presence of additional non-point-like terms of the gravitational potential – for example, disc self-gravity, the presence of a second mass such as a binary star or a planet, oblateness of the central mass or general relativistic corrections – cause the disc to differentially precess at a rate $\omega(a)$, producing a twisted eccentricity pattern and oscillations of the value of $e(a, t)$ and $\varpi(a, t)$. Pressure, disc self-gravity, and viscosity transport these oscillations as waves throughout the disc.

For the purpose of reducing the complexity of the equations describing the dynamics of eccentric discs, $e(a, t)$ and $\varpi(a, t)$ can be effectively merged in one complex quantity $E(a, t)$ (e.g. Ogilvie 2001; Goodchild & Ogilvie 2006; Ogilvie 2008):

$$E(a, t) = e(a, t)e^{i\varpi(a, t)}. \qquad (A.1)$$

We note that the typographical difference between the eccentricity $e$ and the mathematical constant $e$.

Both linear ($|E| \ll 1$ and $|aE_a| \ll 1$, that is $q \ll 1$ in Eq. (6)) and non-linear theory predict the existence of a special set of solutions of the form (e.g. Teyssandier & Ogilvie 2016; Ogilvie & Lynch 2019):

$$E(a, t) = e(a)e^{i\omega t}. \qquad (A.2)$$

In absence of dissipation, this expression implies that the disc have a stationary $e(a)$ profile that coherently precesses at the rate $\omega_i$ throughout the whole disc. If dissipation is present the complex phase will also include a twist $\varpi(a)$ of the pericentre orientation, however the whole profile still precesses rigidly conserving the shape of the twist. These configurations are called ‘eccentric eigenmodes’.

Eccentric eigenmodes are analogous to normal quantum states. In the linear regime, $E(a, t)$ can be described as a linear superposition of eccentric eigenmodes

$$E(a, t) = \sum_{i} c_{i} E_{i}(a, t), \qquad (A.3)$$

where $c_i$ are complex coefficients. Out of the linear regime, non-linear eigenmodes can be always identified even though they do not superimpose in the same way.

The coexistence of multiple eccentric eigenmodes precessing at different rates overall produces some peculiar secular oscillations in the shape of the eccentricity profile and precession rate (e.g. see the evolution of the disc eccentricity in Thun et al. 2017 or Ragusa et al. 2018).

## A.2. Linear theory and eccentric eigenmodes

In this section we provide an explanation of the origin of the functional form of eccentric eigenmodes.

In the approximation of linear perturbations from circular orbits, the unperturbed state has $e = 0$ implying $a = R$, making it reasonable to develop the equations in standard cylindrical coordinates.

The evolution of $E(R, t)$ in a 3D locally isothermal disc, due to pressure effects only (i.e. no self-gravity, no companion mass) is governed by (Ogilvie 2008; Ogilvie & Barker 2014; Teyssandier & Ogilvie 2016):

$$\Sigma R^2 \Omega \frac{\partial E}{\partial t} = \frac{i}{R} \frac{\partial}{\partial R} \left( \frac{1}{2} \Sigma c_s^2 R^3 \frac{\partial E}{\partial R} \right) + \frac{i r}{2} \frac{d}{dR} \left( \Sigma c_s^2 \right) E - \frac{i}{2R} \frac{\partial}{\partial R} \left( \Sigma \frac{dc_s^2}{dr} R^3 E \right) + \frac{3i}{2R} \Sigma \frac{d}{dR} \left( c_s^2 R^2 \right) E. \qquad (A.4)$$

We note that eccentric eigenmodes profiles of the form $E_i(R, t) = e_i(R)e^{i\omega_i t}$ cause Eq. (A.4) to reduce to the form of:

$$\omega_i \mathcal{A}(R) e_i(R) = \frac{d}{dR} \left( \mathcal{B}(R) \frac{de_i}{dR}(R) \right) + C(R) e_i(R). \qquad (A.5)$$

where $\mathcal{A}(R)$, $\mathcal{B}(R)$ and $C(R)$ are real functions of the disc parameters.

It can be shown that, after performing some change of variables, Eq. (A.5) can be rewritten in the form of the Shrödinger equation:

$$-\frac{d^2 \Xi_i(x)}{dx^2} + V \Xi_i(x) = \mathcal{E}_i \Xi_i(x), \qquad (A.6)$$

where $x$ and $a$ are related by:

$$\frac{dx}{dR} = \left[ \frac{\omega_0 \mathcal{A}(R)}{\mathcal{B}(R)} \right]^{1/2}. \qquad (A.7)$$

Functions $e_i(R)$ and $\Xi_i(x)$ are related by:

$$\Xi_i(x) = e_i(R)[\mathcal{A}(R)\mathcal{B}(R)]^{1/4}. \qquad (A.8)$$

The effective energy eigenvalue $\mathcal{E}_i$ is given by:

$$\mathcal{E}_i = -\frac{\omega_i}{\omega_0}, \qquad (A.9)$$

where $\omega_0$ relates with the disc thickness ($H/R)_0$ and mean motion $\Omega_0$ at the inner disc edge as follows:

$$\omega_0 = \frac{1}{2} \left(\frac{H}{R}\right)_0 \Omega_0. \qquad (A.10)$$

The effective potential $V(x)$ reads:

$$V = -\frac{C(R)}{\omega_0 \mathcal{A}(R)} + (\mathcal{A}(R)\mathcal{B}(R))^{-1/4} \frac{d^2}{dx^2} [\mathcal{A}(R)\mathcal{B}(R)]^{1/4}. \qquad (A.11)$$

As a consequence, in analogy to quantum mechanics, there exist a countable infinity of eccentric eigenmodes $E_i(R, t)$. Prograde ($\omega_i$) modes correspond to the finite number of negative energy bound states of the potential. These solutions have the form of Eq. (A.2), that is, they feature a stationary eccentricity profile precessing at a rate $\omega_i$.
### A.3. Functional form of the eccentricity profile in circumbinary discs

Following on from the previous section, we discuss the radial eccentricity profile of eccentric eigenmodes. We solve Eq. (A.4) assuming radial power-law profiles of surface density $\Sigma \propto R^{-\sigma}$ and soundspeed $c_s^2 \propto R^{-q}$. Under this assumption, it can be shown that $E(R, t)$ can be expanded in series as (Goodchild & Ogilvie 2006; Ogilvie 2008)

$$E(R, t) = \sum_{i} a_{i}(t) \Psi_{i}(R), \qquad (A.12)$$

where $\{a_i(t)\}$ are complex coefficients, which depend on time only. The functions $\Psi_i(R)$ are functions of radius only and read

$$\Psi_i(R) = R^{\sigma/2+q-1} Z_v \left( (R/R_i)^\beta \right), \qquad (A.13)$$

where $Z_v$ is a form of Bessel function that can be written as

$$Z_v(x) = AJ_v(x) + BY_v(x), \qquad (A.14)$$

where $J_v$ and $Y_v$ are Bessel functions of the first and second kind, respectively, with order $v$ which is related to $q$, $\sigma$, and $\beta$ by $v^2 = \beta^{-2}(2q + q\sigma + q^2 + \sigma^2/4 - 5)$. A, B, and $\{R_i\}$ are set by the (zero gradient) boundary conditions, $\partial_R \Psi_i|_{a_{\text{in}}} = \partial_R \Psi_i|_{a_{\text{out}}} = 0$ and an appropriate normalisation condition.

In a locally isothermal, 3D, powerlaw, disc with only pressure forces then $\beta = (1 + 2q)/4$ and $\Psi_i$ correspond to the linear eccentric eigenmodes (with $\Psi_i = e_i(R)$ as defined in Eq. (A.2)). The fundamental mode is then given by

$$E_0(R, t) = a_0(t)\Psi_0(R). \qquad (A.15)$$

When a quadrupole potential (e.g. from a binary), and nonlinear effects are included the fundamental mode will contain contributions from the other $\Psi_i$'s. However, at order zero, the radial part of the fundamental mode can be written in the form

$$e(R) = \mathcal{A}R^{-\gamma} \left\{ J_v \left( \left[ \frac{R}{R_0} \right]^\beta \right) + C Y_v \left( \left[ \frac{R}{R_0} \right]^\beta \right) \right\}, \qquad (A.16)$$

where $C, \beta, \gamma$, and $R_0$ are constants whose value depend on the powers of the profiles of $\Sigma$ and $c_s$, as discussed above, as well as other disc properties and other assumptions (e.g. boundary conditions), while $\mathcal{A}$ sets the amplitude of the mode.

A simpler functional form can be found by solving Eq. (A.4) by using the WKB approximation (Shi et al. 2012; Lee et al. 2019b; Muñoz & Lithwick 2020). Again, assuming radial power-law profiles of surface density and soundspeed, it can be shown that the eccentricity profile in a disc has the form

$$e(R) = \mathcal{A}R^{-\gamma} \exp \left[ (R/r_0)^b \right]. \qquad (A.17)$$

Similarly to Eq. (A.16), $b, g$, and $r_0$ are constants whose value depend on the powers of the profiles of $\Sigma$ and $c_s$ as well as other disc properties, while $\mathcal{A}$ sets the amplitude of the mode.

### Appendix B: Vertical Structure including bulk viscosity

It is possible to include bulk viscosity for calculating the disc vertical structure, as follows: we introduce the $T^{zz}$ bulk viscous stress tensor term Eq. (20), which becomes:

$$\frac{\partial v^z}{\partial t} + v^a \frac{\partial v^z}{\partial a} + v^\phi \frac{\partial v^z}{\partial \phi} + v^z \frac{\partial v^z}{\partial z} = -\frac{GM}{R^2} \frac{z}{R} - \frac{1}{\rho} \frac{\partial}{\partial z} (p - T^{zz}). \qquad (B.1)$$

where $T^{zz}$ is the bulk viscous stress tensor

$$T^{zz} = \mu_b \left( \frac{1}{J} \frac{\partial}{\partial \phi} (J\Omega) + \frac{\partial v^z}{\partial z} \right), \qquad (B.2)$$

parametrised by a $\mu_b$ given by (Ogilvie 2001; Ogilvie & Barker 2014; Lynch & Ogilvie 2021):

$$\mu_b = \alpha_b p \Omega_0^{-1}. \qquad (B.3)$$

The addition of the $T^{zz}$ term causes Eq. (29) to become

$$\Omega^2 \frac{\partial^2 H}{\partial \phi^2} = -\Omega \frac{\partial H}{\partial \phi} \frac{\partial \Omega}{\partial \phi} - \frac{GM}{R^3} H + \frac{p_0}{\rho_0 H} \left[ 1 - \frac{\alpha_b}{\Omega_0} \left( \frac{1}{J} \frac{\partial}{\partial \phi} (J\Omega) + \frac{\Omega}{H} \frac{\partial H}{\partial \phi} \right) \right]. \qquad (B.4)$$

We note that in order to produce the same plots for $\alpha_b \neq 0$ one needs specify $q$ and $\alpha$ in Eq. (7) and (8) making the parameter space larger. Since the solution would not be general, we do not provide an example plot.

In practice, the main effect of the bulk viscosity term is an azimuthal shift of the location of the minimum and maximum of $H$ with respect to the pericentre and apocentre of the orbit, where minimum and maximum of $H$ are located for $\alpha_b = 0$, respectively. In addition, bulk viscosity produces also small changes in the maximum and minimum vertical velocity reached along the orbit.

## Appendix C: A more general choice of the form of $\rho(a, \phi, z)$ and $p(a, \phi, z)$

It can be shown that the final form of Eq. (28) and (29) does not change for different assumptions on $\rho(a, \phi, z)$ and $p(a, \phi, z)$, with the only caveat that the form of $\rho$ and $p$ must guarantee the separability of the solution in the variables in the midplane $(a, \phi)$, or $(R, \phi)$, and $z$ (see e.g. Ogilvie & Barker 2014).

$$\rho(a, \phi, z) = \rho_0(a, \phi) F_{\rho}(z/H), \qquad (C.1)$$

$$p(a, \phi, z) = p_0(a, \phi) F_{p}(z/H), \qquad (C.2)$$

$$\frac{\partial}{\partial z} F_{\rho}(z/H) = -\frac{z}{H^2} F_{\rho}(z/H), \qquad (C.3)$$

which are fully satisfied when the disc is assumed to be locally isothermal as we did in Eqs. (23) and (24).

## Appendix D: Channel maps and moment 1 maps

For completeness, we present in Fig. D.1 a selection of channel maps obtained from the RT simulation for both the $i = 0^{\circ}$ and $i = 30^{\circ}$. Both channel maps are obtained smoothing the ‘infinitesimal’ resolution synthetic image from MCFOST using a circular Gaussian beam with $\Delta x = 0.05''$ and $\Delta v = 0.05 \text{ km s}^{-1}$.
<figure><img src="image_17.png" /><figcaption>**Fig. D.1.** Selection of channel maps from the RT model for the inclination *i* = 0° (top panels) and *i* = 30° (bottom panels), after smoothing the synthetic image from MCFOST using a circular Gaussian beam with Δ*x* = 0.05" and Δ*v* = 0.05 kms⁻¹. The velocity of each channel is marked in the top left corner of each panel. The grey circle in the bottom left corner represents the corresponding beam size.</figcaption></figure>