Sampling distribution from two independent Bernoulli populations Let's assume that we have samples of two independent Bernoulli random variables, $\mathrm{Ber}(\theta_1)$ and $\mathrm{Ber}(\theta_2)$.
How do we prove that $$\frac{(\bar X_1-\bar X_2)-(\theta_1-\theta_2)}{\sqrt{\frac{\theta_1(1-\theta_1)}{n_1}+\frac{\theta_2(1-\theta_2)}{n_2}}}\xrightarrow{d} \mathcal N(0,1)$$?
Assume that $n_1\neq n_2$.
 A: Put 
$a=\frac{\sqrt{\theta_1(1-\theta_1)}}{\sqrt{n_1}}$, $b=\frac{\sqrt{\theta_2(1-\theta_2)}}{\sqrt{n_2}}$,
$A=(\bar{X}_1-\theta_1)/a$, 
$B=(\bar{X}_2-\theta_2)/b$. We have
$A\to_d N(0,1),\ B\to_d N(0,1)$. 
In terms of characteristic functions it means
$$ \phi_A(t)\equiv {\bf E}e^{itA}\to e^{-t^2/2},\ \phi_B(t)\to e^{-t^2/2}.$$
We want to prove that 
$$D:=\frac{a}{\sqrt{a^2+b^2}}A-\frac{b}{\sqrt{a^2+b^2}}B\to_d N(0,1) $$
Since $A$ and $B$ are independent,
$$\phi_D(t)=\phi_A\left(\frac{a}{\sqrt{a^2+b^2}}t\right)\phi_B\left(-\frac{b}{\sqrt{a^2+b^2}}t\right)\to e^{-t^2/2},
$$
as we wish it to be.
This proof is incomplete. Here we need some estimates for uniform convergence of characteristic functions. However in the case under consideration we can do explicit calculations. Put $p=\theta_1,\ m=n_1$.
\begin{align}
\phi_{X_{1,1}}(t) &= 1+p(e^{it}-1),  \\
\phi_{\bar X_{1}}(t) &= (1+p(e^{it/m}-1))^m,  \\
\phi_{\bar X_{1}-\theta_1}(t) &= (1+p(e^{it/m}-1))^m e^{-ipt},  \\
\phi_{A}(t) &= (1+p(e^{it/\sqrt{mp(1-p)}}-1))^m e^{-ipt\sqrt{m}/\sqrt{p(1-p)}}  \\[5pt]
 &= \left( \left(1+p(e^{it/\sqrt{mp(1-p)}}-1)\right)e^{-ipt/\sqrt{mp(1-p)}}\right)^m \\[5pt]
&=\left( 1-\frac{t^2}{2m}+O(t^3m^{-3/2}) \right)^m 
\end{align}
as $t^3m^{-3/2}\to 0$. Thus, for a fixed $t$,
$$\phi_D(t)=\left( 1-\frac{a^2t^2}{2(a^2+b^2)n_1}+O(n_1^{-3/2}) \right)^{n_1} 
\left( 1-\frac{b^2t^2}{2(a^2+b^2)n_2}+O(n_2^{-3/2}) \right)^{n_2} \to e^{-t^2/2}
$$
(even if $a\to 0$ or $b\to 0$), since $\left|e^{-y}-(1-y/m)^m\right|\le {y^2}/{2m}\ $ when $\ y/m<1/2$ (see https://math.stackexchange.com/questions/2566469/uniform-bounds-for-1-y-nn-exp-y/ ).
Note that similar calculations may be done for arbitrary (not necessarily Bernoulli) distributions with finite second moments, using the expansion of characteristic function in terms of the first two moments.
A: Proving your statement is equivalent to proving the (Levy-Lindenberg) Central Limit Theorem which states
If $\{Z_i\}_{i=1}^n$ is a sequence of i.i.d random variable with finite mean $\mathbb{E}(Z_i) = \mu $ and finite variance $\mathbb{V}(Z_i) = \sigma^2$ then $$ \sqrt{n}(\bar{Z} - \mu) \to^d N(0,\sigma^2)$$
Here $\bar{Z} = \sum_i Z_i/n$ that is the sample variance.
Then it is easy to see that if we put 
$$ Z_i = X_1i - X_2i $$ with $X_{1i}, X_{2i}$ following a $Ber(\theta_1)$ and $Ber(\theta_2)$  respectively the conditions for the theorem are satisfied, in particular
$$ \mathbb{E}(Z_i) = \theta_1 - \theta_2 = \mu $$
and
$$ \mathbb{V}(Z_i)= \theta_1(1-\theta_1) +\theta_2(1-\theta_2)= \sigma^2 $$
(There's a last passage, and you have to adjust this a bit for the general case where $n_1 \neq n_2$ but I have to go now, will finish tomorrow or you can edit the question with the final passage as an exercise )
