?This is a Mathematica notebook file associated with B. Pittel and D. Poole's paper "Asymptotic distribution of the numbers of vertices and arcs of the giant strong component in sparse random digraphs" ArXiv 1405.4022.

?In this file, you can evaluate the covariance matrix B=B(c) of the centered and properly scales pair of the numbers of vertices and arcs of the largest strong component of D(n,m=cn), c>1.

?First, let's find your theta value for your desired arc density c

THETA[c_]:=FindRoot[1-x-E^(-c*x),{x,0.5}][[1]][[2]]

?The following is a long list of functions that we'll need to evaluate B=PSI+EXTRA. Note that PSI is 2x2 matrix whose entries are given in integral form and EXTRA can be solved explicitly for (albeit a little lengthly)

dzia[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=-(((E^(-zi)) (((-1+(E^zi)))^2) g)/((((a-bi))^2) (-1+(E^zi)-zi)))
dzibi[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((E^(-zi)) (((-1+(E^zi)))^2) g)/((((a-bi))^2) (-1+(E^zi)-zi)))
dzibo[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=0
dzig[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((E^(-zi)) (((-1+(E^zi)))^2))/((a-bi) (-1+(E^zi)-zi)))
dzoa[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=-(((E^(-zo)) (((-1+(E^zo)))^2) g)/((((a-bo))^2) (-1+(E^zo)-zo)))
dzobi[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=0
dzobo[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((E^(-zo)) (((-1+(E^zo)))^2) g)/((((a-bo))^2) (-1+(E^zo)-zo)))
dzog[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((E^(-zo)) (((-1+(E^zo)))^2))/((a-bo) (-1+(E^zo)-zo)))
dzsa[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((g/((a-bi) (a-bo)))-(((a-bi-bo) g)/((a-bi) (((a-bo))^2)))-(((a-bi-bo) g)/((((a-bi))^2) (a-bo))))/(((E^zs)/(-1+(E^zs)))-(((E^(2 zs)) zs)/(((-1+(E^zs)))^2))+(((E^zs) zs)/(-1+(E^zs)))))
dzsbi[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((-(g/((a-bi) (a-bo)))+(((a-bi-bo) g)/((((a-bi))^2) (a-bo))))/(((E^zs)/(-1+(E^zs)))-(((E^(2 zs)) zs)/(((-1+(E^zs)))^2))+(((E^zs) zs)/(-1+(E^zs)))))
dzsbo[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((-(g/((a-bi) (a-bo)))+(((a-bi-bo) g)/((a-bi) (((a-bo))^2))))/(((E^zs)/(-1+(E^zs)))-(((E^(2 zs)) zs)/(((-1+(E^zs)))^2))+(((E^zs) zs)/(-1+(E^zs)))))
dzsg[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(((a-bi-bo) (E^(-zs)) (((-1+(E^zs)))^2))/((a-bi) (a-bo) (-1+(E^zs)-zs)))
Exa2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi/(bi+bo)+bo(E^(zi))(bi*zi*zo)(1+(bi*zi*zo)/((E^(zo)-1)g))/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Exb2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo/(bi+bo)+bi(E^(zo))(bo*zi*zo)(1+(bo*zi*zo)/((E^(zi)-1)g))/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)

Exab[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*E^(zo)bo*zi*zo/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)+bo*E^(zi)bi*zi*zo/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Exar[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*E^(zo)*(a-bi-bo)*zi*zo/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)

Exat[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo*E^(zi)bi*zi*zo(a-bi-bo)zi*zo/((bi+bo)(E^(zi)-1)(E^(zo)-1)g(E^(zo)-1)g)

Exak[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*zo*E^(zo)/((bi+bo)(E^(zo)-1))+bo*E^(zi)bi*zi*zo(1+zi)/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Exbr[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*E^(zo)bo*zi*zo(a-bi-bo)zi*zo/((bi+bo)(E^(zo)-1)(E^(zi)-1)g(E^(zi)-1)g)

Exbk[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo*zi*E^(zi)/((bi+bo)(E^(zi)-1))+bi*E^(zo)bo*zi*zo(1+zo)/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)

Exbt[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo*E^(zi)*(a-bi-bo)*zi*zo/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Exr2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*E^(zo)(a-bi-bo)zi*zo(1+(a-bi-bo)zi*zo/((E^(zi)-1)g))/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)

Exrt[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=0

Exrk[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*E^(zo)(a-bi-bo)zi*zo(1+zo)/((bi+bo)(E^(zo)-1)(E^(zi)-1)g)

Ext2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo*E^(zi)(a-bi-bo)zi*zo(1+(a-bi-bo)zi*zo/((E^(zo)-1)g))/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Extk[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bo*E^(zi)(a-bi-bo)zi*zo(1+zi)/((bi+bo)(E^(zi)-1)(E^(zo)-1)g)

Exk2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=bi*zo*E^(zo)(1+zo)/((bi+bo)(E^(zo)-1))+bo*zi*E^(zi)(1+zi)/((bi+bo)(E^(zi)-1))
f2[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=g(zs)^2/(zi*zo)
df2a[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=2g(zs)/(zi*zo)dzsa[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi^2*zo)dzia[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi*zo^2)dzoa[{a,bi,bo,g,zi,zo,zs}]

df2bi[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=2g*zs/(zi*zo)dzsbi[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi^2*zo)dzibi[{a,bi,bo,g,zi,zo,zs}]

df2bo[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=2g*zs/(zi*zo)dzsbo[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi*zo^2)dzobo[{a,bi,bo,g,zi,zo,zs}]

df2g[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=(zs)^2/(zi*zo)+2g(zs)/(zi*zo)dzsg[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi^2*zo)dzig[{a,bi,bo,g,zi,zo,zs}]-g(zs)^2/(zi*zo^2)dzog[{a,bi,bo,g,zi,zo,zs}]

df1a[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((a-bi)(a-bo))/(g(a-bi-bo))*df2a[{a,bi,bo,g,zi,zo,zs}]+f2[{a,bi,bo,g,zi,zo,zs}]*(-(((a-bi) (a-bo))/((((a-bi-bo))^2) g))+((a-bi)/((a-bi-bo) g))+((a-bo)/((a-bi-bo) g)))

df1bi[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((a-bi)(a-bo))/(g(a-bi-bo))*df2bi[{a,bi,bo,g,zi,zo,zs}]+f2[{a,bi,bo,g,zi,zo,zs}]*((((a-bi) (a-bo))/((((a-bi-bo))^2) g))-((a-bo)/((a-bi-bo) g)))

df1bo[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((a-bi)(a-bo))/(g(a-bi-bo))*df2bo[{a,bi,bo,g,zi,zo,zs}]+f2[{a,bi,bo,g,zi,zo,zs}]*((((a-bi) (a-bo))/((((a-bi-bo))^2) g))-((a-bi)/((a-bi-bo) g)))

df1g[{a_,bi_,bo_,g_,zi_,zo_,zs_}]:=((a-bi)(a-bo))/(g(a-bi-bo))*df2g[{a,bi,bo,g,zi,zo,zs}]+f2[{a,bi,bo,g,zi,zo,zs}]*(-(((a-bi) (a-bo))/((a-bi-bo) (g^2))))
f2a[c_,th_,z_]:=df2a[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f2bi[c_,th_,z_]:=df2bi[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f2bo[c_,th_,z_]:=df2bo[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f2g[c_,th_,z_]:=df2g[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f1a[c_,th_,z_]:=df1a[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f1bi[c_,th_,z_]:=df1bi[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f1bo[c_,th_,z_]:=df1bo[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

f1g[c_,th_,z_]:=df1g[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]
Ea2[c_,th_,z_]:=Exa2[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Eb2[c_,th_,z_]:=Exb2[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Er2[c_,th_,z_]:=Exr2[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Et2[c_,th_,z_]:=Ext2[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ek2[c_,th_,z_]:=Exk2[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Eab[c_,th_,z_]:=Exab[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ear[c_,th_,z_]:=Exar[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Eat[c_,th_,z_]:=Exat[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Eak[c_,th_,z_]:=Exak[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ebr[c_,th_,z_]:=Exbr[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ebt[c_,th_,z_]:=Exbt[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ebk[c_,th_,z_]:=Exbk[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Ert[c_,th_,z_]:=Exrt[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Erk[c_,th_,z_]:=Exrk[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]

Etk[c_,th_,z_]:=Extk[{(E^z-1)/(c*E^(2z))(2z*E^z-c(E^z-1)), (E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)),(E^z-1)/(c*E^(2z))(z*E^z-c(E^z-1)), z^2/c,z, z,c*th}]
EDf1Df1[c_,th_,z_]:=(Ea2[c,th,z]+2 Eab[c,th,z]+Eb2[c,th,z])(f1a[c,th,z])(f1a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1bi[c,th,z])(f1a[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f1bo[c,th,z])(f1a[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f1g[c,th,z])(f1a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1a[c,th,z])(f1bi[c,th,z])+(Er2[c,th,z]-2Ear[c,th,z]+Ea2[c,th,z])(f1bi[c,th,z])(f1bi[c,th,z])+(Ert[c,th,z]-Ebr[c,th,z]-Eat[c,th,z]+Eab[c,th,z])(f1bo[c,th,z])(f1bi[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f1g[c,th,z])(f1bi[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f1a[c,th,z])(f1bo[c,th,z])+(Ert[c,th,z]-Eat[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1bi[c,th,z])(f1bo[c,th,z])+(Et2[c,th,z]-2Ebt[c,th,z]+Eb2[c,th,z])(f1bo[c,th,z])(f1bo[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f1g[c,th,z])(f1bo[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f1a[c,th,z])(f1g[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f1bi[c,th,z])(f1g[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f1bo[c,th,z])(f1g[c,th,z])+(Ek2[c,th,z])(f1g[c,th,z])(f1g[c,th,z])
EDf2Df2[c_,th_,z_]:=(Ea2[c,th,z]+2 Eab[c,th,z]+Eb2[c,th,z])(f2a[c,th,z])(f2a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f2bi[c,th,z])(f2a[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f2bo[c,th,z])(f2a[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f2g[c,th,z])(f2a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f2a[c,th,z])(f2bi[c,th,z])+(Er2[c,th,z]-2Ear[c,th,z]+Ea2[c,th,z])(f2bi[c,th,z])(f2bi[c,th,z])+(Ert[c,th,z]-Ebr[c,th,z]-Eat[c,th,z]+Eab[c,th,z])(f2bo[c,th,z])(f2bi[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f2g[c,th,z])(f2bi[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f2a[c,th,z])(f2bo[c,th,z])+(Ert[c,th,z]-Eat[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f2bi[c,th,z])(f2bo[c,th,z])+(Et2[c,th,z]-2Ebt[c,th,z]+Eb2[c,th,z])(f2bo[c,th,z])(f2bo[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f2g[c,th,z])(f2bo[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f2a[c,th,z])(f2g[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f2bi[c,th,z])(f2g[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f2bo[c,th,z])(f2g[c,th,z])+(Ek2[c,th,z])(f2g[c,th,z])(f2g[c,th,z])
EDf1Df2[c_,th_,z_]:=(Ea2[c,th,z]+2 Eab[c,th,z]+Eb2[c,th,z])(f1a[c,th,z])(f2a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1bi[c,th,z])(f2a[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f1bo[c,th,z])(f2a[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f1g[c,th,z])(f2a[c,th,z])+(-Ear[c,th,z]+Ea2[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1a[c,th,z])(f2bi[c,th,z])+(Er2[c,th,z]-2Ear[c,th,z]+Ea2[c,th,z])(f1bi[c,th,z])(f2bi[c,th,z])+(Ert[c,th,z]-Ebr[c,th,z]-Eat[c,th,z]+Eab[c,th,z])(f1bo[c,th,z])(f2bi[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f1g[c,th,z])(f2bi[c,th,z])+(-Eat[c,th,z]+Eab[c,th,z]-Ebt[c,th,z]+Eb2[c,th,z])(f1a[c,th,z])(f2bo[c,th,z])+(Ert[c,th,z]-Eat[c,th,z]-Ebr[c,th,z]+Eab[c,th,z])(f1bi[c,th,z])(f2bo[c,th,z])+(Et2[c,th,z]-2Ebt[c,th,z]+Eb2[c,th,z])(f1bo[c,th,z])(f2bo[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f1g[c,th,z])(f2bo[c,th,z])+(Eak[c,th,z]+Ebk[c,th,z])(f1a[c,th,z])(f2g[c,th,z])+(-Erk[c,th,z]+Eak[c,th,z])(f1bi[c,th,z])(f2g[c,th,z])+(-Etk[c,th,z]+Ebk[c,th,z])(f1bo[c,th,z])(f2g[c,th,z])+(Ek2[c,th,z])(f1g[c,th,z])(f2g[c,th,z])
psi[c_,th_,z_]:={{EDf1Df1[c,th,z]2(1-E^(-z))/c, EDf1Df2[c,th,z]2(1-E^(-z))/c},{EDf1Df2[c,th,z]2(1-E^(-z))/c,EDf2Df2[c,th,z]2(1-E^(-z))/c}}

PSI[c_]:=NIntegrate[psi[c,THETA[c],z],{z,c*THETA[c],c}]
R[c_]:={{c*E^(2c)/((E^c-1)(E^c-1-c))+E^(2c)+E^(2c)/((E^c-1)^2), -c*E^(2c)/(2(E^c-1)(E^c-1-c))-E^(2c)/(E^c-1)^2, -c*E^(2c)/(2(E^c-1)(E^c-1-c))-E^(2c)/(E^c-1)^2},{-c*E^(2c)/(2(E^c-1)(E^c-1-c))-E^(2c)/(E^c-1)^2,c*E^(2c)/(2(E^c-1)(E^c-1-c))+E^(3c)/(E^c-1)^2,E^(2c)/(E^c-1)^2},{-c*E^(2c)/(2(E^c-1)(E^c-1-c))-E^(2c)/(E^c-1)^2,E^(2c)/(E^c-1)^2,c*E^(2c)/(2(E^c-1)(E^c-1-c))+E^(3c)/(E^c-1)^2}}
K1[c_]:=Inverse[R[c]]
vf1a[c_,th_]:=df1a[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf1bi[c_,th_]:=df1bi[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf1bo[c_,th_]:=df1bo[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf1g[c_,th_]:=df1g[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf2a[c_,th_]:=df2a[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf2bi[c_,th_]:=df2bi[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf2bo[c_,th_]:=df2bo[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
vf2g[c_,th_]:=df2g[{1-E^(-2c),E^(-c)-E^(-2c),E^(-c)-E^(-2c),c,c,c,c*th}]
M[c_,th_]:={{vf1a[c,th],vf1bi[c,th],vf1bo[c,th]},{vf2a[c,th],vf2bi[c,th],vf2bo[c,th]}}
Mt[c_,th_]:=M[c,th].K1[c].Transpose[M[c,th]]
EXTRA[c_]:=Mt[c,THETA[c]]

B[c_]:=PSI[c]+EXTRA[c]

?Phew! Now we can finally compute the covariance matrix B. For example, I find the covariance matrix for c=ln(4) (in this case theta=1/2).

B[Log[4.]]

MatrixForm[%]

?Let's see how close this comes to simulation data from Mathematica. The following two functions are used to generate pairs of these parameters for many instances of D(n,m=cn).

GetLargestComponent[gra_]:=Module[{list,size,result},
list=ConnectedComponents[gra];
size=0;
For[i=1,i≤Length[list],i++,
If[size<Length[list[[i]]],
size=Length[list[[i]]];
result=list[[i]]]
];
Return[result];
];

ListofSizesandOrders[n_,m_,trials_]:=Module[{graph, verts,result},
result=Table[0,{j,1,trials},{l,1,2}];
For[k=1,k≤trials,k++,
graph=RandomGraph[{n,m},DirectedEdges->True];
verts=GetLargestComponent[graph];
result[[k]][[1]]=Length[verts];
result[[k]][[2]]=EdgeCount[Subgraph[graph,verts]];
];
Return[result];
];

?Having defined those, let's look at sample covariance for D(10000, ln(4) 10000 approx 13863) with 10000 trials. This simulation takes about 4 minutes on an iMac with 3.2 GHz processor

LIST=ListofSizesandOrders[10000,13863,10000]
MatrixForm[Covariance[LIST]/10000.]

?We used Mathematica to determine the asymptotics of the entries of PSI(c) as c decreases to 1. Alternatively, we can consider finding the asymptotics as theta decreases to zero. Again, theta and c are related by
?? 1-theta = e^(-c theta) or equivalently
??c=ln(1-theta) / (-theta)

?The simplest way to view the asymptotics of PSI(c) are through a related covariance matrix A(c) determined by PSI(c). For instance,
??A_1,1 = B_1,1
?? A_1,2 = A_2,1 = B_1,2 - B_1,1
?? A_2,2 = B_1,1 + B_2,2 - 2B_1,2.
?? A is the covariance matrix of (number of vertices of giant , number of arcs - number of vertices of giant). This latter term is often called the excess. The reason for this is that the leading order terms of number of vertices and number of arcs of the giant are asymptotically the same as c decreases to 1.

?One can integrate the following three series over z in [c theta, c] to get the 1st three orders (theta, theta^2, theta^3). You will need the integrals of the subsequent three functions as well.

Series[EDf1Df1[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,8},{z,0,1}]
Series[EDf1Df2[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,8},{z,0,1}]
Series[EDf2Df2[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,8},{z,0,1}]

Series[EDf1Df1[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,3}]
Series[EDf1Df2[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,3}]
Series[EDf2Df2[Log[1-th]/(-th),th,z]2(1-E^(-z))/(Log[1-th]/(-th)),{th,0,3}]

?Note the following relation between these functions:

h[z_]:=((32 (-1+3 (E^z)-3 (E^(2 z))+(E^(3 z))-z+2 (E^z) z-(E^(2 z)) z+(E^z) (z^2)-(E^(2 z)) (z^2)+(E^(2 z)) (z^3)) )/((((-1+(E^z)))^3) (z^3)))
g11[z_]:=-((16 (-7+21 (E^z)-21 (E^(2 z))+7 (E^(3 z))-4 z+8 (E^z) z-4 (E^(2 z)) z+4 (E^z) (z^2)-4 (E^(2 z)) (z^2)+7 (E^(2 z)) (z^3)) )/(3 ((((-1+(E^z)))^3) (z^3))))
g12[z_]:=-((8 (-5+15 (E^z)-15 (E^(2 z))+5 (E^(3 z))+z-2 (E^z) z+(E^(2 z)) z-(E^z) (z^2)+(E^(2 z)) (z^2)+5 (E^(2 z)) (z^3)) )/(3 ((((-1+(E^z)))^3) (z^3))))
g22[z_]:=((16 (-2+6 (E^z)-6 (E^(2 z))+2 (E^(3 z))-5 z+10 (E^z) z-5 (E^(2 z)) z+5 (E^z) (z^2)-5 (E^(2 z)) (z^2)+2 (E^(2 z)) (z^3)) )/(3 (((-1+(E^z)))^3) (z^3)))
Simplify[g11[z]+g22[z]-2g12[z]]

?The leading order terms are the following:
??A_1,1 = 20 (THETA)+ O(THETA^2)
?? A_1,2 = 15 (THETA)^2 + O(THETA^3)
??A_2,2 = 34/3 (THETA)^3 + O(THETA^4)

? The leading order terms arising from EXTRA will be drowned out by these leading order terms.

Series[Mt[Log[1-th]/(-th),th][[1]][[1]],{th,0,3}]
Series[Mt[Log[1-th]/(-th),th][[1]][[2]]-Mt[Log[1-th]/(-th),th][[1]][[1]],{th,0,3}]
Series[Mt[Log[1-th]/(-th),th][[1]][[1]]+Mt[Log[1-th]/(-th),th][[2]][[2]]-2*Mt[Log[1-th]/(-th),th][[1]][[2]],{th,0,4}]