Python 3, 83 bytes

f=lambda b,k=1,p=1:(k<b)|p%k|math.gcd(b**3-b,k)-1and f(b,k+1,p*k*k)or k
import math

Try it online!

Uses xnor's corollary to Wilson's theorem. I still find it unbelievable that the shortest way to test for coprimality is import math;math.gcd(x,y)<2, but I couldn't find a shorter way.

Python 3, 81 bytes

f=lambda b,k=1,p=1:(k<b)|p%k|math.gcd(b**3-b,k)-1and f(b,k+1,p*k)or k
import math

Try it online!

Uses xnor's corollary to Wilson's theorem. I still find it unbelievable that the shortest way to test for coprimality is import math;math.gcd(x,y)<2, but I couldn't find a shorter way.

-2 due to xnor, since we only need the regular factorial here (instead of factorial squared). The special case k=4 is taken care of by the other conditions.

Python 3, 83 bytes

f=lambda b,k=1,p=1:(k<b)|p%k|math.gcd(b**3-b,k)-1and f(b,k+1,p*k*k)or k
import math

Try it online!

Python 3, 83 bytes

f=lambda b,k=1,p=1:(k<b)|p%k|math.gcd(b**3-b,k)-1and f(b,k+1,p*k*k)or k
import math

Try it online!

Uses xnor's corollary to Wilson's theorem. I still find it unbelievable that the shortest way to test for coprimality is import math;math.gcd(x,y)<2, but I couldn't find a shorter way.