The circumstellar disk response to the motion of the host star
^{1} Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, 1121 Budapest, Konkoly Thege Miklós út 15–17, Hungary
email: regaly@konkoly.hu
^{2} Department of Astrophysics, University of Vienna, 1180 Vienna, Austria
^{3} Research Institute of Physics, Southern Federal University, Stachki Ave. 194, 344090 RostovonDon, Russia
Received: 20 June 2016
Accepted: 11 January 2017
Context. Gridbased hydrodynamics simulations of circumstellar disks are often performed in the curvilinear coordinate system, in which the center of the computational domain coincides with the motionless star. However, the center of mass may be shifted from the star due to the presence of any nonaxisymmetric mass distribution. As a result, the system exerts a nonzero gravity force on the star, causing the star to move in response, which can in turn affect the evolution of the circumstellar disk.
Aims. We aim at studying the effects of stellar motion on the evolution of protostellar and protoplanetary disks. In protostellar disks, a nonaxisymmetric distribution of matter in the form of spiral arms and/or massive clumps can form due to gravitational instability. Protoplanetary disks can also feature nonaxisymmetric structures caused by an embedded highmass planet or a largescale vortex formed at viscosity transitions.
Methods. We use 2D gridbased numerical hydrodynamic simulations to explore the effect of stellar motion. We adopt a noninertial polar coordinate system centered on the star, in which the stellar motion is taken into account by calculating the indirect potential caused by the nonaxisymmetric disk, a highmass planet, or a largescale vortex. We compare the results of numerical simulations with and without stellar motion.
Results. We found that the stellar motion has a moderate effect on the evolution history and the mass accretion rate in protostellar disks, reducing somewhat the disk size and mass, while having a profound effect on the collapsing envelope, changing its inner shape from an initially axisymmetric to a nonaxisymmetric configuration. Protoplanetary disk simulations show that the stellar motion slightly reduces the width of the gap opened by a highmass planet, decreases the planet migration rate, and strengthens the largescale vortices formed at the viscosity transition.
Conclusions. We conclude that the inclusion of the indirect potential is recommended in gridbased hydrodynamics simulations of circumstellar disks which use the curvilinear coordinate system.
Key words: accretion, accretion disks / stars: formation / stars: protostars / planetdisk interactions / protoplanetary disks / methods: numerical
© ESO, 2017
1. Introduction
The formation mechanisms of stars and planets are intrinsically connected. The stars form owing to the gravitational collapse of molecular clouds and circumstellar disks develop around the forming stars thanks to the conservation of angular momentum. The circumstellar disk, also known as the protostellar disk, is embedded in the parental nebula in the early phase, in which the mass accretion from the parental nebula onto the disk exceeds the mass loss from the disk to the star. If the protostellar disk grows sufficiently large (about ≳ 0.07 M_{⊙}) it becomes gravitationally unstable, often resulting in the formation of gravitationally bound clumps. As the parental nebula is consumed, the nebular accretion ceases and the protostellar disk evolves into the protoplanetary disk, which is thought to be the birth place for planets.
The evolution of gaseous protostellar and protoplanetary disks and their interactions with the newly born planets can be studied by means of hydrodynamical models. In hydrodynamical simulations, one solves the fluid dynamics equations describing the time evolution of mass density, momentum, and energy. Since the advection term in the momentum equation makes the system of equations nonlinear, the analytical investigations are applicable only for highly idealized flows. Therefore, numerical methods have been developed since the mid1960s to solve the system of hydrodynamical equations numerically.
Nowadays, there are two widely used different approaches to solving those differential equations numerically: smoothedparticle hydrodynamics (SPH) and gridbased hydrodynamics. In the SPH method (see a review of Springel 2010) the fluid is handled as a set of discrete elements, referred to as particles, while in the gridbased approach (see a review of Teyssier 2015) the flow is handled with discrete values of fluid properties on a numerical grid.
In astrophysical models that are designed to study the formation of stars and planets, gravitational interaction naturally enters the system of hydrodynamical equations. If any nonaxisymmetric mass distributions are present (such as massive clumps and spirals in gravitationally unstable protostellar disks or massive planets and largescale vortices in protoplanetary disks), the center of mass of the star plus disk system (or its barycenter) does not coincide with the position of the star. This leads to the motion of the star in response to the net gravitational force of the disk and planets.
In the SPH method, the star is usually described as an additional pointsized particle and its motion in the Lagrangian approach does not require special attention. Conversely, the gridbased methods often use cylindrical or spherical coordinate systems, which are naturally centered on the star. If the barycenter of the star plus disk system is offset from the center of the numerical grid, the star will move in response to the net nonzero gravity force exerted on it from the rest of the system. However, the use of curvilinear coordinate systems with a singularity point in the coordinate center makes it extremely difficult to calculate directly the stellar motion and its gravitational response on the disk and planet system. To circumvent this problem, a noninertial frame of reference is often introduced, which moves with the star (positioned in the coordinate center) in response to the net gravitational force of the disk and planets. The resulting acceleration of the star is often described by the socalled indirect potential (see details in the next section).
Gridbased numerical models dealing with the gravitational instability and fragmentation of protostellar disks usually neglect the stellar motion (see, e.g., Pickett et al. 2003; Mejí a et al. 2005; Boley et al. 2006, 2007; Boley & Durisen 2008; Cai et al. 2008; Tsukamoto & Machida 2011; Boss 2012; SteimanCameron et al. 2013; Vorobyov 2013; Lin 2015). Michael & Durisen (2010) made for the first time a focused numerical investigation of the effect of stellar motion induced by nonaxisymmetry developed in a gravitationally unstable circumstellar disk. They concluded that previous findings regarding disk evolution remain valid, while the disk gravitational instability weakens slightly if the stellar motion is taken into account. Their study is, however, limited to isolated disks and neglects the effect of the infalling envelope. Therefore, their results are more applicable to protoplanetary rather than to protostellar disks, although their adopted disktostar mass ratio (≲0.1) is typical for the latter. To develop further this research, we reinvestigate the effect of stellar motion on the formation and evolution of protostellar disk.
The majority of gridbased numerical investigations dealing with the gravitational interactions in protoplanetary disks take into account stellar motion by means of the indirect potential. However, there are several works which neglect stellar motion (some examples include, but are not limited to Lubowetal 1999; Tanaka et al. 2002; Ou et al. 2007; Lyra et al 2009a,b; Zhu et al. 2011, 2014; Fung et al. 2014; Szulagyi et al. 2014; Zhu & Stone 2014). Therefore, it is also important to explore the effect of stellar motion caused by mass asymmetry in protoplanetary disks.
Mittal & Chiang (2015) found that if an artificial displacement of the barycenter of the protoplanetary diskstar system is strong enough, a largescale vortex can form via the Rossby wave instability (RWI, Lovelace et al. 1999). Recently, Zhu & Baruteau (2016) showed that such a displacement can be sustained by the vortex itself, forcing the star to wobble. They also found that neglecting the stellar motion significantly weakens the vortex. We note that Zhu & Baruteau (2016) prescribed an initially RWunstable disk by artificially setting an axisymmetric sharp Gaussian density bump in the disk. We, on the other hand, intend to investigate the effect of stellar motion on the formation of a RWunstable disk and the evolution of largescale vortices in a more selfconsistent manner.
The outline of the paper is the following. In Sect. 2, the methodologies used for calculating the indirect potential in the gridbased hydrodynamic codes are presented. Section 3 presents the effect of stellar motions on the evolution of a protostellar disk. Section 4 presents the effect of stellar motion on the evolution of a protoplanetary disk hosting a massive embedded planet and a largescale vortex formed at the viscosity transition. The paper ends with our conclusions on the importance of stellar motion in Sect. 5.
2. Simulation methodology
To investigate the effect of stellar motion on distinct evolutionary stages of accrection disks formed around young stellar objects, we employ two independent numerical hydrodynamics codes. The reason for using two codes (FEoSaD and GFARGO) is that they are specifically designed to study the protostellar (early) and protoplanetary (late) stages of disk evolution, respectively. On one hand, the FEoSaD code (Formation and evolution of a star and its circumstellar disk) is best suited for describing the early evolution of young stellar objects, including important interconnections and feedback mechanisms between the central star, its disk and infalling envelope, such as the mass infall onto the disk, disk selfgravity, stellar irradiation, and disk radiative cooling. On the other hand, the GFARGO code, the Graphical Processor Unit (GPU) accelereated version of the FARGO code, is specifically designed for simulating the planetdisk interactions, utilizing the Fast advection in rotating gaseous objects (FARGO) algorithm (Masset 2000). If the radial velocities are small compared to the Keplerian velocity (which is true for protoplanetary disks), the FARGO algorithm is about an order of magnitude faster than traditional advection schemes. Nevertheless, both codes bear some similarity as they both employ the thindisk approximation and a similar finitedifference solution procedure on a cylindrical numerical grid. In the following, detailed descriptions of the applied protostellar and protoplanetary disk models are given.
2.1. Protostellar disk model
The numerical model for the formation and evolution of a star and its circumstellar disk (FEoSaD) is described in detail in Vorobyov & Basu (2010, 2015) and is briefly reviewed here for the reader’s convenience. Numerical simulations start from a collapsing prestellar core of a certain mass, angular momentum, and temperature. The forming star is described by the Lyon stellar evolution code (Chabrier & Baraffe 1997; Baraffe et al. 2012), while the formation and longterm evolution of the circumstellar disk are described using numerical hydrodynamics simulations in the thindisk limit. The evolution of the star and the disk are interconnected: the stellar evolution code takes the mass accretion rate onto the star provided by hydrodynamic simulations and returns the stellar photospheric and accretion luminosities, which are used to calculate the disk heating rate by stellar irradiation.
The main physical processes taken into account when modeling the disk formation and evolution include viscous and shock heating, irradiation by the forming star, background irradiation (T_{bg} = 10 K), radiative cooling from the disk surface, and selfgravity. The code is written in the thindisk limit, complemented by a calculation of the vertical scale height Z using an assumption of local hydrostatic equilibrium. The resulting model has a flared structure, which guaranties that both the disk and envelope receive a fraction of the irradiation energy from the central protostar. The pertinent equations of mass, momentum, and energy transport are (1)
where subscripts p and p′ refer to the planar components (r,φ) in polar coordinates, Σ is the mass surface density, e is the internal energy per surface area, is the vertically integrated gas pressure calculated via the ideal equation of state as with γ = 7/5, is the velocity in the disk plane, and is the gradient along the planar coordinates of the disk. The gravitational acceleration in the disk plane, , takes into account selfgravity of the disk, which can be calculated by solving the Poisson integral (see details in Vorobyov & Basu 2010), and the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor Π, which is provided in Vorobyov & Basu (2010). We parameterize the magnitude of kinematic viscosity ν using the alpha prescription of Shakura & Sunyaev (1973) with a spatially and temporally uniform α_{visc} = 5 × 10^{3}. The cooling and heating rates, Λ and Γ, take the blackbody cooling and heating of the disk due to stellar and background irradiation into account (see details in Vorobyov & Basu 2010).
To explore the effects of stellar motion, Eq. (2) is formulated in the noninertial frame of reference moving with the star in response to the gravity force of both the disk and envelope. The acceleration of the star’s frame of reference g_{∗} can be expressed as (4)where dm is the mass in a grid cell with position vector r^{′}. In practice, we find g_{∗} by first calculating its Cartesian components g_{∗ ,x} and g_{∗ ,y} as where φ_{k} is the azimuthal angle of the grid cell (j,k). The summation is performed over all grid zones and the force (per unit stellar mass) acting from the grid cell (j,k) onto the star can be expressed in the following form (7)where m_{j,k} is the mass in the grid cell (j,k) and r_{j} is the radial distance to the grid cell (j,k).
Once the Cartesian components g_{∗ ,x} and g_{∗ ,y} of the stellar acceleration are known in every grid cell, the corresponding polargrid components can be found using the standard coordinate transformation formula where φ is the azimuthal coordinate of a given grid cell.
When the gravitational potential Φ is used in Eq. (2), rather than the gravitational acceleration g, it is convenient to introduce the socalled indirect potential as Φ_{ind} = r·g_{∗}. The gradient of Φ_{ind} equals g_{∗} because the latter depends on r′ rather than on r (see Eq. (4)). The last term in the right hand side of Eq. (2) can then be substituted with the following term: −Σ∇(Φ + Φ_{ind}), and the indirect potential in a grid cell (j,k) can be calculated as: (10)where summation is performed over all grid zones. In this paper, we use −∇(Φ + Φ_{ind}) instead of the acceleration g−g_{∗} in Eq. (2).
Equations (1)–(3) are solved using the method of finite differences with a timeexplicit, operatorsplit solution procedure in polar coordinates (r,φ) on a numerical grid with 512 × 512 grid zones. The radial grid is logarithmically spaced, while the azimuthal grid is equispaced. Advection is treated using the thirdorder piecewise parabolic interpolation scheme (Colella & Woodward 1984). The update of the internal energy per surface area due to cooling and heating is done implicitly using the NewtonRaphson method of root finding, complemented by the bisection method where the NewtonRaphson iterations fail to converge.
To avoid too small time steps, we introduce a “sink cell” at r_{sc} = 6.0 au and impose a free outflow inner boundary condition so that the matter is allowed to flow from the inner active computational region into the sink cell, but is prevented from flowing in the opposite direction. At the outer computational boundary a free inflowoutflow boundary condition is imposed, so that the matter is allowed to flow in and out of the active computational domain. This is done to prevent the artificial accumulation of matter near the outer computational boundary in the noninertial frame of reference (but see discussion in Sect. 3.1). More details on the code and the solution procedure can be found in Vorobyov & Basu (2010, 2015). The parameters of the initial model for studying the evolution of protostellar disks (model 1) are given in Sect. 3.1.
2.2. Protoplanetary disk model
The effect of stellar motion in a protoplanetary disk can be important if the latter has a largescale asymmetry. The sources of these asymmetries can be eccentric gaps opened by a giant planet (Kley & Dirksen 2006; Regály et al. 2010), global disk eccentricities excited due to the gravitational perturbation of a massive perturber (Lubow 1991; Kley et al. 2008) or a giant planet (Regály et al. 2014), and largescale vortices excited by the Rossby wave instability triggered by the gapopening planet (de ValBorro et al. 2007) or by the viscosity transition (Li et al. 2005). Migration of planets can also be affected by stellar motion as it is governed by the angular momentum exchange between the planet and the disk.
We investigate the gap and disk eccentricity formation caused by an embedded nonmigrating giant planet on a circular orbit (models 2 and 3), the type II migration of highmass planets (models 4 and 5), and the formation of a largescale vortex in the vicinity of a sharp viscosity transition (model 6), where the outer edge of the socalled dead zone is assumed to be located. To this end, we employ the twodimensional locally isothermal hydrodynamical simulations with and without taking the effect of the indirect potential term into account.
For the hydrodynamical modeling we use the GPU based version of the FARGO code (Masset 2000), which solves the vertically integrated continuity and NavierStokes equations (Eqs. (1) and (2)) on the cylindrical coordinate system using a locally isothermal approximation. In our simulations, the star is located at the center of the numerical domain. The indirect potential has two sources: 1) the indirect potential of the star due to the shift of the barycenter by the disk material (11)similar to what is described in Sect. 2.2; and 2) the indirect potential caused by the shift of the barycenter due to a massive planet (12)In Eqs. (11) and (12) m_{j,k} and x_{j,k}, y_{j,k}, are the mass and Cartesian coordinates of the grid cell j,k, while M_{p} and x_{p}, y_{p} are the mass and Cartesian coordinates of the planet, respectively. In simulations which do not take the effect of the indirect potential into account, we simply neglect the calculation of and .
In the local isothermal approximation, the radial temperature profile is T(R) ~ R^{1}, and the equation of state of the gas, , depends only on the density and the local sound speed c_{s}(R) = Ω(R)H(R), where H(R) is the local pressure scale height. We note that both c_{s}(R) and H(R) remain constant in time due to the locally isothermal assumption. The disk is assumed to be flat, for which case H(R) = hR, where h is the disk aspect ratio, which is set to h = 0.05 for all models. The disk’s viscosity is approximated by the αprescription of Shakura & Sunyaev (1973) with 5 × 10^{5} ≤ α_{visc} ≤5 × 10^{3}.
We adopt 1 au for the unit of length and 1.0 M_{⊙} for the unit of mass. The unit of time is such that the orbital period of a planet at 1 au is 2π. With this scaling the gravitational constant equal unity. For models 2–5 we use a frame that corotates with the planet.
We use logarithmically distributed N_{R} = 256 radial and N_{φ} = 512 equidistant azimuthal grid cells for all protoplanetary disk models. To model the eccentricity excitation of the inner disk, the relatively close inner boundary is required, thus our computational domain extends from R_{in} = 0.3 au to R_{out} = 30 au in models 2–5 (see, e.g., Regály et al. 2014). For model 6, we use a slightly larger disk inner and outer domain boundary, that is, R_{in} = 3 au and R_{out} = 50 au. The larger inner boundary is required in order to model the longterm evolution of largescale vortices. The larger outer disk boundary is necessary to avoid artificial boundary effects due to the vicinity of the disk’s outer edge to the largescale vortex.
To take the disk thickness into account, we use ϵH(R_{p}) as the smoothing of the gravitational potential of the planet with ϵ = 0.6 in protoplanetary disk simulations following the receipt in Kley et al. (2012) and Müller et al. (2012). In the planetbearing disk models, a circumplanetary disk forms around the giant planet, which affects its migration history (models 4 and 5). To circumvent this effect we remove a certain amount of the torque exerted on the planet by the disk according to Crida el al. (2009) in models 4 and 5. Thus, the force exerted on the planet by the material inside the planetary Hill sphere is multiplied by 1−exp(−(d/R_{H})^{2}), where R_{H} = a_{p}(q/ 3)^{1/3} is the planetary Hill radius and d is the distance between the given grid cell center and the planet.
To speedup the simulations, the damping boundary condition (de ValBorro et al. 2006) is applied at the inner boundary. The damping criteria at the inner domain edge is set such that the waves are killed close to the inner boundary R < 1.2 R_{in}. With this setting, the time step (which is calculated according to the CourantFriedrichLevi criteria described in Masset 2000) remains considerably large because sharp density enhancements (i.e., shock waves) cannot appear at the boundary. The outflow boundary is used at the disk’s outer edge in all protoplanetary disk models.
3. Early disk evolution (protostellar disk)
In this section, we investigate the effect of stellar motion on the early evolution of protostellar disks, when they are still actively accreting material from the parental cores. For this purpose, we use the FEoSaD code described in Sect. 2.1. We start our numerical simulations from the gravitational collapse of a starless cloud core, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations when the age of the star exceeds 0.4 Myr.
3.1. Initial conditions
For the initial conditions, we consider a gravitationally unstable prestellar core submerged into a lowdensity environment, both described by the following radial gas surface density distribution (15)where Σ_{0} is the gas surface density at the center of the core and the radius of the central plateau, c_{s} the initial sound speed in the core, R_{core} the radius of the core, and Σ_{ext} the gas surface density in the external environment. The gas surface density distribution of the core can be obtained (to within a factor of unity) by integrating the 3D gas density distribution characteristic of BonnorEbert spheres with a positive densityperturbation amplitude A (Dapp & Basu 2009). The value of A is set to 1.2. The core has a fixed ratio R_{core}/r_{0} = 6.0, implying that it is initially unstable to gravitational collapse. The initial temperature is set everywhere throughout the core to 10 K. For the chosen core radius of R_{core} = 0.09 pc and central surface density of Σ_{0} = 2 × 10^{2} g cm^{2}, the resulting initial core mass is M_{core} = 1.18M_{⊙}. The total computational domain extends to 0.18 pc and, for the value of Σ_{ext} = 6.4 × 10^{5} g cm^{2}, the total mass contained in the external environment is M_{ext} = 0.02 M_{⊙}, which is only about 1% that of the core mass.
The initial angular velocity is described by the following equation (16)as expected for slowly contracting prestellar cores with conservation of angular momentum (Basu 1997). Here, Ω_{0} = 0.49 km s^{1} pc^{1} is the angular velocity in the center of the core. For the chosen parameters of our model, the ratio of rotational to gravitational energy is β = 5.5 × 10^{3}. We hereafter refer to this model as model 1.
The reason for adding an external constantdensity environment to the prestellar core (rather than just considering an isolated core) is that the shape of the core near its outer edge can be artificially distorted in the noninertial, accelerated frame of reference. We employed the logarithmically spaced grid in the radial direction, which resulted in a good numerical resolution in the inner regions where the protostellar disk forms and evolves, but in a poor numerical resolution near the outer boundary. The imperfections in the implementation of the outer computational boundary lead to artificial accumulation or depression in the gas density near the boundary when the material in the computational domain is accelerated due to the indirect potential. This artificial distortion of the core shape often has an undesirable feedback effect on the disk dynamics. We experimented with different outer boundary conditions (free boundary, extrapolated density and velocity profiles), but the result was not satisfactory.
The development of special numerical schemes that allow for the moving boundaries (so that the flow of material through the outer boundary caused by the acceleration of the coordinate system is minimized) could potentially help to alleviate the problem. Our temporary solution to this problem was to extend the outer boundary sufficiently far from the disk and fill in the rest of the numerical grid with a lowdensity gas so that its gravitational potential would have little effect on the dynamics of the system, even if this external environment becomes distorted due to the boundary effects. As we show later in the paper, this helps to mitigate the problem for the initial several hundred thousand years of protostellar disk evolution. However, because of the pressure disbalance between the core and the external environment, the outer regions of the core gradually expand (while the inner part continues collapsing onto the star) and ultimately reach the outer computational boundary. At this time instance, we have to stop the simulations.
3.2. The effect of stellar motion on the protostellar disk and envelope
Fig. 1 Gas surface density distribution in model 1 in the inner 1400 × 1400 au box focusing on the disk. Six evolution times counted from the formation of the star are shown. The scale bar is in g cm^{2} (log units). 

Open with DEXTER 
In this section, we focus on possible differences in the properties of the disk and infalling envelope that may arise when stellar motion is taken into account. We differentiate between models with and without stellar motion by adding the plus sign to the corresponding model number, that is, model 1+ means that the stellar motion is taken into account, while model 1 means the opposite. The chosen initial parameters of the prestellar core (M_{core} = 1.18 M_{⊙} and β = 5.5 × 10^{3}) imply that the core collapse will lead to the formation of a gravitationally unstable disk (Vorobyov 2013).
Figure 1 presents the gas surface density in model 1 in the inner 1400 × 1400 au^{2} box at several time instances after the formation of the central star. We note that the total computational domain is much greater, but for now we focus only on the innermost regions occupied by the disk. Gravitational instability develops in the disk soon after its formation due to the continuing mass loading from the parental core. The spiral structure has a complicated pattern and dense gaseous clumps form within the spiral arms. The vigorous gravitational instability and fragmentation continues for at least 0.3 Myr until the gas reservoir in the parental core exhausts and the infalling envelope (that continuously replenishes the disk material lost via accretion onto the star) dissipates.
Fig. 2 Similar to Fig. 1, but for model 1+ for which case stellar motion is included. 

Open with DEXTER 
Figure 2 shows the gas surface density in model 1+ at the same time instances as in model 1. The stellar motion is introduced 2 × 10^{4} yr after the instance of disk formation^{1}, or 2.8 × 10^{4} yr after the formation of the central star; the corresponding gas surface density distribution is shown in the upperleft panel in Fig. 2. This delay is introduced to let the disk grow and acquire some material. Otherwise, we found that the disk becomes significantly distorted in the initial formation stage, leading to the termination of numerical simulations. This problem is likely to be related to the imperfections of the inner computational boundary (near to which the disk naturally forms) and requires further investigation.
A visual inspection of Figs. 1 and 2 indicates that the stellar motion has a notable effect on the disk radius and the strength of gravitational instability: the disk radius is smaller and the disk fragmentation is less frequent in model 1+ than in model 1. To quantify these differences, we calculated the number of fragments in the disk every 2000 yr using the fragment tracking mechanism explained in detail in Vorobyov (2013). This algorithm requires that the fragment be pressuresupported, having a negative pressure gradient with respect to the center of the fragment, and that the fragment be kept together by gravity, having the deepest potential well at the center of the fragment. The results are shown in Fig. 3 in the top panel (model 1) and bottom panel (model 1+). While model 1 forms fragments almost continuously starting from t = 0.03 Myr until t = 0.27 Myr, model 1+ exhibits only sporadic fragmentation episodes. The maximum number of the fragments in the disk at any given time is also greater in model 1 (N_{fr,max} = 7) than in model 1+ (N_{fr,max} = 4).
To understand the reason for this difference, we plot the total disk masses, stellar masses, and disk radii as a function of time in Fig. 4. In particular, the red and black lines show the data for models with and without stellar motion, respectively. The black dashdotted lines correspond to the envelope mass, which is similar in both models. Evidently, model 1+ has a somewhat smaller and less massive disk than model 1. We think this is caused by the fact that, in the model with stellar motions, part of the total angular momentum of the collapsing core goes into the orbital motion of the star around the center of mass of the whole system. On the other hand, in the model without stellar motion, all angular momentum of the collapsing core is transferred to the disk. As a result, the disk in model 1 has more angular momentum and a greater support against gravity than in model 1+, helping the former to sustain a more massive and extended disk. This conjecture is supported by direct calculation of the specific angular momentum of the disk: it is equal to J_{disk} = 2.4 × 10^{25} cm^{2} s^{1} in model 1+ and J_{disk} = 3.3 × 10^{25} cm^{2} s^{1} in model 1.
Fig. 3 Number of fragments in the disk at a given time elapsed since the formation of the star. The top and bottom panels correspond to models 1 and 1+ (without and with stellar motion, respectively). 

Open with DEXTER 
Fig. 4 Top panel: stellar, disk, and envelope masses in model 1 (blackcolored lines) and model 1+ (redcolored lines). The envelope mass is the same in both models. Bottom panel: disk radius in model 1 (blacksolid line) and model 1+ (red solid line). 

Open with DEXTER 
The spiral structure in Figs. 1 and 2 is complicated and does not show any regular, longlasting pattern. This is typical for unstable protostellar disks subject to fragmentation caused by the continuous massloading from parental cores (e.g., Vorobyov & Basu 2010; Tsukamoto et al. 2013). Nevertheless, it is possible to analyze the spiral pattern by calculating the global Fourier amplitudes defined as (17)where M_{disk} is the disk mass and R_{disk} is the disk’s physical outer radius. When the disk surface density is axisymmetric, the amplitudes of all modes are equal to zero. When, say, C_{m}(t) = 0.1, the perturbation amplitude of spiral density waves in the disk is 10% that of the underlying axisymmetric density distribution.
Figure 5 presents the time evolution of the global Fourier amplitudes for the first four modes in model 1+ (top panel) and model 1 (bottom panel). The m = 1 spiral mode clearly dominates in model 1+, except for a short time period around t = 0.14–0.2 Myr, when all the four modes have similar amplitudes. In model 1, on the other hand, the dominance of the m = 1 mode is only marginal. At the later times, the amplitudes in model 1+ are somewhat smaller than in model 1, which is consistent with a lower disk mass in model 1+. The dominance of the m = 1 mode in the model with stellar motion may be related to the SLING (stimulation by the longrange interaction of Newtonian gravity) effect discussed in, for example, Shu et al. (1990). We note that Michael & Durisen (2010) found no SLING effect in their numerical hydrodynamics simulations with stellar motion, possibly due to the fact that they considered an isolated disk with a smaller disktostar mass ratio.
Fig. 5 Global Fourier amplitudes of the disk as a function of time in model 1 (top panel) and in model 1+ (bottom panel) calculated for the first four spiral modes. 

Open with DEXTER 
Stellar motion has a notable effect not only on the disk, but also on the structure of the infalling envelope – the remnant of the parental prestellar core. Figure 6 presents a zoomedout view of the computational domain, this time covering the 16 000 × 16 000 au^{2} area which includes both the disk and infalling envelope. The disk occupies the inner several hundred au, while the rest is occupied by the envelope. The isosurface density contour lines are added to the gas surface density maps to better characterize the density structure of the envelope. At the very early stage (t = 0.03 Myr), just after the stellar motion has been turned on, the envelope has an axisymmetric shape consistent with the initial configuration. However, already at t = 0.08 Myr the deviation of the inner envelope structure from a purely axisymmetric shape becomes evident. This is because the stellar motion has a feedback effect not only on the disk, but also on the envelope, at least on its inner part. However, after t ≈ 0.2 Myr the expansion wave caused by the pressure imbalance at the interface between the core and the external environment reaches the outer computational boundary and the boundary effects discussed in Sect. 3.1 start to artificially distort the shape of the whole envelope, as is seen in the t = 0.28 Myr panel. At this stage, we terminate the simulations.
Fig. 6 Gas surface density distribution in the inner 16 000 × 16 000 au^{2} box covering both the disk and the infalling envelope in model 1+ at six times after the formation of the central star. The yellow isosurface density contour lines are added to highlight the inner envelope structure. The scale bar is in log g cm^{2}. 

Open with DEXTER 
To check if the deviation of the inner envelope from a purely axisymmetric configuration is indeed caused by stellar motion, we plot in Fig. 7 the gas surface density distribution in model 1 (without stellar motion) for the inner 8000 × 8000 au^{2} box. The yellow isosurface density contour lines highlight the envelope structure. Evidently, the envelope is axisymmetric, except for very inner regions where it meets with the nonaxisymmetric disk.
Fig. 7 Similar to Fig. 6, only for model 1 without stellar motion. Only four evolutionary times are shown. 

Open with DEXTER 
Figure 8 gives the position of the center of mass (CoM) of the whole system (star+disk+envelope+external environment) on the polar grid (r,φ). The CoM shows large excursions from the initial position (r = 0,φ = 0), as well as smallscale wobbling. The large excursions are the result of the gradual transformation in the structure of the infalling envelope from a purely axisymmetric to largely nonaxisymmetric shape, while the small wobbling is the reaction to the nonaxisymmetric gravitational perturbation of the disk. Since the mass contained in the envelope is significant (it is always greater than the disk mass and is also greater than the stellar mass in the early evolution phase, see Fig. 4) the effect of the nonaxisymmetric inner envelope on the stellar motion is stronger than that of the disk. The spiral arms and fragments in the disk often tend to cancel the effect of each other on the star, which also acts to reduce the stellar wobbling due to the disk.
Fig. 8 Position of the center of mass in model 1+ at the polar coordinate system. The center of mass is initially at (r = 0,φ = 0), but shows significant excursions during the subsequent evolution due to the nonaxisymmetric structure of the infalling envelope. Smallscale wobbling in response to the nonaxisymmetric disk is also evident. 

Open with DEXTER 
Finally, we note that in our adopted thindisk approximation all the disk and envelope mass is concentrated in the midplane. The immediate consequence of this assumption is that any nonaxisymmetry in the disk or envelope has a stronger effect on the star as compared to the real 3D case. In the 2D case, the gravitational force has only two plane components, while in the 3D case there is also a vertical component. The latter would contribute to the vertical excursions of the star (rather than radial ones due to the plane components) and would vanish in systems with equatorial symmetry. This problem can be partly solved in the 2D case by splitting the disk vertical column into i layers of equal width and calculating the input from every vertical layer individually. The resulting force acting upon the star from the grid cell (j,k) can than be expressed as (18)where the vertical height of each layer is denoted by z_{i} and the summation is performed over all vertical layers. The mass of each layer m_{i,j,k} can be found by assuming a constant or Gaussian distribution of the gas volume density within the vertical column. A factor of two expresses the fact that the disk has an equatorial symmetry. Our test calculations indeed indicate that the resulting envelope structure is less distorted and the problem of the outer computational boundary becomes less severe, though not solved completely. From a general point of view, we note that the gravitational instability in the 3D disks is expected to be weaker than in the 2D ones, which may also reduce the gravitational effect of the disk upon the star. Other 3D effects, such as disk warps, may reduce the radial excursions of the star, but introduce notable vertical motions.
3.3. Protostellar accretion rates
Fig. 9 Mass accretion rates as a function of time in model 1 (top panel) and model 1+ (bottom panel). 

Open with DEXTER 
In this section, we investigate the effect of stellar motion on the mass accretion rate defined as the mass passing through the inner computational boundary per hydrodynamical time step, Ṁ = −2πr_{sc}Σv_{r}. Previous numerical simulations without stellar motion demonstrated that Ṁ often exhibits timevariable behavior with episodic bursts, if the disk mass is sufficient to trigger gravitational instability and fragmentation (e.g., Vorobyov & Basu 2006, 2010; Machida et al. 2011; Vorobyov & Basu 2015). This accretion pattern is part of the emerging episodic accretion paradigm, wherein the mass accretion onto the star is not constant or gradually declining in time, as predicted by the classic spherical collapse models (e.g., Shu 1977), but is timevariable with strong accretion bursts separated by longer periods of quiescent, lowrate accretion (see, e.g., a recent review by Audard et al. 2014). It is therefore important to assess the effect of stellar motion on the mass accretion rate.
Figure 9 presents the mass accretion rates as a function of time in models 1+ (top panel) and model 1 (bottom panel). In the early evolution t ≲ 0.1 Myr, the accretion rate in models with and without stellar motion exhibits a similar behavior: it quickly rises to a peak value at t = 0 Myr (the instance of protostar formation), declines gradually for a short period of time (≈0.01 Myr), during which the matter is directly accreted from the infalling envelope onto the protostar, and then develops orderofmagnitude variations when the disk around the protostar forms and gains mass. Accretion bursts amounting to Ṁ = a few × 10^{5}M_{⊙} yr^{1} are evident in both models.
In the later evolution (t ≳ 0.1 Myr), the mass accretion rate in model 1+ is characterized by a somewhat reduced variability in comparison to model 1, though still showing several orderofmagnitude accretion bursts. This reduction in the amplitude of accretion variations is consistent with a lower disk mass in model 1+. Overall, the stellar motion somewhat reduces the accretion variability, particularly in the later evolution stages, mainly due to a reduced disk mass. This reduction can, however, be offset by a small increase in the initial core mass and/or angular momentum (which would result in the formation of a more massive and gravitationally unstable disk).
4. The late disk evolution (protoplanetary disk)
In this section, we investigate how stellar motion affects the evolution of protoplanetary disks perturbed by a massive planet, by the migration of planets, and by the formation of largescale vortices in protoplanetary disks. We assume that the protoplanetary disk has already been formed and an embedded highmass planet revolves on a circular orbit (models 2–5). In simulations dealing with the formation of vortices at the dead zone outer edge, we impose a sharp viscosity transition (model 7). Similar to the previous section, we run protoplanetary disk simulations with and without stellar motion. Models in which we take the effect of the indirect potential into account are labeled with an additional plus sign.
The simulations span 4 × 10^{4} yr (corresponding to about 3500 orbits at 5 au) in models 2–5. To study the vortex formation, much longer simulations are required. Therefore, we run simulations up to 2 × 10^{5} yr in model 6, corresponding to about 4811 orbits at 12 au – a distance where the largescale vortex develops.
4.1. Initial conditions
Fig. 10 Density snapshots in models 2 and 3 for which case a nonmigrating giant planet is embedded in a disk with α_{visc} = 5 × 10^{3} (panel a)) and 5 × 10^{5} (panel b)). In each panel the left and right columns correspond to simulations in which the indirect potential is taken into account (model 2+ and 3+) and neglected (model 2 and 3), respectively. Snapshots are taken after the gap is opened at t = 10^{3} yr corresponding to about 100th orbits (upper row), while at the end of simulation at t = 4 × 10^{4} yr corresponding to about 3600th orbits (lower row). 

Open with DEXTER 
The initial surface density profile of the gas has the following form: Σ(R,φ) = Σ_{0}R^{− p}, where p is set to 1.0 and Σ_{0} is set to 1.06815 × 10^{4}, which corresponds to a disk of 0.02 M_{⊙} inside the computational domain. Since the Toomre Q parameter (Toomre 1964) is hM_{∗}/πR^{2}Σ(R,φ) > 1 throughout the disk in all models, it seems plausible to neglect the disk selfgravity^{2}. We note, however, that the disk selfgravity affects the migration rate of planets; it slightly slows down the migration (see, e.g., Crida el al. 2009) even for lowmass disks where selfgravity is weak. The vortex formation at the planetary gap edge is also known to be affected by the disk selfgravity, which shifts the most unstable modes to higher mode numbers, or even suppresses the RWI for sufficiently high disk masses with Q ≤ 2.5 (see, e.g., Mamatsashvili & Rice 2009; Lin & Papaloizou 2011). Recently, Zhu & Baruteau (2016) showed that including selfgravity results in significantly weakened largescale vortices if Q ≤ π/ (2h). Here we want to investigate the effect of the indirect term only, which requires longterm simulations. Therefore, the computationallydemanding selfgravity calculation is neglected for simplicity.
In models 2 and 3, we study the disk response to a nonmigrating highmass planet (q = 4.7 × 10^{3}, corresponding to M_{p} = 5 M_{Jup} for a solar mass star). In model 2 the viscosity parameter is set to a canonical value of α = 5 × 10^{3}. Model 3 is a nearly inviscid version of model 2 with α = 5 × 10^{5}. With these settings, we expect the formation of various types of disk asymmetry, e.g., a local disk eccentricity at the gap edge (see, e.g., Kley & Dirksen 2006; DAngelo et al. 2006; Regály et al. 2010), global disk eccentricity inside the planetary orbit (see, e.g., Regály et al. 2012, 2014), and largescale vortex formation at the outer gap edge (see, e.g., Regály et al. 2012). We note that the low viscosity assumption is necessary to maintain a longlasting vortex at the planetary gap edges as was shown earlier by de ValBorro et al. (2007), Ataiee et al. (2013), Fu et al. (2014).
In models 4 and 5, the migration of highmass planets (q = 9.54 × 10^{4} and 2.38 × 10^{3}, corresponding to M_{p} = 1.0 and 2.5 M_{Jup}) is investigated. With these planetary masses, no excitation of disk eccentricity is observed, for which case the migration history might be more complicated. The planets are freely migrating in the disk (after the gap is fully fledged), having a canonical viscosity of α = 5 × 10^{3}. With these models, we investigate the effect of stellar motion on the migration history of highmass planets.
In models 2–5, initially the planet orbits at a_{p} = 5 au and is not allowed to accrete gas, that is, its mass is constant. For nonmigrating cases (models 2 and 3), the planetary orbit is kept fixed. For migrating planets (models 4 and 5), the planet is allowed to feel the disk potential after the gap is fully fledged.
To model the formation and evolution of the largescale vortex developed due to the RWI in model 6, we introduce a sharp viscosity transition at the dead zone outer edge. We assume that the disk has an accretionally inactive dead zone, where the value of α is smoothly reduced such that α(R) = αδ_{α}(R). The viscosity reduction is given by (19)where α_{mod} = 0.01 is the depth of the turbulent viscosity reduction. To quantify the radius of viscosity reduction, R_{dze} = 12 au is used, where we adopt the results of Matsumura & Pudritz (2005), who found that R_{dze} lies between 12 au and 36 au, depending on the density of the disk. We assume ΔR_{dze} = 1H_{dze}, where H_{dze} = R_{dze}h is the disk scaleheight at the viscosity reduction, which corresponds to ΔR_{dze} = 0.6 au. We note that the total width of the viscosity transition given by Eq. (19) is about 2ΔR_{dze}, which corresponds to 1.2 au.
4.2. Disk response to a giant planet (models 2 and 3)
We investigate the effect of stellar motion on the disk response to the gravitational perturbation of a giant planet embedded in the disk in models 2 and 3. The viscous gap opening criterium is satisfied (i.e., ) in models 2 (viscous) and 3 (nearly inviscid) for a five Jupiter mass planet (Lin & Papaloizou 1993). As a result, a gap forms after a couple of orbits, independent of whether stellar motion is taken into account or not.
The upper and lower rows of Fig. 10 show the gap structure at t = 10^{3} yr and at the end of simulations t = 4 × 10^{4} yr, respectively. In the early phases, the depletion of the horseshoe region is significantly slower for models without indirect potential (compare upper rows in panels a and b of Fig. 10). In the later phases, the gap width is larger in models where the stellar motion is neglected. From the azimuthally averaged radial density profiles (Fig. 11) it is clearly visible that the inner gap edge forms nearly at similar distance, while the outer gap edge is always at larger distances in models without stellar motion.
Fig. 11 Azimuthally averaged radial density profiles in models 2 and 3 at the end of simulations (corresponding to the lower row of Fig. 10). Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER 
Another notable difference can be identified regarding the largescale vortex formation at the gap edges. The largescale vortex develops due to the excitation of the RWI at the gap edges thanks to the development of a sharp density gradient (Li et al. 2005; Lyra et al 2009b; Regály et al. 2012). However, due to the viscous evolution of the gap edge, the density jumps smear out and the gap edge is not RWunstable any more. As a result, the vortex dissolves. We observed this phenomenon regardless of whether the stellar motion is taken into account or not. We emphasize that the density gradient at the outer gap edge smears out faster in model 2 without stellar motion. This explains why no vortex is present after ten orbits in the high viscosity models 2+ and 2 (upper figures in panel a of Fig. 10) and is present only in model 3+ (upper figures in panel b of Fig. 10). We note, however, that the vortex lifetime is observed to be shorter than 100 orbits even in the lowviscosity regime, independent of whether the stellar motion is on or off.
As can be seen in Fig. 10, a significant level of disk eccentricity develops by the end of the simulations in both models. Not only the gap edge but also the inner disk become globally eccentric, as was shown in Regály et al. (2014). One can calculate the disk eccentricity at each grid cell using Eqs. (5)–(7) of Regály et al. (2010). After azimuthal averaging, this results in the radial eccentricity profile of the disk shown in Fig. 12. Evidently, the eccentricity in the inner disk (with respect to the planetary orbit) is nearly independent whether we take the stellar motion into account or not. However, the outer disk eccentricity, especially close to the gap outer edge, tends to be larger for models without stellar motion. Since the gap width is inferred from the azimuthally averaged radial density profile (Fig. 11), the somewhat larger disk eccentricity (apparent in the vicinity of the gap) can explain the larger gap.
Radial averaging of the eccentricity profiles results in a sole number which represents the disk global eccentricity. It is clearly visible on the plot of the disk global eccentricity as a function of time (Fig. 13) that a quasistatic eccentric disk state with e_{disk} ≃ 0.15 is developed by the end of the simulations, independent of whether the stellar motion is on or off. However, the eccentricity growth timescale is smaller for models with stellar motion (models 2+ and 3+).
Fig. 12 Azimuthally averaged eccentricity profile of the disk at the end of the simulations in models 2 and 3. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER 
Fig. 13 Evolution of averaged disk eccentricity in models 2 and 3. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER 
4.3. Migration of giant planets (models 4 and 5)
4.3.1. On migration timescales
The effect of stellar motion on the radial migration of highmass planets is investigated in models 4 and 5 assuming α = 5 × 10^{3}. In these models, the embedded planets with q = 9.54 × 10^{4} and 2.38 × 10^{3} (corresponding to M_{p} = 1 M_{Jup} and 2.5 M_{Jup}, respectively) are initially on circular orbits. According to our simulations, no significant planetary and disk (local or global) eccentricities develop for these planetary masses throughout the simulations. The planets are allowed to migrate only after the formation of a fully fledged gap (after about 100 orbits) in order to model type II migration exclusively. Since the planets are not allowed to accrete gas, their mass remains constant throughout the simulation. The migration history is shown in Fig. 14 for models 4+/4, and 5+/5.
After opening the gap, planets start migrating inwards in all models. The ratio of the local disk mass (which is the mass of the disk inside the planetary orbit) to the planetary mass is , where R_{0} = a_{p} + 2.5R_{H}(a_{p}). Since M_{disk}(R_{0}) /M_{p}> 1 as long as a_{p}> 1.5, the planets should migrate in the diskdominated type II regime according to Crida & Morbidelli (2007). Assuming that there is now gas flow across the planetary orbit, the change in the planetary semimajor axis as a function of time can be expressed as (20)where ν(R) = αc_{s}(R)^{2}/ Ω(R). Therefore, the migration rate is independent of the planet mass and the planet migrates on the viscous timescale. However, this is not the case and Fig. 14 demonstrates a clear dependence of the migration rate on the planetary mass. Duffell et al. (2014) have shown that the gas can cross the planetary orbit, for which case the above simple description of type II migration is not valid. In agreement with Dürmann & Kley (2015), we also found that the gapopening planet does not migrate at the viscous timescale and its migration rate is inversely proportional to its mass.
Moreover, a significant effect of stellar motion on planetary migration can be observed. If the indirect potential is included, the migration rate decreases significantly. However, as the planet’s orbital distance decreases the migration rate becomes less dependent on whether stellar motion is taken into account or not.
4.3.2. On migration torques
Fig. 14 Migration history of highmass planets in models 4 and 5. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER 
Let us analyze the torques exerted on the migrating planet. The change in the semimajor axis a_{p} of the planet having a mass q due to a small change in its angular momentum can be expressed as (e.g., Vorobyov 2013) (21)where is the angular velocity of the planet and Γ is the torque exerted on the planet. The torque Γ has two sources: one coming from the nonaxisymmetric disk and the other from the star itself if stellar motion is allowed.
The disk torque can be calculated numerically by summing the torques exerted on the planet by the individual grid cells, that is, (22)where Fy_{i,j} and Fx_{i,j} are the x and ycomponents of the gravity force felt by the planet from the grid cells with coordinates x_{i,j} and y_{i,j}. Here, x_{p}, y_{p} and x_{bc}, y_{bc} are the coordinates of the planet and the barycenter of the starplanetdisk system, respectively. The components of the gravitational force of the disk felt by the planet can be written as: (23)and (24)where d_{i,j} = [(x_{i,j}−x_{p})^{2} + (x_{i,j}−x_{p})^{2}] ^{1/2} is the distance of the center of grid cell (i,j) to the planet. The term ϵH(d_{p}) appears due to the applied smoothing of the planetary potential. Taking into account that we also use a torque cutoff (see details in Sect. 2.2), the force exerted by the individual disk cells, F_{i,j}, can be given as (25)where A_{i,j} and Σ_{i,j} are the surface area and surface mass density of the grid cell (i,j), respectively.
If the barycenter of the stardiskplanet system is not centered on the star, a nonzero torque is exerted on the planet by the star itself (see, e.g., Ataiee et al. 2013; Regály et al. 2013). The stellar torque, Γ_{∗}, can be given as (26)where (27)and (28)Since we use a frame of reference that corotates with the planet, y_{p} becomes equal to zero and Eq. (28) vanishes, resulting in (29)Now, let us compare the total torque Γ_{tot} = Γ_{disk} + Γ_{∗} felt by the planet in models with and without stellar motion. In models without stellar motion, x_{bc} = 0 and y_{bc} = 0. In the corotating frame y_{p} = 0, thus the second term of Eq. (22) and obviously Γ_{∗}, vanish, resulting in the following expression for the total torque (30)where the index “NI” means that the indirect potential is not taken into account. However, if stellar motion is taken into account, the second term in Eq. (22) and the stellar torques are nonzero. In this case, the total torque can be given as (31)We emphasize that here we implicitly assume that the gas density distribution is independent of whether the stellar motion is on or off. With a plausible assumption of x_{bc} > 0 in our models^{3}, the second term inside the first bracket in Eq. (31) acts to slow down the planetary migration. On the contrary, the second term and the stellar torque in Eq. (31), both being dependent on y_{bc} which can be either positive or negative, can either increase or decrease the planetary migration rate. We note that in the corotating frame the value of y_{bc} is determined solely by the barycenter of the disk itself.
Fig. 15 Barycenter shift (in au) of the starplanet (x_{bcsp}, x_{bcsp}) and starplantdisk systems (x_{bc}, x_{bc}) as a function of time for models 4+ and 5+. Since the shifts are highly oscillating, the running average over 5 snapshots is also shown. 

Open with DEXTER 
Figure 15 shows the barycenter coordinates of the starplanetdisk system (x_{bc}, y_{bc}) and the barycenter coordinates of the starplanet system (x_{bcsp}, y_{bcsp}) as a function of time in models 4+ and 5+. As one can see, both the x_{bc} and y_{bc} oscillate around the barycenter of the starplanet system. The amplitude of the oscillation decreases with time and practically vanishes by t = 10^{4} yr (corresponding to about 1000 orbits of the planet at 5 au). This means that the barycenter of the disk slowly converges to the barycenter of the starplanet system, even if the planet is migrating. We emphasize that y_{bc} converges to zero, therefore both Γ_{∗} and the term proportional to y_{bc} in Eq. (31) vanish. Since x_{bc} ≃ x_{bcsp} = qx_{p}/ (1 + q), the total torque becomes , which is too small to explain the observed slowdown of migration for models 4+ and 5+. Thus, the assumption that the disk surface density distribution is the same in models 4 and 4+ and also in models 5 and 5+ must be wrong.
Fig. 16 Torques exerted on the planet in models 4/4+ (upper panel) and 5/5+ (lower panel). Significant departures can be observed only for the outer disk torque (Γ_{outer disk}), which eventually explains the decline of the total torque, Γ_{tot}, (see orange solid and dashed curves for models 4/5 and 4+/5+, respectively) in the models where stellar motion is taken into account. 

Open with DEXTER 
According to Eq. (21), the change in the planet’s semimajor axis is proportional to the normalized torque, Γ/Γ_{0}, where Γ_{0} = qa_{p}Ω(a_{p}) is the normalization factor. Figure 16 shows Γ/Γ_{0} against time in models 4/4+ and 5/5+. It is notable that Γ/Γ_{0} has a significantly larger magnitude in models 4 and 5 than in models 4+ and 5+ (see the solid and dashed running average curves), which explains the observed differences in the migration rates. Taking into account that y_{bc} ≃ 0 and Γ_{∗} ≃ 0, and also that the magnitude of the inner disk torque is not sensitive to the indirect potential, the difference in the migration rates can be explained by the difference in the outer disk torque alone. Since the magnitude of the difference in the outer disk torque decreases with time, the migration rates of models 4/4+ and 5/5+ slowly decrease.
To explain the decline of the disk torque, we plot in Fig. 17 the ratio of the surface densities in model 5+ and 5 () at the beginning of the planetary migration. The inner disk has a slight mass enhancement (~ 10%), while beyond the gap outer wall there is a significant decline (≃30%) in density ratio in models 4+/5+ compared to models 4/5. Although the difference in the density distributions weakens with time, it is still present at the end of the simulations. The outer disk removes angular moment from the planet, with the magnitude being proportional to the density. Therefore, the angular moment removal is smaller in model 4+ and 5+ than in model 4 and 5, respectively, which explains the observed migration rates.
In summary, it is neither the stellar torque nor the barycenter shift that causes a slowdown of the migration rate in models with stellar motion. Actually, the planetary migration slows down because the outer disk, with respect to the planetary orbit, contains less mass, which results in a less efficient angular moment removal from the planet.
Fig. 17 Ratio of the surface mass densities in models 4+/4 (left panel), and models 5+/5 (right panel) at the beginning of the planetary migration. The inner and outer gap edges (defined by half maximum density) are found to be nearly identical in models 4/4+ and 5/5+ shown with the dashed circles. 

Open with DEXTER 
4.4. Largescale vortex formation (model 6)
Fig. 18 Evolution of largescale vortex at viscosity transition in model 6. Left and right columns correspond to simulations in which the indirect potential is taken into account (model 6+) and neglected (model 6), respectively. The density distribution (Σ) is normalized by the initial one (Σ_{0}). 

Open with DEXTER 
The formation and evolution of a largescale vortex is investigated in model 6, in which a sharp viscosity transition is set by Eq. (19) at R_{dze} = 12 au. Due to the sudden change in the accretion rate at the viscosity transition, a sufficiently sharp density enhancement forms, where the disk becomes unstable to the RWI (Varniere & Tagger 2006; Terquem 2008). The fastest growing mode is found to be m = 5 in our models, thus initially five anticyclonic vortices form which later coalesce (Lovelace et al. 1999; Li et al. 2000). As a result, a single, largescale vortex develops at the viscosity transition.
Because the density bump formed at the viscosity transition strengthens with time, the RWI excitation criterium remains fulfilled until the end of simulation. As a result, the vortex is longlived (may be present for the disk lifetime) as was shown in Regály et al. (2012) and later confirmed by Lin (2014).
Figure 18 shows four snapshots of the disk inner region (density distribution normalized by the initial one) for model 6+ (left column) and model 6 (right column). The excitation of the RWI can be observed in both models by about 1.5 × 10^{3} yr. The vortex coagulation, that is, the formation of the m = 1 largescale vortex (upper row of Fig. 18) is faster in model 6+ than in model 6. However, a fully fledged single vortex develops in both simulations by ~ 10 × 10^{3} yr (corresponding to about 75 orbits at 12 au). Generally, the largescale vortex is stronger in model 6+ than in model 6 throughout the simulation. By the end of the simulation, the largescale vortex practically transformed to a ringlike density structure in model 6. In contrast, the vortex lasts to the end of the simulation in model 6+.
Figure 19 shows the gradient of the azimuthally averaged radial density profile (∂Σ(R) /∂R) at the end of simulation for models 6+ and 6. It is evident that the density bump is not sharp enough at R ≃ 10 au (at a radial distance of the vortex center) to maintain the vortex for model 6, which eventually causes the slow dissolution of the vortex.
Figure 20 shows the distance between the barycenter of the system and the star as a function of time. Evidently, the barycenter shift is significant in model 6+, whereas it is much smaller in model 6, reaching its maximum (~ 0.014 au) at 5 × 10^{4} yr. It is also evident that the magnitude of the barycenter shift slowly decreases after reaching its maximum in model 6, whereas it continuously increases in model 6+. We note, however, that the effect of the barycenter shift in model 6 has no effect on the evolution of the vortex; it only reflects the asymmetry in the density distribution.
We note that Mittal & Chiang (2015) proposed that a largescale vortex can be sustained in massive disks due to an artificial displacement of the star from the barycenter of the system. However, they applied an indirect potential artificially. Recently, Zhu & Baruteau (2016) found that the indirect potential indeed promotes the formation of a largescale vortex, however, in their simulation an artificially imposed axisymmetric Gaussian density bump served as the background of the density gradient that excites the RWI. Here we not only confirm the above findings but also observe that the stellar motion promotes the excitation of the RWI.
Fig. 19 Density gradient profile at the end of simulations for models 6+ and 6. 

Open with DEXTER 
Fig. 20 Barycenter shift against time in models 6+ and 6. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER 
5. Conclusions
In this paper, we presented a detailed investigation of the effect of stellar motion on the evolution of protostellar and protoplanetary disks by means of 2D numerical hydrodynamic simulations. The stellar motion is caused by various nonaxisymmetric structures in the disk and infalling envelope, which emerge during the early and late evolutionary phases. Two distinct codes were used: FEoSaD for protostellar and GFARGO for protoplanetary disks. Both codes employ the gridbased cylindrical coordinate system centered on the star. The motion of the central star is taken into account by calculating the indirect potential, which occurs when the equations of motion for the disk are formulated in the noninertial frame of reference moving with the star.
In the protostellar disks, asymmetries can develop due to disk gravitational fragmentation resulting in the formation of massive clumps and spiral arms. In the protoplanetary disks, asymmetries can develop due to an embedded planet or largescale vortex forming, for example, at viscosity transitions. In both cases, significant effects of the stellar motion were found. Our key results for the early disk evolution are as follows:

1)
When the stellar motion is taken into account, part of theenvelope angular momentum is transferred to the orbital motionof the protostar around the center of mass, reducing the net angularmomentum of the forming disk. This leads to a somewhat smallerdisk mass and radius, and, as a consequence, to a weakergravitational instability and fragmentation as compared to thecase without stellar motion.

2)
The character of mass accretion rate is also moderately affected by the stellar motion, leading to a mild reduction in the accretion variability inherent to gravitationally unstable disks. This is consistent with the overall reduction in the disk mass and size, and in the strength of the gravitational instability of the disk. A Fourier analysis of the spiral structure shows that models with stellar motion are dominated by onearmed spiral modes, while other loworder modes have roughly a similar strength.

3)
Stellar motion has a profound effect on the collapsing envelope, changing its shape from an initially axisymmteric state to a profoundly nonaxisymmetric configuration. This causes a significant excursion of the center of mass of the system, complemented with smallscale wobbling due to the gravitational perturbation of the disk.
Our key results for the protoplanetary disk evolution are as follows:

4)
The gap formation is significantly slower if the stellar motion isneglected: after about 100 orbits of a nonmigrating5 M_{Jup} planet, the gap is only partially depleted (the planetary horseshoe region is filled with gas). On the contrary, the gap is fully depleted if the stellar motion is taken into account. However, the gap is completely devoid of gas after about several hundred orbits independent of whether the stellar motion is on or off. The time required for the horseshoe gas to vanish is inversely proportional to the disk viscosity.

5)
Independent of the applied viscosity, the gap carved by a nonmigrating 5 M_{Jup} planet is significantly narrower, if the stellar motion is taken into account. If the density gradient at the outer gap edge is large enough for a long time (which requires low viscosity, e.g., α ≤ 5 × 10^{4}) a longterm largescale vortex forms. However, we observed the vortex formation only if the stellar motion is taken into account.

6)
The azimuthally averaged disk eccentricity (excited by a nonmigrating 5 M_{Jup} planet) is only slightly modified by the effect of stellar motion. However, the time required for the formation of a quasisteady eccentric disk is shorter in models with stellar motion.

7)
The radial migration of highmass planets (investigated in the range of 1.0−2.5 M_{Jup}) is also affected by the stellar motion: the type II migration rate is significantly smaller if the effect of stellar motion is taken into account. This phenomenon can be explained by the fact that the gap outer edge is depleted in gas if the stellar motion is taken into account, resulting in a smaller magnitude of the outer disk torque.

8)
Regarding the largescale vortex formation at a sharp viscosity transition, we found that the formation of the RWI unstable density bump and the vortex coagulation are stimulated by the stellar motion. If the stellar motion is included, the fully fledged largescale vortex is stronger (the azimuthal asymmetry is larger). On the contrary, the largescale vortex dissolves (only a ringlike perturbation remains) within about 5000 orbits at the distance of the vortex center in models where the effect of the stellar motion is neglected.
Our general conclusion is that the stellar motion has a notable effect on the evolution of gravitationally unstable protostellar disks and a strong effect on the planet or largescale vortexbearing protoplanetary disks. Therefore, inclusion of the indirect potential is recommended in gridbased hydrodynamics simulations of circumstellar disks.
Finally, we note that special care with the computational boundaries needs to be taken when implementing the indirect potential, especially when the mass of the circumstellar disk is commensurable to that of the central star. The acceleration of the noninertial frame of reference may lead to an artificial distortion of the disk or/and the infalling envelope, especially near the outer computational boundary if the numerical resolution is insufficient, as is the case for logarithmically spaced grids. To mitigate this problem, we extended the computational domain and filled it with a low density external environment. However, a better solution needs to be found in the future, probably involving the introduction of moving grids or boundaries.
Acknowledgments
The authors are thankful to the referee for insightful comments. The project was partly supported by the joint OeADOMAA program through project 90öu25 and by the Hungarian OTKA Grant No. 119993. Zs.R. acknowledges support from Momentum grant of the MTA CSFK Lendület Disk Research Group K101393 and NVIDIA Corporation for donating of the Tesla 2075 and K40 GPUs used for this research. E. I. Vorobyov acknowledges support from the Austrian Science Fund (FWF) under research grant I2549N27. The simulations were performed on the Vienna Scientific Cluster (VSC2 and VSC3) and on the Shared Hierarchical Academic Research Computing Network (SHARCNET).
References
 Ataiee, S., Pinilla, P., Zsom, A., et al. 2013, A&A, 553, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: Univ. Arizona Press), 387 [Google Scholar]
 Baraffe, I., Vorobyov, E. I., & Chabrier, G. 2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Basu, S. 1997, ApJ, 485, 240 [CrossRef] [Google Scholar]
 Boley, A. C., & Durisen, R. H. 2008, ApJ, 685, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C., Mej, A. C., Durisen, R. H., et al. 2006, ApJ, 651, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C., Hartquist, T. W., Durisen, R. H., & Michael, S. 2007, ApJ, 656, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P. 2012, MNRAS, 419, 1930 [NASA ADS] [CrossRef] [Google Scholar]
 Cai, K., Durisen, R. H., Boley, A. C., Pickett, M. K., & Mejía, A. C. 2008, ApJ, 673, 1138 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039 [NASA ADS] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 D’Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698 [NASA ADS] [CrossRef] [Google Scholar]
 Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duffell, P. C., Haiman, Z., MacFadyen, A. I., D’Orazio, D. J., & Farris, B. D. 2014, ApJ, 792, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fu, W., Li, H., Lubow, S., & Li, S. 2014, ApJ, 788, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Fung, J, Shi, J., & Ciang, E. 2014, ApJ, 782, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Müller, T. W. A., Kolb, S. M., BenitezLlambay, P., & Masset, F. 2012, A&A, 546, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, H. F. J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. 2014, MNRAS, 437, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M. K. 2015, MNRAS, 448, 3806 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 749 [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. C. B. 2011, MNRAS, 415, 1426 [NASA ADS] [CrossRef] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H. 1991, ApJ, 381, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S H., Seibert, M., & Artymowicz, P. 1999, 526, 1001 [Google Scholar]
 Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009a, A&A, 493, 1125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009b, A&A, 497, 869 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Machida, M. N., Inutsuka, S., & Matsumoto, T. 2011, ApJ, 729, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Mamatsashvili, G. R., & Rice, W. K. M. 2009, MNRAS, 394, 2153 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., & Pudritz, R. E. 2005, ApJ, 618, L137 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Mejía, A. C., Durisen, R. H., Pickett, M. K., & Cai, K. 2005, ApJ, 619, 1098 [NASA ADS] [CrossRef] [Google Scholar]
 Michael, S., & Durisen, R. H. 2010, MNRAS, 406, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Mittal, T., & Chiang, E. 2015, ApJ, 798, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ou, S., Ji, J., Liu, L., & Peng, X. 2007, ApJ, 667, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Pickett, B. K., Mejía, A. C., Durisen, R. H., et al. 2003, ApJ, 590, 1060 [NASA ADS] [CrossRef] [Google Scholar]
 Regály, Zs., Sándor, Zs., Dullemond, C. P., & van Boekel, R. 2010, A&A, 523, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Regály, Zs., Juhász, A., Sándor, Zs., & Dullemond, C. P. 2012, MNRAS, 419, 1701 [NASA ADS] [CrossRef] [Google Scholar]
 Regály, Zs., Sándor, Zs., Csomós, P., & Ataiee, S. 2013, MNRAS, 433, 2626 [NASA ADS] [CrossRef] [Google Scholar]
 Regály, Zs., Király, S., & Kiss, L. 2014, ApJ, 785, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shu, F. S. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. S., Tremaine, S., Adams, F. C., & Ruden, S. P. 1990, ApJ, 358, 495 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010, ARA&A, 48, 391 [NASA ADS] [CrossRef] [Google Scholar]
 SteimanCameron, T. Y., Durisen, R. H., Boley, A. C., Michael, S., & McConnell, C. R. 2013 ApJ, 768, 192S [NASA ADS] [CrossRef] [Google Scholar]
 Szulagyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Takeuchi, T., & Ward, W. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C. E. J. M. L. J. 2008, ApJ, 689, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2015, ARA&A 53, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Tsukamoto, Y., & Machida, M. N. 2011, MNRAS, 416, 591 [NASA ADS] [Google Scholar]
 Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2013, MNRAS, 436, 1667 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Vorobyov, E. I. 2013, A&A, 552, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Varniere, P., & Tagger, M. 2006, A&A, 446, 13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956 [NASA ADS] [CrossRef] [Google Scholar]
 Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
 Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., & Baruteau, C. 2016, MNRAS 458, 3918 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., & Stone, J. M. 2014, ApJ, 795, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X. 2014, ApJ, 785, 122 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Gas surface density distribution in model 1 in the inner 1400 × 1400 au box focusing on the disk. Six evolution times counted from the formation of the star are shown. The scale bar is in g cm^{2} (log units). 

Open with DEXTER  
In the text 
Fig. 2 Similar to Fig. 1, but for model 1+ for which case stellar motion is included. 

Open with DEXTER  
In the text 
Fig. 3 Number of fragments in the disk at a given time elapsed since the formation of the star. The top and bottom panels correspond to models 1 and 1+ (without and with stellar motion, respectively). 

Open with DEXTER  
In the text 
Fig. 4 Top panel: stellar, disk, and envelope masses in model 1 (blackcolored lines) and model 1+ (redcolored lines). The envelope mass is the same in both models. Bottom panel: disk radius in model 1 (blacksolid line) and model 1+ (red solid line). 

Open with DEXTER  
In the text 
Fig. 5 Global Fourier amplitudes of the disk as a function of time in model 1 (top panel) and in model 1+ (bottom panel) calculated for the first four spiral modes. 

Open with DEXTER  
In the text 
Fig. 6 Gas surface density distribution in the inner 16 000 × 16 000 au^{2} box covering both the disk and the infalling envelope in model 1+ at six times after the formation of the central star. The yellow isosurface density contour lines are added to highlight the inner envelope structure. The scale bar is in log g cm^{2}. 

Open with DEXTER  
In the text 
Fig. 7 Similar to Fig. 6, only for model 1 without stellar motion. Only four evolutionary times are shown. 

Open with DEXTER  
In the text 
Fig. 8 Position of the center of mass in model 1+ at the polar coordinate system. The center of mass is initially at (r = 0,φ = 0), but shows significant excursions during the subsequent evolution due to the nonaxisymmetric structure of the infalling envelope. Smallscale wobbling in response to the nonaxisymmetric disk is also evident. 

Open with DEXTER  
In the text 
Fig. 9 Mass accretion rates as a function of time in model 1 (top panel) and model 1+ (bottom panel). 

Open with DEXTER  
In the text 
Fig. 10 Density snapshots in models 2 and 3 for which case a nonmigrating giant planet is embedded in a disk with α_{visc} = 5 × 10^{3} (panel a)) and 5 × 10^{5} (panel b)). In each panel the left and right columns correspond to simulations in which the indirect potential is taken into account (model 2+ and 3+) and neglected (model 2 and 3), respectively. Snapshots are taken after the gap is opened at t = 10^{3} yr corresponding to about 100th orbits (upper row), while at the end of simulation at t = 4 × 10^{4} yr corresponding to about 3600th orbits (lower row). 

Open with DEXTER  
In the text 
Fig. 11 Azimuthally averaged radial density profiles in models 2 and 3 at the end of simulations (corresponding to the lower row of Fig. 10). Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER  
In the text 
Fig. 12 Azimuthally averaged eccentricity profile of the disk at the end of the simulations in models 2 and 3. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER  
In the text 
Fig. 13 Evolution of averaged disk eccentricity in models 2 and 3. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER  
In the text 
Fig. 14 Migration history of highmass planets in models 4 and 5. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER  
In the text 
Fig. 15 Barycenter shift (in au) of the starplanet (x_{bcsp}, x_{bcsp}) and starplantdisk systems (x_{bc}, x_{bc}) as a function of time for models 4+ and 5+. Since the shifts are highly oscillating, the running average over 5 snapshots is also shown. 

Open with DEXTER  
In the text 
Fig. 16 Torques exerted on the planet in models 4/4+ (upper panel) and 5/5+ (lower panel). Significant departures can be observed only for the outer disk torque (Γ_{outer disk}), which eventually explains the decline of the total torque, Γ_{tot}, (see orange solid and dashed curves for models 4/5 and 4+/5+, respectively) in the models where stellar motion is taken into account. 

Open with DEXTER  
In the text 
Fig. 17 Ratio of the surface mass densities in models 4+/4 (left panel), and models 5+/5 (right panel) at the beginning of the planetary migration. The inner and outer gap edges (defined by half maximum density) are found to be nearly identical in models 4/4+ and 5/5+ shown with the dashed circles. 

Open with DEXTER  
In the text 
Fig. 18 Evolution of largescale vortex at viscosity transition in model 6. Left and right columns correspond to simulations in which the indirect potential is taken into account (model 6+) and neglected (model 6), respectively. The density distribution (Σ) is normalized by the initial one (Σ_{0}). 

Open with DEXTER  
In the text 
Fig. 19 Density gradient profile at the end of simulations for models 6+ and 6. 

Open with DEXTER  
In the text 
Fig. 20 Barycenter shift against time in models 6+ and 6. Solid and dashed lines correspond to models in which the indirect potential is taken into account and neglected, respectively. 

Open with DEXTER  
In the text 