Why is this geospatial search so slow?

Maybe try this too:

WITH -140.0312186 AS Xcenter, 39.33233141 AS Ycenter, 
pu,
0.5 AS Xslop, // where slop is some bounding box value is equivalent to +/-10000
0.5 AS Yslop // where slop is some bounding box value is equivalent to +/-10000
WHERE 
  pu.location.x < Xcenter) + Xslop AND 
  pu.location.x > Xcenter) - Xslop AND 
  pu.location.y < Ycenter) + Yslop AND 
  pu.location.y > Ycenter) - Yslop AND
  distance(point({srid:4326, x:Xcenter, y:Ycenter}), pu.location) < maxdistance

The four pu.location tests are short circuit booleans, which will avoid the expensive (unindexed) distance calculations if the pu.location is outside the bounding box.
I'm not sure if Neo4j indexes the individual values of location.x and location.y, but you could create indexes for them. With this index, it will go faster.

The Xslop and Yslop may be trickier to calculate as they are in degrees and distance will vary by what location they are at.

See this:

and