would be the case if the disc were truncated by a co-rotating stellar magnetosphere (shearless boundary layer). At the outer boundary we assume the torque gradient $\partial z / \partial r \,(r_{\rm out}) = 0$. Neither boundary condition is critical as we are interested in the flow near the planet, away from either boundary.

The surface density of the disc is initialized with the similarity solution for an isolated viscous accretion disc with $\nu \propto r^1$ (@lp_1974; @hcg_1998): $$\begin{aligned}
 \label{eqn:sim}
    \Sigma(r,0) = \frac{M_{\rm disc}}{2 \pi r_1^2} \frac{r_1}{r}  e^{-r/r_1} \,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\, 
    ({\rm simulation})\end{aligned}$$ where $M_{\rm disc}  = 15.5 \, M_{\rm J} \simeq 0.015  \, M_\odot$ is the initial mass of the disc and $r_1 = 30$ au is a characteristic disc radius (where the diffusion time is $r_1^2/\nu \simeq 5$ Myr). We consider two fixed planet masses, $M_{\rm p} = 0.3 \,M_{\rm J} < M_{\rm thermal}$ and $M_{\rm p} = 10 \,M_{\rm J} > M_{\rm thermal}$. Planet masses that freely grow are modeled in section \[sec:mass\].

At every timestep, we first advance $\Sigma(r,t)\rightarrow \Sigma(r,t+\Delta t)$ for all $r$ according to (\[eq:kep\_simp\]) using the implicit solver, and then we advance $\Sigma(r_{\rm p},t)\rightarrow \Sigma(r_{\rm p},t+\Delta t)$ using (\[eqn:euler\]) and (\[eqn:subgrid\]). 