The only way I see it possible to do the matrix-multiplication in the MixColumns operation of AES is by shifting the bits in the multiplied number, and then reduce with the polynomial if needed.

This can be done in constant time for a bitsliced implementation as shown below in C-like pseudo-code (the example is for multiplication by $3$ in $GF(2^8)$):

vector bit[8] = "bits from bitsliced bytes of plaintext";
vector new_bit[8];

//First multiply with 2 in GF(2^8) with polynomial x^8+x^4+x^3+x+1:
new_bit[0] = {0} ^ bit[7];
new_bit[1] = bit[0] ^ bit[7];
new_bit[2] = bit[1];
new_bit[3] = bit[2] ^ bit[7];
new_bit[4] = bit[3] ^ bit[7];
new_bit[5] = bit[4];
new_bit[6] = bit[5];
new_bit[7] = bit[6];

//XOR M2 with original bits to get the result for multiplying with 3
for (int bit = 0; bit < 8; bit++){
 new_bit[i] = new_bit[i] ^ bit[i];
}

However, this does not seem to be the way that Käsper and Schwabe do it in their bit-sliced AES MixColumns operation. Instead they are doing rotates. I have read their paper and looked at the layout of their bit-sliced state, but I'm unable to understand, how they arrive at those equations. If someone understands it, could they perhaps help a poor soul out there to understand it too? Just sending me in the direction for a paper that explains it more in depth would mean the world to me as well!

I'm currently trying to implement a bit-sliced mix columns implementation for a cipher other than AES (I'm trying to bit-slice the primates cipher: http://primates.ae/wp-content/uploads/primatesv1.02.pdf), and have currently made the "simple" MixColumns that I described. I'm curious in whether it can be done better, after having seen that they do it differently in their AES implementation. There's not the same symmetry in this MixColumns matrice, so I doubt a bit that it is possible.