\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}  
