Very interesting. Following Sean's suggestion, here are a few ideas in an attached notebook.
DensityPlot3D[
Log10[density[x, y, z]], {x, y, z} \[Element]
ImplicitRegion[-rx <= x <= rx && -rx <= z <= rx && 0 <= y <= rx, {x,
y, z}], AxesLabel -> {"X", "Y", "Z"},
ColorFunction -> "TemperatureMap"]
Attachments: