As currently written, the new systems that are being created are non-minimal. For example:

s = ct.tf('s')
ct.config.defaults['xferfcn.display_format'] = 'zpk'  # simpler output format
siso = (s+1) / ((s+2) * (s+3))
mimo = ct.combine_tf([
  [1/s,         1/(s+4)], 
  [1/(s+5), (s+6)/(s**2 + 2*s + 2)]
])
print(siso * mimo)

<TransferFunction>: sys[219]
Inputs (2): ['u[0]', 'u[1]']
Outputs (2): ['y[0]', 'y[1]']

Input 1 to output 1:
      (s + 1) (s + 5)
---------------------------
(s) (s + 2) (s + 3) (s + 5)

Input 1 to output 2:
        (s) (s + 1)
---------------------------
(s) (s + 2) (s + 3) (s + 5)

Input 2 to output 1:
        (s + 1) (s + (1-1j)) (s + (1+1j))
-------------------------------------------------
(s + (1-1j)) (s + (1+1j)) (s + 2) (s + 3) (s + 4)

Input 2 to output 2:
             (s + 1) (s + 4) (s + 6)
-------------------------------------------------
(s + (1-1j)) (s + (1+1j)) (s + 2) (s + 3) (s + 4)

Note the various pole/zero cancellations that appear. It seems to me that these are spurious and shouldn't show up (eg, (s+4) in the [2, 2] transfer function).

This same type of situation occurs in state space and is consistent with the results that occur if you convert the SISO system to a matrix system by hand. I'm not sure (yet) what is gong on, but this seems like something we should fix, either for this specific case or for the more general case of multiplying MIMO systems.

I need to think a bit more about what is causing this, but wanted to flag it here in case someone else has some insights. It seems to me that perfect pole/zero cancellations (and the equivalent in state space) should get sorted out ahead of time, either by not creating them or by removing them before we send things back to the user.