Representing a continuous LTI system as a discrete one I am aware that there are different ways to represent a continuous time system in discrete domain (e.g. bilinear transform, impulse invariance transform).

But my problem is as follows: Given an arbitrary LTI system $h_c(t)$, its output is given as $y(t) = (x * h_c)(t)$.

Now suppose $x(t)$ is given as:

$$
x(t) = \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc}\left(\frac{t-n T}{T}\right)
$$

i.e., Whittaker interpolation formula. I can then write a discrete-time system:

$$
y[n] = \sum_{q=-\infty}^{\infty} h[q] x[n-q] ,
$$

such that $y[n]$ is identical to $y(n T)$.

My question: What is the relation between $h[n]$ and $h_c(t)$?

There must be an exact one and it is neither given by impulse invariance nor bilinear transform.
I am aware that when $H_c(f)$ is bandlimited to $1/(2B)$, the answer is just impulse invariance transform, i.e., $h[n] = h_c(nT)$.
If this is not satisfied, it will result in aliasing but this should not prevent such a relationship from existing.
The sequence $h[n]$ consists of samples of a low pass filtered version of $h_c(t)$. Define that low pass filtered impulse response as

$$\tilde{h}_c(t)=\int_{-\infty}^{\infty}h_c(\tau)\frac{\sin(\pi(t-\tau)/T)}{\pi(t-\tau)/T}d\tau\tag{1}$$

and

$$h[n]=\tilde{h}_c(nT)\tag{2}$$

The output signal is

$$\begin{align}y(t)&amp;=(x\star h_c)(t)\\&amp;=\int_{-\infty}^{\infty}x(t-\tau)h_c(\tau)d\tau\\&amp;=\int_{-\infty}^{\infty}\sum_{k=-\infty}^{\infty}x[k]\frac{\sin(\pi(t-kT-\tau)/T)}{\pi(t-kT-\tau)/T}h_c(\tau)d\tau\\&amp;=\sum_{k=-\infty}^{\infty}x[k]\int_{-\infty}^{\infty}h_c(\tau)\frac{\sin(\pi(t-kT-\tau)/T)}{\pi(t-kT-\tau)/T}d\tau\\&amp;=\sum_{k=-\infty}^{\infty}x[k]\tilde{h}_c(t-kT)\tag{3}\end{align}$$

Consequently, we have

$$\begin{align}y(nT)&amp;=\sum_{k=-\infty}^{\infty}x[k]\tilde{h}_c((n-k)T)\\&amp;=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\end{align}\tag{4}$$

This results is also obvious from the fact that $(x\star h_c)(t)=(x\star\tilde{h}_c)(t)$ must hold, because $x(t)$ is band-limited. And since $\tilde{h}_c(t)$ is also band-limited, we can simply take samples from $x(t)$ and $\tilde{h}_c(t)$ to obtain the result.