© The Authors 2024

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; 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 ϖ(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 〈c_s〉(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), ϖ(a), M_a(a), and 〈c_s〉(a), predicts the following: the radial motion υ_R, the orbital motion υ_θ, the vertical motion υ_z, and disc morphology, in the form of H/R and Σ(a, ϕ) – 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 ϖ(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 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, hydro-dynamic 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, ϕ, 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.

Fig. 1 Example of the morphology of a disc composed by a set of nested ellipses with both eccentricity e(a) and pericentre phase ϖ(a) varying for different values of the semi-major axis.

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, θ, z) when dealing with eccentric disc geometry is very inconvenient: a single value of R spans through multiple elliptical orbits for different values of θ. As mentioned above, since the orbital structure of an eccentric disc is defined by e(a) and ϖ(a), a much better choice is an elliptical coordinate system using (a, ϕ), where a is the semimajor axis of the ellipse, while the azimuthal coordinate ϕ is related to the true anomaly by f = ϕ − ϖ (a). Since we assume that the focus of the ellipses is fixed at the origin, ϕ is coincident with the polar azimuthal coordinate θ.

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, ϕ, z) → (x, y, z) coordinates are

$$x(a,\varphi )=R(a,\varphi )\text{\hspace{0.17em}cos}[\varphi ]$$(1) $$y(a,\varphi )=R(a,\varphi )\text{\hspace{0.17em}sin}[\varphi ],$$(2) $$z=z,$$(3)

where R(a, ϕ) is

$$R(a,\varphi )=\frac{a\left[1-e{(a)}^2\right]}{1+e(a)\text{\hspace{0.17em}cos}[\varphi -\varpi (a)]}.$$(4)

The Jacobian matrix from this transformation is

$$\mathcal{𝒥}=\frac{\partial (x,y,z)}{\partial (a,\varphi ,z)}.$$(5)

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

$$\begin{array}{l}J(a,\varphi )=R{(a,\varphi )}^3\frac{1-e(a)\left[e(a)+2a{e}_a(a)\right]}{a^2{\left[1-e{(a)}^2\right]}^2}\hfill \\ \cdot \{1-q\cos [\varphi -\varpi (a)-\alpha ]\},\hfill \end{array}$$(6)

where the angle α 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ϖ_a by

$$q\cos (\alpha )=\frac{\left[1+e{(a)}^2\right]a{e}_a(a)-\left[1-e{(a)}^2\right]e(a)}{1-e(a)\left[e(a)+2a{e}_a(a)\right]},$$(7) $$q\sin (\alpha )=\frac{e(a)\left[1-e{(a)}^2\right]a{\varpi }_a(a)}{1-e(a)\left[e(a)+2a{e}_a(a)\right]}.$$(8)

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

$$\tan (\alpha )=\frac{ae(a)\left[1-e{(a)}^2\right]{\varpi }_a(a)}{\left[1+e{(a)}^2\right]a{e}_a(a)-\left[1-e{(a)}^2\right]e(a)}.$$(9)

We note that since the coordinate z does not depend on a and 0, the expression of J is the same both when we deal with a 2D disc using the (a, ϕ) coordinates only, and for a full 3D disc (a, ϕ, 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

$$\text{d}V=J\text{d}a\text{d}\varphi \text{d}z,$$(10) $$\text{d}A=J\text{d}a\text{d}\varphi .$$(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 ≥ 1. Similarly α represents the true anomaly, that is the angle from the pericentre longitude ϖ, where the maximum orbital compression occurs.

2.2 Orbital velocity field: Keplerian motion

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

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

$$v_R(a,\varphi )=\dot{R}=a{\text{Ω}}_0\frac{e(a)}{\sqrt{1-e^2(a)}}\sin [\varphi -\varpi (a)],$$(12) $$v_{\theta ,\text{K}}(a,\varphi )=R\dot{\theta }=a{\text{Ω}}_0\frac{1+e(a)\text{\hspace{0.17em}cos}[\varphi -\varpi (a)]}{\sqrt{1-e^2(a)}},$$(13)

where the dot notation indicates total time derivative, d/dt, and Ω₀ is the mean motion

$${\text{Ω}}_0=\sqrt{\frac{GM}{a^3}}\text{.}$$(14)

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

$$\text{Ω}(a,\varphi )\equiv \dot{\theta }={\text{Ω}}_0\frac{{\{1+e(a)\text{\hspace{0.17em}cos}[\varphi -\varpi (a)]\}}^2}{{\left[1-e^2(a)\right]}^{3/2}}.$$(15)

The standard analytical treatment of differential operators in the eccentric coordinate system requires the definition of contravariant velocities as υ^a ≡ ȧ and ${\upsilon }^{\varphi }\equiv \dot{\varphi }$. With these definitions in mind we can define υ^a and υ^ϕ 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,$$(16)

while, given the coincidence between ϕ and θ, υ^ϕ is

$$v^{\varphi }=\text{Ω}.$$(17)

We remark that υ^ϕ is a contravariant variable, with the dimension of an angular velocity. While υ_θ 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, υ_θ and υ^ϕ 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}\left(J{v}^a\rho \right)-\frac{1}{J}\frac{\partial }{\partial \varphi }\left(J{v}^{\varphi }\rho \right)-\frac{\partial }{\partial z}\left(v^z\rho \right),$$(18)

obtained from the expression of the divergence of a vector field u in a generic coordinate system {xⁱ},

$${\nabla }_i\left(u^i\right)=\frac{1}{J}\frac{\partial }{\partial x_i}J{u}^i.$$(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^{\varphi }\frac{\partial v^z}{\partial \varphi }+v^z\frac{\partial v^z}{\partial z}=-\frac{GM}{R^2}\frac{z}{R}-\frac{1}{\rho }\frac{\partial p}{\partial z}.$$(20)

We now look for steady solutions and substitute υ^a = 0 and υ^ϕ = Ω(a, ϕ) – as we prescribed in Sect. 2.2. After these operations, Eqs. (18) and (20) reduce to

$$\frac{\partial }{\partial z}\left(v^z\rho \right)=-\frac{1}{J}\frac{\partial }{\partial \varphi }(J\text{Ω}\rho ),$$(21) $$\text{Ω}\frac{\partial v^z}{\partial \varphi }+v^2\frac{\partial v^z}{\partial z}=-\frac{GM}{R^2}\frac{z}{R}-\frac{1}{\rho }\frac{\partial p}{\partial z}.$$(22)

We assume the disc to be locally isothermal so that ρ and p have Gaussian form (e.g. Lodato 2008 for a review of classical disc dynamics)

$$\rho =\frac{{\rho }_0}{\sqrt{2\pi }}\text{\hspace{0.17em}exp\hspace{0.17em}}\left(-\frac{z^2}{2H^2}\right),$$(23) $$p=\frac{p_0}{{\rho }_0}\rho ,$$(24)

where p₀ and ρ₀ represent the density value in the midplane. We note that $p_0/{\rho }_0=c_{\text{s}}^2$ is the squared local isothermal sound speed at the disc midplane; H(a, ϕ), is the standard deviation of the vertical density distribution, here representing the disc local vertical scale-height¹. The form of the final equations does not change for different assumptions on ρ 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 υ^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{\text{d}}{\text{d}t}\frac{z}{H}=\frac{\partial }{\partial t}\frac{z}{H}+\text{Ω}\frac{\partial }{\partial \varphi }\frac{z}{H}+v^z\frac{\partial }{\partial z}\frac{z}{H}=0.$$(25)

This ansatz is physically the homogeneous vertical expansion of the gas column. When the disc is stationary (∂/∂t = 0) such a condition is satisfied if υ^z has the form²

$$v^z=\text{Ω}\frac{\partial H}{\partial \varphi }\frac{z}{H}$$(27)

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

$$\text{Ω}\frac{\partial {\rho }_0}{\partial \varphi }=-{\rho }_0\left[\frac{1}{J}\frac{\partial }{\partial \varphi }(J\text{Ω})+\frac{\text{Ω}}{H}\frac{\partial H}{\partial \varphi }\right],$$(28) $${\text{Ω}}^2\frac{{\partial }^{2}H}{\partial {\varphi }^2}=-\text{Ω}\frac{\partial H}{\partial \varphi }\frac{\partial \text{Ω}}{\partial \varphi }-\frac{GM}{R^3}H+\frac{c_{\text{s}}^2}{H}$$(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 ϖ = 0, and the resulting υ^z from Eq. (27). The value of ρ₀ solving Eq. (28) can then be easily obtained by noting that Eq. (28) can be rewritten as

$$\frac{1}{J}\frac{\partial }{\partial \varphi }\left(J\text{Ω}H{\rho }_0\right)=0,$$(30)

implying that JΩHρ₀ = ℱ (a) depends on the semi-major axis only and is a conserved mass flux along the orbit. Further considerations on ℱ(a) will be used in the next section to discuss the disc morphology.

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_apo/H_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 ≈ 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.

Fig. 2 Vertical height H and vertical velocity υz as a function of ϕ Top panel: H/H0 as a function of ϕ, where H0 = cs/Ω0, obtained solving Eq. (29). Bottom panel: υz/cs as a function of ϕ, obtained from Eq. (27) after solving Eq. (29), calculated at an elevation z = H above the midplane − cs is the most natural normalisation2 for υz. Solution to Eq. (29) was obtained solving for periodic solutions across one orbit, i.e. H(a, 0) = H(a, 2π), using a shooting method for solving the boundary value problem.

2.4 Morphology

Eccentric discs might present steady over-densities, caused by orbital compression (i.e. e_a ≠ 0 and ϖ_a ≠ 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 Σ(a, ϕ).

We obtained in Eq. (30) conservation of mass flux along elliptical orbits. We note that $\text{Σ}={\int }_z\rho \text{d}z=H{\rho }_0$, thus we can rewrite ℱ(a) as³

$$ℱ(a)=J(a,\varphi )\text{Ω}(a,\varphi )\text{Σ}(a,\varphi ).$$(31)

We note that integrating the r.h.s. of Eq. (31) over an orbital time t_orb implies

$${\displaystyle {\int }_0^{t_{\text{ort}}}J(a,\varphi )\text{Σ}(a,\varphi )\text{Ω}(a,\varphi )\text{d}t}={\displaystyle {\int }_0^{2\pi }J(a,\varphi )\text{Σ}(a,\varphi )\text{d}\varphi }=M_a(a)$$(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_ada is the mass of an infinitesimal eccentric annulus of width da.

Since ℱ(a) does not depend on ϕ, we have t_orbℱ(a) = M_a(a), which can be rewritten as ℱ(a) = M_a(a)Ω₀/2π, where we recall Ω₀ is the mean motion (Eq. (14)). In the light of these considerations, Eq. (31) can be rewritten as (Ogilvie 2001; Statler 2001)

$$\text{Σ}(a,\varphi )=\frac{M_a(a){\text{Ω}}_0}{2\pi J(a,\varphi )\text{Ω}(a,\varphi )}.$$(33)

If e_a(a) = 0 and ϖ_a(a) = 0 – that is, the disc has vanishing eccentricity gradient and twist – the expression of Σ(a, ϕ) simplifies to an expression that does not depend on ϕ.

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 (ϖ_a(a) + 0), the overdensity is located at angle 0 satisfying ϕ − ϖ(a) = α, where α is presented in Eq. (9).

It is important to note that any e(a) = const, ϖ(a) = 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 δυ_θ,P required for a circular accretion disc, that is, (e.g. Lodato 2008 for a review of classical disc dynamics):

$$v_{\theta }=v_{\theta ,\text{K}}\sqrt{1+\eta },$$(34) $$\eta ={\left({\text{Ω}}_{0}a\right)}^{-2}\left(\frac{a}{{\text{Σ}}_0}\frac{\text{d}P}{\text{d}a}\right)$$(35) $$P={\langle c_{\text{s}}\rangle }^2{\text{Σ}}_0,$$(36)

where υ_θ,K is the Keplerian velocity (Eq. (13)), Σ₀ constitutes an averaged value along the orbit defined as

$${\text{Σ}}_0=\frac{M_a}{2\pi a},$$(37)

and 〈c_s〉 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 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⁴, above a certain mass-ratio threshold (q ≳ 10⁻³, 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 encouter between the two stars (Cuello et al. 2019). Furthermore numerous numerical works on circumbinary discs, 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 ≲ 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 υ_θ,_circ = const and υ_R = 0 along an orbit, an eccentric disc has υ_θ and υ_R varying around the orbit. The azimuthal velocity υ_θ (Eqs. (13), (34) for the pressure corrected) reaches its maximum at the phase of pericentre ϖ and minimum at the phase of the orbit apocentre. Similarly, υ_R (Eq. (12)) oscillates around υ_R = 0, reaching its absoulute maxima (|υ_R|) at true anomaly f = π/2 and f = 3π/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 ≳ 6 for e ≳ 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 ∝ H.

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

The amplitude of these vertical υ_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 υ_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, υ_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 peri-centre phase variations across the disc result in overdensities and depletions of the disc surface density, Σ (Eq. (33)). However, one should always keep in mind that a change in the disc thickness might affect the volume density, ρ, at the midplane without affecting the surface density: for example, a disc with e = const and ϖ = const across its entire domain, produces no variations in the surface density morphology, that is, Σ(θ) = const. But, since the volume density at the midplane is ρ = Σ/H, ρ changes as the inverse of H – that is, ρ is maximum at pericentre and minimum at apocentre even when Σ 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 υ_R, υ_θ, υ_z and to cartesian velocities as υ_x, υ_y, υ_z.

We approached the problem as follows (Sect. 3.3 for more details): we first extracted from the simulations the profiles of e(a), ϖ(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 (υ_x, υ_y, υ_z). We finally compared the predictions of the theoretical model with the simulations.

We note that quantities generated using the theoretical model are indicated with a subscript ‘th’ (e.g. x_th) while those calculated directly from the simulation with subscript ‘sim’ (e.g. x_sim).

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 M2/M1 = 0.1, H/R = 0.05) in Ragusa et al. (2020) after t = 500tbin. The dump features a cavity with an eccentricity of ecav ≈ 0.2 and a cavity with a semi-major axis acav ≈ 2.7abin.

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₂/M₁ > 0.05 presented in Ragusa et al. (2020) show a rapid growth of the disc eccentricity after t ≈ 400–500t_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₂/M₁ = 0.1). In particular, its t = 500t_bin dump (shown in Fig. 3), where t_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_cav ≈ 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₂/M₁ > 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⁶ particles, and consisted of a circular, live binary with mass-ratio M₂/M₁ = 0.1 surrounded by a circumbinary disc with initial mass M_d = 5 · 10⁻³ M_bin (M_d = 4.8 · 10⁻³ M_bin after t = 500 t_bin), where M_bin = M₁ + M₂. The equation of state was chosen to be locally isothermal – with radial temperature profile producing c_s ∝ R^−0.25 and H/R = 0.05 at R_in = 2a_bin. A circular cavity with R = R_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) α_ss ≈ 5 × 10⁻³ (obtained using α_AV = 0.2, β = 2 Price et al. 2018a)⁵ at the beginning of the simulation.

3.2 Obtaining maps of v_R, v_θ, v_z, H, Σ from the simulation

For the comparison with the analytical model, we start generating projected x − y maps of υ_x,_sim, υ_y,sim, υ_z,sim, H_sim, and Σ_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 υ_x,sim = 〈υ_x〉_z and υ_y,sim = 〈υ_y〉_z, where 〈·〉_z indicates density weighted integration along the vertical direction and υ_x and υ_y the velocities of individual SPH particles in the density weighted integration. Then, υ_R,sim and υ_θ,sim are obtained projecting υ_2D = {υ_x,sim, υ_y,sim} on radial and azimuthal unit vectors, respectively.

The map for H_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_sim as the second moment of the particles vertical density distribution. Again, here z² indicates that the value of z² 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 σ = H_sim.

Similarly, the map for υ_z,sim is obtained performing a density weighted integration of the quantity ${\upsilon }_{\text{z,sim}}=\sqrt{\pi /2}{\langle {\upsilon }_z\text{sign}\left(z\right)\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 υ_z is anti-symmetric with respect to the midplane, the quantity υ_z sign(z) has the nice property of providing a υ_z,sim maintaining the sign of the velocities of the z > 0 half of the disc during the average across the mid-plane.| In contrast, ${\langle |{\upsilon }_z^{}|\rangle }_z\text{\hspace{0.17em}\hspace{0.17em}or\hspace{0.17em}\hspace{0.17em}}{\langle {\upsilon }_z^2\rangle }_z$ would produce a map that is everywhere positive, while the vertical average of 〈υ_z〉_z, would produce υ_z,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}\left(z\right)\rangle }_z=H\sqrt{2/\pi }$. Under the assumption that υ_z ∝ z/H, as derived in Eq. (27), the introduction of the pre-factor makes sure that υ_z,sim traces the vertical velocity of particles at z ~ H_sim above the midplane. We anticipate that the results presented in the following sections confirm the validity of this assumption.

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), ϖ(a), their derivatives, e_a and ϖ_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 (υ_x, υ_y, υ_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), ϖ(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), ϖ_a(a) with respect to the semi-major axis. We set the filter window to be ω_sg = 50 profile points (corresponding to a Δa ~ 1.6a_bin) and polynomial order p_sg = 2. We note that filtering is a very important step, since noisy datasets of e(a), ϖ(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 υ_R,th and υ_θ,th. We then solve Eq. (29) to obtain the vertical scale height H_th and use it for calculating υ_z,th from Eq. (27). We finally use Eq. (33) to calculate the surface density morphology Σ_th. Results are plotted as 2D maps where the region with a < a_cav has been excised – where a_cav = 2.7a_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 υ_R, υ_θ, and υ_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 υ_R and υ_θ. In the same figure, the residuals obtained subtracting υ_R,circ = 0 and ${\upsilon }_{\theta ,\text{circ}}=\sqrt{G{M}_{\text{bin}}/R}$. We normalise the residuals for υ_θ with υ_θ,circ, that is, the circular Keplerian velocity at the x − y pixel radius; we normalise υ_R using e_cavυ_θ,circ where e_cav = 0.2, which represent the maximum vR 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 7 shows a comparison between the vertical velocity field in the simulation υ_z,sim against the one obtained from the theoretical model υ_z,th. In Fig. 7, we also show the residuals from the subtraction of the two, and one where υ_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 υ_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 Σ for both the simulation and the theoretical model, obtained from Eq. (33).

Fig. 4 Profiles of e(a) (top panel), ϖ(a) (middle panel), and Ma (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.

Fig. 5 Maps of υR and υθ. Top panels: comparison of the velocity fields υR and υθ in the form of residuals between the simulation and the theoretical predictions in Eqs. (12) and (34): (υR,sim − υR,th)/υθ,circ (top left panel) and (υθ,sim − υθ,th)/υθ,circ (top right panel), respectively. Bottom panels: similar to top panels but subtracting a circular Keplerian velocity profile (υθ,circ) from the υR,sim and υθ,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; (υR,sim − υR,circ)/υθ,circ (bottom left panel), where obviously υR,circ = 0, and (υθ,sim − υθ,circ)/υθ,circ (bottom right panel). In the bottom right panel it can be clearly seen that the simulation azimuthal velocity υθ,sim at large radii is sub-Keplerian, as υθ,circ does not account for the pressure support term.

Fig. 6 Maps of disc aspect-ratio H/R. Left panel: disc Hsim/R from simulation, calculated as Hsim/R = 〈z2(R)〉1/2/R; we note that 〈z2(R)〉 = Hsim by definition, with it being the vertical second moment of the density distribution. Right panel: disc Hth/R from solving Eq. (29) using e(a), ϖ(a), and cs(R) from the simulation as input.

Fig. 7 Maps of disc vertical velocity υz. Top panels: υz,sim/cs velocity field from the simulation (top left panel) compared to the theoretical model υz,th/cs obtained solving Eq. (29) (top right panel). Bottom panels: residuals of the simulation against the theoretical model (υz,sim − υz,th)/cs (bottom left panel), i.e. the top left panel minus the top right panel of this figure, and with respect to a circular disc (υz,sim − υz,circ)/cs (bottom right panel), i.e. υz,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.

Fig. 8 Σ density field from the simulation (left panel) and relative analytical model (right panel) 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 υ_R field show a maximum residual value of |Δυ_R|_max/(e_cavυ_θ,circ) ~ 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 |Δυ_R|_avg ≲ 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 υ_R, that is, in a circular disc, produces |Δυ_R,circ|_max/(e_cavυ_θ,circ) ~ 100%. The comparison with the circular case, shows as expected an increase of the radial velocity reaching its maximum in absolute value |υ_R| at true anomaly f = π/2 and f = 3π/2.

The υ_θ field similarly has a maximum residual value of Δυ_θ/υ_θ,circ ~ 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 |Δυ_θ,circ|/υ_θ,circ ~ 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 precence of a ‘breathing’ mode that causes the disc to compress towards the midplane at peri-centre 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 ~ 3, as visible in Fig. 2 for eccentricity e ≈ e_cav = 0.2.

The vertical velocity field υ_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 υ_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 υ_z field in the numerical simulation shows a small level of asymmetry with respect to the pericentre longitude: that is, the positive υ_z lobe looks slightly different in shape with respect to the negative one. Finally, the simulation υ_z map appears to show a slight shift with respect to the pericentre and apocentre of the location where υ_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 υ_z vanishes, consistently with the one observed in the simulation. However, the new theoretical prescription for υ_z fails in any case to reproduce the m = 1 spiral pattern visible in the residuals of υ_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 υ_z, which is beyond the scope of the current paper.

In conclusion, we cannot identify the precise origin of the discrepancies in the predicted υ_z and the simulated one. Such a spiral feature in the υ_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 υ_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 ≳ 0.1 mass-ratio regime. To our knowledge, such spiral features in the vertical υ_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 υ_z 3D field produces υ_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 ¹²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 ¹²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₁ = 0.91 M_⊙ and M₂ = 0.09 M_⊙ and to radiate as two black bodies with temperature T₁ = 4200 K and T₂ = 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 ${\upsilon }_{\text{sc}}=29.7/\sqrt{15}$ to have the velocities in km s⁻¹. The total gas mass is 5 M_J, 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⁶ ρ¹²_CO/ρ_H = 5 × 10⁻⁵.

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 ∝ s^7/2.

We use Nγ = 10⁷ 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_γ = 10⁷ 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 ‘infinite’ spatial resolution (set by the hydrodynamical simulation), while spectral resolution bins in the cube are Δυ = 0.004 km s⁻¹ for the face-on case and Δυ = 0.01 km s⁻¹ 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 Δs = 0.05 arcsec and the velocity with Δυ = 0.05 km s⁻¹ (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 υ_z ≠ 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⁻¹ 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 υ_z,max ≈ 0.25 km s⁻¹, 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 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 𝒮 = {x, y, z} defining the 3D disc upper emission surface, and velocities υ = {υ_x, υ_y, υ_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 τ = 1 surface above the midplane H_τ=1 output from the Monte Carlo radiative transfer simulation. The value of τ 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 ={\displaystyle {\int }_{+\infty }^{z(x,y)}\rho (x,y,z)\underset{v}{\text{\hspace{0.17em}max}}\left[k_v(x,y,z)\right]\text{dz}}=1$$(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_τ=1/H_sim for a comparison with the vertical scale height of the simulation in Fig. 11. The average value 〈H_τ=1 /H_sim〉 ≈ 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_τ=1 is not visible, as that would be the case only if υ_z = const throughout the entire gas column. Indeed, τ defined in Eq. (38) should not be confused with τ_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 τ). Since the opacity k_v of the material changes with the material velocity, τ_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_τ=1.

The disc has a top and a bottom surface, thus for our 3D eccentric and circular model we define them to be 𝒮_top = {x, y, 2.5H_th} for the top surface and 𝒮_bottom = {x, y, −2.5H_th}.

We define the velocity field on the disc surfaces as υ_top = {υ_x, υ_y, υ_z} for the top layer and υ_bottom = {υ_x, υ_y, −υ_z} for the bottom layer, where υ_x and υ_y are obtained projecting υ_R and υ_θ from Eqs. (12) and (34) on to the x and y unit vectors, while υ_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 υ_circ on to Cartesian unit vectors.

Similarly, after rescaling the numerical simulation as described in Sect. 5.1 (a_bin = 15 au, M_bin = 1M_⊙, d = 130 pc, that is, x_sc = 15, ${\upsilon }_{\text{sc}}=29.7/\sqrt{15}$), we define the upper surface in the simulation to be 𝒮_sim = x_sc{x, y, 2.5H_sim} and υ = υ_sc{υ_x,sim, υ_y,sim, 2.5υ_z,sim}, where, we note, we applied an additional scaling of a factor 2.5 to the original z_sim and υ_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 υ_z ∝ z/H_th.

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 Δυ = 0.05 km s−1. 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 υz 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.

Fig. 10 Brightness map of the 0.72 km s−1 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−1 for vertical velocities, a bright contribution in higher velocities channels is also expected to be observable.

Fig. 11 Map of the ratio between Hτ=1, obtained from radiative transfer, tracing the upper emission surface and Hsim, representing the simulation scale height. The ratio Hτ=1/Hsim ≈ 2–3 across the disc. We consider Hτ=1 ≈ 2.5Hsim a conservative estimate of the location of the upper emission surface – ‘conservative’ since the higher in the disc atmosphere, the larger the vertical velocity.

5.4 Morphology of projected velocity residuals in circumbinary discs

Consider an observer located at z = + ∞, and define the projected velocity υ_proj = −υ_z (following the convention that redshifted velocities are positive). At this point, the observer, looking down on the x − y plane observing υ_proj will see the kinematics and geometry of a face-on disc. We do this for both the eccentric and 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 𝒮 and υ), 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 argument 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 simulation, the other two for the eccentric and circular disc theoretical models.

For the analysis described in this section, we consider only 𝒮_top, 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 υ_proj 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 circularmodel.

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 Δυ > 0.3 km s⁻¹. 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).

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τ=1 = 2.5H (as discussed in Sect. 5.3), and rescaling the simulations to have a binary with abin = 15 au, Mbin = 1 M⊙, 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 υz 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 υz,circ = 0 (as the circular case does not have vertical motion) from υz, 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 υz ~ 0.4 km s−1 in the upper emission layer.

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 ϖ(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 υ_θ υ_R, along with the surface density Σ 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 υ_z ∝ 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 υ_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 Δυ_z ≈ 0.15 km s⁻¹.

• (iv)

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

• (v)

  Both i = 0° and i = 30° 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 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°), the M9 map reveals the vertical velocity reaches a maximum of υ_z,max ~ 0.25 km s⁻¹. However, channel maps show significant signal at higher velocities from the upper layers of the disc, making υ_z eccentric signatures potentially detectable in face-on discs even for eccentricities e ≪ 0.2 (e.g. υ_z,max ~ 0.4 km s⁻¹ on the τ = 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).

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 ϖ(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 ϖ(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 ω(a), producing a twisted eccentricity pattern and oscillations of the value of e(a, t) and ϖ(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 ϖ(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){\text{e}}^{\text{i}\omega (a,t)}.$$(A.1)

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

Both linear (|E| ≪ 1 and |aE_a| ≪ 1, that is q ≪ 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){\text{e}}^{\text{i}\omega t}.$$(A.2)

In absence of dissipation, this expression implies that the disc have a stationary e(a) profile that coherently precesses at the rate ω_i throughout the whole disc. If dissipation is present the complex phase will also include a twist ϖ(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)={\displaystyle \sum _i}{c}_i{E}_i(a,t),$$(A.3)

where c_i are complex coefficients. Out of the linear regime, nonlinear 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):

$$\begin{array}{l}\text{Σ}{R}^2\text{Ω}\frac{\partial E}{\partial t}=\frac{\text{i}}{R}\frac{\partial }{\partial R}\left(\frac{1}{2}\text{Σ}{c}_{\text{s}}^2{R}^3\frac{\partial E}{\partial R}\right)+\frac{\text{i}r}{2}\frac{\text{d}}{\text{d}R}\left(\text{Σ}{c}_{\text{s}}^2\right)E\hfill \\ -\frac{\text{i}}{2R}\frac{\partial }{\partial R}\left(\text{Σ}\frac{\text{d}{c}_{\text{s}}^2}{\text{d}r}{R}^{3}E\right)+\frac{3\text{i}}{2R}\text{Σ}\frac{\text{d}}{\text{d}R}\left(c_{\text{s}}^2{R}^2\right)E.\hfill \end{array}$$(A.4)

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

$${\omega }_i\mathcal{𝒜}(R)e_i(R)=\frac{\text{d}}{\text{d}R}\left(ℬ(R)\frac{\text{d}{e}_i}{\text{d}R}(R)\right)+C(R)e_i(R).$$(A.5)

where 𝒜(R), ℬ(R) and 𝒞(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{{\text{d}}^2{\text{Ξ}}_i(x)}{\text{d}{x}^2}+V{\text{Ξ}}_i(x)={ℰ}_i{\text{Ξ}}_i(x),$$(A.6)

where x and a are related by:

$$\frac{\text{d}x}{\text{d}R}={\left[\frac{{\omega }_0\mathcal{𝒜}(R)}{ℬ(R)}\right]}^{1/2}.$$(A.7)

Functions e_i(R) and Ξ_i(x) are related by:

$${\text{Ξ}}_i(x)=e_i(R){[\mathcal{𝒜}(R)ℬ(R)]}^{1/4}.$$(A.8)

The effective energy eigenvalue E_i is given by:

$${ℰ}_i=-\frac{{\omega }_i}{{\omega }_0},$$(A.9)

where ω₀ relates with the disc thickness (H/R)₀ and mean motion Ω₀ at the inner disc edge as follows:

$${\omega }_0=\frac{1}{2}{\left(\frac{H}{R}\right)}_0{\text{Ω}}_0.$$(A.10)

The effective potential V(x) reads:

$$V=-\frac{C(R)}{{\omega }_0\mathcal{𝒜}(R)}+{(\mathcal{𝒜}(R)ℬ(R))}^{-1/4}\frac{{\text{d}}^2}{\text{d}{x}^2}{[\mathcal{𝒜}(R)ℬ(R)]}^{1/4}.$$(A.11)

As a consequence, in analogy to quantum mechanics, there exist a countable infinity of eccentric eigenmodes E_i(R, t). Pro-grade (ω_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 ω_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 Σ ∝ R^−σ and soundspeed $c_{\text{s}}^2\propto R^{-q}$. Under this assumption, it can be shown that E(R, f) can be expanded in series as (Goodchild & Ogilvie 2006; Ogilvie 2008)

$$E(R,t)={\displaystyle \sum _i}{a}_i(t){\text{Ψ}}_i(R),$$(A.12)

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

$${\text{Ψ}}_i(R)=R^{\sigma /2+q-1}{Z}_v\left({\left(R/R_i\right)}^{\beta }\right),$$(A.13)

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

$$Z_v(x)=A{J}_v(x)+B{Y}_v(x)$$(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, σ, and ß by v² = β⁻²(2q + qσ + q² + σ²/4 – 5). A, B, and {R_i} are set by the (zero gradient) boundary conditions, ${\partial }_R{\Psi }_i{|}_{a_{in}}={\partial }_R{\Psi }_i{|}_{a_{out}}=0$ and an appropriate normalisation condition.

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

$$E_0(R,t)=a_0(t){\text{Ψ}}_0(R).$$(A.15)

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

$$e(R)=\mathcal{𝒜}{R}^{-\gamma }\left\{J_v\left({\left[R/R_0\right]}^{\beta }\right)+C{Y}_v\left({\left[R/R_0\right]}^{\beta }\right)\right\},$$(A.16)

where 𝒞, ß, γ, and R₀ are constants whose value depend on the powers of the profiles of Σ and c_s, as discussed above, as well as other disc properties and other assumptions (e.g. boundary conditions), while 𝒜 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{𝒜}{R}^{-g}\text{\hspace{0.17em}exp\hspace{0.17em}}\left[{\left(R/r_0\right)}^b\right].$$(A.17)

Similarly to Eq. (A.16), b, g, and r₀ are constants whose value depend on the powers of the profiles of Σ and c_s as well as other disc properties, while 𝒜 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^{\varphi }\frac{\partial v^z}{\partial \varphi }+v^z\frac{\partial v^z}{\partial z}=-\frac{GM}{R^2}\frac{z}{R}-\frac{1}{\rho }\frac{\partial }{\partial z}\left(p-T^{zz}\right).$$(B.1)

where T^zz is the bulk viscous stress tensor

$$T^{zz}={\mu }_b\left(\frac{1}{J}\frac{\partial }{\partial \varphi }(J\text{Ω})+\frac{\partial v^z}{\partial z}\right),$$(B.2)

parametrised by a μ_b given by (Ogilvie 2001; Ogilvie & Barker 2014; Lynch & Ogilvie 2021):

$${\mu }_b={\alpha }_{b}p{\text{Ω}}_0^{-1}.$$(B.3)

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

$$\begin{array}{l}{\text{Ω}}^2\frac{{\partial }^{2}H}{\partial {\varphi }^2}=-\text{Ω}\frac{\partial H}{\partial \varphi }\frac{\partial \text{Ω}}{\partial \varphi }-\frac{GM}{R^3}H+\hfill \\ +\frac{p_0}{{\rho }_{0}H}\left[1-\frac{{\alpha }_b}{{\text{Ω}}_0}\left(\frac{1}{J}\frac{\partial }{\partial \varphi }(J\text{Ω})+\frac{\text{Ω}}{H}\frac{\partial H}{\partial \varphi }\right)\right].\hfill \end{array}$$(B.4)

We note that in order to produce the same plots for α_b ≠ 0 one needs specify q and a 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 α_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 ρ(a,φ,z) and p(a,φ,z)

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

$$\rho (a,\varphi ,z)={\rho }_0(a,\varphi )F_{\rho }(z/H)$$(C.1) $$p(a,\varphi ,z)=p_0(a,\varphi )F_p(z/H),$$(C.2) $$\frac{\partial }{\partial z}{F}_p(z/H)=-\frac{z}{H^2}{F}_{\rho }(z/H),$$(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° and i = 30°. Both channel maps are obtained smoothing the ‘infinite’ resolution synthetic image from MCFOST using a circular Gaussian beam with Δx = 0.05″ and Δv = 0.05 kms⁻¹.

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−1. 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.

References

All Figures

	Fig. 1 Example of the morphology of a disc composed by a set of nested ellipses with both eccentricity e(a) and pericentre phase ϖ(a) varying for different values of the semi-major axis.
In the text

Fig. 2 Vertical height H and vertical velocity υz as a function of ϕ Top panel: H/H0 as a function of ϕ, where H0 = cs/Ω0, obtained solving Eq. (29). Bottom panel: υz/cs as a function of ϕ, obtained from Eq. (27) after solving Eq. (29), calculated at an elevation z = H above the midplane − cs is the most natural normalisation2 for υz. Solution to Eq. (29) was obtained solving for periodic solutions across one orbit, i.e. H(a, 0) = H(a, 2π), using a shooting method for solving the boundary value problem.

In the text

	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 M2/M1 = 0.1, H/R = 0.05) in Ragusa et al. (2020) after t = 500tbin. The dump features a cavity with an eccentricity of ecav ≈ 0.2 and a cavity with a semi-major axis acav ≈ 2.7abin.
In the text

Fig. 4 Profiles of e(a) (top panel), ϖ(a) (middle panel), and Ma (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.

In the text

Fig. 5 Maps of υR and υθ. Top panels: comparison of the velocity fields υR and υθ in the form of residuals between the simulation and the theoretical predictions in Eqs. (12) and (34): (υR,sim − υR,th)/υθ,circ (top left panel) and (υθ,sim − υθ,th)/υθ,circ (top right panel), respectively. Bottom panels: similar to top panels but subtracting a circular Keplerian velocity profile (υθ,circ) from the υR,sim and υθ,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; (υR,sim − υR,circ)/υθ,circ (bottom left panel), where obviously υR,circ = 0, and (υθ,sim − υθ,circ)/υθ,circ (bottom right panel). In the bottom right panel it can be clearly seen that the simulation azimuthal velocity υθ,sim at large radii is sub-Keplerian, as υθ,circ does not account for the pressure support term.

In the text

	Fig. 6 Maps of disc aspect-ratio H/R. Left panel: disc Hsim/R from simulation, calculated as Hsim/R = 〈z2(R)〉1/2/R; we note that 〈z2(R)〉 = Hsim by definition, with it being the vertical second moment of the density distribution. Right panel: disc Hth/R from solving Eq. (29) using e(a), ϖ(a), and cs(R) from the simulation as input.
In the text

Fig. 7 Maps of disc vertical velocity υz. Top panels: υz,sim/cs velocity field from the simulation (top left panel) compared to the theoretical model υz,th/cs obtained solving Eq. (29) (top right panel). Bottom panels: residuals of the simulation against the theoretical model (υz,sim − υz,th)/cs (bottom left panel), i.e. the top left panel minus the top right panel of this figure, and with respect to a circular disc (υz,sim − υz,circ)/cs (bottom right panel), i.e. υz,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.

In the text

	Fig. 8 Σ density field from the simulation (left panel) and relative analytical model (right panel) from Eq. (33).
In the text

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 Δυ = 0.05 km s−1. 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 υz 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.

In the text

	Fig. 10 Brightness map of the 0.72 km s−1 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−1 for vertical velocities, a bright contribution in higher velocities channels is also expected to be observable.
In the text

	Fig. 11 Map of the ratio between Hτ=1, obtained from radiative transfer, tracing the upper emission surface and Hsim, representing the simulation scale height. The ratio Hτ=1/Hsim ≈ 2–3 across the disc. We consider Hτ=1 ≈ 2.5Hsim a conservative estimate of the location of the upper emission surface – ‘conservative’ since the higher in the disc atmosphere, the larger the vertical velocity.
In the text

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τ=1 = 2.5H (as discussed in Sect. 5.3), and rescaling the simulations to have a binary with abin = 15 au, Mbin = 1 M⊙, 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 υz 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 υz,circ = 0 (as the circular case does not have vertical motion) from υz, 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 υz ~ 0.4 km s−1 in the upper emission layer.

In the text

	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−1. 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.
In the text