Nonstationary solutions for stationary ARMA equations By "stationary" I mean "weakly stationary".
Consider a "stationary" AR(1) equation:
$$X_t=\varphi X_{t-1}+\varepsilon_t,$$
where $t\in\mathbb{Z}$ are discrete time moments, $\varepsilon_t$ a zero-mean white noise (just some iid sequence), $\varphi\in(-1,1)$. It is well known that there is a stationary solution (that is, a discrete time series satisfying the equation). Denote it by $X_t.$ However, we can introduce another time series $Y_t=X_t+\varphi^t$, which appears to be a nonstationary solution for the "stationary" equation (clearly, $\mathbb{E}[Y_t]$ is not free of $t$, since $X_t$ is evidently zero-mean).
Given more general stationary AR($p$) process, is it possible to somehow damage the weak stationarity property? Or, in general, is it true that any stationary discrete time AR (or even ARMA) equation has a nonstationary solution?
 A: Terminology being used in the question is not quite correct. You're mixing up the model (or equations) and solution to the model. 
It does not make sense to speak of an equation (in this case, a system of stochastic difference equations) being stationary or non-stationary. Stationarity, of lack thereof, is a property of a solution. An equation can have stationary or non-stationary solutions.
What you have found are two solutions, one stationary and one non-stationary, to the AR(1) equation when AR parameter $|\phi| \neq 1 $. (If $|\phi| > 1$, substitute $-t$ for $t$ in your example.) In contrast, when $|\phi| = 1$, there are only non-stationary solutions.
The answer to your question is, yes, this generalizes to the AR(p) case. The AR(p) equation(s)
$$
\Phi(L)X_t = \epsilon_t, \; t = \cdots -1, 0, 1, \cdots
$$
has both stationary and non-stationary solutions if the polynomial $\Phi(z^{-p})$ has no roots on the unit circle and all roots are real.
For example, suppose the AR(2) model
$$
X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \epsilon_t
$$
has a stationary solution $(X_t)$ and $z^2 - \phi_1 z - \phi_2$ has two real roots $a$ and $b$, then
$$
X_t + a^t + b^{t-1}
$$
is a non-stationary solution. 
Notice setting $\phi_2 = 0$ and considering $z - \phi_1 = 0$ recovers your AR(1) example.
A: You are correct that $Y_t = X_t + \varphi^t$ is a non-stationary solution to the $\text{AR}(1)$ equation.  Thus, contra your terminology, it is not a "stationary equation" and it is not "damaged" by having a non-stationary solution.  The non-stationary solution is simply a non-stationary model that is a form of $\text{AR}(1)$ process, consistent with the recursive equations.  There are also non-stationary solutions to the $\text{AR}(p)$ equations.
A: If you let your process go on, then you'll notice how the term $\varphi^t$ disappears:
$$\lim_{t\to\infty}\varphi^t=0$$
So, $$E[Y_t]=E[X_t]+E[\varphi^t]=_{t\to\infty}0$$
this is despite the right hand side having being dependent on finite $t$.
So the answer to your question is that your process $Y_t$ is not non-stationary. Hence, it doesn't serve as a counter-example.
Additional thoughts.
You formulated your question in terms of solutions of the stochastic processes. Look at what is the solution of the AR(1) process.
For instance, if you forecast $\tau$ steps ahead you get: $$X_{t+\tau}=\varphi^\tau (X_t+\sum_{s=1}^\tau \varepsilon_{t+s}\varphi^{-s})$$
You can see how it simply collapses to the noise around zero as $\tau$ grows, no matter what was the initial $X_t$. When you add your $\varphi^\tau$ term it also disappears, so the stable solution is the same: noise around zero:
$$X_{t+\tau}+\varphi^\tau=\varphi^\tau (X_t+\sum_{s=1}^\tau \varepsilon_{t+s}\varphi^{-s}+\varphi^s)$$
