Symbolic conversion of transfer function to state-space model using Sympy looks incorrect

Question

I'm trying out the new Sympy control module and a bit puzzled by the output of this symbolic conversion of a Laplace transfer function to a continuous-time state-space model:

import sympy
assert sympy.__version__ >= '1.13.0'
from sympy.physics.control.lti import StateSpace, TransferFunction

s, K, T1 = sympy.symbols('s, K, T1')
G = TransferFunction(K, T1 * s + 1, s)
print(G)
Gss = G.rewrite(StateSpace)
print(Gss)

Output:

TransferFunction(K, T1*s + 1, s)
StateSpace(Matrix([[-1/T1]]), Matrix([[1]]), Matrix([[K]]), Matrix([[0]]))

When I simulated this I got a steady-state gain 3 times higher than K. So I would have expected for example, the C matrix to be Matrix([[K/T1]]) not Matrix([[K]]), or some similar adjustment to the B matrix.

If I plug some numbers into Python Control I get a result consistent with my expectation:

import control as con

K = 2
T1 = 3
G = con.tf(K, [T1, 1])
print(G)
Gss = con.ss(G)
print(Gss)

Output:

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


   2
-------
3 s + 1

<StateSpace>: sys[14]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
States (1): ['x[0]']

A = [[-0.33333333]]

B = [[1.]]

C = [[0.66666667]]

D = [[0.]]