computer code from R. M. Feenstra, “Path distributions for describing eigenstates of the harmonic oscillator and other 1-dimensional problems”,  arXiv:2306.11155 [quant-ph].

computations for Fig. 1(a) and Fig. 2

examine values of R^1 propagator times wavefunction, Eq. (14), as plotted in Fig. 2(a)

Clear["Global`*"];
k=1; T=10000;np=200000; pmin=-6;pmax=4;delp=(pmax-pmin)/np;
Plot[{Im[Exp[I(((p-k))^2)T/2 ]Sqrt[T/(2π I)]],Re[Exp[I(((p-k))^2)T/2 ]Sqrt[T/(2π I)]]},{p,0.7,1.3},PlotRange->{{.7,1.3},{-70,70}},AspectRatio->1/5,PlotStyle->{{Thickness[0.008],ColorData[97,"ColorList"][[2]]},{Thickness[0.008],ColorData[97,"ColorList"][[1]]}},AxesOrigin->{0.7,-70}]

examine phasor values (integral of R^1 propagator times wavefunction), Eq. (13), as plotted in Fig. 1(a)

Clear["Global`*"];
k=1; T=10000;np=200000; pmin=-6;pmax=4;delp=(pmax-pmin)/np;
tab=N[delp (((T/(2π I)))^(1/2))]Table[
N[Exp[I (((p-k))^2) T/2]],{p,pmin,pmax,delp}];
list1=Accumulate[Re[tab]];
list2=Accumulate[Im[tab]];
listplot=Transpose[{list1,list2}];
ListLinePlot[listplot,Joined->True,AspectRatio->1/2,PlotRange->{{-.5,1.5},{-.5,.5}}]

obtain analytic form for integral above (not actually used for plotting)

Clear["Global`*"];
Integrate[(((T/(2π I)))^(1/2))Exp[I (((p-k))^2) T/2],{p,pc-Δp,pc+Δp}]

examine the *difference* between phasor values, Eq. (15), as plotted in Fig. 2(b)

nskip=Round[2 /(Sqrt[T]delp)];
Print[nskip];
list1p=Sqrt[T]Drop[list1,nskip]/2;
list1n=Sqrt[T]Drop[list1,-nskip]/2;
list1d=list1p-list1n;
list2p=Sqrt[T]Drop[list2,nskip]/2;
list2n=Sqrt[T]Drop[list2,-nskip]/2;
list2d=list2p-list2n;
ListLinePlot[{Table[{pmin+(i+nskip/2) delp,list2d[[i]]},{i,Length[list1d]}],Table[{pmin+(i+nskip/2) delp,list1d[[i]]},{i,Length[list1d]}]},PlotRange->{{.7,1.3},{-25,35}},AspectRatio->1/2,PlotStyle->{{Thickness[0.008],ColorData[97,"ColorList"][[2]]},{Thickness[0.008],ColorData[97,"ColorList"][[1]]}},AxesOrigin->{0.7,-25}]