(*23 May 2026*)
(*Quantum statistical mechanical gauge invariance*)
(*Johanna Müller and Matthias Schmidt*)
(*arXiv:2509.20494*)

(*The first list below contains the following entries.
{1: normalization, 2: density profile n, 3: n beta H,
 4: beta H, 5: F_id beta H, 6: F_id}
The entry 1 is the normalization that the remaining entries still need.
The second list “res” below contains the following results, expressed in units of length a.
{1: n, 2: covariance(beta F_id|beta H), 3: covariance(beta F_ext|beta H), 4: x beta d n/d beta}*)

repl={a->Sqrt[(hbar/(m omega))],hbar->1,m->1,omega->1};

((HermiteH[n,x/a])/(Sqrt[(2^n)n!a Sqrt[Pi]]))Exp[-((x^2)/(2(a^2)))];
Sum[{1,((%)^2),((%)^2)beta en,beta en,(-((hbar^2)/m)D[((D[%,x])^2),x]+((hbar^2)/(4m))D[((%)^2),{x,3}])beta en,(-((hbar^2)/m)D[((D[%,x])^2),x]+((hbar^2)/(4m))D[((%)^2),{x,3}])}Exp[-beta  en]/.en->hbar omega((1/2)+n),{n,0,25}];
%/.a->1/.repl;
res=beta{((%[[2]])/(%[[1]]))/beta, ((%[[5]])/(%[[1]]))-((%[[6]])/(%[[1]]))((%[[4]])/(%[[1]])),-x (((%[[3]])/(%[[1]]))-((%[[2]])/(%[[1]]))((%[[4]])/(%[[1]]))),x  beta D[((%[[2]])/(%[[1]])),beta]};
Plot[res[[2(*1,2,3,4*)]]/.{{beta->2},{beta->1},{beta->.5},{beta->4}}//Evaluate,{x,-5,5},PlotRange->All]

(*density profile*)
Table[{beta->bb},{bb,.5,6,.5}];
Plot[res[[1]]/.%//Evaluate,{x,-4,4},PlotStyle->Evaluate[{RGBColor[1-(beta-1)/6.,(beta-1)/6,1-0beta/5.]}/.%],PlotRange->All]

(*covariance of Hamiltonian and ideal force density*)
Table[{beta->bb},{bb,.5,6,.5}];
Plot[res[[2]]/.%//Evaluate,{x,-4,4},PlotStyle->Evaluate[{RGBColor[1-(beta-1)/6.,(beta-1)/6,1-0beta/5.]}/.%],PlotRange->All]

(*covariance of Hamiltonian and external force density*)
Table[{beta->bb},{bb,.5,6,.5}];
Plot[res[[3]]/.%//Evaluate,{x,-4,4},PlotStyle->Evaluate[{RGBColor[1-(beta-1)/6.,(beta-1)/6,1-0beta/5.]}/.%],PlotRange->All]

(*covariance of Hamiltonian and the external force density, alternative route as a consistency check*)
Table[{beta->bb},{bb,.5,6,.5}];
Plot[res[[4]]/.%//Evaluate,{x,-4,4},PlotStyle->Evaluate[{RGBColor[1-(beta-1)/6.,(beta-1)/6,1-0beta/5.]}/.%],PlotRange->All]

(*version with physical units as a check for correct scaling*)
((HermiteH[n,x/a])/(Sqrt[(2^n)n!a Sqrt[Pi]]))Exp[-((x^2)/(2(a^2)))];
Sum[{1,((%)^2),((%)^2)beta en,beta en,(-((hbar^2)/m)D[((D[%,x])^2),x]+((hbar^2)/(4m))D[((%)^2),{x,3}])beta en,(-((hbar^2)/m)D[((D[%,x])^2),x]+((hbar^2)/(4m))D[((%)^2),{x,3}])}Exp[-beta  en]/.en->hbar omega((1/2)+n),{n,0,1(*truncate upper limit to keep expressions lean for further symbolic manipulations*)}];
beta{((%[[2]])/(%[[1]]))/beta, ((%[[5]])/(%[[1]]))-((%[[6]])/(%[[1]]))((%[[4]])/(%[[1]])),-m (omega^2)x(((%[[3]])/(%[[1]]))-((%[[2]])/(%[[1]]))((%[[4]])/(%[[1]]))),m (omega^2) x beta D[((%[[2]])/(%[[1]])),beta]};
{%[[1]]a,%[[2]] (a^2),%[[3]](a^2),%[[4]](a^2)(*scaling to make quantities dimensionless*)}//FullSimplify ;
%/.x->xs a/.a->Sqrt[(hbar/(m omega))]/.beta->bho /(hbar omega)//FullSimplify
(*The symbolic results below indeed have no units and depend merely on the scaled lenth xs = x/a and on the energy ratio bho = beta hbar omega. Hence the scaling is correct.*)
{%[[2]]+%[[3]]==0(*sum rule*),%[[3]]==%[[4]](*consistency check*)}//Simplify