I am interested in this question and since no-one has answered yet I will offer my further thoughts (up for debate or confirmation):
I agree with the OP's perspective that the Discrete Fourier Transform (DFT) fundamentally represents when cast to domains that extend to $\pm \infty$ periodic sequences in time (and in frequency) and therefore is fundamentally non-causal due to the time periodicity specifically.
However, to utilize these transforms with context of causality, we can consider the cases that best match an approximation of the Continuous Time Fourier Transform for the same causal representation of our sampled system. Such cases are the ones with no time aliasing, which would occur when $x[n]=0$ for $n\ge N/2$ given $N$ samples with $n \in [0, N-1]$. When we can restrict ourselves to causal sequences specifically (meaning when we know the underlying process is causal) then this restriction can be extended to be $x[n]=0$ for $n\ge N$ given $N$ samples with $n \in [0, N-1]$. This means basically that we know the time duration of the data capture exceeds the expected response time of the underlying continuous time process.
Further in frequency for such an equivalently causal time domain waveform the phase would be going increasingly negative versus frequency indicating the delay that is indeed properly being modeled.
It is the counter case, when for a known causal sequence $x[n]$ has significant non-zero values in the upper-half of the time domain samples right up to the $N-1$ boundary that we cannot be reasonably assured we are not observing the effects of time domain aliasing. We are no longer able to uniquely distinguish between causal and non-causal sequences, or more specifically and practically when we know the sequence is causal; we are no longer able to distinguish from underlying continuous-time responses that went beyond the sample duration of the sequence due to time domain aliasing (so specifically in the OP's case there is no way to know if the solution is representing longer response times than the samples are providing).
This is not very different from what I would typically do when evaluating results in the frequency domain: if I observe strong spectral content right up to the Nyquist boundary I am not assured that my sampling rate is high enough (or filtering is not tight enough) as frequency domain aliasing is likely occurring. While if the spectrum rolls off to sufficiently low levels prior to the Nyquist boundary (and if I can assume the system in question would not have higher spectral content at further offsets) then I am reasonably assured that the digital representation accurately represents the continuous time spectrum. This obviously requires knowledge of the system and what the likely spectral content is beyond the sampling rate, but is very analogous to the time domain challenge presented here.
For further details on matching an underlying CTFT with the DFT scroll down to the "Exact CTFT Result from DFT" section of my answer at this post: Why do we have to rearrange a vector and shift the zero point to the first index, in preparation for an FFT?