Bence Mélykúti, 30/10/2012-4/3/2013

(* Constants: k1=κ_1, km1=κ_-1, kb=κ_b, ku=κ_u;
Variables dependent on tbig=T: mbig=M, dbig=D, gbigst=G^* *)

(* Computing τ and σ *)

eqa= k1 mbig^2-km1 dbig
eqb=ku gbigst - kb dbig gbig
eqaindbiggst=eqa/.{mbig->tbig-2*dbig-2*gbigst} (* computing in terms of D and G^* *)
eqaindbiggst2=Expand[eqaindbiggst];
eqaindbiggst3=Collect[eqaindbiggst2,dbig];
eqaindbiggst4=Simplify[eqaindbiggst3/k1]
eqaindbiggst5=eqaindbiggst3/.{gbigst->((kb dbig)/(kb dbig+ku))gbigtotal};
eqaindbiggst6=Collect[Expand[Simplify[Expand[eqaindbiggst5*(kb dbig+ku)^2]]],dbig] (* this is the quartic equation *)

(* elements of the quadratic formula *)
eqaindbiggst4=eqaindbiggst4/.{gbigst->gbigsthat}
afordbig=Coefficient[eqaindbiggst4,dbig,2]
bfordbig=Coefficient[eqaindbiggst4,dbig,1]
cfordbig=Coefficient[eqaindbiggst4,dbig,0]
bp2a=Expand[(bfordbig/(2afordbig))]
Clear[discr,sqrtdiscrp2a]
discr[tbig_]=Expand[bfordbig^2-4 afordbig cfordbig] (* discriminant: b^2-4ac *)
sqrtdiscrp2a[tbig_]=((Sqrt[discr[tbig]])/(2afordbig))

lengthofexpansion=2

Clear[gbigsthat]
Series[sqrtdiscrp2a[tbig],{tbig,Infinity,lengthofexpansion}]
dbighat=-bp2a-Series[sqrtdiscrp2a[tbig],{tbig,Infinity,lengthofexpansion}]
τ=dbighat;
σ=Simplify[tbig-2 τ-2 gbigsthat];
(* cτ=Sqrt[(km1/(8 k1))] => km1 -> cτ^2*8*k1 *)
$Assumptions={cτ>0}
τ=Simplify[τ/. km1 -> cτ^2*8*k1]
σ=Simplify[σ/. km1 -> cτ^2*8*k1]

This is τ and σ.

(* Computation of the coefficients by Cauchy's coefficient formula and the saddle point I: figures *)

(* the substitution of floating point values for constants *) 
(* values used for figures in paper: *)
(*tbigf=100;k1f=1;km1f=6;kbf=0.1;kuf=2;nf=tbigf+0(* or +2 *);gbigstf=0.5;*)
(* values used for simulation in GNU Octave: *)
tbigf=100;k1f=5/60;km1f=1;kbf=1/6;kuf=5;nf=tbigf+0(* or +2 *);gbigstf=0.0001;
cτf=Sqrt[(km1f/(8 k1f))];
flvalues={tbig->tbigf,k1->k1f,km1->km1f,kb->kbf,ku->kuf,cτ->cτf,gbigst->gbigstf}
τf=(tbig/2)-cτ (Sqrt[tbig])+(((cτ)^2)-gbigst)/.flvalues
σf=2 cτ (Sqrt[tbig])-2 ((cτ)^2)/.flvalues
LogPlot[Abs[Exp[τf*x^2+σf*x]/x^(nf+1)],{x,-2.4,1.9},PlotRange->Automatic]
Plot3D[{Log10[Abs[Exp[τf*(x+I y)^2+σf*(x+I y)]/(x+I y)^(nf+1)]]},{x,-2.4,1.9},{y,-2,2},PlotRange->{-10,99}, ClippingStyle->None,AxesLabel->{Re,Im}]
Plot3D[Log10[Abs[Exp[τf*r^2*Exp[2 I θ]+σf*r*Exp[I θ]-nf*I θ]/r^nf]],{r,10^-7,2.5},{θ,-Pi,Pi},PlotRange->{-10,99}, ClippingStyle->None,AxesLabel->{r,ϑ}]

(* Computation of the coefficients by Cauchy's coefficient formula and the saddle point II: computing the saddle point *)

(* Note: we are solving quadratic equations, therefore there are two solutions for each
r_{(T,1),T-2,1} is rtbigtbigm2[[1]] and r_{(T,1),T-2,2} is rtbigtbigm2[[2]],
 r_{(T,1),T,1} is rtbigtbig[[1]] and r_{(T,1),T,2} is rtbigtbig[[2]] *)
rbp2a=Simplify[(σ/(4 τ))];
Simplify[(1/(4 τ))Sqrt[σ^2 + 8 τ tbig]];
(* r_{(T,1),T-2} *)rtbigtbigm2={Simplify[-rbp2a+Simplify[(1/(4 τ))Sqrt[σ^2 + 8 τ (tbig-2)]]],Simplify[-rbp2a-Simplify[(1/(4 τ))Sqrt[σ^2 + 8 τ (tbig-2)]]]}
(* r_{(T,1),T} *)rtbigtbig={Simplify[-rbp2a+Simplify[(1/(4 τ))Sqrt[σ^2 + 8 τ tbig]]],Simplify[-rbp2a-Simplify[(1/(4 τ))Sqrt[σ^2 + 8 τ tbig]]]}

(* Decomposition of the integral *)

Simplify[τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]])]
Simplify[tbig*Log[rtbigtbig[[1]]]-(tbig-2)*Log[rtbigtbigm2[[1]]]]
Simplify[τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]]-(τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]])]
factor1[tbig_]=Exp[Simplify[τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]]-(τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]])]]

This last expression is factor1.

(* Computation of the integral -- Error terms *)

(* Error term 1: The original tails of the integrals ([-π,-θ_0] and [θ_0,π]) are exponentially small *)
Series[Cos[2θ]-1,{θ,0,3}]
cos2θ=Series[Cos[2tbig^(-2/5)]-1,{tbig,Infinity,1}]
cos2θm2=Series[Cos[2(tbig-2)^(-2/5)]-1,{tbig,Infinity,1}]
cosθ=Series[Cos[tbig^(-2/5)]-1,{tbig,Infinity,1}]
cosθm2=Series[Cos[(tbig-2)^(-2/5)]-1,{tbig,Infinity,1}]

Clear[linttailstbigm2,linttailstbig,inttailstbigm2,inttailstbig]
linttailstbigm2[tbig_]=Simplify[τ*rtbigtbigm2[[1]]^2*cos2θm2+σ*rtbigtbigm2[[1]]*cosθm2]
linttailstbig[tbig_]=Simplify[τ*rtbigtbig[[1]]^2*cos2θ+σ*rtbigtbig[[1]]*cosθ]
inttailstbigm2=Exp[linttailstbigm2[tbig]]
inttailstbig=Exp[linttailstbig[tbig]]
err1=inttailstbig;err1m2=inttailstbigm2;

Simplify[-τ*rtbigtbig[[1]]^2-σ*rtbigtbig[[1]]]
Simplify[-(((Sqrt[2]+2))/2)σ*rtbigtbig[[1]]]
Simplify[-σ*rtbigtbig[[1]]]

(* Error term 2: The difference bw original and Gaussian central regions is small *)
(* We show that the term in the exponent with order θ^3 is O[tbig^(-1/5)] *)
{Series[Exp[I θ]-1,{θ,0,3}], Series[Exp[I θ]-1,{θ,0,3}]^2, Series[Exp[I θ]-1,{θ,0,3}]^3}
exp=Collect[(n-((σT)/2)r)Series[Exp[I θ]-1,{θ,0,3}]^2-(n/3)Series[Exp[I θ]-1,{θ,0,3}]^3,θ]
{Coefficient[exp,θ,2],Coefficient[exp,θ,3]}
Simplify[{Coefficient[exp,θ,2],Coefficient[exp,θ,3]}/.{n->tbig,r->rtbigtbig[[1]],σT->σ}]
err2=Simplify[Coefficient[exp,θ,3]*n^(-6/5)/.{n->tbig,r->rtbigtbig[[1]],σT->σ}]
Simplify[{Coefficient[exp,θ,2],Coefficient[exp,θ,3]}/.{n->tbig-2,r->rtbigtbigm2[[1]],σT->σ}]
err2m2=Simplify[Coefficient[exp,θ,3]*n^(-6/5)/.{n->tbig-2,r->rtbigtbigm2[[1]],σT->σ}]

(* Error term 3: The superfluous tails of the Gaussian are exponentially small *)
err3=((Sqrt[2])/(Sqrt[tbig-(( σ)/2)*rtbigtbig[[1]]]))Exp[-(tbig-(( σ)/2)*rtbigtbig[[1]])*(tbig^(-2/5))^2]
err3m2=((Sqrt[2])/(Sqrt[tbig-2-(( σ)/2)*rtbigtbigm2[[1]]]))Exp[-(tbig-2-(( σ)/2)*rtbigtbigm2[[1]])*((tbig-2)^(-2/5))^2]

These were the error terms.

(* factor2 and putting it together with factor1 *)
Simplify[tbig-(( σ)/2)*rtbigtbig[[1]]]
Simplify[tbig-2-(( σ)/2)*rtbigtbigm2[[1]]]
fromGaussint[tbig_]=Simplify[Sqrt[Simplify[tbig-(( σ)/2)*rtbigtbig[[1]]]]]
(* this gives the same as above: *)Series[fromGaussint[tbig],{tbig,Infinity,lengthofexpansion}]
fromGaussintm2[tbig_]=Simplify[Sqrt[Simplify[tbig-2-(( σ)/2)*rtbigtbigm2[[1]]]]]
(* this gives the same as above: *)Series[fromGaussintm2[tbig],{tbig,Infinity,lengthofexpansion}]
factor2[tbig_]=Simplify[((fromGaussint[tbig])/(fromGaussintm2[tbig]))]

This last expression is factor2.

(* Taster: the final, crucial calculation without the error terms *)

ftbigtbigm2patbigtbig=factor1[tbig]*factor2[tbig]/.cτ->Sqrt[(km1/(8 k1))]
Simplify[{ku,ku+kb*τ}/.cτ->Sqrt[(km1/(8 k1))]]
Simplify[{ku,ku+kb*τ*ftbigtbigm2patbigtbig}/.cτ->Sqrt[(km1/(8 k1))]]
Simplify[{kb*τ,kb*τ+ku/ftbigtbigm2patbigtbig}/.cτ->Sqrt[(km1/(8 k1))]] (* this is used in paper *)
prelimPgbigst0=Together[(ku/(ku+kb*τ*ftbigtbigm2patbigtbig))]/.cτ->Sqrt[(km1/(8 k1))]
prelimPgbigst1=Together[Simplify[((kb*τ)/(kb*τ+ku/ftbigtbigm2patbigtbig))]]/.cτ->Sqrt[(km1/(8 k1))] (* this is used in paper *)
(* checking P(G*=0)+P(G*=1)=1 *)Simplify[prelimPgbigst0+prelimPgbigst1]

(* The final, crucial calculation repeated with error terms *)

err1
Simplify[err2/I]
err3

(* the integral for f_T,T without the φ(r_T,T)/2π term that is found in factor1:
(((Sqrt[Pi])/(Sqrt[tbig-(( σ)/2)*rtbigtbig[[1]]]))-error3)*(1+Abs[error2])+(2π-2θ_0)*error1  *)

error2=SeriesData[tbig,Infinity,{1},0,1,5]

factor2enumwerrorterms[tbig_]=((Sqrt[Pi])/(Sqrt[Simplify[tbig-2-(( σ)/2)*rtbigtbigm2[[1]]]]))*error2
factor2denomwerrorterms[tbig_]=((Sqrt[Pi])/(Sqrt[Simplify[tbig-(( σ)/2)*rtbigtbig[[1]]]]))*error2

factor2werrorterms[tbig_]=Simplify[((factor2enumwerrorterms[tbig])/(factor2denomwerrorterms[tbig]))]

ftbigtbigm2patbigtbigwerrorterms=factor1[tbig]*factor2werrorterms[tbig]/.cτ->Sqrt[(km1/(8 k1))]
Simplify[{ku,ku+kb*τ*ftbigtbigm2patbigtbigwerrorterms}/.cτ->Sqrt[(km1/(8 k1))]]
Simplify[{kb*τ,kb*τ+ku/ftbigtbigm2patbigtbigwerrorterms}/.cτ->Sqrt[(km1/(8 k1))]]
Pgbigst0werrorterms=Together[(ku/(ku+kb*τ*ftbigtbigm2patbigtbigwerrorterms))]/.cτ->Sqrt[(km1/(8 k1))]
Pgbigst1werrorterms=Together[Simplify[((kb*τ)/(kb*τ+ku/ftbigtbigm2patbigtbigwerrorterms))]]/.cτ->Sqrt[(km1/(8 k1))]

This is the final result.

(* checking P(G*=0)+P(G*=1)=1 *)
Simplify[Pgbigst0werrorterms+Pgbigst1werrorterms]

(* For the peace of mind, here we check that
-the saddle point at r_1 dominates the one at r_2
-it's lower at r_1 than at -r_2 and at 1
-it's lower at r_2 than at -r_1 and at -1
*)

(* φ_TT(r_TT1)/φ_TT(r_TT2) *)Simplify[τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]]-(τ*rtbigtbig[[2]]^2+σ*rtbigtbig[[2]]-tbig*Log[-rtbigtbig[[2]]])]
(* φ_T,T-2(r_T,T-2,1)/φ_T,T-2(r_T,T-2,2) *)Simplify[τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]]-(τ*rtbigtbigm2[[2]]^2+σ*rtbigtbigm2[[2]]-(tbig-2)*Log[-rtbigtbigm2[[2]]])]

(* φ_TT(-r_TT2)/φ_TT(r_TT1) *)Simplify[τ*rtbigtbig[[2]]^2-σ*rtbigtbig[[2]]-tbig*Log[-rtbigtbig[[2]]]-(τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]])]
(* φ_T,T-2(-r_T,T-2,2)/φ_T,T-2(r_T,T-2,1) *)Simplify[τ*rtbigtbigm2[[2]]^2-σ*rtbigtbigm2[[2]]-(tbig-2)*Log[-rtbigtbigm2[[2]]]-(τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]])]
(* φ_TT(1)/φ_TT(r_TT1) *)Simplify[τ+σ-tbig*Log[1]-(τ*rtbigtbig[[1]]^2+σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]])]
(* φ_T,T-2(1)/φ_T,T-2(r_T,T-2,1) *)Simplify[τ+σ-(tbig-2)*Log[1]-(τ*rtbigtbigm2[[1]]^2+σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]])]

(* φ_TT(-r_TT1)/φ_TT(r_TT2) *)Simplify[τ*rtbigtbig[[1]]^2-σ*rtbigtbig[[1]]-tbig*Log[rtbigtbig[[1]]]-(τ*rtbigtbig[[2]]^2+σ*rtbigtbig[[2]]-tbig*Log[-rtbigtbig[[2]]])]
(* φ_T,T-2(-r_T,T-2,1)/φ_T,T-2(r_T,T-2,2) *)Simplify[τ*rtbigtbigm2[[1]]^2-σ*rtbigtbigm2[[1]]-(tbig-2)*Log[rtbigtbigm2[[1]]]-(τ*rtbigtbigm2[[2]]^2+σ*rtbigtbigm2[[2]]-(tbig-2)*Log[-rtbigtbigm2[[2]]])]
(* φ_TT(-1)/φ_TT(r_TT2) *)Simplify[τ-σ-tbig*Log[1]-(τ*rtbigtbig[[2]]^2+σ*rtbigtbig[[2]]-tbig*Log[-rtbigtbig[[2]]])]
(* φ_T,T-2(-1)/φ_T,T-2(r_T,T-2,2) *)Simplify[τ-σ-(tbig-2)*Log[1]-(τ*rtbigtbigm2[[2]]^2+σ*rtbigtbigm2[[2]]-(tbig-2)*Log[-rtbigtbigm2[[2]]])]

τ*rtbigtbigm2[[1]]^2
τ*rtbigtbig[[1]]^2
σ*rtbigtbigm2[[1]]
σ*rtbigtbig[[1]]

rtbigtbig[[2]]-1
Log[rtbigtbig[[2]]]
rtbigtbigm2[[2]]-1
Log[rtbigtbigm2[[2]]]