Solution (source code)

= Solution

Assume a steady spherical outflow of fixed composition. Let $\bar m$ be its mean particle <mass>, and define the constant <isothermal sound speed> by $a^2=k_BT_0/\bar m$. Subsonic here refers to this isothermal closure, so $0<u<a$. <Mass conservation> and radial momentum balance give
$$
\dot M=4\pi r^2\rho u,\qquad
u\frac{du}{dr}=-a^2\frac{d\ln\rho}{dr}-\frac{GM}{r^2}.
$$
Logarithmically differentiating the first equation yields $d\ln\rho/dr=-2/r-u'/u$. Substitution eliminates the density gradient and gives the <isothermal Parker wind> equation
$$
\left(u-\frac{a^2}{u}\right)\frac{du}{dr}
=\frac{2a^2}{r}-\frac{GM}{r^2}.
$$
At a smooth speed maximum, $u'=0$, while the left prefactor is finite and negative. The right side must vanish, hence
$$
\boxed{r_M=\frac{GM}{2a^2}=\frac{GM\bar m}{2k_BT_0}.}
$$
For $r<r_M$ both the right side and the prefactor are negative, so $u'>0$; for $r>r_M$, $u'<0$. This is a true maximum on a <subsonic isothermal breeze>, not the accelerating <transonic branch>.

There is an underdetermination in the requested density estimate. Introduce the peak <Mach number> $\mathcal M_M=u_M/a$, where $0<\mathcal M_M<1$. The <density at the maximum speed of an isothermal breeze> is exactly
$$
\boxed{\rho_M=\frac{\dot M}{4\pi r_M^2u_M}
=\frac{\dot M a^3}{\pi G^2M^2\mathcal M_M}
=\frac{\dot M}{\pi G^2M^2\mathcal M_M}
\left(\frac{k_BT_0}{\bar m}\right)^{3/2}.}
$$
The hypotheses do not fix $\mathcal M_M$. A boundary density or <pressure>, or an equivalent wind integration constant, is needed. The definite statement following from subsonicity is the lower bound
$$
\boxed{\rho_M>\rho_{\rm char},\qquad
\rho_{\rm char}=\frac{\dot M}{\pi G^2M^2}
\left(\frac{k_BT_0}{\bar m}\right)^{3/2}.}
$$
If a maximum speed of order $a$ is intended as an extra estimate, then $\rho_M\sim\rho_{\rm char}\propto\dot M T_0^{3/2}M^{-2}$. A very slow breeze can have a much larger density, so setting $u_M=a$ as an exact equality would contradict the stated always-subsonic assumption. The mean particle mass is another implicit fixed-composition input needed to express a <temperature> as a squared sound speed.

One can see the free parameter directly by integrating the wind equation:
$$
\frac{u^2}{2}-a^2\ln\frac ua
=2a^2\ln\frac r{r_M}+\frac{GM}{r}+C.
$$
Different $C$ give different subsonic maxima at the same radius. The original plot illustrates this nonuniqueness. These breeze branches require nonzero ambient <pressure> at large radius; they are not the usual transonic escape into a low-pressure exterior.

\Image[/past-exam-of-the-astrophysical-fluid-dynamics-course-of-the-university-of-cambridge/2017-subsonic-breezes.png]
{title=Different subsonic breeze maxima at the same critical radius demonstrate the free density parameter}

Maintaining an isothermal expanding wind generally requires heat exchange; isothermal motion is not adiabatic motion.https://nnp.ucsd.edu/the_mhd_muddle/PDF/Parker_1958_Dynamics_Gas_MagneticFields_HL.pdf[PDF p. 4, printed p. 667, section III: uniform temperature maintained by heating]{ref}