lower limits.

So far we have described how we compute the mass sink term, which includes the sub-grid Lindblad torque, in equation (\[eq:dsig\_dt\]). The remaining diffusive term is solved in a standard way. We first change variables to $z \equiv r^{1/2} \nu \Sigma$ and $y \equiv 2 r^{1/2}$ so that the diffusive portion of equation reads $$\begin{aligned}
 \label{eq:kep_simp}
\pd{z}t &= \frac{12 \nu}{y^2} \pd{^2 z}{y^2}\end{aligned}$$ with non-constant diffusion coefficient $12 \nu/y^2$. We solve equation as an initial value problem using an implicit scheme (e.g., @ptv_2007). Our computation grid extends from an inner boundary of $r_\mathrm{in} = 0.01 \, \mathrm{au}$ to an outer boundary of $r_\mathrm{out} = 500 \, \mathrm{au}$, and is divided into 300 cells that are uniform in $\Delta y$. We fix the timestep $\Delta t = 10^{-4} t_{\nu,{\rm p}}$, where $t_{\nu,{\rm p}} \equiv r_{\rm p}^2 / \nu (r_{\rm p}) \simeq 1.7$ Myr is the viscous diffusion timescale at the planet’s orbital radius of $r_{\rm p} = 10$ au (where $h \simeq 0.054$). Recognizing that our transformed variable $z$ is proportional to the viscous torque $2 \pi \nu \Sigma r^3  d \Omega / dr \propto r^{1/2} \nu \Sigma$, we use a torque-free inner boundary condition, $z (r_\mathrm{in}) = 0$, as 