What are the formulas to convert discrete-time state space ($A_d,B_d,C_d,D_d$) to continuous-time state space ($A,B,C,D$)? i.e., Convert the state space
$x_{k+1} = A_d x_{k} + B_d u_k$
$y_k = C_d x_{k} + D_d u_k$
To
$\dot{x}(t) = A x(t) + B u(t)$
$y(t) = C x(t) + D u(t)$
x: state vector
y: output vector
u: input (or control) vector
One can try to do the inverse of zero order hold discretization:
\begin{align} \begin{bmatrix}A_d & B_d \\ 0 & I\end{bmatrix} &= e^{\begin{bmatrix}A & B \\ 0 & 0\end{bmatrix} T}, \tag{1a} \\ C_d &= C, \tag{1b} \\ D_d &= D. \tag{1c} \end{align}
So one has to take the matrix logarithm of the matrix containing $A_d$ and $B_d$ from $(1a)$. However, this is not always well defined, especially in the context of state space models of physical systems where $A$ and $B$ should only contain finite real values. For example if $A_d$ has one eigenvalue at some negative real value then $A$ should have one eigenvalue which has only a nonzero imaginary part and not an associated complex conjugate eigenvalue. Another example would be if $A_d$ has eigenvalues at zero, since the logarithm of zero is not well defined. However, one could also interpret those eigenvalues as a delay. Furthermore, a matrix logarithm it not unique since $e^{M+2\,\pi\,i\,k\,I}=e^{M},\ \forall\,k\in\mathbb{Z}$. However, one could add the assumption that the imaginary parts of the eigenvalues of $A$ have the smallest possible absolute value (which can also be interpret as to saying that the discrete system is not "aliasing" the dynamics of the continuous system).
Another popular discretization method is the bilinear transform, a.k.a. the Tustin's method, which I believe gives a better behaved back and forth mapping between the continuous and discrete time representation. Namely, for discrete to continuous one can use the following relation
$$ z = \frac{2 + s\,T}{2 - s\,T}. \tag{2} $$
Substituting $(2)$ in
$$ z\,x = A_d\,x + B_d\,u, \tag{3} $$
and factoring out $s$ yields
$$ s \left((A_d + I)\,x + B\,u\right) = \frac{2}{T} \left((A_d - I)\,x + B\,u\right). \tag{4} $$
Now by defining a new state vector as $x' = (A_d + I)\,x + B\,u$ and solving it for $x$ yields
$$ x = (A_d + I)^{-1} (x' - B\,u). \tag{5} $$
Substituting $(5)$ in $(4)$ and $y = C_d\,x + D_d\,u$ yields
\begin{align} s\,x' &= A\,x' + B\,u, \tag{6a} \\ y &= C\,x' + D\,u, \tag{6b} \end{align}
with
\begin{align} A &= \frac{2}{T} (A_d - I)\,(A_d + I)^{-1}, \tag{7a} \\ B &= \frac{2}{T} \left(I - (A_d - I)\,(A_d + I)^{-1}\right) B_d, \tag{7b} \\ C &= C_d\,(A_d + I)^{-1}, \tag{7c} \\ D &= D_d - C_d\,(A_d + I)^{-1} B_d, \tag{7d} \\ \end{align}
which should always allow one to define a mapping from discrete to continuous if minus one is not an eigenvalue of $A_d$.