One can use non-commutative probability theory to compute the spectral distributions of $X_\theta$ and $X_\theta - 2 \tan(\theta) P$ with respect to $\delta_e$.  Since $G = \mathbb{Z}_2 * \mathbb{Z}_2$, the operators $T_s$ and $T_t$ are freely independent with respect to the vector state given by $\delta_e$.  Each of them has the spectral distribution given by the Bernoulli measure $(1/2)(\delta_{-1} + \delta_1)$.  Thus, the spectral distribution of $T_s + T_t$ is the free convolution of two Bernoulli distributions.  Using the $R$-transform (see Chapter 4 of https://arxiv.org/pdf/1908.08125.pdf), we can compute that the free convolution is the measure $\mu$ satisfying
$$
F_\mu(z) = \sqrt{z^2 - 4},
$$
where $F_\mu(z) = 1 / \int (z - t)^{-1}\,d\mu(t)$, also known as the *$F$-transform*.  Here we use the square root which is defined on the upper half-plane and close to $z$ when $z \to \infty$.  This means that $\mu$ is an arcsine law
$$
d\mu(x) = \frac{1}{\pi \sqrt{4 - x^2}} \chi_{(-2,2)}(x)\,dx
$$
This computation is known; see the last full paragraph at the bottom of page 2 here: https://arxiv.org/pdf/1008.5205.pdf.

Next, let $\nu$ be the spectral measure of $T_s + T_t - 8 \tan(\theta)$.  This just shifts the measure $\mu$ by $-8 \tan(\theta)$, and we have $F_\nu(z) = F_\mu(z + 8 \tan(\theta))$.  Since the support of $\mu$ is $[-2,2]$, the support of $\nu$ is $[-2 - 8 \tan(\theta), 2 - 8 \tan(\theta)]$ and thus the norm of $X_\theta$ is $2 + 8 |\tan \theta|$, at least on the cyclic subspace generated by $\delta_e$.  However, because $\delta_e$ is cyclic for the right shift operators by $s$ and $t$, it is cyclic for the commutant of $X_\theta$ and hence the norm of $X_\theta$ on the entire space agrees with the norm on the cyclic subspace generated by $\delta_e$.

Finally, let's compute the spectral measure of $X_\theta - 2 \tan(\theta) P$, which we will call $\rho$.  This can be done using rank-one perturbation formulas.  Alternatively, $X_\theta$ and $P$ are Boolean independent (see https://mast.queensu.ca/~speicher/papers/boolean.ps).  The spectral measure of $-2 \tan(\theta) P$ with respect to $\delta_e$ is $\delta_{-2 \tan(\theta)}$.  The theory of Boolean independence tells us that
$$
F_\rho(z) - z = (F_\nu(z) - z) + (F_{\delta_{-2\tan(\theta)}}(z) - z) = F_\nu(z) - z + 2 \tan(\theta).
$$
Hence,
$$
F_\rho(z) = \sqrt{(z + 8 \tan \theta)^2 - 4} + 2 \tan(\theta).
$$
Thus, for $z$ in the upper half-plane,
$$
G_\rho(z) = \frac{1}{\sqrt{(z + 8 \tan \theta)^2 - 4} + 2 \tan(\theta)},
$$
where $G_\rho(z) = \int (z - t)^{-1} \,d\rho(t)$.  Using the fact that this is rank-one perturbation that the norm of the new operator $X_\theta - 2 \tan(\theta) P$ is the maximum of the norm of $X_\theta - 2 \tan(\theta) P$ on the cyclic subspace generated by $\delta_e$ and the norm on the orthogonal complement, which will be bounded by $\lVert X_\theta \rVert$. Thus, to show that the two operators have different norms, it suffices to show that the support radius of $\rho$ is strictly larger than that of $\nu$.

Looking at $G_\rho$, the support of the measure $\rho$ agrees with the support of $\nu$ except with the addition of atoms at the point where $(2 \tan \theta)^2 = (z + 8 \tan \theta)^2 - 4$ according to the version of the square root described above.  Assume $\tan \theta > 0$ since the other case is symmetric (and recall $\theta$ was assumed to be in $(-\pi/2,\pi/2)$).  Then we are looking at a point on the negative side where $$(z + 8 \tan \theta)^2 = 4 + 4 \tan^2 \theta = 4 \sec^2 \theta.$$  Hence, $z = -2 \sec \theta - 8 \tan \theta$ is the location of the atom.  Since the support of $\nu$ was already farther on the negative side, we have
$$
\lVert X_\theta \rVert = 2 + 8 |\tan \theta| < 2 \sec \theta + 8 |\tan \theta| = \lVert X_\theta - (2 \tan \theta) P \rVert.
$$

Similar computations can be done for $\mathbb{Z}_2^{*n}$.  Let $T_1$, \dots, $T_n$ be the operators corresponding to the generators.  Then each $T_1 + \dots + T_n$ as a spectral distribution $\mu_n$ is the free convolution of $n$ Bernoulli laws.  This can be done directly with the $R$-transform, or using Proposition 3.1 of this paper: https://arxiv.org/pdf/math/0703295.pdf.  This results in
$$
F_{\mu_n}(z) = \frac{(n-2)z + n \sqrt{z^2 - 4(n-1)}}{2(n-1)}.
$$
Thus, similar computations can be done as in the $n = 2$ case.

Unfortunately, for general groups, it is not as easy to compute because the operators might not be self-adjoint.  However, the tools of free probability could still be used in theory to compute the spectral distribution of $X^*X$.