Meridional Circulation driven by Planetary Spiral Wakes in Radiative and Magnetized Protoplanetary Discs

In this work, we solve the ideal MHD equations including radiative transfer in a spherical coordinate system ($r$, $\theta$, $\phi$), i.e.

Here $\rho$ is the gas density, $\mathbf{v}$ the gas velocity, and $\mathbf{B}$ the magnetic field. Furthermore, $P_{t}$ is the total pressure, including both thermal and magnetic pressure, giving $P_{t}=P+\frac{B^{2}}{2}$, and $P$ is given by an ideal equation of state

where $u$ is the internal energy, that is given by $u=\frac{k_{B}T}{(\gamma-1)\mu m_{p}}$, with $k_{B}$ being the Boltzmann constant, $\mu$ the mean molecular weight of the gas, $m_{p}$ the proton mass, and $T$ the temperature of the gas. The adiabatic index $\gamma$ is taken to be $1.42$ and the mean molecular weight $\mu$ is taken to be $\mu=2.353$ because the gas is considered a mix of hydrogen and helium in agreement with solar abundance (Flock et al., 2013; Szulágyi et al., 2016; Flock et al., 2017). $\Phi(\mathbf{r})$ is the gravitational potential due to the central star and the orbiting planet, i.e. $\Phi=-\frac{GM_{\star}}{r}-\frac{GM_{p}}{\sqrt{|\mathbf{r}-\mathbf{r_{p}}|^{2}+\epsilon^{2}}}$, where $\mathbf{r_{p}}$ is the position of the planet and $\epsilon$ is the gravitational softening length. $E_{t}=\rho u+\rho\frac{v^{2}}{2}+\frac{B^{2}}{2}$ is the total energy of the gas, including thermal, kinematic, and magnetic energy, but not including radiation energy, which is called $E_{R}$.

The equations for radiative transfer are written following the formalism by Mihalas & Mihalas (1984) and make use of the Planck opacity $\kappa_{P}$ and the Rosseland opacity $\kappa_{R}$ (Kley, 1989; Commerçon et al., 2011). $E_{R}$ is the radiation energy, $a_{R}=\frac{4\sigma_{b}}{c}$ the radiation constant, and $\lambda$ is a flux limiter chosen so that the radiation diffusion approaches the black body radiation in the optically thin scenario, while $\lambda\to 1/3$ in optically thick parts (Levermore & Pomraning, 1981). In particular, in our setups

with $R=\frac{|\nabla E_{R}|}{\rho\kappa_{R}E_{R}}$. This way, $R$ gets very small in the optically thick scenario, giving $\lambda=1/3$, while $R$ gets very large when the opacity is low, meaning that $\lambda\sim\frac{1}{R}\sim\frac{\rho\kappa_{R}E_{R}}{|\nabla E_{R}|}$, and the diffusion term becomes

that approaches the black-body radiation, as in thermal equilibrium $E_{R}=\frac{B(T)}{c}=\frac{4\sigma_{b}T^{4}}{c}$. (Kley, 1989; Commerçon et al., 2011; Flock et al., 2016; Szulágyi et al., 2016).

In our case, non-ideal MHD effects have not been considered as, for example, no resistivity has been included. Nevertheless, these effects may be important when considering low-ionized gas like the one of protoplanetary discs, then we briefly discuss possible implementations and consequences of non-ideal mechanisms in Section 4.

For all runs, the equations were solved via the PLUTO code (Mignone et al., 2007) with the HLLD Riemann solver (Miyoshi & Kusano, 2005), parabolic reconstruction (Mignone, 2014) and second order Runge-Kutta time integrator. In the MHD framework, we used the Constrained Transport module to ensure $\nabla\cdot\mathbf{B}=0$ (Balsara & Spicer, 1999; Mignone et al., 2007). The radiative transfer equations were not solved in the same Godunov formalism, but an implicit solution was implemented in the code in order to be able to use reasonable timestep lengths. In particular, the radiative transfer equations were solved by a module implemented in the PLUTO and described by Flock et al. (2013).

The PLUTO code was modified consistently by adding a module able to solve the implicit radiative transfer equations (Flock et al., 2013, 2016, 2017). A new variable was implemented, i.e. the radiation energy $E_{R}$, that at the initial state follows $E_{R}=\frac{4\sigma_{b}T^{4}}{c}$. For simplicity, the Planck and the Rosseland mean opacities have been considered equal, with their values calculated as a function of pressure and temperature according to the opacity table by Zhu et al. (2009), with a floor value of $10\,cm^{2}g^{-1}$. The radiative transfer solver of PLUTO has been validated by comparing its results to the results of the code JUPITER (Szulágyi et al., 2016), which includes a radiative transfer solver as well. More details and plots are provided in Appendix B. In order to keep the runs stable, a floor density equal to $7\times 10^{-9}M_{\odot}/AU^{3}$ ($10^{-6}$ in code units) was implemented.

The computational grid domain always stretches logarithmically from 2.08 AU to 13 AU in the radial direction, then linearly from $0$ to $2\pi$ in azimuth and linearly from $\sim 0.436\pi$ to $\pi/2$ in co-latitude, so that only the upper half of the disc is simulated and the lower part is just considered to be symmetric and mirrored in the analysis. This gives a grid opening angle of $\sim 0.2$ radians. This range turned out to be large enough in capturing the meridional circulation flow. In fact, in all the runs presented in Section 3, the gas density at upper latitudes is comparable to the floor density, i.e. 3 or 4 orders of magnitude lower than the gas density in the mid-plane.

In the low-resolution runs, we had 210x22x720 ($r$, $\theta$, $\phi$) cells (44 effective cells in the co-latitude direction when mirrored), while high-resolution runs had 420x44x1440 cells. The Courant number has been chosen to be $0.25$. Units have been chosen to have $G=1$, then $M_{\odot}$ has been used for mass, $5.2AU$, i.e. Jupiter’s semi-major axis, for length, and consequently, $\frac{P_{J}}{2\pi}$ for time, where $P_{J}$ is the orbital period of Jupiter.

In the azimuthal direction, the boundary condition is simply periodic. The boundary condition in the co-latitude direction is symmetric at mid-plane, meaning that scalar quantities ($\rho$, $P$, $E$), parallel velocities, and perpendicular magnetic fields simply reflect with the same sign in the ghost cells, while perpendicular velocities and parallel magnetic fields reflect but switching sign. At the inner co-latitude boundary, that in these coordinates corresponds to the higher latitudes, the values in the ghost cells are equal to the value in the first active cell, with the caveat that the velocity in the $\theta$-direction is set to $0$ if positive, i.e. if gas is being injected in the grid domain. The radiation energy is chosen so that the radiation temperature, which is defined as $T_{R}^{4}=\frac{cE_{R}}{4\sigma_{b}}$, is set to a value of $10K$ in the ghost cells to allow cooling of the upper disc layers. The radial and azimuthal components of magnetic fields in ghost cells are set equal to the first active cell, while the $\theta$-component is calculated to keep $\nabla\cdot\mathbf{B}=0$.

Radially, similar prescriptions are applied. All ghost-cell quantities have the same value as the domain cells except the radial velocity, which is set to 0 whenever the gas radial velocity points towards the domain, in order to prevent inflow, and the radial magnetic field, which is automatically calculated to keep it divergence-free. This time though, the radiation temperature is set to $10K$ only if the gas temperature is lower than that value, otherwise, in ghost cells, we set $T_{R}=T$.

First of all, the central star is taken to be a Sun-like star, with $M_{\star}=M_{\odot}$, $R_{\star}=R_{\odot}$, and $T_{\star}=T_{\odot}$ on the surface. The surface density profile of the disc has the form $\Sigma=\Sigma_{0}(R/R_{0})^{-1/2}$, with $\Sigma_{0}$ calculated so that the total mass of the disc at the beginning is $M=0.01M_{\odot}$. The vertical density profile follows

with $H=hR$ and the aspect ratio $h=0.05$ being constant. Here $R=r\sin(\theta)$ and $z=r\cos(\theta)$ are the correspondent cylindrical coordinates, with $R_{0}$ being any reference radius, for example, the semi-major axis of the forming planet.

We initialize our setups with a vertical constant temperature, while the radial profile is calculated to ensure a constant aspect ratio $h$ throughout the disc. This brings to a profile $T\propto R^{-1}$. The initial velocities are set so that the radial and the co-latitude components are set to $0$, while the azimuthal component is set to be equal to the Keplerian velocity $v_{\phi}=\sqrt{\frac{GM_{\odot}}{R}}$.

The initial magnetic field is assigned to be poloidal, i.e. $B_{r}=B_{\phi}=0$, while $B_{\theta}=-B_{0}\left(\frac{R_{0}}{R}\right)$, where $B_{0}$ is calculated so that a given value for the magnetic-$\beta$, i.e. $\beta=\frac{P}{P_{\mathrm{mag}}}=\frac{2P}{B^{2}}$, is obtained at planet location. The components of the magnetic field $\mathbf{B}$ are assigned via a magnetic potential $\mathbf{A}$, defined so that $\mathbf{B}=\nabla\times\mathbf{A}$, that ensures that $\nabla\cdot\mathbf{B}=0$. In particular, in this case, $A_{r}=A_{\theta}=0$, and $A_{\phi}=\frac{B_{0}}{\sin(\theta)}$. This configuration ensures also that the disc does not have magnetic forces at the beginning of the simulation, since $\nabla\times\mathbf{B}=0$.

This setup was run for 50 orbits with no magnetic fields so that a new thermal equilibrium is reached under the heating and cooling effects. Then, every simulation with or without magnetic fields was run starting from this relaxed initial condition, which we identify as time 0, meaning that the first relaxation with no magnetic fields happens between $t=-50$ orbits and $t=0$. In all setups, the simulation ran then for 200 Jupiter-orbits without the planet reaching a steady state with MHD (or HD) + radiative effects. Then, a Jupiter-like planet is gradually grown at $5.2$ AU over 150 orbits, and it was assumed to be in a circular orbit. The planet was a point-mass with the gravitational potential according to Equation 1, with a softening length equal to the cell size at the planet location, meaning $0.04AU$ in low-resolution and $0.02AU$ in high resolution. The mass of the planet is grown via a continuous function as

where the time $t$ is in units of Jupiter-orbits here. Since the planet is implemented as a point mass kept artificially in a circular orbit, it does not migrate nor change its eccentricity or inclination. The simulation was then run for another 250 orbits to let the disc relax and reach a new equilibrium with the planet.

Before analyzing the meridional circulation and gap formation, we first examined the effect of the initial conditions on the final disc structure. In order to do that, we run a set of 4 different low-resolution MHD simulations, each of them starting with a different initial value of the parameter $\beta$, i.e. $\beta=10^{5},10^{6},10^{7},10^{8}$ at planet location. This has been done by choosing the initial value of $B_{0}$, described in Section 2, accordingly. These runs showed that the initial magnetic field magnitude is not relevant in terms of the final outcome of the system. In fact, even though the total magnetic energy of the disc had different initial values in the four cases, it rapidly grew to the same magnetic field strength value over the first 100 orbits, as shown in Figure 1. After that, this value remained roughly constant at a global $\beta\sim 10$ for the entire simulation, even after the planet is injected and fully grown. The only difference we see is the timescale of the magnetic energy growth, which is longer in the cases with a larger initial $\beta$. The fast growth of magnetic fields and magnetic energy is due to the upwind term in the induction equation for $\frac{\partial\mathbf{B}}{\partial t}$ (see Equation 1). In fact, the fast orbital velocities quickly generate a strong toroidal field starting from the vertical-field configuration. In this and the following examples, a global $\beta$ parameter of the disc was defined and calculated as

where the integrals were done over the whole disc volume $V$. This global parameter is the ratio of the total thermal and magnetic energy in the disc, similarly to the local $\beta$, which is the ratio of the local thermal and magnetic pressure. It also happens to be approximately equal to the local $\beta$ at the planet’s location at the beginning of the simulation, while the two differ at the end of the runs when most of the magnetic energy is found at higher latitudes (see plots in Section 3.4).

We found that the longer growing timescale for the magnetic field strength keeps the disc structure stable and prevents large perturbations in the density and temperature profiles, which could affect the results. In fact, the case with the initial $\beta=10^{5}$ proved to be particularly subject to perturbations across the whole disc, that modify the dynamics at the planet location. On the other hand, higher initial $\beta$’s showed to produce much smoother results and to preserve the structure of the disc around the planet much better (Figure 2). The lowest initial $\beta$ that produced unquestionably stable and smooth results is $\beta=10^{7}$ and that is why this value has been chosen as the initial condition for most of the MHD runs. Nevertheless, the case with initial $\beta=10^{5}$ will be mostly analyzed in Appendix A.

Every simulation was run for 200 orbits without the planet, then the planet itself took 150 orbits to grow via a mass-taper function, and another 350 orbits are carried out to let the disc stabilize and reach the steady state. Even though in cases with low viscosity many thousands of orbits are needed in order to reach perfect equilibrium (Hammer et al., 2018), in our case, a few hundred orbits turn out to be enough to reach a sufficiently steady state. In order to check that, we verified that the density and temperature profile did not change significantly over time at the end of the simulation. In particular, Figure 3 shows the average density and temperature profiles of the disc between 560 and 600 orbits and their standard deviation. It is clear that, except for some smaller fluctuations through the disc due to MHD effects, both the density and temperature profiles around the planet are very stable. This suggests that the disc structure after 600 orbits can be taken as a final configuration, and we do not expect it to change significantly even though the same runs were run for longer.

When a Jupiter-like planet forms in a protoplanetary disc, it is expected to open up a gap, while its spiral wakes cause a meridional circulation of material to launch. While the spiral wakes themselves bring material up and down in the entire circumstellar disc, the downward flow is the strongest above the planet causing the vertical accretion stream (Szulágyi et al., 2014; Szulágyi & Garufi, 2021). Figure 4 shows two sections of the protoplanetary disc, one (right side) correspondent to the planet’s location, the other (left side) on the opposite side in the protoplanetary disc after 600 orbits. The color scale corresponds to the vertical velocity component, showing the strong upward and downward flows in the entire circumstellar disc driven by the spiral wakes of the planet, but with the highest magnitude just above the planet. In addition, the gas streamlines are overplotted as well, showing the vertical circular motions, i.e. the meridional circulation.

Figure 5 shows the same as Fig. 4 , but in the MHD case with an initial $\beta=10^{7}$. In this case, too, the meridional circulation pattern and the velocity magnitude around the planet appear to be very similar to the HD case, but in other areas of the circumstellar disc further away from the planet, the vertical velocities and circulation are much stronger than in the HD simulation. In this case, naturally, the entire circumstellar disc is more turbulent.

We followed the accretion flow down to the planet, with the help of 3D velocity streamlines representing the true gas motion (Figure 6). The left panel is the HD simulation, while the right one is the MHD simulation with initial $\beta=10^{7}$. The color map represents the volume density at the mid-plane, where the planet is located on the right of the plot, while the streamline’s colors represent the velocity magnitude in the rotating frame normalized to the Keplerian velocity. The 3D streamlines are above the mid-plane, and show how the flow is going toward the planet from the vertical direction. As the streamlines in the planetary gap encounter the spiral wake of the planet (which is a shock-front), the gas slows down losing angular momentum, Therefore those streamlines spiral down to the circumplanetary disc, some directly hitting the planet’s polar regions and thus being accreted by the planet. Most of the gas entering the Hill-sphere will not be accreted by the planet, the streamlines in the mid-plane region of the CPD spiral outward and bring away gas from the planet region back to the circumstellar disc. Such recycling effect is observed regularly in disc-planet simulations, recently again described by Moldenhauer et al. (2021).

All in all, the accretion flow does not change whether magnetic effects are included or not, highlighting that the flow dynamics are really driven by the planetary spiral arms. Between the two panels in Figure 6, the density values are the main difference and that is the reason why the accretion rate and Hill sphere’s masses are so different, as described in Section 3.6.

We also examined the azimuthal velocity in the mid-plane, normalized by the local Keplerian velocity, to compare the HD and MHD cases, see Figure 7. In the HD case, Figure 7 (left panel) shows clearly the super and sub-Keplerian regions due to the pressure rise and fall at the outer and inner gap edges. These mostly axisymmetric structures are strongly disturbed in the MHD case (see the right panel in Figure 7), which looks much more turbulent and less axisymmetric. Which part of the mid-plane is sub-Keplerian or super-Keplerian is important in the context of dust trapping, because the Keplerian areas are particularly subject to dust accumulation (Kato et al., 2010). From the plots, it is clear that in the absence of magnetic fields, the disc structure is more regular, with rings of sub- and super-Keplerian material easily identified. On the other hand, in the MHD case, the disc structure is again much more turbulent, as expected. Rings are not easy to define since they do not have clear boundaries nor shapes. Furthermore, such a fragmented velocity profile would continuously change with time, while gas is orbiting around the star with different Keplerian velocities. We can infer then that dust accumulation should be much more difficult in the ideal MHD case compared to the HD one and, consequently, formation mechanisms such as the streaming instability (Youdin & Goodman, 2005) are much more difficult to happen.

The mixing due to the meridional circulation and (mainly) the encounter with the spiral shock wake of the planet have an effect also on the viscous behavior of the gas. On one hand, the circular motions generate vortexes on the lower scales, which then cause the growth of turbulence, that can be described by a turbulent, or effective, viscosity in the disc, that can be quantified via an $\alpha$ parameter (Shakura & Sunyaev, 1973). Furthermore, as the gas parcels encounter the spiral shock fronts of the spiral wake of the planet, they lose angular momentum and drive accretion as well. This can be translated into an effective viscosity as well, given the angular momentum loss mechanism.

We estimated the effective viscosity of the circumstellar disc by calculating the Reynolds stress as shown by Flock et al. (2011, 2013). The $r-\phi$ component of the Reynolds tensor is defined as $T_{R}=\rho v^{\prime}_{\phi}v^{\prime}_{r}$, where $v^{\prime}=v-<v>$, and $<>$ is hereafter the average over time. The corresponding $\alpha_{R}$ parameter, i.e. relative to velocity turbulence or to the Reynolds stress, is then calculated as

where the integral is over the whole disc domain. Figure 8 shows the evolution of $\alpha_{R}$ in the HD and MHD cases, with different initial $\beta$’s, for the first 200 orbits, before the planet starts forming, and then once the planet is completely formed from 400 to 600 orbits. The plot shows that the kinematic effective viscosity is lowest for the HD case, i.e. $\alpha_{R}\simeq 10^{-6}$, while the magnetized cases have a higher viscosity, between $10^{-3}$ and $10^{-5}$, as specified in Table 1. This viscosity acts also as a heating mechanism, balancing radiative cooling.

A similar estimate was already done by Stoll et al. (2017), which found that the presence of a planet up to $100M_{\oplus}$ did not affect much the turbulence in the disc in the case of 3D HD inviscid and vertically isothermal simulations. In our case, when a Jupiter-like planet is present in the disc, it clearly dominates turbulent viscosity through the spiral wakes. This causes the effective $\alpha$ to grow up to values between $10^{-3}$ in the HD case and $10^{-2.2}$ in the MHD cases. This means that in the HD case, the spiral wakes contribute to three orders of magnitude higher alpha in the entire disc, meaning that the effects of spirals are the main angular momentum loss mechanism in the disc. In our research field, since the MRI is no longer considered the main angular momentum loss mechanism, there has been an ongoing search for other processes that can remove angular momentum, for example, magnetic winds (Gressel et al., 2015; Béthune et al., 2017) and various instabilities like the vertical shear instability (Stoll & Kley, 2014; Barker & Latter, 2015), the zombie wave instability (Marcus et al., 2015, 2016), Rossby-wave instability (Lovelace et al., 1999) and non-ideal MHD effects (Bai & Stone, 2011, 2013). The spiral shock wakes of planets were not really considered before, or mainly just in the CPD (Szulágyi et al., 2014; Zhu & Baruteau, 2016; Szulágyi & Mordasini, 2017), or without quantification of the effective alpha (Szulágyi & Garufi, 2021). In the MHD simulations (Table 1), the $\alpha_{R}$ values jump to higher final values than in the HD, but in these MHD runs the initial (first 200 orbits without the planet) $\alpha_{R}$ was also much higher than in the HD case (Fig. 9). One can conclude that even in the MHD case also the planet spirals contribute to the effective viscosity.

Velocity fluctuations and material circulation are not the only sources of turbulence in the disc though. In fact, following again Flock et al. (2011, 2013), magnetic fields can give their contribution to the effective viscosity via the Maxwell stress, calculated from the $r-\phi$ component of the Maxwell tensor, i.e. $T_{M}=-\frac{B_{\phi}B_{r}}{4\pi}$, where the full value of the magnetic field is taken into account. In this case, the magnetic equivalent viscosity parameter is calculated as

and the total effective viscosity is estimated as $\alpha=\alpha_{R}+\alpha_{M}$. The evolution of the total viscosity is shown in Figure 9. Now, $\alpha$ is much larger in magnetized discs, reaching values between $10^{-3}$ and $10^{-2}$ before the planet is injected and up to $10^{-1}$ after that, compared to values around $10^{-6}$ and $10^{-3}$ respectively in HD simulations (see Table 2). Once again, the presence of the planet deeply influences the global viscosity of the disc because of the turbulence and magnetic fields developing in the planetary spiral wakes. We should expect a much more viscous behaviour in MHD cases then, that can influence momentum transport and gap opening, even though we do not expect viscous heating to be much larger in the magnetic case in the bulk of the circumstellar disc, except the very inner regions close to the star.

In fact, when the Maxwell stress tensor was defined, the total value of magnetic fields was considered, instead of only its fluctuations, as in the case of the Reynolds stress. As described by Flock et al. (2011, 2013, 2017), the mean magnetic field strength acts as a momentum transport mechanism, equivalently to viscosity. At the same time, fluctuations in the magnetic field structure underline the presence of instabilities, such as the Magneto-Rotational Instability (MRI, Balbus & Hawley 1998), that translates into turbulence and into viscous heating as well. In general, we can then distinguish between a mean Maxwell tensor $T_{M,\mathrm{mean}}=-\frac{<B_{\phi}><B_{r}>}{4\pi}$ and a residual or turbulent tensor $T_{M,\mathrm{res}}=T_{M}-T_{M,\mathrm{mean}}$ where again $T_{M}=-\frac{B_{\phi}B_{r}}{4\pi}$. These tensor components can be used to identify relative viscosity parameters $\alpha_{M,\mathrm{mean}}$ and $\alpha_{M,\mathrm{res}}$, that sum up to retrieve the total equivalent viscosity $\alpha=\alpha_{R}+\alpha_{M,\mathrm{mean}}+\alpha_{M,\mathrm{res}}$. Trivially, $\alpha_{R}$ is the only contribution in HD simulations, while the three components can have different importance in MHD simulations. In particular, in our MHD models, the dominant component is by far the one due to mean magnetic fields. In fact, the Reynolds component is found to be a couple of orders of magnitude lower, and the residual magnetic fields give a negligible equivalent $\alpha$ from the turbulent Maxwell stress tensor. This means that, at least in this resolution, magnetized discs are expected to be much more viscous than HD setups, but, at the same time, viscous heating should not substantially kick in.

After having analyzed the meridional circulation of gas, in this section we address the issue of gap opening, comparing the HD and the MHD cases. First of all, without magnetic fields, a gap is quite clearly opened in the disc. This is observable in Figure 10 where density and temperature profiles of the HD case are shown. In the 2D sections, the gap is visible around the planetary orbit and presents a lower density, other than a lower temperature. Because of the gap formation, gas piles up in the inner and outer disc, increasing the density, especially at the inner boundary. At the planet location, we can also see a concentration of material, in both figures. The resolution is not high enough to determine the structure of such a clump, but that could be in general a trace of the formation of a circumplanetary disc, given the peak in both density and temperature. Since we expect the disc to slowly reach a definitive equilibrium configuration after thousand of orbits (Hammer et al., 2018), it is likely that all the remaining gas in the gap will eventually clear out.

The same section and quantities have been measured in the MHD case and are shown in Figure 11, with a snapshot of the mid-plane and vertical sections taken after 600 orbits, with a fully-formed planet. In this case, the presence of a gap is not clearly spotted. The density and temperature distributions reveal a much more turbulent structure, as already mentioned in Section 3.2, and there is no strong decrease in temperature or density around the planetary orbit, even though circumplanetary material is again visible right around the planet. This again shows that the gap formation and the CPD formation are not necessarily linked, but there can be circumplanetary material (either an envelope or a disc) even when a planetary gap is not present or is yet to form.

More details are shown in Figure 12. These plots compare the mid-plane density (top) and temperature (bottom) profiles at the planet’s location of the HD simulation with the four MHD setups (initial $\beta=10^{5},10^{6},10^{7},10^{8}$). In the HD case, the gap structure is clearly visible in both density and temperature, with a decrease of more than one order of magnitude in the first case. At the same time, the peaks in density and temperature exactly at the planet’s location reveal the presence of circumplanetary material, and the formation of the gap forces most of the gas to pile up in the outer and especially in the inner part of the circumstellar disc, enhancing density close to the inner radial boundary.

On the other hand, the MHD runs show a completely different behavior. First of all, the case with $\beta=10^{5}$ reveals itself again as a more perturbed case, due to the faster growth of magnetic fields and their larger effects on the circumstellar disc structure. Nevertheless, the other runs show, as expected, very similar profiles. The gap is mostly absent, as there is not any clear decrease in the density distribution. The circumplanetary material is again present, with a similar density profile, but higher in magnitude, suggesting more massive CPDs. Even the pile-up of materials in the inner protoplanetary disc is not detected, as a consequence of the lack of a deep gap. At the same time, the temperature profile resembles closely the HD results in the vicinity of the planet, with a small positive offset, while the disc appears to be much colder in the rest of its structure.

The lack of a gap opening in MHD is somehow expected from the results of Section 3.2. In fact, the conditions for a gap to open can be analytically estimated by comparing the timescales of the different mechanisms acting in the gap formation process, i.e. the gravitational torque and the viscous timescales (Lin & Papaloizou, 1980; Papaloizou & Lin, 1984). In this context, Armitage (2010) showed that, following Takeuchi et al. (1996), a gap is expected to open when the following condition is met

where $\alpha$ is the viscosity parameter, $h$ the aspect ratio, and $m=M_{p}/M_{\star}$ the mass ratio between the planet and the star. In our case, with $h=0.05$ and $m=10^{-3}$, the conditions yields to $\alpha<8\times 10^{-3}$. During a simulation $h$ can decrease down to $\sim 0.035$ (see Figure 23), giving the condition $\alpha<1.5\times 10^{-2}$. Since the equivalent viscosity before the planet is injected is about $10^{-6}$ in HD runs, this condition is fully met and a gap is allowed to open. On the other hand, an equivalent $\alpha$ between $10^{-3}$ and $10^{-2}$, such as the one shown by MHD runs, would be exactly on the edge and would not allow the formation of a full gap, and it even grows at the end of the simulations to much higher values. Moreover, the presence of strong magnetic fields in the gap would perturb the dynamics, providing an extra pressure term that would act against the opening of a gap. This effect is not dominant, as the $\beta$ parameter never gets below $10$, as shown in Figure 13, but at the same time not completely negligible.

In fact, an extra pressure term can be quite important in controlling the mechanism of gap formation. The simpler way to study its effect is to enhance the thermal pressure in an HD disc. The faster way to do it is to deactivate radiative transfer and let the disc evolve adiabatically over time. This would prevent the disc to radiate energy away and would produce higher temperatures and, consequently, higher pressures, even though a proper equilibrium may not be reached if the disc keeps getting hotter over time. The result is shown in Figure 14. The adiabatic run is actually much hotter than the others, meaning also that the thermal pressure is much higher in the disc. Consequently, the torque of the planet is not dominant in the gap-opening process anymore. This translates into a lack of gap, with a density profile very similar to the MHD case, even though slightly lower.

To elaborate more on the effective viscosity, even though the gap opening process is prevented by the high effective viscosity, an MHD simulation with an equivalent $\alpha$ is not directly comparable to a viscous HD setup. In fact, in the MHD case, the effective viscosity is dominated by mean magnetic fields, and consequently, we do not expect the temperature to be affected as much as it is in a viscous run. In order to prove that, we ran two HD simulations implementing viscosity, with artificial $\alpha$ parameters set to be $10^{-3}$ and $10^{-2}$, and the results are shown in Figure 15. In the Figure, we observe that the density profile becomes similar to the MHD case (initial $\beta=10^{7}$) while $\alpha$ increases, even though the HD setups show a slightly higher density in the entire circumstellar disc. The thermodynamics is extremely different though. Despite the viscous behavior, the magnetized disc remains much colder than the HD viscous cases, especially far away from the planet. This confirms that the role of MHD can be compared to the effects of viscosity regarding the density distribution, while they are not quite similar with regard to the temperature of the disc.

All the previous results have been obtained in runs with radiative transfer, except the adiabatic comparison shown in Figure 14. We also tested whether similar results would be obtained in case other strategies are implemented together with MHD. For example, using an isothermal equation of state would keep the disc temperature low and possibly damp any effective viscosity, changing in fact the meridional circulation pattern and the gap-opening process.

As shown in Figure 16, in the HD case, an isothermal run produces only smaller differences compared to the radiative tests. In fact, we have a similar gap profile, even though the isothermal case is slightly more massive. On the other hand, the situation changes when comparing magnetized discs. In fact, the isothermal setup appears to have a much clearer gap, almost comparable to the HD cases. This means that the combination of an isothermal treatment with MHD does not produce the same effects of the MHD + radiative transfer implementation. The main difference is to be found in how the equivalent $\alpha$ evolves in the disc. In fact, in this case, the disc is not dominated by the mean Maxwell contribution anymore. This is found to be actually comparable to the contribution of the turbulent velocities (Reynolds stress) for most of the disc evolution. The total viscosity parameter is itself always around $10^{-5}$ before the planet is injected. This means that viscosity is much lower and produces the same behavior as in the HD case, i.e. a gap opening in the disc, a similar density profile, and a meridional circulation pattern.

The different gap-opening mechanism in isothermal and radiative cases comes from the actual magnitude of magnetic fields that develop in the disc. In fact, in the radiative case, the global $\beta$ parameter decreases down to $\sim 10$, while in the isothermal run, $\beta$ never decreases lower than $\sim 10^{3}$, as shown in Figure 17. The turbulent motion, caused by the fluctuations of radiative transport, has proven to be more important in radiative HD simulations, with an equivalent $\alpha_{R}$ that can get almost one order of magnitude above the one shown in the isothermal run, and this causes the magnetic field growth to be much more sustained in the disc. As a consequence, the equivalent $\alpha_{M}$ is also much stronger and prevents a gap from opening.

The different behaviors in the gap-opening process and in the meridional circulation are expected to have consequences on the accretion rates onto the circumplanetary disc and hence onto the planet as well. In order to check that, we followed the amount of mass contained in the Hill sphere around the planet (hereafter defined as $M_{\mathrm{Hill}}$ or Hill mass), with the Hill radius defined as

where $a$ is the semi-major axis of the planet, $M_{p}$ is the mass of the planet, and $M_{\star}$ the mass of the star. We used the Hill-sphere because of the low resolution in the planet’s vicinity, which prevents to calculate the accretion rate directly on the planet’s surface. However, of course, what material enters the Hill-sphere might not end up on the planet, in fact, it is known that there is a significant outflow from the Hill-sphere based on many previous works (Tanigawa et al., 2012; Szulágyi et al., 2014; Ormel et al., 2015a, b; Cimerman et al., 2017). The value of the Hill mass is of course variable in time as the mass of the planet is also changing during a single simulation and the size of the Hill sphere changes accordingly. For the sake of our analysis, we took into account only the final value, after the disc has reached stability with the fully grown planet.

We studied the evolution of $M_{\mathrm{Hill}}$, comparing the final Hill mass in the HD case with the Hill mass in the four MHD cases (initial $\beta=10^{5},10^{6},10^{7},10^{8}$). Firstly, once again, the MHD case with initial $\beta=10^{5}$ proves to be different compared to the other MHD cases. In fact, the perturbations at the planet location do not allow the planet to accrete the same amount of mass and a total $M_{\mathrm{Hill}}\simeq 2\times 10^{-6}M_{\odot}\sim 2\times 10^{-3}M_{J}$ has been measured. Secondly, the final Hill mass, after the disc has relaxed, is quite different in MHD and HD cases. In fact, a magnetized disc (except of course the case with initial $\beta=10^{5}$) produces a Hill mass of about $8\times 10^{-6}M_{\odot}\sim 8\times 10^{-3}M_{J}$, where $M_{J}$ is the Jupiter mass, while the HD case shows a Hill mass that is a factor 4 lower ($2\times 10^{-6}M_{\odot}\sim 2\times 10^{-3}M_{J}$). This is mainly due to the lack of a gap. In fact, the minimum density in the planetary gap is very different in the two cases, as shown in Figure 12, and even though the circumplanetary density has similar profiles in shape, this offset causes the final value of $M_{\mathrm{Hill}}$ to be higher in the MHD case.

The same analysis can be done when an isothermal treatment is implemented. Since the MHD isothermal case has already shown to give different results compared to its radiative counterpart, we expect the Hill mass to decrease. In fact, the two isothermal cases (MHD and HD) turn out to be indeed extremely similar, not only in the gap opening process but also in the accretion onto the planet, reaching both a Hill mass value slightly below $2\times 10^{-6}M_{\odot}$. This means that the lower magnetic energy developed in the disc, as shown by Figure 17, is responsible for the different behaviour of the disc. In general, isothermal runs show a slightly lower Hill mass compared to the radiative HD one, mainly because the disc is slightly warmer around the planet since it cannot radiate energy away (Figure 16), and this causes the gas to settle less towards the mid-plane.

Lastly, we compared the evolution of $M_{\mathrm{Hill}}$ in the MHD case ($\beta=10^{7}$) to a variety of HD runs. We have already seen that the isothermal HD run produces much lower values for the Hill mass. At the same time, also the adiabatic setup, despite being able to reproduce the absence of the gap, gives low and comparable values, again because the high temperatures do not allow gas to settle around the planet. On the other hand, radiative and viscous HD setups are much more interesting. In fact, even though the case with $\alpha=10^{-3}$ shows a low mass, at about $4\times 10^{-6}M_{\odot}$, the $\alpha=10^{-2}$ shows a much higher value, at around $12\times 10^{-6}M_{\odot}$. This again confirms that effective viscosity is the key factor in the gap-opening process and in planetary accretion as well. Having higher $\alpha$’s not only prevents a gap from fully opening but also fosters higher accretion rates onto the planet. In particular, we can infer that a value for $\alpha$ between $10^{-3}$ and $10^{-2}$ would produce a Hill mass similar to the MHD case.

Other than the mass in the Hill sphere, one could measure the mass accretion rate, i.e. the time derivative of the Hill mass itself $\frac{dM_{\mathrm{Hill}}}{dt}$. In general, these values will not be stable in time, nor easily measurable. In fact, shocks and heating/cooling processes around the planet cause the temperature in the Hill sphere to oscillate, and, consequently, also the Hill mass will not have a completely smooth evolution, but will rather show oscillations around a mean value. In fact, when the Hill sphere is colder, more mass can be compressed in it, while when it is warmer, the gas expands and the outer layer is pushed away by thermal pressure (Szulágyi et al., 2016), decreasing the overall density of the Hill sphere. This is also why radiative simulations are so important in determining realistic accretion rates.

In particular, in both HD and MHD radiative cases, the accretion rate values initially span from $-10^{-6}\frac{M_{\odot}}{\mathrm{orbit}}$ to $+10^{-6}\frac{M_{\odot}}{\mathrm{orbit}}$ with an average value around $+10^{-8}\frac{M_{\odot}}{\mathrm{orbit}}$ (Figure 18). For comparison, this would translate into an accretion rate of $10^{-6}\frac{M_{J}}{\mathrm{year}}$, in agreement with the assumed formation timescale of Jupiter. Towards the end of the simulation, after 600 orbits, the oscillation amplitude of the HD cases decreases by about one order of magnitude, reaching the same amplitude values we observed in isothermal cases, both HD and MHD, as shown again in Figure 18. This proves that the combined effect of magnetic fields and radiative transfer makes the meridional circulation, hence accretion, more turbulent and less smooth than all the other combinations.

In order to verify the effect of higher resolutions, we also run some selected configurations in a 420x44x1440 setup (8 times the number of cells as in the low-resolution case). In particular, three different radiative cases were tested. One HD case and two MHD cases with $\beta=10^{5},10^{7}$. Again, an HD case was run for about 50 orbits, and the three setups were initiated from there, obviously adding magnetic fields when necessary. The first, and not surprising result, is that the two HD simulation results (low- and high-resolution) are almost identical. In practice, increasing the resolution does not produce relevant effects on the hydrodynamics of the disc. Even the effective viscosity in the disc turns out to be very similar, as shown in Figure 19, with values oscillating around $10^{-3}$ in both resolutions after the planet is formed.

Comparing various resolutions is especially important in MHD since magnetic effects are known to be sensitive to resolution (Fromang & Nelson, 2006). Starting from the equivalent viscosity, we found that in the case with initial $\beta=10^{5}$, the value for $\alpha$ was roughly the same in the two resolutions (Figure 19) with some more oscillations in the high-resolution run. The latter happens because the residual magnetic fields are the main source of viscosity, in this case. In fact, even though residuals were found not to be significant in the low-resolution runs, in high-resolution we start to resolve the magnetic instabilities, such as the MRI. Consequently, that causes the residual $\alpha$ to be much more important compared to the total viscosity value. This means, in particular, that one should start to see the heating effect of viscosity, which was not important in lower resolution. In the case with initial $\beta=10^{7}$ though, this is even more relevant. In fact, in this case, the total viscosity was found to be initially lower in the high-resolution run and to reach roughly the same value at the end of the simulation. Given that the total magnetic $\beta$ has found to be even lower than $\sim 10$ at the end of the simulation despite the higher temperature (and thermal pressure), we already have a hint that MHD-driven instabilities could play a major role in this case.

In fact, as shown in Figure 20, the density and temperature profiles in the $\beta=10^{5}$ cases are similar enough, with some differences due to a slight overheating in the high-resolution run. On the other hand, heating is much more prominent in the case with initial $\beta=10^{7}$ because of the residual viscosity, causing the disc to be much hotter, thus the density profile also to be different. This highlights that going to high resolutions is a powerful tool to resolve and detect magnetic instabilities that can modify the outcome of the model in certain circumstances, at least in an ideal-MHD scenario. It is known, that more accurate ionization models predict a low coupling between the gas and the magnetic fields, that damps any MHD-driven instability in the so-called dead zone, where planets are forming (Gammie, 1996). In this more realistic case, increasing the resolution would not resolve any instabilities better and it is expected to give roughly the same results.