(* Mathematica code used to compute concurrence in the presence of white noise in "Entanglement generation and scaling from noisy quenches across a quantum critical point" by R. Jafari et al. For questions, please contact R. Jafari at raadmehr.jafari@gmail.com *) (*$PreRead=(#/.s_String/;StringMatchQ[s,NumberString]&&Precision@ToExpression@s\[Equal]MachinePrecision\[RuleDelayed]s<>"`50."&);*) SetDirectory[NotebookDirectory[]]; Clear["Global`*"]; L=200; \[Gamma]=1.; J=1.; h1=-30.; h2=30; r=1.; lstA={0.001,0.002,0.003,0.004}; Do[A=lstA[[z]]; Concfile2=OpenWrite["Numerical_NNConcurrence_VS_Inverse_V,gamma="<>ToString[N[\[Gamma]]]<>",xi="<>ToString[N[A]]<>",h1="<>ToString[N[h1]]<>",h2="<>ToString[N[h2]]<>",N="<>ToString[N[L]]<>".txt",FormatType->ScientificForm]; lstV={100.0,80.0,40.0,20.0,18.0,16.0,14.0,12.0,10.0,8.0,6.0,4.0,2.0,1.0,1.0/2,1.0/3,1.0/4,1.0/5,1.0/6,1.0/7,1.0/8,1.0/9,1.0/10,1.0/11,1.0/12.0,1.0/13.0,1.0/14.0,1.0/15.0,1.0/16.0,1.0/18.0,1.0/20.0,1.0/22.0,1.0/24.0,1.0/26.0,1.0/28.0,1.0/30.0,1.0/32.0,1.0/34.0,1.0/36.0,1.0/38.0,1.0/40.0,1.0/42.0,1.0/44.0,1.0/46.0,1.0/48.0,1.0/50.0,1.0/52.0,1.0/54.0,1.0/56.0,1.0/58.0,1.0/60.0,1.0/62.0,1.0/64.0,1.0/66.0,1.0/68.0,1.0/70.0,1.0/72.0,1.0/74.0,1.0/76.0,1.0/78.0,1.0/80.0,1.0/82.0,1.0/84.0,1.0/86.0,1.0/88.0,1.0/90.0,1.0/92.0,1.0/94.0,1.0/96.0,1.0/98.0,1.0/100.0,1.0/110.0,1.0/120.0,1.0/130.0,1.0/140.0,1.0/150.0,1.0/160.0,1.0/170.0,1.0/180.0,1.0/190.0,1.0/200.0,1.0/210.0,1.0/220.0,1.0/230.0,1.0/240.0,1.0/250.0,1.0/260.0,1.0/270.0,1.0/280.0,1.0/290.0,1.0/300.0}; Do[v=lstV[[u]]; Table[f[i],{i,1,L/2,1}]; Table[g[x],{x,1,L/2,1}]; For[m=1,m<=L/2,m=m+1, k=N[((2m-1)*\[Pi])/L]; hz=(h2-Cos[k]); hy=Sin[k]; R={{(Sqrt[((hz)^2+(hy)^2)]-hz)/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)],-I*hy/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)]},{-I*hy/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)],(Sqrt[((hz)^2+(hy)^2)]-hz)/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)]}}; Rinv={{(Sqrt[((hz)^2+(hy)^2)]-hz)/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)],I*hy/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)]},{I*hy/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)],(Sqrt[((hz)^2+(hy)^2)]-hz)/Sqrt[2*Sqrt[((hz)^2+(hy)^2)](Sqrt[((hz)^2+(hy)^2)]-hz)]}}; Clear[\[Rho]11]; Clear[\[Rho]12]; Clear[\[CapitalDelta]11]; Clear[\[CapitalDelta]12]; {Sol1,Sol2}=NDSolveValue[{v*(\[Rho]11^\[Prime])[t]==-2*(Sin[k]*Re[\[Rho]12[t]]),v*(\[Rho]12^\[Prime])[t]==-Sin[k]+2 Sin[k] \[Rho]11[t]+2 I (I A^2-t+Cos[k]) \[Rho]12[t],\[Rho]11[h1]==1,\[Rho]12[h1]==0},{\[Rho]11,\[Rho]12},{t,h1,h2}]; \[Rho]11k=Chop[Sol1[h2]]; \[Rho]12k=Chop[Sol2[h2]]; \[Rho]21k=Conjugate[\[Rho]12k]; \[Rho]22k=1-\[Rho]11k; \[Rho]k=(\[Rho]11k \[Rho]12k \[Rho]21k \[Rho]22k ); rhoa=Rinv.\[Rho]k.R; f[m]=Chop[rhoa[[2,2]]]; g[m]=Chop[rhoa[[1,2]]]; ]; AA[s_,p_]=(2*I)/L*Sum[Re[g[q]]*Sin[((2q-1)*\[Pi])/L*(p-s)],{q,1,L/2}]+KroneckerDelta[s,p]; AB[s_,p_]=1/L*Sum[2*Chop[1-2*f[q]]*(Cos[((2q-1)*\[Pi])/L*(p-s)]*(h2-J*Cos[((2q-1)*\[Pi])/L])/Sqrt[(h2-J*Cos[((2q-1)*\[Pi])/L])^2+(\[Gamma]*Sin[((2q-1)*\[Pi])/L])^2]-Sin[((2q-1)*\[Pi])/L*(p-s)]*(\[Gamma]*Sin[((2q-1)*\[Pi])/L])/Sqrt[(h2-J*Cos[((2q-1)*\[Pi])/L])^2+(\[Gamma]*Sin[((2q-1)*\[Pi])/L])^2])- (4*Im[g[q]]*(Sin[((2q-1)*\[Pi])/L*(p-s)]*(h2-J*Cos[((2q-1)*\[Pi])/L])/Sqrt[(h2-J*Cos[((2q-1)*\[Pi])/L])^2+(\[Gamma]*Sin[((2q-1)*\[Pi])/L])^2]+Cos[((2q-1)*\[Pi])/L*(p-s)]*(\[Gamma]*Sin[((2q-1)*\[Pi])/L])/Sqrt[(h2-J*Cos[((2q-1)*\[Pi])/L])^2+(\[Gamma]*Sin[((2q-1)*\[Pi])/L])^2])),{q,1,L/2}]; BB[s_,p_]=(2*I)/L*Sum[Re[g[q]]*Sin[((2q-1)*\[Pi])/L*(p-s)],{q,1,L/2}]-KroneckerDelta[p,s]; Sxx=Chop[-(1/4)*AB[r+1,r]]; Syy=Chop[-(1/4)*AB[r,r+1]]; Szz=Chop[1/4*(AB[r,r]*AB[r+1,r+1]-AA[r,r+1]*BB[r,r+1]-AB[r+1,r]*AB[r,r+1])]; Mz=Chop[-(1/2)*AB[r,r]]; Sxy=Chop[-(I/4)*BB[r,r+1]]; Syx=Chop[-(I/4)*AA[r,r+1]]; \[Rho]11=Chop[Mz+Szz+1/4]; \[Rho]11c=Conjugate[\[Rho]11]; \[Rho]22=\[Rho]33=Chop[-Szz+1/4]; \[Rho]22c=\[Rho]33c=Conjugate[\[Rho]22]; \[Rho]44=Chop[-Mz+Szz+1/4]; \[Rho]44c=Conjugate[\[Rho]44]; \[Rho]23=Chop[Sxx+Syy+I*(Sxy-Syx)]; \[Rho]23c=Conjugate[\[Rho]23]; \[Rho]14=Chop[Sxx-Syy-I*(Sxy+Syx)]; \[Rho]14c=Conjugate[\[Rho]14]; \[Lambda]a=Chop[1/2 (2 \[Rho]23 \[Rho]23c+\[Rho]22c \[Rho]33+\[Rho]22 \[Rho]33c-\[Sqrt](4 \[Rho]22 \[Rho]23 \[Rho]23c \[Rho]33+4 \[Rho]22c \[Rho]23 \[Rho]23c \[Rho]33+\[Rho]22c^2 \[Rho]33^2+4 \[Rho]22 \[Rho]23 \[Rho]23c \[Rho]33c+4 \[Rho]22c \[Rho]23 \[Rho]23c \[Rho]33c-2 \[Rho]22 \[Rho]22c \[Rho]33 \[Rho]33c+\[Rho]22^2 \[Rho]33c^2))]; \[Lambda]b=Chop[1/2 (2 \[Rho]23 \[Rho]23c+\[Rho]22c \[Rho]33+\[Rho]22 \[Rho]33c+\[Sqrt](4 \[Rho]22 \[Rho]23 \[Rho]23c \[Rho]33+4 \[Rho]22c \[Rho]23 \[Rho]23c \[Rho]33+\[Rho]22c^2 \[Rho]33^2+4 \[Rho]22 \[Rho]23 \[Rho]23c \[Rho]33c+4 \[Rho]22c \[Rho]23 \[Rho]23c \[Rho]33c-2 \[Rho]22 \[Rho]22c \[Rho]33 \[Rho]33c+\[Rho]22^2 \[Rho]33c^2))]; \[Lambda]c=Chop[1/2 (2 \[Rho]14 \[Rho]14c+\[Rho]11c \[Rho]44+\[Rho]11 \[Rho]44c-\[Sqrt](4 \[Rho]11 \[Rho]14 \[Rho]14c \[Rho]44+4 \[Rho]11c \[Rho]14 \[Rho]14c \[Rho]44+\[Rho]11c^2 \[Rho]44^2+4 \[Rho]11 \[Rho]14 \[Rho]14c \[Rho]44c+4 \[Rho]11c \[Rho]14 \[Rho]14c \[Rho]44c-2 \[Rho]11 \[Rho]11c \[Rho]44 \[Rho]44c+\[Rho]11^2 \[Rho]44c^2))]; \[Lambda]d=Chop[1/2 (2 \[Rho]14 \[Rho]14c+\[Rho]11c \[Rho]44+\[Rho]11 \[Rho]44c+\[Sqrt](4 \[Rho]11 \[Rho]14 \[Rho]14c \[Rho]44+4 \[Rho]11c \[Rho]14 \[Rho]14c \[Rho]44+\[Rho]11c^2 \[Rho]44^2+4 \[Rho]11 \[Rho]14 \[Rho]14c \[Rho]44c+4 \[Rho]11c \[Rho]14 \[Rho]14c \[Rho]44c-2 \[Rho]11 \[Rho]11c \[Rho]44 \[Rho]44c+\[Rho]11^2 \[Rho]44c^2))]; \[Lambda]max=Max[\[Lambda]a,\[Lambda]b,\[Lambda]c,\[Lambda]d]; Concurrence1=Max[0,2*Sqrt[\[Lambda]max]-Sqrt[\[Lambda]a]-Sqrt[\[Lambda]b]-Sqrt[\[Lambda]c]-Sqrt[\[Lambda]d]]; (*\[Lambda]a=Chop[Sqrt[\[Rho]11*\[Rho]44]+Abs[\[Rho]14],10^-6]; \[Lambda]b=Chop[Sqrt[\[Rho]22*\[Rho]33]+Abs[\[Rho]23],10^-6]; \[Lambda]c=Chop[Abs[Sqrt[\[Rho]11*\[Rho]44]-Abs[\[Rho]14]],10^-6]; \[Lambda]d=Chop[Abs[Sqrt[\[Rho]22*\[Rho]33]-Abs[\[Rho]23]],10^-6]; Concurrence1=Max[0,\[Lambda]a-\[Lambda]b-\[Lambda]c-\[Lambda]d];*) WriteString[Concfile2,FortranForm@N[1/v]," ",FortranForm@N[Concurrence1],"\n"]; ,{u,1,Length[lstV]}]; Close[Concfile2]; ,{z,1,Length[lstA]}];