Here are some results for running the BE problem past sink formation:

From the BP collapse (quarter of the grid):

On the left side is a plot of a zoomed in section of the BE sphere (sink formation is at frame = 96). You can see the density compression wave converging on the center. In the right panel are the sinks. At first one forms, then a next, and a next, to a multitude. This is not what we'd expect. What we should expect is exact symmetry: at the center of the collapsing sphere a sink should form, stay stationary, and accrete the surrounding envelope. My assumption is that this is not happening because we are modeling a quadrant of the sphere— to be at the center of the sphere, a sink would have to form along the vertex of the center-most cell.

Except I realized this wasn't a great way to see what was happening… Here is a better view —

Checking the sink file data, I see that the first sink formed at r~(0.004, 0.004, 0.004). The smallest dx=0.0075, so this is approximately at the center of the center smallest cube on the grid. Ideally the sink would form at the origin, corresponding with the peak in phi. The sink then begins to migrate to the center… but then off the domain. I checked this in the sink data file, and indeed the sink acquired a position with negative coords. At frame 103 or so, a second sink was at position r=(0.001) and the first was at r=(-0.02, -0.02, -0.02).

The initial sink being attracted to the origin makes sense. It does not feel gravitational attraction from material at r>rsink by Gauss's law, and so is attracted to the center of the sphere as it is the location of the deepest part of the potential well. However, how is the sink position mapped to the domain, such that it can be off the grid? Does the sink always form at the center of a cell, and then can travel around subsequently non-restricted to cell centers? If so, why can't it form at any given r, not just (n-½)dx? And lastly, if we want to calculate the change in rho® over time after sink formation, what is best way to model the problem? Fixed grid, entire box, odd number of cells would insure a sink to form in the center (and remain in center), but the effective resolution would be extremely computationally heavy. Should this be a 1D calculation?

From the Matched collapse:

Here is early (pre sink) zoomed out evolution:

Here is late zoomed in (sink and post sink formation):