Variance of Cohen's $d$ statistic Cohen’s $d$ is one of the most common ways we measure the size of an effect (see Wikipedia). It simply measures the distance between two means in terms of the pooled standard deviation. How can we derive the mathematical formula of variance estimate of Cohen's $d$?  
December 2015 edit: Related to this question is the idea of calculating confidence intervals around $d$. This article states that 
$$\sigma_{d}^2 = \dfrac{n_{+}}{n_{\times}} + \dfrac{d^2}{2n_{+}} $$
where $n_{+}$ is the sum of the two sample sizes and $n_{\times}$ is the product of the two sample sizes. 
How is this formula derived?
 A: Note that the variance expression in the question is an approximation.  Hedges (1981) derived the large sample variance of $d$ and approximation in a general setting (i.e. multiple experiments/studies), and my answer pretty much walks through the derivations in the paper. 
First, the assumptions we'll utilize are the following:
Let's assume we have two independent treatment groups, $T$ (treatment) and $C$ (control). Let $Y_{Ti}$ and $Y_{Cj}$ be the scores/responses/whatever from subject $i$ in group $T$ and subject $j$ in group $C$, respectively.
We assume the responses are normally distributed and the treatment and control groups share a common variance, i.e.
\begin{align*}
Y_{Ti} &\sim N(\mu_T, \sigma^2), \quad i = 1, \dots n_T \\
Y_{Cj} &\sim N(\mu_C, \sigma^2), \quad j = 1, \dots n_C
\end{align*}
The effect size we're interested in estimating in each study is $\delta = \frac{\mu_T - \mu_C}{\sigma}$.  The estimator of the effect size we'll use is
\begin{equation*}
d = \frac{\bar{Y}_T - \bar{Y}_C}{\sqrt{\frac{(n_T - 1)S_T^2 + (n_C - 1)S_C^2}{n_T + n_C - 2}}}
\end{equation*}
where $S_k^2$ is the unbiased sample variance for group $k$.  
Let's consider the large-sample properties of $d$.  
First, note that:
\begin{equation*}
\bar{Y}_T - \bar{Y}_C \sim N \Bigg( \mu_T - \mu_C, \,\sigma^2\frac{n_T + n_C}{n_T n_C} \Bigg)
\end{equation*} 
and (being loose with my notation):
\begin{equation}
\frac{(n_T - 1)S_T^{2}}{\sigma^2(n_T + n_C - 2)} = \frac{1}{n_T + n_C - 2}\frac{(n_T - 1)S_T^{2}}{\sigma^2} \sim \frac{1}{n_T + n_C- 2}\chi_{n_T - 1}^2 \tag{1}
\end{equation}
and
\begin{equation}
\frac{(n_C - 1)S_C^{2}}{\sigma^2(n_T + n_C - 2)} = \frac{1}{n_T + n_C - 2}\frac{(n_C - 1)S_C^{2}}{\sigma^2} \sim \frac{1}{n_T + n_C- 2}\chi_{n_C - 1}^2 \tag{2}
\end{equation}
Equations (1) and (2) lead to the fact that (again, being loose with my notation):
\begin{equation*}
\frac{1}{\sigma^2}\frac{(n_T - 1)S_T^{2} + (n_C - 1)S_C^{2}}{n_T + n_C - 2} \sim \frac{1}{n_T + n_C - 2}\chi_{n_T + n_C - 2}^2
\end{equation*}
Now, some clever algebra:
\begin{align*}
d &= \frac{\bar{Y}_T - \bar{Y}_C}{\sqrt{\frac{(n_T - 1)S_T^2 + (n_C - 1)S_C^2}{n_T + n_C - 2}}} \\\\
&= \frac{\left(\sigma\sqrt{\frac{n_T + n_C}{n_T n_C}}\right)^{-1}(\bar{Y}_T - \bar{Y}_C)}{\left(\sigma\sqrt{\frac{n_T + n_C}{n_T n_C}}\right)^{-1}\sqrt{\frac{(n_T - 1)S_T^2 + (n_C - 1)S_C^2}{n_T + n_C - 2}}} \\\\
&= \frac{\frac{(\bar{Y}_T - \bar{Y}_C) - (\mu_T - \mu_C)}{\sigma\sqrt{\frac{n_T + n_C}{n_T n_C}}} + \frac{\mu_T - \mu_C}{\sigma\sqrt{\frac{n_T + n_C}{n_T n_C}}}}{\left(\sqrt{\frac{n_T + n_C}{n_T n_C}}\right)^{-1}\sqrt{\frac{(n_T - 1)S_T^2 + (n_C - 1)S_C^2}{\sigma^2(n_T + n_C - 2)}}} \\\\
&= \sqrt{\frac{n_T + n_C}{n_T n_C}}\left(\frac{\theta + \delta\sqrt{\frac{n_T n_C}{n_T + n_C}}}{\sqrt{\frac{V}{\nu}}}\right)
\end{align*}
where $\theta \sim N(0,1)$, $V \sim \chi^2_{\nu}$, and $\nu = n_T+n_C-2$. Thus, $d$ is $\sqrt{\frac{n_T + n_C}{n_T n_C}}$ times a variable which follows a non-central t-distribution with $n_T + n_C - 2$ degrees of freedom and non-centrality parameter of $\delta\sqrt{\frac{n_T n_C}{n_T + n_C}}$.
Using the moment properties of the non-central $t$ distribution, it follows that:
\begin{equation*}
\mathrm{Var}(d) = \frac{(n_T + n_C - 2)}{(n_T + n_C - 4)}\frac{(n_T + n_C)}{n_T n_C}(1+ \delta^2\frac{n_T n_C}{n_T + n_C}) - \frac{\delta^2}{b^2} \tag{3}
\end{equation*}
where 
\begin{equation*}
b = \frac{\Gamma\left(\frac{n_T + n_C - 2}{2}\right)}{\sqrt{\frac{n_T+n_C-2}{2}}\Gamma\left(\frac{n_T+n_C-3}{2}\right)} \approx 1 - \frac{3}{4(n_T+n_C-2)-1}
\end{equation*}
So Equation (3) provides the exact large sample variance.  Note that an unbiased estimator for $\delta$ is $b d$, with variance:
\begin{equation*}
\mathrm{Var}(bd) = b^2\frac{(n_T + n_C - 2)}{(n_T + n_C - 4)}\frac{(n_T + n_C)}{n_T n_C}(1+ \delta^2\frac{n_T n_C}{n_T + n_C}) - \delta^2 
\end{equation*}
For large degrees of freedom (i.e. large $n_T+n_C-2$), the variance of a non-central $t$ variate with $\nu$ degrees of freedom and non-centrality parameter $p$ can be approximated by $1 + \frac{p^2}{2\nu}$ (Johnson, Kotz, Balakrishnan, 1995).  Thus, we have:
\begin{align*}
\mathrm{Var}(d) &\approx \frac{n_T + n_C}{n_T n_C}\left(1 + \frac{\delta^2\left(\frac{n_T n_C}{n_T + n_C}\right)}{2(n_T+n_C-2)}\right)
 \\\\
&= \frac{n_T + n_C}{n_T n_C} + \frac{\delta^2}{2(n_T+n_C-2)}
\end{align*}
Plug in our estimator for $\delta$ and we're done.
