I'm looking at the following definition of a causal AR(p) (autoregressive) model:
An AR(P) model $\phi(B)x_t= \epsilon_t$ is said to be causal if it has a stationary solution $$x_t=\epsilon_t +\sum^{\infty}_{i=1}b_i\epsilon_{t-i}$$ where the $\sum^{\infty}_{i=1}|{b_i}|<\infty$ and $\{\epsilon_t\}$ is a white noise process.
How do we know that this solution is stationary?
There are technical issues associated with the infinite sum. Identifying and resolving them requires some attention to the definitions, so let's begin there.
A white noise process $\epsilon_t$ (in discrete time indexed by the integers) when (a) its mean $E[\epsilon_t]=0$ for all $t$ and (b) its autocorrelation function $$R_\epsilon(h,t)=E[\epsilon_t\epsilon_{t+h}]$$ is finite for all integral $h$ and zero when $h\ne 0.$ See https://en.wikipedia.org/wiki/White_noise#Discrete-time_white_noise. Accordingly, we may drop references to $t$ in the notation for the autocorrelation.
A stationary process $X_t$ is one for which given any dimension $n\ge 1$ and any vector $(h_1,\ldots, h_n)$ of integers, the joint distribution of $(X_{t+h_1}, X_{t+h_2}, \ldots, X_{t+h_n})$ does not vary with $t.$ See https://en.wikipedia.org/wiki/Stationary_process#Definition.
(There are weaker definitions of stationarity. You will have no trouble establishing the result for any of them by emulating the argument used here for this strong stationarity condition.)
A countable linear combination of random variables $X_0, X_1, X_2, \ldots, X_n, \ldots$ given by coefficients $b_0, b_1, b_2, \ldots$ is the limit (if it exists) of the finite partial sums, $$\sum_{i=0}^\infty b_i X_i = \lim_{n\to\infty} \sum_{i=1}^n b_i X_i.\tag{*}$$
But limit in what sense? The strictest sense is the pointwise limit of the random variables; but the one most germane to models is the limit in distribution. It doesn't matter which definition we adopt, as you will soon see.
First notice that the definition of a white noise process is so general that it doesn't even guarantee a white noise process is stationary! When we begin with a non-stationary process $\epsilon_t,$ the stated result will generally be false. As a counterexample, let $b_i=0$ for all $i:$ it asserts $\epsilon_t + \sum b_i \epsilon_{t-i} = \epsilon_t$ is stationary, flatly contradicting the assumption.
To make any progress, then, we must assume that $\epsilon_t$ is a stationary white noise process.
With this assumption, the derivation I gave at https://stats.stackexchange.com/a/566582/919 shows that any finite linear combination of a stationary process is stationary. Applying this result to the situation in the question with the coefficients $(1, b_1, \ldots, b_{n})$ shows
The process $X^{(n)}_t = \epsilon_t + \sum_{i=1}^n b_i \epsilon_{t-i}$ is stationary for $n=1, 2, \ldots.$
Thus, every partial sum on the right hand side of $(*)$ is a stationary process.
We need to show the limit of a sequence of stationary processes is itself stationary. But this is trivial: the definition of stationarity means all finite $n$-variate marginal distributions are unchanged by time translation, so upon taking limits we deduce all finite $n$-variate marginal distributions of the limit are unchanged, whence (by definition of stationarity) the limiting process is stationary. Applying this to the sequence $n\to X^{(n)}$ of processes implies
$$X_t = \lim_{n\to\infty} X_t^{(n)} = \epsilon_t + \sum_{i=1}^\infty b_i \epsilon_{t-i}$$
is stationary, QED.
