The paper does not specify the exact algorithm or class of distributions $\mathcal{D}$
for which such algorithm fails to refute XQUATH, and some classes of distributions
$\mathcal{D}$ do not satisfy XQUATH, so some additional assumptions about $\mathcal{D}$ are needed.
In this answer we will try to follow
the idea of the paper and try to show that it fails for a wide range of distributions $\mathcal{D}$.
This answer is currently incomplete though (see the bottom).
We assume that circuits $C \in \mathcal{D}$ are expressed as
$C = \pi_m U_m \pi_{m-1} U_{m-1} \cdots \pi_1 U_1 \pi_0$ where
- $\pi_l$ are "classical" in a sence that their matrix elements satisfy $\left|(\pi_l)_{ij}\right| \in \{0, 1\}$.
- Operators $U_{l}$ are random unitaries operating on a small number of qubits.
- For each structure of $C$, which includes $m$, $\{\pi_l\}_{l=0..m}$,
and qubits on which $U_l$ operate, unitaries $U_l$ are independent and each of them is picked
from a distribution $\mathcal{D}_l$.
- $\mathcal{D}_l$ is closed under multiplication by diagonal matrices with
$\pm 1$ or $\pm i$ on the diagonal.
- For each fixed structure of $C$ the expectation of $p_0$ over the choice of $\{U_l\}$ is $2^{-n}$.
- Some additional assumptions may be needed to show $\beta_2 = O(1/N^\epsilon)$ for $\epsilon > 0$ below.
The assumption 5 is to exclude circuits where some bits of the output can be
recovered using a classical simulation ignoring quantum effects (i.e. ignoring $U_l$-s).
E.g. consider a circuit where qubit 1 is only used by the following operations: (1) NOT gates;
(2) control gates where qubit 1 is used as a control for some classical operation on other qubits;
(3) a single CCNOT operation on qubit 1 controlled by conjunction of at least 2 other qubits.
In this case we can predict the state of qubit 1 at the end of the circuit with probability of
at least 3/4 assuming control qubits for the operation (3) are as likely to be in
$\left|11\right>$ state as in any other of 4 states, and, thus,
refute XQUATH for that particular $\mathcal{D}$.
We want to maximize $I = 2^{2n}\mathbb{E}\left((p_0 - 2^{-n})^2 - (p_0 - p)^2\right)$. In particular,
if we manage to find an algorithm providing $p$ in polynomial time and satisfying $I = \Omega(2^{-n})$
then we will refute XQUATH. Let's assume for now that the
structure of $C$ is fixed and the expectation in the formula for $I$ is with respect to the choice
of $U_l$. Let $p_0 = 2^{-n} (1 + \tilde p_0)$ and let's
search for $p$ in the form of $p = 2^{-n} (1 + \alpha_1 p_1)$ where $\alpha_1$ depends only on
the structure of $C$ and not on $U_l$ or the random choices within our algorithm. Then
$$I = 2 \alpha_1 \mathbb{E}(\tilde p_0 p_1) - \alpha_1^2 \mathbb{E}(p_1^2).$$
The optimal choice of $\alpha_1$ is $\alpha_1 = \mathbb{E}(\tilde p_0 p_1) / \mathbb{E}(p_1^2)$ which gives
$$I = \left.\left(\mathbb{E}(\tilde p_0 p_1)\right)^2\middle/\mathbb{E}(p_1^2)\right..$$
In order to come up with the estimator $p_1$, following the idea from the paper,
we consider the Feynman paths $z = (z_0, \dots, z_{2m+1})$. We can write
$p_0 = \left|\sum_z \left<0\right|C_z\left|0\right>\right|^2$ where
$$C_z = \left|z_{2m+1}\right>\left<z_{2m+1}\right|\pi_m \left|z_{2m}\right>
\left<z_{2m}\right| U_{m} \left|z_{2m-1}\right>... \left<z_0\right|.$$
Some of these paths trivially result in $\left<0\right|C_z\left|0\right> = 0$,
so we can restrict our attention to nontrivial ones, i.e. where $z_{2m+1} = 0$, $z_0 = 0$,
and $\left|\left<z_{2l+1}\right|\pi_l \left|z_{2l}\right>\right| = 1$. Assume there are $N$ nontrivial Feynman paths and pick random $k$ of them to compute $p_1$. In order to simplify the notation, we assume that the picking procedure works in the following way: first, we pick a random permutation $\sigma$,
apply it to the list $\{2^{n/2}\left<0\right|C_z\left|0\right>\}_{z}$
of coefficients resulting from nontrivial Feynman paths.
Then denote the resulting list of coefficients as $a_1,\dots,a_N$,
and take $a_1,\dots,a_k$ as our sample. Given that
$$\tilde p_0 = \left|\sum_{i=1}^N a_i\right|^2 - 1 = \sum_{i=1..N} \left(\left|a_i\right|^2 - \frac1N\right) + \sum_{i,j=1..N;i\neq j} a_i \bar a_j,$$
a reasonable choice of $p_1$ is
$$p_1 = \alpha_2 p_2 + \alpha_3 p_3\textrm{, where}\; p_2 = \sum_{i=1..k} \left(\left|a_i\right|^2 - \frac1N\right)\!,\; p_3 = \sum_{i,j=1..k;i\neq j} a_i \bar a_j.$$
Let
$$
\mathbb{E}\left(\left|a_1\right|^2 - \frac1N\right)^2 = \beta_1/N^2,\qquad
\mathbb{E}\left(\left(\left|a_1\right|^2 - \frac1N\right)\left(\left|a_2\right|^2 - \frac1N\right)\right) = \beta_2/N^2.$$
Then if we can show that $\beta_2 = O(1/N)$, then by computing $\mathbb{E}(\tilde p_0 p_s)$ and $\mathbb{E}(p_s^2)$ for $s=2,3$ we can show that $I = O(k^2/N^2)$ and, hence, show that Feynman path averaging algorithm above does not refute XQUATH (since $\log_2 N = \sum_{l=1}^{m} \log_2 \dim U_l$). It is not clear, though, which additional assumptions are needed to show $\beta_2 = O(1/N)$ or, at least, $\beta_2 = O(1/N^\epsilon)$ for some $\epsilon > 0$.
