Here, "simultaneous measurement" of an observable does not refer to simultaneously in time, but rather the ability to estimate several observables from the same measurement setting acting on $\rho\otimes \rho$ (and therefore using only a single copy of $\rho \otimes \rho$). Basically, you may estimate both $|\text{Tr}(\rho P_i)|$ and $|\text{Tr}(\rho P_j)|$ for $P_i \neq P_j$ using $O(1)$ copies$^1$ of $\rho \otimes \rho$. However, because estimation involves statistical fluctuations, you can't naively scale this up for an arbitrary number of Pauli observables.
It may make more sense to think of the protocol like this: The paper says you can attempt to estimate $\lvert\text{Tr}(P \rho)\rvert^2$ within $\epsilon$ error for all $n$-qubit Pauli operators $P$ using just $O(1)$ copies of $\rho \otimes \rho$ - instead of , say, $O(n)$. But if you do this, there is a decent chance that a few of your estimates are off by more than $\epsilon$. So you scale up the number of measurements for each observable by $\log (4^n)$ to guarantee that all observables are estimated to $\epsilon$ accuracy with high enough probability, recovering the $O(n)$ scaling. This is the meaning of the discussion about the "union bound" in Appendix E.2.c.
Example of simultaneous measurement with Bell basis
For completeness, I will show how a single Bell measurement configuration may be used on $\rho \otimes \rho$ to estimate $|\text{Tr}(\rho P_i)|$ "simultaneously". Define the Bell state
\begin{equation}
|\Phi_{00}\rangle := (|00\rangle + |11\rangle)/\sqrt{2}, \tag{1}
\end{equation}
and then define the rest of a Bell basis according to
\begin{equation}
|\Phi_{ij}\rangle := (Z^i \otimes X^j)|\Phi_{00}\rangle. \tag{2}
\end{equation}
Measuring $\rho\otimes \rho$ in the Bell basis means using a PVM $\{|\Phi_{ij}\rangle \langle \Phi_{ij}|\, \forall \, i,j\in\{0,1\}\}$, and observing outcome $(i,j)$ with probability $p_{ij} =\langle \Phi_{ij} |(\rho\otimes \rho)|\Phi_{ij}\rangle$. We repeat this many times (see footnote 1) to produce an estimate $\hat{p}_{ij}$ for the probability distribution $p_{ij}$.
Now that we have a way to estimate $p_{ij}$, we would like to express each $|\text{tr}(\rho P_j)|$ in terms of these probabilities. Since the Pauli operators form a orthonormal basis for Hermitian matrices, we decompose $\rho \otimes \rho$ as
\begin{align}
\rho \otimes \rho &= \sum_{ij} \text{Tr}\left[(\rho \otimes \rho)(P_i \otimes P_j)\right](P_i \otimes P_j) \tag{3}
\\&:= \sum_{i}c_{ii}^2 (P_i \otimes P_i) + \sum_{i\neq j} c_{ij}(P_i \otimes P_j), \tag{4}
\end{align}
where we identify $c_{ii}^2 = |\text{Tr}(\rho P_i)|^2$ as the quantities we're looking for. Since $n=1$, we can compute by hand the probability of observing each outcome $(i,j)$. A useful identity will be
\begin{align}
|\Phi_{00}\rangle \langle \Phi_{00}| = \frac{1}{4} (I\otimes I + X\otimes X - Y\otimes Y + Z\otimes Z). \tag{5}
\end{align}
Then, compute:
\begin{align}
p_{00} &= \text{Tr}\left( |\Phi_{00}\rangle\langle \Phi_{00} | \rho \otimes \rho \right) \tag{6a-c}
\\&=\frac{1}{4} \text{Tr}\left[\left( I\otimes I + X\otimes X - Y\otimes Y + Z\otimes Z \right) \left(\sum_{i}c_{ii}^2 (P_i \otimes P_i) + \sum_{i\neq j} c_{ij}(P_i \otimes P_j) \right)\right]
\\&= \frac{1}{4} \left(c_{00}^2 + c_{11}^2 - c_{22}^2 + c_{33}^2\right) .
\end{align}
Similarly, using Eq. (2) and some Pauli conjugation identities, one finds
\begin{align}
p_{01} &= \text{Tr}\left( (I\otimes X)|\Phi_{00}\rangle\langle \Phi_{00} | (I \otimes X)\rho \otimes \rho \right) \tag{7}
\\&= \frac{1}{4}\left(c_{00}^2 + c_{11}^2 + c_{22}^2 - c_{33}^2\right) \tag{8}\\
p_{10} &= \text{Tr}\left( (Z\otimes I)|\Phi_{00}\rangle\langle \Phi_{00} | (Z \otimes I)\rho \otimes \rho \right) \tag{9}
\\&= \frac{1}{4}\left(c_{00}^2 - c_{11}^2 + c_{22}^2 + c_{33}^2\right) \tag{12}
\\p_{11} &= \text{Tr}\left( (Z\otimes X)|\Phi_{00}\rangle\langle \Phi_{00} | (Z \otimes X)\rho \otimes \rho \right) \tag{10}
\\&= \frac{1}{4}\left(c_{00}^2 - c_{11}^2 - c_{22}^2 - c_{33}^2\right) \tag{11}.
\end{align}
Rearrange these equations to solve for the desired expectation values in terms of $p_{ij}$:
\begin{align}
c_{00}^2 = |\text{Tr}(\rho I)|^2 &= p_{00} + p_{01} + p_{10} + p_{11} = 1 \tag{12}\\
c_{11}^2 = |\text{Tr}(\rho X)|^2 &= p_{00} + p_{01} - p_{10} - p_{11}\tag{13}\\
c_{22}^2 = |\text{Tr}(\rho Y)|^2 &= -p_{00} + p_{01} + p_{10} - p_{11}\tag{14}\\
c_{33}^2 = |\text{Tr}(\rho Z)|^2 &= p_{00} - p_{01} + p_{10} - p_{11} \tag{15}\\
\end{align}
Then, just define estimates for the expectation values $|\text{Tr}(\rho P_i)|:=\hat{c}_{ii}$ that are computed from $\hat{p}_{ij}$ in an analogous manner to the above system. Thus, for $n=1$ we see how using the estimated probability distribution $\{\hat{p}_{ij}\}$ output by a single measurement setting in the Bell basis , we may output estimates for all $|\text{Tr}(\rho P_i)|$.
$^1$ Here the $O(1)$ includes repeating the measurement to get an estimate of the average and hides terms in $\epsilon$ and $\delta$ that control the accuracy of this estimate.
