Consider two $n$ qubit quantum states $| \psi_1 \rangle$ and $| \psi_2 \rangle$, where $$ |\langle\psi_2|\psi_1\rangle|^2 \geq 1 - \epsilon. $$
Let $p_x = |\langle x| \psi_1\rangle|^2$ and $q_x = |\langle x| \psi_2\rangle|^2$, for any $x, y \in \{0, 1\}^n$.
I am trying to get upper and lower bounds on the following two quantities:
$$\sum_{x=1}^{2^n} q_x p_x,$$ and $$\sum_{x=1}^{2^n} q_x p^2_x,$$ in terms of known quantities, like $\epsilon$, $p_x$ and $q_x$.
Given you only have a lower bound for $|\langle\psi_2|\psi_1\rangle|^2$ you cannot expect to non-trivially upper bound your quantity: Indeed, for $\psi_1\approx|0\ldots0\rangle\approx\psi_2$ you have $\sum_xq_xp_x,\sum_xq_xp_x^2\approx 1$ (which is as big as these can get). That being said, note the following
Result. The lower bound $$ \boxed{\sum_{x=1}^{2^n}q_xp_x\geq \frac{1-\varepsilon}{2^n}}\tag1 $$ holds and is tight, that is, for all $n$ and all $\varepsilon\in[0,1]$, (1) is valid and there exist pure states $\psi_1,\psi_2\in\mathbb C^{2^n}$ satisfying $|\langle\psi_2|\psi_1\rangle|^2\geq 1-\varepsilon$ such that equality holds in (1).
For the other quantity $\sum_xq_xp_x^2$, numerics suggest that the lower bound is in the realm of $\frac{1-\varepsilon}{2^{2n}}$ although for larger $\varepsilon$ this seems to deviate a bit so there seems to be a bit more to it. In any case one way to show (1) is through the following basic lemma:
Lemma. Given any $A\in\mathbb C^{d\times d}$ define $\mathcal D(A):={\rm diag}(A_{11},\ldots,A_{dd})$ as the usual dephasing channel. Then for all $A\in\mathbb C^{d\times d}$ it holds that $|{\rm tr}(A)|\leq \sqrt d\|\mathcal D(A)\|_2$; equivalently, $$ \Big|\sum_i A_{ii}\Big|^2\leq d\Big(\sum_i|A_{ii}|^2\Big)\,. $$
Proof. The trick is to re-write $\sum_i A_{ii}$ as an inner product $\langle {\bf e}|a\rangle$ where $a:=(A_{ii})_i\in\mathbb C^d$ and ${\bf e}:=(1,\ldots,1)^\top\in\mathbb C^d$. Then we can use Cauchy–Schwarz to estimate \begin{align*} \Big|\sum_i A_{ii}\Big|^2=|\langle {\bf e}|a\rangle|^2\leq\langle{\bf e}|{\bf e}\rangle\langle a|a\rangle=d\Big(\sum_i|A_{ii}|^2\Big)\,. \end{align*} Taking the square root then concludes the proof. $\ \square$
With this in mind, (1) is recovered by setting $d=2^n$ together with \begin{align*} \sum_{x=1}^{2^n}q_xp_x=\sum_{x=1}^{2^n}\big|\big\langle x\big| \; |\psi_1\rangle\langle\psi_2| \; \big|x\big\rangle\big|^2&\geq \frac1{2^n}\Big| \sum_x \big\langle x\big| \; |\psi_1\rangle\langle\psi_2| \; \big|x\big\rangle\Big|^2 \\ &= \frac1{2^n}\big|{\rm tr}\big( |\psi_1\rangle\langle\psi_2| \big)\big|^2\\ &=\frac1{2^n}|\langle\psi_2|\psi_1\rangle|^2\geq \frac{1-\varepsilon}{2^n} \end{align*}
To see that the bound (1) is optimal and cannot be improved further, consider the following example. For an $n$-qubit system and any $\varepsilon\in[0,1]$ define $$ \mu_{n,\pm}:=\sqrt{\frac{1\pm\sqrt{ 1-\sqrt[n]{1-\varepsilon} }}2}\in[0,1]\,. $$ It is easy to see that $|\mu_{n,+}|^2+|\mu_{n,-}|^2=1$ so both $\psi_{1,n}:=(\mu_{n,+},\mu_{n,-})^\top\in\mathbb C^2$ and $\psi_{2,n}:=(\mu_{n,-},\mu_{n,+})^\top\in\mathbb C^2$ are valid 1-qubit pure states. From this we can build our saturating vectors: $$ \psi_1:=\bigotimes_{j=1}^n\psi_{1,n}=\psi_{1,n}\otimes\psi_{1,n}\otimes\ldots\otimes\psi_{1,n}\in\mathbb C^{2^n} $$ as well as $\psi_2:=\bigotimes_{j=1}^n\psi_{2,n}$. We make the following
Claim.
- $|\langle\psi_2|\psi_1\rangle|^2=1-\varepsilon$
- $\sum_{x=1}^{2^n}q_xp_x=\frac{1-\varepsilon}{2^n}$
which would obviously show optimality of (1) for all $n,\varepsilon$. To see 1. note that $2\mu_{n,+}\mu_{n,-}=\sqrt[2n]{1-\varepsilon}$ which implies \begin{align*} |\langle\psi_2|\psi_1\rangle|^2&=\Big|\Big\langle\bigotimes_{j=1}^n\psi_{2,n}\Big|\bigotimes_{j=1}^n\psi_{1,n}\Big\rangle\Big|^2\\ &=\prod_{j=1}^n|\langle\psi_{2,n}|\psi_{1,n}\rangle|^2=|\langle\psi_{2,n}|\psi_{1,n}\rangle|^{2n}\\ &=|2\mu_{n,+}\mu_{n,-}|^{2n}=1-\varepsilon\,, \end{align*} as claimed. For 2. one computes similarly \begin{align*} \sum_{x=1}^{2^n}q_xp_x&=\sum_{j_1,\ldots,j_n\in\{0,1\}}\Big|\Big\langle j_1\ldots j_n\Big|\bigotimes_{j=1}^n\psi_{2,n}\Big\rangle\Big\langle j_1\ldots j_n\Big|\bigotimes_{j=1}^n\psi_{1,n}\Big\rangle\Big|^2\\ &=\sum_{j_1,\ldots,j_n\in\{0,1\}}\prod_{k=1}^n|\langle j_k|\psi_{2,n}\rangle\langle j_k|\psi_{1,n}\rangle|^2\\ &=\Big(\sum_{j\in\{0,1\}}|\langle j|\psi_{2,n}\rangle\langle j|\psi_{1,n}\rangle|^2\Big)^n\\ &=\big(2\mu_{n,+}^2\mu_{n,-}^2)^n\\ &=\Big(\frac{(2\mu_{n,+}\mu_{n,-})^2}2\Big)^n=\frac{(\sqrt[2n]{1-\varepsilon})^{2n}}{2^n}=\frac{1-\varepsilon}{2^n} \end{align*} which is what we had to show.
NB: Not only is the bound (1) valid for arbitrary qudit systems (because the lemma we used to prove it works in arbitrary dimensions), but numerics strongly suggest that the generalized bound $\frac{1-\varepsilon}{d}$ is tight for all dimensions $d$. However, for $d\neq 2^n$ writing down analytic expressions for the optimal $\psi_1,\psi_2$ seems to be more difficult; this is evidenced by the fact that already for one qutrit ($d=3$) and any $\varepsilon\in[0,1]$, the pure states $$ \psi_1=\begin{pmatrix}a\\\sqrt{ac}\\c\end{pmatrix}\qquad\psi_2=\begin{pmatrix}c\\\sqrt{ac}\\a\end{pmatrix} $$ which satisfy $|\langle\psi_2|\psi_1\rangle|^2=1-\varepsilon$ and $\sum_{x=1}^3|\langle x|\psi_1\rangle\langle\psi_2|x\rangle|^2=\frac{1-\varepsilon}3$ come about via the cumbersome expressions $$ a=\frac{\sqrt{-\sqrt{1-\varepsilon }+\sqrt{3} \sqrt{\varepsilon -2 \sqrt{1-\varepsilon }+2}+3}}{\sqrt{6}} $$ and $$ c=\tfrac{3 \sqrt{2} \sqrt{\sqrt{1-\varepsilon }-\sqrt{3} \sqrt{\varepsilon -2 \sqrt{1-\varepsilon }+2}+5}-\sqrt{6} \sqrt{-\sqrt{1-\varepsilon }+\sqrt{3} \sqrt{\varepsilon -2 \sqrt{1-\varepsilon }+2}+3}}{12} $$ This solution as well as numerics suggest that for all $d$ there exist saturating $\psi_1,\psi_2$ in $[0,1]^d$ (i.e., one can choose the entries to be non-negative), and one can even choose $\psi_1$ to have decreasing entries and $\psi_2$ to have increasing entries. However, since OP asked about the qubit case in particular I'll leave it at that.
