Such an approximation follows by the Cornish--Fisher asymptotic expansion for the quantiles of a probability distribution.
The error estimate can be obtained by noting that here the Cornish--Fisher asymptotic expansion can be obtained by inverting the Edgeworth asymptotic expansion for the distribution function of the sum of i.i.d. random variables.
I have finally got around to supplying the straightforward but tedious details.
Recall that $\chi^2_n$ is the distribution of the sum of $n$ i.i.d. copies of $Z^2$, where $Z\sim N(0,1)$. So, if $S_n\sim\chi^2_n$,
we can use the Edgeworth asymptotic expansion for the distribution function, say $F_n$, of $S_n$ as e.g. presented (with an explicit expression for the remainder) in Theorem 1 of Section VI.3 of Petrov's book, so that
\begin{equation*}
F_n(x)=\Phi(z)+\phi(z)\sum_{j=1}^{k-2}P_j(z)n^{-j/2}+O(n^{-(k-1)/2}) \tag{10}\label{10}
\end{equation*}
uniformly in real $x$, where $z:=\dfrac{x-n}{\sqrt{2n}}$;
$\Phi$ and $\phi$ are the c.d.f. and p.d.f. of $N(0,1)$; $k$ is any integer $\ge3$; and the $P_j$'s are certain polynomials with coefficients depending on the moments of the underlying distribution (in our case, the moments of the distribution of $Z^2$).
In our case, after long calculations, we have \eqref{10} with $k=7$ and
\begin{equation*}
\begin{aligned}
P_1(z)&=\frac{\sqrt 2(1-z^2)}3, \\
P_2(z)&=-\frac{ z (3 - 11 z^2 + 2 z^4)}{18}, \\
P_3(z)&=\frac{\sqrt{2} \left(3+69 z^2-399 z^4+145 z^6-10 z^8\right)}{810}, \\
P_4(z)&=\frac{z \left(135-1035 z^2+7947 z^4-4167 z^6+560 z^8-20 z^{10}\right)}{9720}, \\
P_5(z)&=-\frac{\sqrt{2}}{204120}\,(3375+3105 z^2-16470 z^4+142218 z^6 \\
&\qquad\qquad\quad-94941 z^8+18501 z^{10}-1288 z^{12}+28
z^{14}),
\end{aligned}
\tag{20}\label{20}
\end{equation*}
with a remainder $O(n^{-(7-1)/2})=O(n^{-3})$.
The median, say $m_n$, of $\chi^2_n$ is the unique root $x$ of the equation $F_n(x)=1/2$. Using now \eqref{10} and \eqref{20}, expanding in powers of $n^{-1/2}$, and letting
\begin{equation*}
m_{n;4}:=n+\sum_{j=-1}^4 c_j n^{-j/2}
\end{equation*}
for some real constants $c_j$, we see that
\begin{equation*}
m_n=m_{n;4}+O(n^{-5/2})
\end{equation*}
iff
\begin{equation*}
(c_{-1},\dots,c_4)=\Big(0,-\frac23,0,\frac{32}{405},0,\frac{1472}{25515}\Big),
\end{equation*}
so that
\begin{equation*}
m_{n;4}=n-\frac{2}{3}+\frac{32}{405 n}+\frac{1472}{25515 n^2}.
\end{equation*}
The approximation of the median $m_n$ in question was
\begin{equation*}
\tilde m_{n;4}:=n\Big(1-\frac 2{9n}\Big)^3=n-\frac{2}{3}+\frac{4}{27 n}-\frac{8}{729 n^2}.
\end{equation*}
We see that already the coefficient of $n^{-1}$ in $\tilde m_{n;4}$ is not the best possible one. So, there can hardly be a solid reason for the estimate $\tilde m_{n;4}$ of the median $m_n$.
The approximation $m_{n;4}$ of the median $m_n$ is much better than $\tilde m_{n;4}$ even for $n$ as small as $5$ or $10$. E.g.,
\begin{equation*}
m_{5;4}-m_5\approx-0.0000167\quad\text{versus}\quad \tilde m_{5;4}-m_5\approx0.0111,
\end{equation*}
\begin{equation*}
m_{10;4}-m_{10}\approx-6.28\times10^{-6}\quad\text{versus}\quad \tilde m_{10;4}-m_{10}\approx0.00622.
\end{equation*}
For somewhat related results, see this paper and references there. (Recall that the chi-squared distribution is a subspecies of the gamma family of distributions).
