Performance of Spatial Intersections

Hello all
I'm evaluating the performance of Neo4J with doing a geospatial intersection between points. The intent is as follows: Take 300,000 Address points that have a street address, latitude, and longitude and see which addresses are within 50m of each other. The below query works fine on very small datasets, however, it is extremely unwieldly when dealing with a dataset of 300,000 points. Due to the cartesian products, this tends to take quite a long time.

Is there something I can do to optimize this further? I can achieve much faster results using GIS tools such as QGIS, however, I'd prefer to keep a work pipeline within Neo4j if the functionality can support it.

CREATE index for (n:Address) on (n.address) ;
CREATE index for (n:Address) on (n.location) ;

MATCH (p:Address)
WHERE p.location is NULL
AND p.latitude is not NULL
SET p.location = Point({latitude:tofloat(p.latitude), longitude:tofloat(p.longitude)})
RETURN count(*);

CALL apoc.periodic.iterate( 
MATCH (p1:Address) 
WHERE exists(p1.location)
MATCH (p2:Address) 
WHERE exists(p2.location) 
AND distance(p1.location, p2.location) < 50
AND <>
RETURN p1.address, p2.address
{batchSize:10000, iterateList:true, parallel:false}  

OK, so modifying:

RETURN p1.address, p2.address


RETURN p1,p2

greatly increases performance. I had inadvertently left in a debugging return statement in the production query statement.