How do subsequent convolution layers work?
I am not sure about the alternatives described above, but the commonly used methodology is:

Before the application of the non-linearity, each filter output depends linearly on all of the feature maps before within the patch, so you end up with $k_2$ filters after the second layers. The overall number of parameters is $3 \dot{} 3\dot{}k_1 + k_1\dot{} 5 \dot{} 5 \dot{} k_2$.

Bonus 1: Pooling is done per feature map, separately.

Bonus 2: The order of "sliding" does not matter. In fact, each output is computed based on the previous layer, so the output filter responses do not depend on each other. They can be computed in parallel.