\documentclass[11pt]{article} \usepackage{graphicx,natbib,amsmath} \setlength{\textheight}{24.5cm} \setlength{\textwidth}{16.5cm} \setlength{\topmargin}{-20mm} \setlength{\evensidemargin}{-2mm} \setlength{\oddsidemargin}{-2mm} \setlength{\headsep}{4mm} \def\nH{n_{\langle\rm H\rangle}} \def\rhod{\rho_{\rm d}} \def\rhom{\rho_{\rm m}} \def\pabl#1#2{\frac{\partial #1}{\partial #2}} \def\vth{v_{\rm th}} \def\cs{c_{\rm s}} \def\vdr{{v_{\rm dr}}} \def\vdreq{{v^{\hspace{-0.8ex}^{\circ}}_{\rm dr}}} \def\Dd{D_{\rm\hspace{0.05ex}d}} \def\Dgas{D_{\rm gas}} \def\vz{\langle v_z\rangle} \def\tstop{\tau_{\rm stop}} \def\tedd{\tau_{\rm eddy}} \def\St{{S\hspace*{0.1ex}t}} \def\apj{ApJ} \def\pasp{PASP} \def\pasj{PASJ} \def\araa{ARA\&A} \def\aap{A\&A} \def\aaps{A\&AS} \def\apjl{ApJL} \def\apjs{ApJS} \def\mnras{MNRAS} \def\aj{AJ} \def\nat{Nature} \def\icarus{Icarus} \def\pasa{PASA} \DeclareRobustCommand{\Eins}{\text{\usefont{U}{bbold}{m}{n}1}} \def\la{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil $\displaystyle##$\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr <\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr <\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr <\cr\sim\cr}}}}} \def\ga{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil $\displaystyle##$\hfil\cr>\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr >\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr >\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr >\cr\sim\cr}}}}} \begin{document} \centerline{\huge\bf Mixing and Settling in Disks} \vspace*{5mm} \centerline{(c) Peter Woitke, November 2019} \section{Vertical equilibrium between mixing and settling} \citet{Mulders2012} write down the continuity equation for the size-dependent dust mass density $\rhod(a)$ as \begin{equation} \pabl{\rhod}{t} = \pabl{}{z}\bigg(\underbrace{\rho\,\Dd\, \pabl{}{z}\Big(\frac{\rhod}{\rho}\Big)}_{j_{\rm diff}}\bigg) - \pabl{}{z}\Big(\underbrace{\rhod\vdreq}_{j_{\rm set}}\Big) \end{equation} where $\rho$ is the gas mass density, $t$ the time and $z$ the vertical coordinate. It is assumed that \begin{itemize} \item there is no bulk motion of the gas $\vec{v}=0$, \item grains are transported downwards by gravitational settling with their size-dependent equilibrium settling velocity $\vdreq(a)$, \item grains are transported upwards by eddy diffusion according to the size-dependent dust diffusion coefficient $\Dd(a)$, driven by dust particle concentration gradients $\pabl{}{z}\big(\frac{\rhod(a)}{\rho}\big)$, \item grains are passive tracers in this model, their size and mass is assumed not to change, e.g.\ due to coagulation, shattering or ice formation. \end{itemize} The latter assumption is not quite so obvious. Grains could, for example, be shattered in the midplane, then mixed up by turbulence where they would re-grow and settle back to the midplane. Assuming that the grains have relaxed towards a steady state, we have for each size $\pabl{\rhod}{t}\to 0$ and zero net flux $j_{\rm diff}-j_{\rm set}\to 0$ \begin{equation} \rho\,\Dd\,\pabl{}{z}\Big(\frac{\rhod}{\rho}\Big) ~=~ \rhod\vdreq \ , \end{equation} from which we find \begin{equation} \boxed{ \pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ \frac{\vdreq}{\Dd}} \label{basic1} \ . \end{equation} Equation (\ref{basic1}) addresses the speed at which the \ {\it dust densities falls off quicker than the gas} \ with increasing height. Integration from 0 to $z$ results in \begin{equation} \ln\frac{\rhod(z)/\rhod(0)}{\rho(z)/\rho(0)} = \int_0^z \frac{\vdreq(z')}{\Dd(z')}\;dz' \end{equation} As long as the integral on the right side is small ($\ll 1$), the dust densities follow the gas densities. Otherwise (keeping in mind that $\vdreq$ is downward, hence negative), the dust densities will subside sooner than the gas. Larger particles have larger $\vdreq$, so they are more affected than small grains. \section{Basic processes} \subsection{Gravitational settling} For large Knudsen numbers and subsonic particle velocities (Epstein regime), the equilibrium settling velocity, also called terminal fall speed, is given by \citet{Schaaf1963} \begin{equation} \vdreq = -\frac{g_z\,\rhom\,a}{\rho\,\vth} = -g_z \tstop \label{vdr} \end{equation} where $g_z$ is the vertical gravitational acceleration, $a$ is the particle radius, $\rhom$ the dust material density, and where the thermal velocity is given by \begin{equation} \vth = \sqrt{\frac{8kT}{\pi\mu}} = \sqrt{\frac{8}{\pi}}\,\cs \ . \label{vth} \end{equation} $k$ the Boltzman constant, $T$ the local gas temperature, $\mu$ the mean molecular mass, $\cs$ the local isothermal speed of sound, and $p=\cs^2\rho$ the local gas pressure. Equation~(\ref{vdr}) seems widely incorrect in the literature, where often $\cs$ is sloppily used instead of the correct mean of the magnitude of the thermal velocities of gas particles $\vth$, as required in Schaaf's formula, see \citet{Woitke2003}. The second part of Eq.\,(\ref{vdr}) follows from a general relaxation ansatz \begin{equation} \tstop = m\,\left(\pabl{F_{\rm fric}}{\vdr}\bigg|_{\,\vdreq}\right)^{-1} =~ \frac{\rhom\,a}{\rho\,\vth} \ , \end{equation} where $\tstop$ is the stopping or frictional coupling timescale, $m$ the mass of a dust particle, $F_{\rm fric}(\rho,a,\vdr)$ the frictional force for large Knudsen numbers and subsonic particle drift $\vdr\ll\cs$, see Eq.\,(21) in \citet{Woitke2003}. \subsection{Gas diffusion} Eddy diffusion (``turbulent mixing'') is assumed to be dominant over gas-kinetic diffusion. The eddy diffusion coefficient in general is given by a characteristic velocity times a characteristic length \begin{equation} \Dgas = \vz\,L \end{equation} where $L$ is the size of the largest eddies (``mixing length'') and $\vz$ the root-mean square average of the fluctuating part of the vertical gas velocities. \subsection{Dust diffusion} In a Kolmogorov type of power spectrum $P(k)\!\propto\!k^{-5/3}$, there are different turbulent modes associated with different wave-numbers $k$ or different spatial eddy sizes $l$. Any given particle of size $a$ tends to co-move with all sufficiently large and slow turbulent eddies whereas its inertia prevents following the short-term, small turbulence modes. In order to arrive at an effective particle diffusion coefficient, the advective effect of all individual turbulent eddies has to be averaged, and thereby transformed into a collective particle diffusion coefficient. This procedure is carried out with different methods and approximations, see \citet{Dubrulle1995}, \citet{2004ApJ...614..960S} and \citet{Youdin2007}. The result of \citet[][see their Eq.\,27]{2004ApJ...614..960S}, reads \begin{equation} \Dd = \frac{\Dgas}{1+\St} \end{equation} where $\St$ is the Stokes number given by \begin{equation} \St = \frac{\tstop}{\tedd} \ . \end{equation} $\tedd$ is the eddy turnover or turbulence correlation timescale in consideration of the largest eddy, which is assumed to equal the mixing length $L$. \begin{equation} \tedd = \frac{L}{\vz} \ . \end{equation} Our basic equation (\ref{basic1}) for the balance of upward mixing against gravitational settling is hence \begin{equation} \boxed{\pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ -\frac{g_z\,\tstop}{\Dgas}\,(1+\St)} \ . \label{basic2} \end{equation} \section{Disk physics} So far, our approach is quite general, describing an equilibrium between settling and mixing in a relaxed 1D structure, as long as our dust particles are in the Epstein regime and chemically/physically passive. In the following, we introduce a bit of disk jargon which will enable us to make the connections to the disk literature. On Keplerian orbits, at pressure-supported height $z$ above the midplane, the zentrifugal force caused by the rotation around the symmetry axis is balanced by the radial component of gravity \begin{equation} \frac{v_K^2}{r} ~=~ \Omega^2 r ~=~ g_r \ , \end{equation} where $v_K$ is the Keplerian velocity and the gravity has components \begin{equation} g = \frac{GM_\star}{r^2+z^2} \quad\mbox{with components}\quad \frac{g_r}{g} = \frac{r}{\sqrt{r^2+z^2}} \quad\mbox{and}\quad \frac{g_z}{g} = \frac{z}{\sqrt{r^2+z^2}} \end{equation} $r$ is the distance from the symmetry axis and $z$ the height above the midplane. We find \begin{eqnarray} \Omega^2 &=& \frac{GM_\star}{(r^2+z^2)^{3/2}} ~\approx~ \frac{GM_\star}{r^3}\\ v_K^2 &=& \frac{GM_\star\,r^2}{(r^2+z^2)^{3/2}} ~\approx~ \frac{GM_\star}{r}\\ g_z &=& \frac{GM_\star\,z}{(r^2+z^2)^{3/2}} ~=~ z\,\Omega^2 \end{eqnarray} The approximations on the right side are called the ``thin disk'' approximation. $\Omega$ is the Keplerian circular frequency. The condition of vertical hydrostatic equilibrium can be expressed by \begin{equation} \pabl{p}{z} = -\rho\,g_z = -\rho\,z\,\Omega^2 \ . \label{hydrostat} \end{equation} For thin disks with vertically constant temperature $T$ and vertically constant mean molecular weight $\mu$, we can write $\pabl{p}{z}=\cs^2\pabl{\rho}{z}$ which allows us to integrate Eq.\,(\ref{hydrostat}) as \begin{eqnarray} \pabl{}{z} \ln\rho = -z\frac{\Omega^2}{\cs^2} \quad\Rightarrow\quad \rho(z) &\!\!\!=\!\!\!& \rho_0\exp\Big(-\frac{z^2}{2\,H_p^2}\Big)\\ \mbox{with pressure scale height}\quad H_p &\!\!\!=\!\!\!& \cs/\Omega \ , \end{eqnarray} and our basic equation (\ref{basic2}) reads \begin{equation} \boxed{\pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ -\frac{z\,\Omega^2\,\tstop}{\Dgas}\,(1+\St)} \ . \label{basic3} \end{equation} According to an $\alpha$ disk model \citep{Shakura1973}, the typical vertical length scale is $L=H_p$ and the turbulent velocity is $\vz=\alpha\,\cs$, so the eddy diffusion coefficient is given by \begin{eqnarray} \Dgas &=& \alpha\,\cs\,H_p \\ \tedd &=& \frac{H_p}{\alpha\,\cs} ~=~ \frac{1}{\alpha\,\Omega} \\ \tstop &=& \St\;\tedd ~=~ \frac{\St}{\alpha\,\Omega} \ . \label{tstop} \end{eqnarray} Something strange happens here, as \citet{Dubrulle1995}, \citet{2004ApJ...614..960S} and \citet{Riols2018} all have $\tedd = H_p/\cs = 1/\Omega$ which I don't understand. The eddy turnover timescale must be much larger than this, because typical turbulent velocities are $\alpha\cs$, i.e. much smaller than $\cs$, because $\alpha\ll 1$. According to Eq.\,(\ref{tstop}) we then have \begin{equation} \pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ -\frac{z\,\Omega\,\St}{\alpha\,\Dgas}\,(1+\St) \ , \label{basic4} \end{equation} which is not quite as general as Eq.\,(\ref{basic3}), because $\cs=H_p\,\Omega$ has been used, which requires $T\!=\,$const and $\mu\!=\,$const. All quantities on the right hand side of Eq.\,(\ref{basic4}), $\Omega$, $\St$, $\alpha$ and $\Dgas$, as well as $T$, $\cs$, $\mu$ and $\rhom$, will in general depend on height $z$. In addition, $\St$ and $\tstop$ depend on grain size $a$. \section{Application in disks} \subsection{Dubrulle et al.\,(1995)} \citet{Dubrulle1995} ignore all $z$-dependences on the r.h.s.\ of Eq.\,(\ref{basic3}) and replace all those variables by their mid-plane values. They apparently drop the $(1+\St)$ term and assume $\St=\tstop(0)\,\Omega$, i.e.\ without the $\alpha$. In that case, we carry out the $z$-integration over \begin{equation} \pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ -\frac{z\,\Omega^2\,\tstop(0)}{\alpha\,\cs(0)\,H_p} \end{equation} which results in \begin{equation} \ln\frac{\rhod/\rhod(0)}{\rho/\rho(0)} ~=~ -\frac{\Omega^2\,\tstop(0)}{\alpha\,\cs(0)\,H_p}\,\int_0^z z'\,dz' ~=~ -\frac{z^2}{2\,H_p^2}\,\frac{\Omega^2\,\tstop(0)\,H_p}{\alpha\,\cs(0)} ~=~ -\frac{z^2}{2\,H_p^2}\,\frac{\Omega\,\tstop(0)}{\alpha} \ . \end{equation} We introduce the size-dependent dust scale height $H_a$ via $\rhod(a,z)=\rhod(a,0)\exp(-\frac{z^2}{2\,H_a^2})$ to find \begin{eqnarray} \Rightarrow\quad -\frac{z^2}{2\,H_a^2} + \frac{z^2}{2\,H_p^2} &=& -\frac{z^2}{2\,H_p^2}\,\frac{\Omega\,\tstop(0)}{\alpha} \end{eqnarray} \begin{equation} \Rightarrow\quad \boxed{\frac{H_p^2}{H_a^2} ~=~ 1 + \frac{\Omega\,\tstop(0)}{\alpha}} \label{Dubrulle} \end{equation} \paragraph{Notes about the implementation in ProDiMo:}\ \\[2mm] In ProDiMo there are variables called $$ \St' = \frac{\Omega\rhom\,a}{\rho(0)\,\cs(0)} \quad\quad\mbox{and}\quad\quad {\rm Hh2} = \sqrt{3}\,\frac{\St'}{\alpha} \quad\quad\mbox{and}\quad\quad {\rm reduce} = \frac{H_a^2}{H_p^2} = \frac{1}{1+{\rm Hh2}} $$ so the Dubrulle midplane-settling is currently implemented as \begin{equation} \frac{H_p^2}{H_a^2} ~=~ 1 + \frac{\sqrt{3}\,\Omega\,\rhom\,a}{\alpha\,\rho(0)\,\cs(0)} ~=~ 1 + \sqrt{3}\,\frac{\vth}{\cs(0)}\frac{\Omega\,\tstop(0)}{\alpha} \end{equation} Comparison to Eq.\,(\ref{Dubrulle}) reveals a mismatch of a factor $\sqrt{\frac{3\times 8}{\pi}}\approx 2.76$ for Hh2, which is substantial as it enters in the exponential, equivalent to the viscosity parameter $\alpha$ being off by that factor. For hydrostatic models, an additional complication is that one should interpret $H_p$ as an actual pressure scaleheight rather than a density scaleheight, so considering (partial) pressures rather than densities, we write \begin{equation} \frac{f(a,z)\,\cs^2(z)}{f(a,0)\,\cs^2(0)} ~=~ \exp\Big(\frac{\ln(p(z)/p(0))}{\rm reduce}\Big) ~=~ \exp\Big(\frac{\ln\,p(z)}{\rm reduce}-\frac{\ln\,p(0)}{\rm reduce}\Big) \end{equation} This is done to ensure a monotonic decrease of $f(a,z)$ with height in the hydrostatic mode, where $p(z)$ decreases monotonically, but $\rho(z)$ may not, because of sudden changes in $T(z)$ or $\mu(z)$. Since $C=f(a,0)\,\cs^2(0)$ is a constant normalisation factor, and ${\rm const}=\ln\,p(0)/{\rm reduce}$~ is constant, too, we calculate the dust size distribution function $f(a)$ from \begin{equation} f(a,z) ~\propto~ \frac{1}{\cs^2(z)} \exp\Big(\frac{\ln\,p(z)}{\rm reduce}-{\rm const}\Big) \ . \label{Pro1} \end{equation} In the mostly used MCFOST-like mode, where $\rho(r,z)$ is set parametrically, $H_p$ is constant and hence well-defined, $\cs(z)=\cs(0)$ is assumed, and $\cs(0)=H_p \Omega$ is used to determine $\cs(0)$, the ``pressure'' is computed as $p(z)=\cs^2(0)\rho(z)$, and we simply have \begin{equation} \frac{f(a,z)}{f(a,0)} = \exp\Big(-\frac{z^2}{2\,H_p^2\,{\rm reduce}}\Big) \label{Pro2} \end{equation} but we use Eq.\,(\ref{Pro1}) instead, also in the MCFOST-like mode, where it is equally valid. Equation\,(\ref{Pro2}) explains the meaning of ``reduce'', namely the factor by which the gas scale height$^2$ is reduced for size $a$. Eventually, $f(a,0)$ is fixed by the nomalisation condition to have the same dust column density before/after settling in each size bin. \subsection{Riols \& Lesur\,(2018)} \citet{Riols2018} also drop the $(1+\St)$ term and start with Eq.\,(\ref{basic4}), taking into account that $\St$ is, and $\Omega$ and $\Dgas$ might be, $z$-dependent. They furthermore explicitly assume a Gaussian for the vertical gas density distribution, $\rho = \rho(0) \exp\big(-\frac{z^2}{2\,H_p^2}\big)$, and since $\St\,\propto 1/\rho$ they write \begin{eqnarray} \St &=& \St(0)\,\exp\Big(\frac{z^2}{2\,H_p^2}\Big) ~=~ \alpha\,\Omega\,\tstop(0) \,\exp\Big(\frac{z^2}{2\,H_p^2}\Big) \ , \end{eqnarray} Equation\,(\ref{basic4}) hence becomes \begin{equation} \pabl{}{z}\Big(\ln\frac{\rhod}{\rho}\Big) ~=~ -\frac{z\,\Omega\,\St}{\alpha\,\Dgas} ~=~ -\frac{z\,\Omega\,\St(0)}{\alpha\,\Dgas} \,\exp\Big(\frac{z^2}{2\,H_p^2}\Big) \end{equation} and integration yields \begin{equation} \ln\frac{\rhod/\rhod(0)}{\rho/\rho(0)} ~=~ -\int_0^z \frac{z'\,\Omega\,\St(0)}{\alpha\,\Dgas} \,\exp\Big(\frac{z'^2}{2\,H_p^2}\Big)\,dz' \ , \end{equation} which reproduces Eq.\,(33) in \citet{Riols2018} if one keeps in mind that their definition of $\St$ misses the $\alpha$. \begin{equation} \frac{\rhod}{\rhod(0)} ~=~ \exp\Big(-\frac{z^2}{2\,H_p^2}\Big) \exp\left(\int_0^z \frac{z'\,\Omega\,\St(0)}{\alpha\,\Dgas} \,\exp\Big(\frac{z'^2}{2\,H_p^2}\Big)\,dz'\right) \ . \end{equation} With the additional assumption that $\Omega$, $\Dgas$ and $\alpha$ are $z$-independent we find \begin{eqnarray} \ln\frac{\rhod}{\rhod(0)} &=& -\frac{z^2}{2\,H_p^2} - \frac{\Omega\,\St(0)}{\alpha\,\Dgas}\int_0^z z'\,\exp\Big(\frac{z'^2}{2\,H_p^2}\Big)\,dz' \\ -\frac{z^2}{2\,H_a^2} &=& -\frac{z^2}{2\,H_p^2} - \frac{\Omega\,\St(0)}{\alpha\,\Dgas} H_p^2\bigg(\exp\Big(\frac{z^2}{2\,H_p^2}\Big)-1\bigg) \\ \frac{H_p^2}{H_a^2} &=& 1 + \frac{\Omega\,\St(0)\,H_p^2}{\alpha\,\Dgas}\; \Big(\frac{2\,H_p^2}{z^2}\Big) \bigg(\exp\Big(\frac{z^2}{2\,H_p^2}\Big)-1\bigg) \end{eqnarray} \begin{equation} \Rightarrow\quad \boxed{\frac{H_p^2}{H_a^2} ~=~ 1 + \frac{\tstop(0)\Omega}{\alpha}\; \Big(\frac{2\,H_p^2}{z^2}\Big) \bigg(\exp\Big(\frac{z^2}{2\,H_p^2}\Big)-1\bigg) } \end{equation} which clearly shows the improvement with respect to the \citet{Dubrulle1995} forumula (Eq.\,\ref{Dubrulle}). The decrease of the dust scaleheight $H_a$ with respect to Eq.\,(\ref{Dubrulle}) is remarkable, the function $(2/x^2)\big(\exp(x^2/2)-1\big)$ is 3.19 at $x=2$, 20 at $x=3$, 370 at $x=4$ and 21000 at $x=5$. Infrared emission lines from disks, for example, are quite typically emitted from layers located at about $3-4$ scaleheights, see e.g. \citet{Woitke2018b}, their Fig.\,8, where lines from OH, H$_2$O, CO$_2$ and HCN are emitted from $z/r\approx 0.15-0.2$ and the scale height at 1\,au is $H_p=\rm 10\,au\,\big(\frac{1\,au}{100\,au}\big)^{1.15}= 0.05\,au$. \paragraph{Notes about the implementation in ProDiMo:}\ \\[2mm] In ProDiMo the following variables are used $$ \St'(0) = \frac{\Omega\,\rhom\,a}{\rho(0)\,\cs(0)} \quad\quad\mbox{and}\quad\quad {\rm Hh2} = \sqrt{3}\,\frac{\Omega\,\St'(0)\,H_p^2}{\Dgas} \quad\quad\mbox{and}\quad\quad \Dgas = \alpha\,\cs\,H_p $$ $$ {\rm fak} = \Big(\frac{2\,H_p^2}{z^2}\Big) \bigg(\exp\Big(\frac{z^2}{2\,H_p^2}\Big)-1\bigg) \quad\quad\mbox{and}\quad\quad {\rm reduce} = \frac{H_a^2}{H_p^2} = \frac{1}{1+{\rm Hh2\;fak}} $$ The same inconsistency of a factor $\sqrt{\frac{3\times 8}{\pi}}\approx 2.76$ occurs. It remains unclear to me why both \citet{Dubrulle1995} and \citet{Riols2018} have not taken into account the $(\St+1)$ factor in the dust diffusion coefficient. \subsection*{Riols \& Lesur in hydrostatic models} In the hydrostatic case, where $T_{\rm gas}(z)$, $\mu(z)$ and $p(z)$ are consistently calculated from the chemistry and gas energy balance, and where consequently the scale height $H_p$ is not constant and ill-defined, we follow the same concept \begin{equation} \pabl{}{z}\Big(\ln\frac{f(a,z)\cs^2(z)}{p(z)}\Big) ~=~ -\frac{z\,\Omega^2\tstop(a,z)}{D_{\rm gas}(z)}\left(1+\St(a,z)\right) \ . \end{equation} With the same simplifications used above, i.e.\ drop $(1+\St)$, use $\tau_{\rm eddy}=1/\Omega$, use $\tstop=\sqrt{3}\,\rho_{\rm m} a /(\rho \cs)$, use $D_{\rm gas}=\alpha \cs H_p=\alpha \cs^2/\Omega$, the result is \begin{equation} \ln\frac{f(a,z)\cs^2(z)}{f(a,0)\cs^2(0)} ~=~ \ln\frac{p(z)}{p(0)} - \underbrace{\frac{\sqrt{3}\,a}{\alpha} \int_0^z z'\,\frac{\Omega^3\,\rho_{\rm m}(z')}{\cs^3(z')\,\rho(z')}\,dz'}_{F(z)} \label{Pro3} \end{equation} \vspace*{-2mm} \begin{equation} f(a,z) ~\propto~ \frac{1}{\cs^2(z)} \exp\!\Big(\ln\frac{p(z)}{p(0)}-{F(z)}\Big) \label{Pro4} \end{equation} Since $\rho(z)$ becomes small as $z$ increases, $F(z)$ will reach the order of unity at some point, and then soon becomes $\gg1$, which means that $f(a,z)$ rapidly declines, by orders of magnitude, beyond that point. Equations (\ref{Pro3}) and (\ref{Pro4}) are implemented in ProDiMo for the case of hydrostatic models with \citet{Riols2018}-settling. The integral in Eq.\,(\ref{Pro3}) is solved numerically from bottom to top in each column for each dust size bin. If we rather believe in the equations explained in this article, we would have \begin{equation} \ln\frac{f(a,z)\cs^2(z)}{f(a,0)\cs^2(0)} ~=~ \ln\frac{p(z)}{p(0)} - \frac{a}{\alpha} \int_0^z z'\,\frac{\Omega^3\,\rho_{\rm m}(z')} {v_{\rm th}(z')\,\cs^2(z')\,\rho(z')} \big(1+\St(a,z')\big)\,dz' \end{equation} with $\St=\alpha\Omega\tstop$, which, because of the additional $(1+St)$-factor, would lead to even more extreme settling results. \begin{figure*} \vspace*{-3mm} \begin{tabular}{cc} Dubrulle $\alpha=0.01$ & Riols\,\&\,Lesur $\alpha=0.01$\\[-1mm] \hspace*{-2mm}\includegraphics[width=80mm]{epsfig/nH_settle2.pdf} & \hspace*{-7mm}\includegraphics[width=80mm]{epsfig/nH_settle3.pdf} \\[-3mm] \hspace*{-2mm}\includegraphics[width=80mm]{epsfig/dust_settle2.pdf} & \hspace*{-7mm}\includegraphics[width=80mm]{epsfig/dust_settle3.pdf} \\[-3mm] \hspace*{-2mm}\includegraphics[width=80mm]{epsfig/dg_settle2.pdf} & \hspace*{-7mm}\includegraphics[width=80mm]{epsfig/dg_settle3.pdf} \\[-3mm] \end{tabular} \caption{\footnotesize Results for the \citet{Dubrulle1995} and the \citet{Riols2018} settling, both with $\alpha=0.01$. The model is for a MMSN-type hydrostatic disk around an 1\,Myr old T\,Tauri star ($M=1\,M_\odot$, $L_\star=2.24\,L_\odot$, $T_{\rm eff}=4280$\,K) with $\Sigma_{\rm dust}(1\,{\rm au})=1700\rm\,g/cm^2$, $\epsilon=1.5$ and $\dot{M}_{\rm acc}=10^{-8}\rm\,M_\odot/yr$ with viscous heating. The gas density, temperature, and pressure distributions are the same in both models, they are taken from the consistent model with Dubrulle-settling, but then frozen for the right model with Riols\,\&\,Lesur settling.} \end{figure*} \bigskip \bibliographystyle{chicagon} \bibliography{references} \end{document} In comparison to our definitions, $1/\alpha$ is missing in their definition of $\tedd'$ and hence the factor $1/\alpha$ is missing in their definition of the Stokes number $\St'$. Also they use $\cs$ instead of $\vth$ in their $\tstop'$. \begin{eqnarray} \tedd &=& \frac{1}{\alpha}\tedd' \\ \tstop(0) &=& \sqrt{\frac{\pi}{8}}\,\tstop'(0) \\ \St &=& \alpha\sqrt{\frac{\pi}{8}}\,\St' \end{eqnarray}