Given the constraints of the problem, we have for $U$ unitary,
$$
\newcommand{\ketbra}[1]{\lvert#1\rangle\!\langle #1\rvert}
\newcommand{\ket}[1]{\lvert#1\rangle}
\newcommand{\bra}[1]{\langle#1\rvert}
\newcommand{\sqoverlap}[2]{|\langle #1|#2\rangle|^2}
U=\ketbra{\lambda}+\sum_k e^{i\varphi_k}\ketbra{\lambda_k},$$
for some orthonormal basis $\{\ket\lambda\}\cup\{\ket{\lambda_k}\}_k$ and phases $\varphi_k\in\mathbb R$.
Given a target $\ket\mu$, the expectation value of $U^r$ equals
$$
\bra\mu U^r\ket\mu=\sqoverlap{\mu}{\lambda} +
\underbrace{
\sum_k e^{ir\varphi_k} \sqoverlap{\mu}{\lambda_k}
}_{\equiv F_r}.
$$
Note also that the normalisation of $\ket\mu$ implies
$$
\sqoverlap\mu\lambda +
\sum_k |\langle\mu|\lambda_k\rangle|^2 = 1
$$
The problem here is clearly in how to handle $F_r\in\mathbb C$, which has a pretty nasty behaviour: as a function of $r$, it oscillates wildly in the complex plane. Without any constraint on the possible $\ket{\lambda_k}$ and $\varphi_k$, or on the total dimension of $U$, pretty much any shape is possible (and here I cannot not link this beautiful video by 3Blue1Brown). Just to give a mental picture of what I mean, here are two examples of what you can get when you plot $r\mapsto F_r$ for randomly sampled phases and overlaps, for a six-dimensional $U$:

The black dots correspond to the integer values of $r$, the blue line gives the value of $F_r$ for the real values in-between, and the red circle has radius $1-\sqoverlap{\mu}{\lambda}$. Only values of $r$ between $1$ and $100$ are shown.
Plotting even more points (all $r$ between $1$ and $10^3$), we see that the circle is essentially filled:

This is telling us that, in the general case, all we can say about $F_r$ is that it satisfies
$$
|F_r|\le1-\sqoverlap\mu\lambda.
$$
Now, to get the probability of finding the initial state $\ket\mu$ unchanged after $r$ applications of $U$, we need to study the quantity
$|\langle\mu|U^r|\mu\rangle|=\big|\sqoverlap\mu\lambda+F_r\big|$.
When $\sqoverlap\mu\lambda\ge1/2$ we know that $|F_r|\le\sqoverlap\mu\lambda$, and we thus have the bound
$$
(2\sqoverlap{\mu}{\lambda}-1)^2 \le |\langle\mu|U^r|\mu\rangle|^2 \le 1,\quad \forall r.
$$
On the other hand, when the initial squared overlap is smaller than $1/2$, we can have $F_r=-\sqoverlap\mu\lambda$, and thus we get the trivial bound $|\langle\mu|U^r|\mu\rangle|\in[0,1]$.
Putting the two cases together, if we plot the possible values of $|\langle\mu|U^r|\mu\rangle|^2$ as a function of $p\equiv\sqoverlap\mu\lambda$ we get the following:

The interesting conclusion is that we can say something about the possible probabilities of a state remaining unchanged after a number of applications of the same unitary $U$, provided that the initial overlap between the state and one of the eigenvectors of $U$ is big enough.
Let me remark a few other things here:
There is no loss of generality in our initial assumption that the eigenvalue of some $\ket\lambda$ is $1$. We can repeat the same identical argument by just taking an arbitrary eigenvector of $U$ with eigenvalue $e^{i\varphi}$, by only replacing each instance of $U$ with $U e^{-i\varphi}$. Most notably, this does not affect the final result on $|\langle\mu|U^r|\mu\rangle|$.
If $\sqoverlap\mu\lambda=1-\delta$ with $\delta\ge0$ small, then the bound can be simplified to $|\langle\mu|U^r|\mu\rangle|^2\ge 1-4\delta$.
Because we are not putting any constraint on $U$ or $r$ other than the initial overlap $\sqoverlap\mu\lambda$, we can also understand this result as telling us about the probability of $\ket\mu$ being unchanged upon a single application of a unitary $U$, or more generally about the probability of $\ket\mu$ being unchanged upon arbitrary applications of arbitrary unitaries all having an eigenvector close enough to $\ket\mu$.
Mathematica code to generate the plots of $F_r$:
dim = 5;
probs = 0.8 * # / Total @ # & @ RandomReal[{0, 1}, dim];
phases = RandomReal[{0, 2 Pi}, dim];
otherFactor[probs_, phases_, r_] := Total[Exp[I r phases] probs];
Show[
ParametricPlot[ReIm@otherFactor[probs, phases, r], {r, 0, 100},
PlotPoints -> 60, MaxRecursion -> 10,
PlotRange -> ConstantArray[{-.9, .9}, 2]],
Graphics[{
Point@Table[ReIm@otherFactor[probs, phases, r], {r, 1, 100}],
Thick, Red, Circle[{0, 0}, 0.8]
}]
]