Exciting spiral arms in protoplanetary discs from flybys

ABSTRACT

Spiral arms are observed in numerous protoplanetary discs. These spiral arms can be excited by companions, either on bound or unbound orbits. We simulate a scenario where an unbound perturber, i.e. a flyby, excites spiral arms during a periastron passage. We run three-dimensional hydrodynamical simulations of a parabolic flyby encountering a gaseous protoplanetary disc. The perturber mass ranges from |$10\, \rm M_J$| to |$1\, \rm {\rm M}_{\odot }$|⁠. The perturber excites a two-armed spiral structure, with a more prominent spiral feature for higher mass perturbers. The two arms evolve over time, eventually winding up, consistent with previous works. We focus on analysing the pattern speed and pitch angle of these spirals during the whole process. The initial pattern speed of the two arms are close to the angular velocity of the perturber at periastron, and then it decreases over time. The pitch angle also decreases over time as the spiral winds up. The spirals disappear after several local orbital times. An inclined prograde orbit flyby induces similar disc substructures as a coplanar flyby. A solar-mass flyby event causes increased eccentricity growth in the protoplanetary disc, leading to an eccentric disc structure which dampens over time. The spirals’ morphology and the disc eccentricity can be used to search for potential unbound stars or planets around discs where a flyby is suspected. Future disc observations at high resolution and dedicated surveys will help to constrain the frequency of such stellar encounters in nearby star-forming regions.

1 INTRODUCTION

Protoplanetary discs around young stars provide a unique laboratory for studying planet formation and planet–disc interactions. Spiral arms, stretching from tens to hundreds of au, are observed in numerous protoplanetary discs in visible to near-infrared light. Any perturbation in a rotating disc can excite spirals, and the spirals encode the information of the perturbation in the past (Zhu & Zhang 2022). Particularly, the gravitational interaction between a perturber and a disc excites strong spirals, and thus spirals are normally regarded as a signpost for undetected planets (e.g. Grady et al. 1999, 2013; Muto et al. 2012; Wagner et al. 2015; Monnier et al. 2019; Garufi et al. 2020; Muro-Arena et al. 2020). On the other hand, spirals can also be excited through flyby interactions. A flyby is defined as an object undergoing a single periastron passage of a stellar companion within |$100\!-\!1000\, \rm au$|⁠. Stars are formed in relatively dense clusters, and therefore the interaction between stellar counterparts is enhanced (Hillenbrand 1997; Carpenter 2000; Lada & Lada 2003; Porras et al. 2003). From the work of Pfalzner (2013) and Winter et al. (2018a), the probability of stellar flybys are about 30  per cent for solar-type stars during the first million years of stellar evolution, using a background stellar density that is larger than in Taurus. Considering that the spirals excited by an eccentric perturber are dramatically different from those excited by a circular perturber (Zhu & Zhang 2022), studying the spirals excited by an unbound flyby perturber is desired to understand the origin of the observed spirals and how planetary systems are formed.

Flyby encounters have been both simulated numerically and observed in protoplanetary discs. Cuello et al. (2019) ran extensive three-dimensional smoothed particle hydrodynamic simulations of protoplanetary discs (consisting of both gas and dust) being perturbed by a stellar flyby on a parabolic orbit. The authors explored a variety of different orbital inclinations, periastron distances, and mass ratios of the flyby event. They find that spiral structures are produced in both gas and dust, which have a lifetime of the order of thousands of years. The perturber truncates and warps the primary disc for a range of perturber inclinations and periastron distances (Clarke & Pringle 1993; Ostriker 1994; Terquem & Bertout 1996; Bhandare, Breslau & Pfalzner 2016; Xiang-Gruess 2016), leaving the dust disc more compact than the gas disc. Furthermore, there are multiple observational signatures that are present during a recent flyby encounter. Aside from spiral formation, long bridges of material are linked from the primary disc to the intruding flyby (Cuello et al. 2020). Clarke & Pringle (1993) demonstrated that a prograde, coplanar parabolic flyby encounter stripped material off the protoplanetary disc, and the perturber captured a portion of the stripped material. Warps and misalignments are common in the primary disc and are observable in moment one maps. Broken protoplanetary discs can have large mutual misalignments between the inner and outer gas rings that may be generated by a flyby scenario (Nealon, Cuello & Alexander 2020). There have been several cases of flyby candidates in RW Aur (Dai et al. 2015), HV Tau and DO Tau (Winter, Booth & Clarke 2018b), AS 205 (Kurtovic et al. 2018), UX Tauri (Ménard et al. 2020), Z CMa (Dong et al. 2022), and FU Ori (Borchert et al. 2022a). For a recent review on flyby’s shaping protoplanetary discs, see Cuello, Ménard & Price (2023).

Stellar flybys may also be responsible for shaping wide planetary systems (Bailey & Fabrycky 2019; Cai et al. 2019; Batygin et al. 2020; Li, Mustill & Davies 2020; Wang, Perna & Leigh 2020a,b), including our Solar system (Kobayashi, Ida & Tanaka 2005; Pfalzner et al. 2018). Stellar flybys can strongly perturb protoplanetary discs during their periastron passage, which can impact the total size and frequency of planetary systems (Scally & Clarke 2001; Kenyon & Bromley 2002; Olczak, Pfalzner & Spurzem 2006; Rosotti et al. 2014; Steinhausen & Pfalzner 2014; Portegies Zwart 2016; Vincke & Pfalzner 2016; Concha-Ramírez et al. 2019; Jiménez-Torres ). Most of these stellar flyby encounters are especially weak due to the fast relative velocity (Veras & Moeckel 2012). Still, low relative velocity encounters occurring within dense environments have a more dominant gravitational interaction (Li & Adams 2015). Flyby encounters are also hypothesized to excite the migration of outer giant planets, forming hot Jupiters (Rodet, Su & Lai 2021). The origin of the highly inclined Trans-Neptunian objects may also be associated with stellar flybys (Moore, Li & Adams 2020).

Free-floating (rogue) planets may also perturb discs by flyby interactions (Gladman & Chan 2006). A planet can become ejected from its parent system by planet–planet scattering (Rasio & Ford 1996; Lin & Ida 1997; Chatterjee et al. 2008). 75 per cent of exoplanetary system’s with giant planets are expected to undergo planet–planet scattering (e.g. Raymond & Morbidelli 2022). It is found to be easier to eject rocky planets rather than giant planets (Veras & Armitage 2005; Matsumura, Ida & Nagasawa 2013; Carrera, Davies & Johansen 2016), however ejection of giant planets can still occur. Rogue planets are more common to form around high-mass stars (Ma et al. 2016; Barclay et al. 2017). Recently, Mróz et al. (2020) found a free-floating planet candidate from microlensing events. Other mechanism for producing rogue planets include interactions in multistar systems (Spurzem et al. 2009; Kaib, Raymond & Duncan 2013), stellar flybys (Malmberg, Davies & Heggie 2011), or post-main-sequence evolution (Veras et al. 2011).

Protoplanetary discs around single stars are observed to have a lifetime of around |$1\!-\!10\, \rm Myr$| (Haisch, Lada & Lada 2001; Hernández et al. 2007, 2008; Mamajek 2009; Ribas, Bouy & Merín 2015). These lifetimes are long enough for discs to be perturbed by stellar flyby encounters (Cuello et al. 2019; Jiménez-Torres 2020). As the perturber approaches periastron passage, tidal effects prompt the formation of spirals (Ostriker 1994; Pfalzner 2003; Shen et al. 2010; Thies et al. 2010). External companions will excite spiral density waves at Lindblad resonances, which has been investigated analytically (Goldreich & Tremaine 1978, 1979; Lin & Papaloizou 1993; Lubow & Ogilvie 1998; Ogilvie & Lubow 2002) and numerically (Kley 1999; Dong et al. 2011; Zhu et al. 2015; Dong & Fung 2017; Hord et al. 2017; Bae & Zhu 2018; Zhang & Zhu 2020; Ziampras et al. 2020; Muley, Dong & Fung 2021). Spiral are launched at these resonances because the local orbital frequency is a multiple of the Doppler-shifted forcing frequency of the companion. If the companion is an external star, it exerts a strong tidal force where its Roche lobe can reach beyond the location of most of these resonances. The superposition of density waves excited at Lindblad resonances can constructively interfere, leading to the production of multiple spiral arms (Bae & Zhu 2018; Miranda & Rafikov 2019). In the case of a companion on an unbound orbit, density waves are not repeatedly excited, leading to the arms winding up overtime (e.g. Semczuk, Łokas & del Pino 2017; Pringle & Dobbs 2019; Kumar et al. 2022).

Spiral density waves can be described by two main parameters, their pattern speed and their pitch angle. The pattern speed is how quickly a spiral arm appears rotating azimuthally, and the pitch angle is the angle between the tangent to a spiral arm and a perfect circle, which measures how tightly the spiral arms are wound. The pattern speed and pitch angle of spiral arms in protoplanetary discs can be measured from observations (Kraus et al. 2017; Huang et al. 2018; Dong et al. 2018a; Casassus et al. 2021).

This work further examines the formation and evolution of spiral arms in protoplanetary discs that are produced by flyby encounters. We systemically measure the pattern speed and pitch angle of the excited spiral arms in discs around a solar mass star for different perturber masses in the range of |$M_2 = 10\!-\!1000\, \rm M_{\rm J}$|⁠, where M_J is Jupiter mass. The evolution of the pattern speed and pitch angle can be compared to observations in order to infer the presence of unbound perturbers. The origin of spirals has direct implications for both planet formation and disc evolution (Dong, Najita & Brittain 2018b; Brittain et al. 2020). The layout of the paper is as follows. In Section 2, we explain the set-up of the hydrodynamical flyby simulations and showcase the results in Section 3. In Section 4, we discuss how our results are related to the observed properties of the spiral arms seen in various systems. Lastly, we draw our conclusions in Section 5.

2 METHODS

To model a flyby of a massive object perturbing a protoplanetary disc, we use the three-dimensional smoothed particle hydrodynamical (SPH) code phantom (Price & Federrath 2010; Price et al. 2018). The code can model a wide variety of geometries of parabolic orbits. On top of this, the system’s angular momentum is conserved with the same order of accuracy as that of the time-stepping scheme. In phantom, linear and angular momentum are both conserved to round-off error, typically ∼10⁻¹⁶ in double precision. Recently, phantom has been well tested for parabolic flyby encounters interacting with gaseous and dusty discs (Cuello et al. 2019, 2020; Ménard et al. 2020; Nealon et al. 2020). The reference frame within our simulations is centred on the centre of mass of the system. A summary of the simulations is given in Table 1. The simulation ID is given in the first column, where the ‘M’ corresponds to the mass of the perturber and the numbers following give the mass value in units of Jupiter mass (M_J). For example, M100 is a |$100\, \rm M_J$| perturber. We also consider a higher resolution case, M100h, and inclined encounters, M100iy and M100ix.

2.1 Circumprimary disc set-up

We set up a coplanar gaseous protoplanetary disc around a generic solar-type star. The hydrodynamical disc is in the bending waves regime where the disc aspect ratio H/r is larger than the viscosity coefficient α (Shakura & Sunyaev 1973). Within this regime, the warp induced by the interaction of the flyby propagates as a pressure wave with speed c_s/2 (Papaloizou & Pringle 1983; Papaloizou & Lin 1995; Lubow & Ogilvie 2000), where c_s is the sound speed. We model an initially flat disc with 10⁶ Lagrangian particles with a total disc mass of |$0.001\, \rm {\rm M}_{\odot }$| We do, however, include one higher resolution simulation with 10⁷ particles (see Table 1). The low disc mass has a negligible dynamical effect on the orbit of the flyby. The mass of the primary star is |$M_1 = 1\, \rm {\rm M}_{\odot }$|⁠. The inner disc radius is initially set to |$r_{\rm in} = 10\, \rm au$| and the outer radius is |$r_{\rm out} = 100\, \rm au$|⁠. The disc lies initially in the x–y plane. The accretion radius for the central star is set to |$10\, \rm au$|⁠. The accretion radius is a hard boundary. Any Lagrangian particle that passes across the boundary is considered accreted, and the particle’s mass and angular momentum are deposited on to the star. Our primary sink accretion radius is comparable to the initial inner edge of the disc to reduce the computational time by neglecting to resolve close-in particle orbits. The surface density profile is initially a power-law distribution given by

where |$\Sigma _0 = 7.00\, \rm g\,cm^{-2}$| is the density normalization and p is the power-law index. We set p = +3/2 and the density normalization is defined by the total disc mass. We use a locally isothermal equation of state with a disc thickness that is scaled with radius as

where |$\Omega = \sqrt{GM_1/r^3}$| and c_s is proportional to r^−q. We use q = +1/2 so that the initial disc aspect ratio is constant at H/r = 0.05.

The Shakura & Sunyaev (1973) viscosity, α_SS, prescription is given by

where ν is the kinematic viscosity. To calculate α_SS, we use the prescription in Lodato & Price (2010) given as

where 〈h〉 is the mean smoothing length on particles in a cylindrical ring at a given radius (Price et al. 2018). We take the Shakura & Sunyaev (1973) viscosity parameter to be α_SS = 0.005, which sets the artificial viscosity to α_AV = 0.1260 (0.2713 for the high-resolution simulation). Note that a value of α_AV = 0.1 represents the lower limit, below which a physical viscosity is not well resolved in SPH, meaning that viscosities smaller than this lower limit produce disc spreading independent of the value of α_AV (see Meru & Bate 2012, for details). In hydrodynamics, we also include a parameter, β_AV, which provides a non-linear term that was originally introduced to prevent particle penetration in high Mach number shocks (e.g. Monaghan 1989). In general, β_AV = 2.0. In phantom, the smoothing length is specified as a function of the particle number density n given as

where h_fact is the proportionality factor specifying the smoothing length in terms of the mean local particle spacing. The disc is resolved with a shell-averaged smoothing length per scale height of 〈h〉/H ≈ 0.39 and 〈h〉/H ≈ 0.18 for our high-resolution simulation. Lastly, we ignore the effect of disc self-gravity.

2.2 Flyby set-up

We consider an unbound perturber that gravitationally influences a gaseous protoplanetary disc around the primary star. We simulate only parabolic encounters (e = 1), which induce the greatest star-to-disc angular momentum transfer and therefore engender the most prominent spiral features in the disc (Vincke & Pfalzner 2016; Winter et al. 2018b; Cuello et al. 2019). Perturbers on hyperbolic trajectories (e > 1) would have a higher angular velocity at the closest approach, leading to a less significant impact on the generation of structures within the disc (e.g. Winter et al. 2018b). Hyperbolic orbit flybys are also less efficient in capturing material compared to parabolic encounters (e.g. Larwood & Papaloizou 1997; Pfalzner, Umbreit & Henning 2005).

We model flybys on prograde coplanar and inclined trajectories. The centre of mass of the system is set as the centre of the x–y plane. We measure the tilt of the flyby orbit with respect to the z-axis. The distance between the perturber and the central star is defined as r₂. A coplanar perturber initially lies in the x–y plane (perturber’s angular momentum vector is parallel to the z–axis), arrives initially from the negative y direction, and leaves towards the negative y direction. We also consider two prograde inclined trajectories that are tilted by 45°. The inclined models are rotated clockwise with respect to the y-axis (model M100iy) and x-axis (model M100ix), respectively. For the coplanar and y-axis rotation cases, the 3D position of the pericentre does not change (similar to the simulations from Cuello et al. 2019) but that for the x-axis rotation the pericentre position is sent out of the disc plane. This may affect the excitation mechanism and the efficiency in terms of angular momentum transfer.

The probability of an encounter with a periastron distance |$100\!-\!1000\, \rm au$| is roughly 30  per cent but rapidly decreases for increasing stellar cluster age for leaky cluster environments (e.g. Pfalzner 2013). Therefore, we select the minimum separation between the primary and the perturber (periastron distance) to be |$r_{\rm p} = 200\, \rm au$|⁠, and hence the distance between the perturber and the centre of mass at the closest encounter varies with perturber mass. The periastron of the coplanar and y-axis rotated simulations occurs at x = 0, y > 0, and z = 0. For the x-axis rotated simulation, the periastron occurs at x = 0, y > 0, and z < 0. We consider perturber masses of |$10\, \rm M_{J}$|⁠, |$100\, \rm M_{J}$|⁠, and |$1\, \rm M_\odot$|⁠. The accretion radius of the perturber is set to |$1\, \rm au$|⁠.

3 RESULTS

We analyse the primary protoplanetary disc both radially and azimuthally using spherical shells. For our radial analysis, we divide the disc into 300 linear bins in spherical radius, r, which range from |$10\, \rm au$| (the initial inner disc edge) to 100 au (the initial outer disc edge). For our azimuthal analysis, we use 250 bins in the range 0 ≤ θ ≤ 2π at a given radius r, where θ is the azimuthal angle. Within each bin, radial or azimuthal, we calculate the mean properties of the particles, such as the surface density, inclination, longitude of ascending node, and eccentricity.

3.1 Flyby orbit

We now examine the angular velocity and distance of the flyby encounters during the simulations. The parabolic orbit in our model is described by

where r_p is the periastron distance, which occurs at θ = +π/2. The angular speed as a function of r₂ is then

where G is the gravitational constant, M₁ is the mass of the primary star, and M₂ is the mass of the perturber. The relationship between the radial distance, r₂, and time t is given by

where r₂ = r_p when t = t_p. Fig. 1 shows the angular velocity of the perturber for masses M₂ = |$10\, \rm M_J$|⁠, |$100\, \rm M_J$|⁠, and |$1\, \rm M_\odot$|⁠. The peak of the angular velocity occurs when the perturber is at periastron. The more massive the perturber, the higher the angular velocity at periastron, and the time of periastron passage occurs sooner. The derived angular speed for each parabolic companion from equation (8) is shown by the black-dotted curves, which is consistent with the simulation results. Fig. 2 shows the separation of the unbound perturber to the primary star. The troughs in each model represent the time at which the flyby reaches periastron, which is at a distance of |$200\, \rm au$| by construction. The derived radial separation for each parabolic companion from equation (9) is shown by the black-dotted curves, which is consistent with the simulation results.

3.2 Coplanar non-equal mass flybys

In this subsection, we first analyse the simulations with a non-equal mass perturber, models M10 and M100 from Table 1. In Fig. 3, we plot the peak surface density, Σ_peak, at a radius |$r = 65\, \rm au$| as a function of time for |$M_2 = 10,\, 100\, \rm M_J$|⁠. The peak surface density is defined as the maximum surface density value at each time. At |$r = 65\, \rm au$|⁠, we have a relatively stronger density contrast near the arms. The maximum point of the curves corresponds to the strength of the spiral arms that are produced. The density enhancement for the |$10\, \rm M_J$| model is a factor of ∼1.2, while it is a factor of ∼3 for the |$100\, \rm M_J$| model. Next, we visualize the spiral structure that is excited within the disc. Fig. 4 shows the disc evolution at four different times for a |$10\, \rm M_J$| (first row, M10), and |$100\, \rm M_J$| (second row, M100) flyby encounter. For visualization purposes only, the coordinates and velocities are centred on the primary star, marked with a green dot at x = 0 and y = 0. The first column shows the discs at |$t=0\, \rm yr$|⁠, with a trace of the perturber’s parabolic orbit given in the top panel. The second column shows the flyby at periastron (or closest approach), seen by the blue dot to the north. At this time, the perturbers have not excited any substructures in each case. The third and fourth columns show times shortly after periastron passage, where the perturbers excite a two-armed spiral structure in the disc. A lower mass perturber excites weaker spiral arms than a higher mass perturber. The arms begin to wind up on themselves when the perturber moves away from the primary star. The spiral arms generated from a more massive companion is more prominent – leading to a higher surface density concentration within the arms, as shown in Fig. 3. When a two-armed spiral is excited, the arm that is excited initially closest to the perturber is denoted as ‘arm 1’, and the arm excited initially furthest from the perturber is ‘arm 2’. A depiction of arms 1 and 2 are given in the bottom row, third panel in Fig. 4.

We now examine the azimuthal variation of the surface density at a given radius for each coplanar non-equal mass perturber simulation. Fig. 5 shows the azimuthal cut of the disc surface density at a radius of |$r = 65\, \rm au$|⁠. The top panel shows the disc structure under the influence of a |$10\, \rm M_J$| (M10) perturber, and the bottom panel for a |$100\, \rm M_J$| (M100) perturber. The angular position of the flyby is given by the solid blue curve. The time of periastron passage occurs at θ = +π/2 in each simulation. Before the time of the closest approach of the perturber, the protoplanetary disc appears quiescent, with no excited substructures. Soon after periastron passage, the perturber excites a two-armed spiral within the disc, seen by the azimuthal striations in the surface density. From our comparison of different unbound-perturber masses, we find that the more massive perturber excites more prominent arms versus a lower mass perturber. The density enhancement between the two cases is about a factor of 2, which we note can be distinguishable by observations (see Section 4). On average, the spiral arms in each simulation last for approximately |$5000\, \rm yr$| (equivalent to ∼2 Keplerian orbits at r_p) after the closest approach of the flyby, as shown in Fig. 5. When the spiral arms fully wind up, the resulting disc structure restores to the initial disc structure.

In Fig. 6, we display the spiral arm contrast versus azimuthal angle for a |$10\, \rm M_J$| perturber (blue lines) and |$100\, \rm M_J$| perturber (red lines). The solid lines denote the spiral arm at radius |$r = 65\, \rm au$|⁠, and the dotted lines at |$r = 85\, \rm au$|⁠. The two peaks in the surface density represent the individual arms. For example, for |$M_2 = 100\, \rm M_J$|⁠, arm 1 is located just beyond θ ∼ π, and arm 2 is located just below θ ∼ 2π. The density contrast between the arms is identical for each perturber mass at each radius.

In order to map out arm 1 and arm 2, we apply a tracing technique shown in the bottom panel of Fig. 5. At each time-step, we track the highest density located within the spirals. We continue tracing the arms until both arms cannot be fully traced. We then store the azimuthal angle and time that comprise the spiral arms. The tracing procedure is robust, and can be applied to all other simulations. We begin the analysis of the pattern speed and pitch angle at |$t = 8000\, \rm yr$| and track the properties of the spiral arms until |$t = 11500\, \rm yr$|⁠. This selected time range ensures that the spiral arms are better defined, resulting in a high-precision analysis.

The pattern speed of an arm is the derivative of its azimuthal position at a given radius with respect to time. We compute the pattern speeds of the spiral arms by measuring the slope of the curves from our spiral arm tracings. Fig. 7 shows the pattern speed of arm 1 (solid curve) and arm 2 (dotted curve) as a function of time for each perturber mass. We calculate the pattern speed at two different radii: 65 and |$85\, \rm au$|⁠. We smooth the data using a Gaussian filter. The pattern speed of the spiral arms are higher closer to the central star. The angular velocity of the perturbers are shown by the black and grey curves and the Keplerian angular velocity at radii |$r = 65\, \rm au$| and |$85\, \rm au$| are shown by the horizontal dashed lines. The initial pattern speeds of the arms are closely set by the angular velocity of the perturber at periastron. The pattern speed of the spiral arms decreases with time and the arms ultimately dissipate ∼5000 yr after the closest approach of the flyby. For a given radius, the pattern speeds of arm 1 and arm 2 oscillate out of phase with respect to one another. This phenomenon can be seen clearly in the blue curves from Fig. 7, where the trough of one arm occurs at the peak of the other arm. The oscillation amplitude of the arms induced by a lower mass perturber is lower than the oscillations driven by a higher mass perturber. We tested this phenomenon with a bound perturber, finding the pattern speed of the arms does not display such oscillations but instead has a constant pattern speed at a given radius, as expected, thus ruling out numerical causes of the phenomenon. We defer further exploration of this phenomenon to a future work.

For the remainder of this subsection, we analyse only the two-armed spiral excited from the |$100\, \rm M_J$| perturber since the arms are the most conspicuous. Besides the pattern speed, we also calculate the pitch angle of the spirals. For simplicity, we only analyse the pitch angle for arm 1. The pitch angle, β, is defined as in Binney & Tremaine (2008) with

where r is the radial distance, and θ is the azimuthal position. The dr/dθ term is calculated from our simulations by tracing the spiral arms in a (r, θ) grid, using a similar method as in Fig. 5. To visualize the evolution of the pitch angle, we show the azimuthal surface density as a function of logarithmic radius at five different times in Fig. 8. The top panel shows the time at periastron passage of the |$100\, \rm M_J$| perturber (M100). We mark the peak surface density along arm 1 at the four chosen radii, 55, 65, 75, and |$85\, \rm au$|⁠. The calculation of the pitch angles only begins when the spirals are excited. The perturber excites a two-armed spiral structure soon after the closest approach, shown in the second panel (⁠|$t = 8001\, \rm yr$|⁠). As time increases, the spiral arms wind up, decreasing the pitch angle. The lower panels show times at which the spiral structure is winding up. Eventually, the spiral arms dissipate, as shown in the bottom panel. After calculating dr/dθ, we calculate the pitch angle evolution using equation (10) at |$r = 55\, \rm au$|⁠, |$65\, \rm au$|⁠, |$75\, \rm au$|⁠, and |$85\, \rm au$| given in Fig. 9. A Gaussian filter is applied to smooth the data. As the radial distance in the disc increases, the larger the initial pitch angle. At |$r=55\, \rm au$|⁠, the initial pitch angle is β ∼ 4°, while at |$r=85\, \rm au$|⁠, β ∼ 12°. The pitch angle decreases in time as the spiral winds up, regardless of radius, which is consistent with the higher pattern speed closer to the star found in Fig. 7.

3.3 Resolution study

We present a resolution study for the simulation with a |$100\, \rm M_J$| perturber. In SPH simulations, we use an artificial viscosity to model an expected Shakura & Sunyaev (1973) viscosity coefficient. However, the number of Lagrangian particles impacts how close the artificial viscosity is to the actual value. Here, we compare two simulations with 10⁶ particles (M100, α_AV = 0.1260) and 10⁷ particles (M100h, α_AV = 0.2713). Fig. 10 shows the evolution of the disc surface density at 65 au for the two simulations. The top sub-panel is the lower resolved disc, while the bottom sub-panel is the more highly resolved disc. The spirals in the higher resolution are more conspicuous. The lifetime of the spiral arms at higher resolution is longer, but only persists an extra few hundred years. Considering the total lifetime of the spirals for the lower resolution is |$\sim 5000\, \rm yr$|⁠, the higher resolution spirals’ lifetime is increased by less than 10 per cent. We expect this factor in spiral longevity to be greater for lower viscosity. In this work, we select a target viscosity of α_SS ∼ 0.005, while observationally, protoplanetary discs can have viscosities in the range 10⁻³–10⁻⁵ (e.g. Flaherty et al. 2015; Pinte et al. 2016; Flaherty et al. 2020; Villenave et al. 2022). Protoplanetary discs can have higher viscosities depending on the disc size (e.g. Hartmann et al. 1998; Ansdell et al. 2018). Therefore, observed protoplanetary discs perturbed by flybys may have spiral arms that last longer than our numerical simulations presented in this work.

To further analyse simulations with a different initial number of particles, we compare the pattern speed and pitch angle for arm 1 between M100 and M100h in Fig. 11. The more resolved simulation shows that the spiral arm survived for a more extended period of time than the lower resolution spiral. For both the pattern speed and pitch angle, the lower resolution spiral mirrors the higher resolution spiral initially. However, deviations appear near the end of the low-resolution spiral’s lifetime. This is due to the lowly resolved spiral beginning to dissipate. Despite that, the differences between the two simulations are minimal, and modelling the simulations with 10⁶ particles appears sufficient for capturing the dynamics at play.

3.4 Coplanar equal mass flyby

Up to this point, we have only considered sub-stellar flyby perturbers. We now explore a stellar-mass encounter on an unbound parabolic orbit (M1000 from Table 1). Fig. 12 shows snapshots of the disc structure during the stellar-mass flyby. The top-left panel shows the initial disc structure at |$t = 0\, \rm yr$|⁠. The top row panels slowly progress in time as the perturber reaches periastron in the top-right panel. At periastron, there is a tidal tail of material that extends from the primary disc to the perturber. Shortly after the periastron passage (bottom-left panel), the perturber excites two prominent spiral arms with large pitch angles. The arm extending towards the perturber shows substructures. We are uncertain whether these substructures are numerical or physical, i.e. an instability. We ran additional simulations with increased viscosity and particle resolution but the substructures are still present. The arms do not survive for long; however, the inner cavity becomes elongated, which is consistent with the disc being eccentric. The eccentric disc structure survives for the duration of the simulation shown in the bottom-right panel, ∼15,000 yr after periastron passage. The radial separation of the flyby at the end of the simulation is |$\sim 4000\, \rm au$|⁠, meaning the disc is maintaining the eccentric structure even though the perturber is not close to perturb the disc. This eccentricity growth induced in protoplanetary disc by flyby encounters is consistent with fig. 10 from Cuello et al. (2019). We show further analysis of the short-lived spiral arms in the appendix.

We further investigate the eccentricity growth induced by a flyby event by comparing the disc surface density and eccentricity profiles for the |$10\, \rm M_J$|⁠, |$100\, \rm M_J$|⁠, and |$1\, \rm M_\odot$| mass perturbers. The left panel of Fig. 13 shows the azimuthally averaged disc surface density as a function of time, with the colour map denoting the surface density value. The top, middle, and bottom sub-panels are for a |$10\, \rm M_J$|⁠, |$100\, \rm M_J$|⁠, and |$1\, \rm M_\odot$| mass perturber, respectively. The time of periastron passage of the encounter is shown by the vertical magenta line in each sub-panel. The surface density profiles for the sub-stellar mass encounters are similar, with little truncation of the outer radius of the disc. The bottom sub-panel shows the disc surface density for a stellar flyby. After the closest approach, the primary disc is strongly truncated by the perturber. The outer disc edge is truncated to |$\sim 60\, \rm au$|⁠, which is consistent with the tidal truncation radius of an equal-mass circular binary being one-third of the binary separation (Artymowicz & Lubow 1994). The right panel of Fig. 13 shows the azimuthally averaged disc eccentricity as function of time, with the colour map denoting the eccentricity value. The eccentricity is determined by calculating the magnitude of the eccentricity vector for each radial bin. The stellar-mass perturber excites large eccentricity growth in the disc compared to the sub-stellar mass perturbers, especially in the inner regions of the disc. The more massive perturber imparts more energy and angular momentum to the wave that are launched at Lindbald resonances. The disc then undergoes changes in angular momentum where the excited waves dampen, resulting in the disc eccentricity growth (e.g. Lubow 1991). The eccentricity in the disc dampens over time due to pressure forces that recircularize the disc (Pfalzner 2003). The accretion radius of the sink has a minor effect on the eccentricity growth in the outer regions of the disc (see appendix C in Cuello et al. 2019).

3.5 Inclined flyby

We now investigate spiral arms excited by a |$100\, \rm M_J$| flyby on an inclined trajectory with an initial misalignment of 45°. We consider two different orbital orientations, one rotated with respect to the y-axis and the other with respect to the x-axis. Fig. 14 shows the peak surface density, Σ_peak, at a radius |$r = 65\, \rm au$| as a function of time for a coplanar, y-axis rotated, and x-axis rotated flyby. The maximum point of the curves corresponds to the strength of the spiral arms that are produced. A coplanar perturber excites a larger peak surface density. The y-axis rotated flyby excites the second highest peak surface density, followed by the x-axis rotated flyby. In Fig. 15, we show the spiral arm contrast versus azimuthal angle for the coplanar perturber (blue lines) and inclined perturber with a x-axis rotation (red lines). The solid lines denote the spiral arm at radius |$r = 65\, \rm au$|⁠, and the dotted lines at |$r = 85\, \rm au$|⁠. The two peaks in the surface density represent the individual arms. The density contrast between the arms is identical at each radius.

Next, we compare the spiral structure that is excited within the disc between the coplanar and inclined models. Like the coplanar model, we select a time range to conduct the analyse to reduce the noise. For the rotation models, we analyse the spiral arms from |$t = 8500\, \rm yr$| to |$t = 11500\, \rm yr$|⁠. Fig. 16 shows an azimuthal cut of the disc surface density as a function of time. The top panel shows the coplanar flyby, the middle panel shows the y-axis rotation, and the bottom panel shows the x-axis rotation. Soon after periastron passage of the perturber (at θ = +π/2), two spiral arms are excited within the disc for each model. The spiral structure excited by the inclined cases are similar to the spirals excited by a coplanar flyby of the same mass. For both inclined models, arm 1 is first excited at a time of |$\sim 7000\, \rm yr$|⁠, and lasts until |$\sim 15000\, \rm yr$|⁠. The spiral arms generated from the coplanar companion are more prominent – leading to a higher surface density concentration within the arms. Spirals excited by the x-axis rotated flyby are less prominent than spirals excited by the y-axis rotated flyby, as shown back in Fig. 14.

We analyse both the pattern speed and pitch angle of arm 1 for the inclined models and then compare to the coplanar case. The pattern speed and pitch angle calculations are performed using the same techniques presented earlier. Fig. 17 shows the pattern speed of arm 1 at |$65\, \rm au$| for each model, along with the angular velocity of the perturber. The initial pattern speed of arm 1 from each model is slightly higher than the perturber’s angular velocity during periastron passage, but decreases over time. The pattern speed in the inclined cases is similar to the coplanar model but with less oscillations. Fig. 18 shows the pitch angle of arm 1 as |$65\, \rm au$| for the coplanar and inclined models. Similar to the coplanar case, the pitch angle for the inclined models decrease in time as the spiral arms wind up. As the flyby orbit is rotated about the y-axis, the initial pitch angle as a function of radius is lower than the coplanar flyby simulation. The x-axis rotation model has a lower initial pitch angle than the y-axis rotation model. Overall, the inclined scenario shows similar spiral structures as the coplanar model.

4 DISCUSSION

The interaction between the flyby and a protoplanetary disc can heat up the disc through compression and shocks (Lodato et al. 2007). In our phantom simulations, we do not take this heating into account. Protoplanetary discs with excited spirals can also be heated by absorbing stellar radiation, which produces a positive vertical temperature gradient (Chiang & Goldreich 1997). The temperature gradient causes a thermal stratification in the disc, which can alter the morphology of the spiral wake (Lee & Gu 2015). Juhász & Rosotti (2018) found that the pitch angle of the spirals produced by an embedded planet in a thermally stratified disc decreases near the disc mid-plane, and then increases towards the disc surface. These results can be applied to a stellar companion perturbing a disc, where the spiral pitch angles excited by a flyby event will be affected by a thermally stratified disc. Moreover, previous hydrodynamics simulations of flyby–disc interactions used a radial surface density profile of p = 1 to match the observed disc profiles (e.g. Cuello et al. 2019, 2020), which loads more material in the outer disc regions compared to p = 1.5 used in this work. For p = 1, we expect brighter spirals compared to p = 1.5. Given the low disc mass simulated in this work, the dynamical behaviour of the disc during the encounter does not sensitively depend on the initial surface density profile.

The density variation in the spiral arms between our hydrodynamical models with different mass perturbers is about a factor of 2, which may be distinguishable by observations. Cuello et al. (2020) provided observational signatures of flyby encounters. From their results, spiral arms excited by flybys can be observed in the continuum, which means the encounter is ongoing or the perturber is a few thousand au from the host star. Moreover, the spiral arm between the two stars is brighter since it is more illuminated. Inclined encounters will excite spiral arms that are not coplanar with the disc, resulting in unique kinematic signatures that can be resolved with a spectral resolution of the order of 1 km s⁻¹ to accurately map the vertical layers of the disc (Cuello et al. 2019). In Dong & Fung (2017), they showed spiral contrasts in scattered light (see their fig. 5) versus spiral contrast in surface density (see their fig. 1). They found a factor of 2 contrast in surface density roughly translates into a factor of 2 in scattered light surface brightness, and is detectable after convolving the images to achieve a realistic angular resolution. Even if the density enhancement is a factor of 2 at the mid-plane, it can be a lot higher at the surface (Speedie, Booth & Dong 2022). The synthetic observations of flyby-induced spirals may provide direct evidence that they can be visible (e.g. fig. 4 in Cuello et al. 2020). Therefore, the density enhancement of a factor of 2 in the spiral arms between our different mass perturber encounters should be distinguishable by observations.

There are several observational examples of spiral arms in protoplanetary discs with no confirmed companion, with one being Elias 2–27 (e.g. Pérez et al. 2016). As shown in this work and from previous works, a way to produce spirals without a bound companion is a flyby encounter. Another way to produce spirals without any perturber is gravitationally unstable protostellar discs (e.g. Baehr & Zhu 2021). Early evolution of a star–disc system suggests the mass of the central star is comparable to the disc mass, forcing the system to be in the self-gravitating regime (Lin & Pringle 1987, 1990). Numerical simulations of self-gravitating discs with a Toomre parameter |$\mathcal {Q} \sim 1.5\!-\!1.7$| produced non-axisymmetric perturbations that grow into spiral density waves (Durisen et al. 2007). In the context of Elias 2–27, the spiral arms observed in the system are thought to be caused by the gravitational instability (Meru et al. 2017; Tomida et al. 2017). However, Forgan, Ilee & Meru (2018) postulated that a massive stable disc undergoing a stellar companion encounter can produce spirals with pitch angles that increase with radius (e.g. Veronesi et al. 2021). Measuring the pattern speed of the spiral arms in Elias 2–27 offers a more direct assessment of the origin of the spirals. If the pattern speed of the spirals decreasing in time, Elias 2–27 may have undergone a flyby in the past. A single value for each spiral’s pitch angle was measured, β ∼ 15°, by Huang et al. (2018) without a radial dependence. Currently, the measured values do not favour a specific scenario. Shuai et al. (2022) used the Gaia Data Release 3 (DR3) to determine if flyby encounters occurred for 20 systems with spiral arms. For their selected systems, their analysis suggested that an unbound encounter was not the cause for isolated spiral arms in scattered light. We note that the Gaia DR3 does not have radial velocity information for three-dimensional motion tracing and thus only can map two-dimensional on-sky locations of the systems. More observations are needed to ascertain the origin of the spirals.

The spirals excited by a flyby perturber are different from the GI spirals and spirals excited by a bound perturber in several ways. The pitch angle of GI spirals are ∼10° throughout the disc (Baehr & Zhu 2021; Béthune, Latter & Kley 2021), The pitch angle of spirals excited by either a bound or unbound perturber decreases further away from the perturber (Zhu et al. 2015; Bae & Zhu 2018). As shown in this work, the spirals from an unbound perturber would wind up and disappear after several orbits. The spirals from a circularly orbiting perturber have a steady shape, co-orbiting with the perturber. The spirals from an eccentric perturber would also evolve and change its shape with time but normally have multiple complex spiral arms (Zhu & Zhang 2022).

The hydrodynamical simulations conducted in this work only dealt with one specific disc size, which modelled a non-penetrating encounter. Increasing the disc size makes the model approach a grazing encounter. During a grazing encounter, spiral arms should still be produced within the primary disc; however, material may be captured by the perturber, and form a secondary disc (Breslau, Vincke & Pfalzner 2017; Cuello et al. 2019). If the outer disc edge encloses the periastron of the flyby, we then have a penetrating encounter. During a penetrating encounter, the primary disc size should be truncated, and accretion outbursts may occur due to the flyby (Vorobyov et al. 2020; Borchert et al. 2022a, b). As shown by our resolution study (Section 3.3), viscosity affects the lifetime of the spiral arms. Less viscous discs that encounter a flyby should maintain spiral features and disc substructures on a longer time-scale than more viscous discs.

From our hydrodynamical simulations, flyby encounters can significantly affect the structure of protoplanetary discs, which has implications for planet formation. Dust grains can become trapped within spiral arms due to local gas pressure maxima (Haghighipour & Boss 2003a, b; Baehr & Zhu 2021). Dust trapping could promote grain growth, and eventually planetesimal formation (Baehr, Zhu & Yang 2022). Dust particles with Stokes number ∼1 are more prone to become trapped at the local pressure maxima in protoplanetary discs (Paardekooper & Mellema 2006; Rice et al. 2006). Thies et al. (2010) demonstrated that flyby encounters may trigger gravitational instability, which can lead to the formation of brown dwarfs and giant planets in the outer regions of the disc. Hence, massive companions or planets can be formed at great distances from their host stars.

Perturbers with few tens of Jupiter masses are able to excite spiral arms and disc substructures, including eccentricity growth. Another key observational signature in a disc that has been recently perturbed by a flyby is the presence of disc eccentricity that damps over time (see Fig. 13). A significant disc eccentricity can leave an observable imprint in the CO emission-line profiles, which appears as an asymmetry in the line profile (Regály et al. 2011). For example, Flaherty et al. (2015) observed a small asymmetry in the CO line profiles of the HD 163 296 protoplanetary disc, which may be consistent with the disc being eccentric. In the case of disc eccentricity excited by a flyby event, the observed asymmetry in the CO line profiles should slowly dissipate over time as the disc eccentricity dampens. However, given that the eccentric disc around HD 163 296 is on large scale, the time-scale for observable eccentricity decrease should be longer than the detectability time-scale of |$\sim 10\, \rm yr$| (operating time for ALMA).

Based on the results of this work, the changes in the spiral structure occur over time-scales of thousands of years, which are not observable in any practice sense. However, it may be possible to distinguish the origin of spiral arms at a single epoch by measuring the pitch angle of the spirals (see table 1 from Yu, Ho & Zhu 2019). A disc perturbed by a flyby should eventually have more tightly wound spirals than spiral arms excited by bound companions or GI. The pitch angle of spirals excited by both bound companions (the inner arms; Dong et al. 2015) and GI (Baehr & Zhu 2021; Béthune et al. 2021) can be ≳10°.

Is it possible to disentangle the effects of a bound binary from a flyby event? Binary stars will excite eccentricity growth in circumbinary discs (e.g. Dunhill, Cuadra & Dougados 2015; Ragusa et al. 2020; Smallwood, Lubow & Martin 2022; Siwek et al. 2023). Distinguishing these two mechanisms at a single epoch would be rather challenging. However, due to kinematic measurements of stellar masses and recent Gaia data, it is possible to better constrain the orbits of (supposedly) interacting stars. This is done for instance in Weber et al. (2023) (see their fig. 11) for the system AS205N and AS205S. For that specific case, despite the degeneracy, they determine that the orbit is very likely parabolic/hyperbolic, not bound. In addition, for the flyby case, the eccentric disc eventually circularizes once the flyby leaves the system. Moreover, if the dynamical centre of an eccentric disc (i.e. the foci of the disc) can be determined precisely, perhaps in the circumbinary disc cases one can see that the disc is not dynamically centred on the primary star (often the more visible star). The spirals excited by bound stellar companions (e.g. Gonzalez et al. 2020; Poblete et al. 2020; Rosotti et al. 2020; Poblete et al. 2022) evolve differently to spiral excited by stellar flybys. For example, the pattern speed and pitch angle of the spiral arms should remain rather constant in time with a bound companion. This is in direct contrast to flyby events, where the pattern speed and pitch angle evolve in time.

Flybys – although considered as violent and rather destructive events – can generate favourable conditions for planet formation in the long run. After disc truncation, the surface density of the disc is enhanced locally in regions where the eccentricity remains low (see Fig. 13). In any case, the eccentricity excitation would dampen due to pressure forces in the disc that causes the disc to recircularize (Pfalzner 2003). Once the eccentricity excitation dissipates, the disc becomes more compact and denser, which can be a favourable environment for planet formation. Cuello et al. (2019) found that removing gas and small dust grains at large disc radii effectively leads to higher values of dust-to-gas ratios in the inner parts of the disc. Speculatively, this could lead to the formation of self-induced dust traps (Gonzalez, Laibe & Maddison 2017; Vericel et al. 2021) or potentially trigger planetesimal formation (Johansen, Youdin & Mac Low 2009; Bai & Stone 2010; Yang, Johansen & Carrera 2017; Li & Youdin 2021; Schaffer, Johansen & Lambrechts 2021).

5 CONCLUSIONS

We investigate the interaction of a protoplanetary disc with a flyby using three-dimensional SPH simulations. The unbound perturber excites a two-armed spiral structure in the disc. We model both sub-stellar and stellar flybys, as well as coplanar and inclined flyby trajectories. For the sub-stellar flybys, the more massive the perturber, the more prominent the spiral arms. The pattern speed of the spiral arms excited by coplanar and inclined flybys is close to the angular velocity of the perturber at periastron. Overall, spiral arms last for a couple of dynamical time-scales at the radius of the periastron before winding up and disappearing. From our resolution study, we found that disc viscosity plays a role in the lifetime of flyby-included spiral arms. Discs with lower viscosity maintain the spiral arms for a longer period of time versus discs with higher viscosity. The stellar mass flyby strongly perturbs the disc, forcing significant eccentricity growth, which dampens overtime.

From our hydrodynamical simulations of flyby interactions, the pattern speed and pitch angle of the flyby-induced spiral arms evolve in time. Both the pattern speed and pitch angle of spiral arms decrease over time as the spiral arms wind up. Therefore, the pitch angle could be a unique observation signature for flyby-induced spirals. The pattern speed, pitch angle, and general morphology of the spirals excited by an inclined flyby are similar to a coplanar flyby. For typical flyby configurations, the temporal window to search for nearby perturbers is of the order of a few kyr (tens of kyr at most). Here, we postulate that the spirals’ morphology and the disc eccentricity can be used to search for/chase potential unbound stars or planets around discs where a flyby is suspected. Future disc observations at high resolution and dedicated surveys will help to constrain the frequency of such stellar encounters in nearby star-forming regions.

Acknowledgement

We thank the anonymous referee for helpful suggestions that positively impacted the work. We acknowledge the use of SPLASH (Price 2007) for the rendering of the figures. JLS acknowledges funding from the ASIAA Distinguished Postdoctoral Fellowship. CCY and ZZ acknowledge the support by NASA via the Emerging Worlds program (grant number 80NSSC20K0347) and the Astrophysics Theory Program (grant number 80NSSC21K0141). CCY is also grateful for the support by NASA TCAN program (grant number 80NSSC21K0497). RGM acknowledges support from NASA through grant 80NSSC21K0395. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 90783311. This research was funded, in part, by ANR (Agence Nationale de la Recherche) of France under contract number ANR-22-ERCS-0002-01. This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 101042275, project Stellar-MADE).

DATA AVAILABILITY

The data supporting the plots within this article are available on reasonable request to the corresponding author. A public version of the phantom and splash codes are available at https://github.com/danieljprice/phantom and http://users.monash.edu.au/~dprice/splash/download.html, respectively.

References

APPENDIX A: EQUAL-MASS FLYBYS

We examine the azimuthal variation of the surface density at a given radius for the equal-mass perturber simulation. Fig. A1 shows the azimuthal cut of the disc surface density at a radius of |$r = 65\, \rm au$|⁠. The angular position of the flyby is given by the solid blue curve. The time of periastron passage occurs at θ = +π/2 in each simulation. Before the time of the closest approach of the perturber, the protoplanetary disc appears quiescent, with no excited substructures. Soon after periastron passage, the perturber excites a two-armed spiral within the disc, seen by the azimuthal striations in the surface density. Unlike the unequal mass simulations, the two-armed spiral is short-lived, only lasting until |$\sim 6000\, \rm yr$|⁠. Unlike the unequal-mass models, soon after the two-arm spirals are excited, the arm led by the perturber develops wiggles and then merges into the other arm. The remaining one-arm structure is short-lived, lasting for about 6000 yr. In the mean time, the perturber excites eccentricity growth within the disc.

Fig.A2 shows the pattern speed of arm 1 (solid curve) and arm 2 (dotted curve) as a function of time for the equal-mass perturber simulation. We calculate the pattern speed at |$65\, \rm au$|⁠. We smooth the data using a Gaussian filter. The angular velocity of the perturber is shown by the black curve and the Keplerian angular velocity at radii |$r = 65\, \rm au$| is shown by the horizontal dashed line. The initial pattern speeds of the arms are closely set by the angular velocity of the perturber at periastron. The pattern speeds of the arms also display an oscillation in the amplitude, similar to Fig. 7.

To visualize the evolution of the pitch angle, we show the azimuthal surface density as a function of logarithmic radius at five different times in Fig. A3. The top panel shows the time at periastron passage of the |$1\, \rm {\rm M}_{\odot }$| perturber (M1000). We mark the peak surface density along arm 1 at the radii, |$r = 65\, \rm au$|⁠. The calculation of the pitch angles only begins when the spirals are excited. The perturber excites a two-armed spiral structure, shown in the top panel. As time increases, the initial spiral arms are short-lived. The lower panels show times at which the spiral structure is winding up and also begin having a chaotic evolution. We calculate the pitch angle evolution using equation (10) at |$r = 65\, \rm au$|⁠, given in Fig. A4. We smooth the data using a Gaussian filter. The pitch angle first increases in time as the spiral arms are excited and then begins to decrease in time as the spiral winds up.