Use strict inequality for search radius check?

We currently only check if distance2 <= search_radius^2.
However, this causes problems in some computations. For example, this here in the surface tension computation:

julia> support_radius = 0.02; distance = support_radius; (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25
ERROR: DomainError with -6.938893903907228e-18:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
     ⋮ internal @ Base.Math
 [2] ^(x::Float64, y::Float64)
   @ Base.Math ./math.jl:1206
Use `err` to retrieve the full stack trace.

This does not happen when distance < search_radius.
The question is if there are any downsides to this strict inequality check. I don't see any.

Note that distance2 < search_radius^2 is not sufficient:

julia> 0.0004 * (1 - 0.7eps()) < 0.02^2
true

julia> sqrt(0.0004 * (1 - 0.7eps())) < 0.02
false

So we would have to do something like

distance2 < search_radius^2 * (1 - eps())

Start by locating the neighborhood-search code that checks whether distance2 is within search_radius^2, then reproduce the Julia floating-point examples from the issue. Done means resolving whether a strict or tolerance-based boundary check is appropriate and assessing its downsides; no specific file or test is named.