Stepdis[d_]:=Join[Table[UnitVector[d,k],{k,d}],Table[-UnitVector[d,k],{k,d}]];
T[N_]:=Flatten[Table[{i,j,k},{i,0,N-1},{j,0,N-1},{k,0,N-1}],2]

Step distribution of random walk and 3-dimensional torus of size length N.

TraceRW[x_,p_,N_]:=Mod[Accumulate[Join[{x},RandomChoice[Stepdis[3],p]]],N];

Computes the trace associated to a random walk started at x on the torus T after p step.

palpha[alpha_,N_]:=IntegerPart[alpha*1.516386059151978*Log[N^3]*N^3];

Computes the number of iterations to be done at level alpha in a box of size N.

In order to illustrate the role of alpha_*, let us now compute the vacant set at alpha=0.5, 0.55... (recall alpha_*=0.67), and save the plot so that the points which have at least 1 neighbor are highlighted in red.

Nu=400;x={0,0,0};step=0.05;goal=0.45;alpha=0;RW={};Vu=T[Nu];K=1;
Do[alpha=alpha+step;Print[alpha];
RW=TraceRW[x,palpha[step,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
,{k,1,goal/step}]
step=0.05;goal=0.45;
Do[alpha=alpha+step;Print[alpha];
RW=TraceRW[x,palpha[step,Nu],Nu];
Vu=Complement[Vu,RW];
x=RW[[Nu]];
If[Length[Vu]<2,MindisV=Infinity,ne=Nearest[Vu];
Distance=MapThread[ManhattanDistance,{Vu,ne[Vu,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Vu[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Vu,Neighbor];
MindisV=Min[Distance];
Export[StringJoin["Vacant_set_torus_size=",TextString[Nu],"alpha=",TextString[alpha],"red_distance_smaller",TextString[K],".jpg"],
Labeled[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotStyle->{Directive[{Blue,PointSize[0.007]}],Directive[{Red,PointSize[0.007]}]}],Style[StringJoin["alpha=",TextString[alpha],", closest=",TextString[MindisV]],30]]];
Print[Labeled[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotStyle->{Directive[{Blue,PointSize[0.007]}],Directive[{Red,PointSize[0.007]}]}],Style[StringJoin["alpha=",TextString[alpha],", closest=",TextString[MindisV]],30]]]
,{k,1,goal/step}]

For comparison, let us now do the same simulation but for independent Bernoulli random variables

Bernoulli[N_,p_]:=RandomInteger[{1,N},{RandomVariate[BinomialDistribution[N^3,p]],3}]
Nu=400;step=0.05;goal=0.45;alpha=0.45;K=1;
Do[alpha=alpha+step;Print[alpha];
Balpha=Bernoulli[Nu,Nu^(-3*alpha)];
K=1;If[Length[Balpha]<2,MindisV=Infinity,ne=Nearest[Balpha];
Distance=MapThread[ManhattanDistance,{Balpha,ne[Balpha,2][[All,2]]}]];
Neighbor={};
Do[If[Distance[[k]]<=K,Neighbor=Join[Neighbor,{Balpha[[k]]}]]
,{k,1,Length[Distance]}];
NoNeighbor=Complement[Balpha,Neighbor];
MindisV=Min[Distance];
Export[StringJoin["Bernoulli_torus_size=",TextString[Nu],"alpha=",TextString[alpha],"red_distance_smaller",TextString[K],".jpg"],Labeled[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotStyle->{Directive[{Blue,PointSize[0.007]}],Directive[{Red,PointSize[0.007]}]}],Style[StringJoin["alpha=",TextString[alpha],", closest=",TextString[MindisV]],30]]];
Print[Labeled[ListPointPlot3D[{NoNeighbor,Neighbor},ImageSize->1000,BoxRatios->{1,1,1},AxesStyle->Opacity[0],TicksStyle->Opacity[0],PlotStyle->{Directive[{Blue,PointSize[0.007]}],Directive[{Red,PointSize[0.007]}]}],Style[StringJoin["alpha=",TextString[alpha],", closest=",TextString[MindisV]],30]]]
,{k,1,goal/step}]