your approach really physicalizes arrays from 2D to 4D

No, it does not. It uses clump and slice, both of which are affine transformations (as shown by the P2Child=>1 in https://github.com/PDLPorters/pdl/blob/master/Basic/Slices/slices.pd for the relevant operations). Physicalising only happens on output or use in calculations, and when loop fusion comes in, not even then. I'd love to hear of an approach that Kronecker-products m1,n1 and m2,n2 matrices without producing a m1*m2,n1*n2 in-memory object. While you're pondering that, please show some numbers demonstrating these "grave consequences" :-)

Is this a bug in PDL::NiceSlice, documentation flaw, or just me reading it incorrectly?

You're going to need to show some code that you're saying is behaving unexpectedly. This worked as described (with output):

pdl> p sequence(2)->(*1)
[
 [0]
 [1]
]
pdl> p sequence(2)->(*)
[
 [0]
 [1]
]
pdl> p sequence(2)->(*,*)
[
 [
  [0]
 ]
 [
  [1]
 ]
]
[download]

Naturally, you need to make sure a use PDL::NiceSlice is in scope. You will have noticed that generally I avoid using it in core PDL code, because it both slows down module-loading, and has limitations; see PDL::LinearAlgebra, where I had to remove a lot of it in favour of building arguments to slice that took account of adding optional parts for PDL::Complex.