Idea behind Carleson's theorem modern proof "intitial reductions" I'm having troubles to understand the philosophy behind the modern proof of Carleson's theorem. For convenience, let me state precisely what I am asking for. 
For any $f \in L^2(\mathbb{R})$, let $\mathcal{Cf}:=\sup_{N \in \mathbb{Z}} \left\vert P_{-}(e^{iN\cdot}f)\right\vert$ be the maximal Carleson operator; where $P_{-}$ is the projection on negative Fourier spectrum $\{\xi <0\}$.  The Carleson theorem essentially states that:

Theorem (Carleson). $|\{\mathcal{C}f > \lambda\}|_{L^2} \lesssim \lambda^{-2}\|f\|^2_{L^2}$. 

In the modern proof of this theorem [I'm reading this], one instead studies the operators defined by 
$Q_{\xi}f:= \sum_{s \in T} \mathbf{1}_{\omega_s^+}(\xi) \langle f, \varphi _s\rangle \varphi _s$, where the notation is as follows:


*

*$T$ denotes the set of all tiles $I_s \times \omega_s$ such that $I_s, \omega_s$ are dyadic intervals, such that the area of the tile $I_s \times \omega_s$ is one. 

*$\omega_s^+$ stands for the upper half of the interval, and the $\varphi_s$ are functions such that $\hat \phi_s$ have Fourier support inside $\omega_s^{-}$ (the lower half of the interval)


It is not difficult to pass from $Q_{\xi}$ to the Carleson operator, one can take averages and get that the operator:
$$Q:=\lim_{Y \to \infty} \frac{1}{Y^2}\int_{[1,2] \times [0,Y]^2} Dil^2_{2^{-\lambda}}Tr_{-y}Mod_{-\xi}Q_{\xi} Mod_{\xi}Tr_yDil_{2^{\lambda}}^2d\lambda dy d\xi\,,$$
commutes with translations and dilations and its Kernel is made of functions with Fourier support lying on $\{\xi >0\}$, thus this operator is $P_{-}$. 
My question is then:

How is it somewhat 'natural' to come up with the operators $Q_{\xi}$? How one can guess that such an operator has a similar behavior to that of $\mathcal C$? 

I think there is some 'discretization' idea behind but I do not see how is it natural in any sense. Put in another way my question is: starting from $\mathcal C$ and $P_{-}$ how does one introduces the operators $Q_{\xi}$?
Does anybody has some good insights?
 A: I have my own confusions here, but let me share my thoughts.
As you mention, there is a discretization here. If you want to decompose the operator $P_-$, you use the standard decomposition $\sum_k\hat{\varphi}_k = 1_{(-\infty,0]}$, where $\hat{\varphi}_k(\xi) := \hat{\varphi}(\xi/2^k)$ is supported at frequencies $\vert\xi\vert\sim 2^k$. You can think of $\varphi_k$ as attached to the tile $I_s\times\omega_s = [-2^{-k-1},2^{-k-1}]\times[-2^k,0]$ ---this tile doesn't belong to the mesh $\mathcal{D}$, but let's ignore these "unfortunate technicalities", as Fefferman put it.
Since $P_-^2f = P_-f$, then we get the decomposition 
$$
P_- f = \sum_{k,k'} \varphi_k*\varphi_{k'}*f.
$$
We may assume that $\varphi_k*\varphi_{k'} = 0$ unless $k=k'$; you can take Fourier transform to see that this is morally true. For each term in the series we get
$$
\begin{align}
\varphi_k*\varphi_k*f(x) &= \int f(z)\varphi(y-z)\varphi(x-y)\,dydz \\
&= \int \varphi_k(y)\int f(z)\varphi_k(x-z-y)\,dzdy \\
&= \int \textrm{Tr}_y\varphi_k(x)\langle f,\textrm{Tr}_y\varphi_k\rangle\,dy \\
&= \sum_{\vert I\vert = 2^{-k}}\frac{1}{2^{-k}}\int_{-2^{-k-1}}^{2^{-k-1}}2^{-\frac{k}{2}}\textrm{Tr}_{y+c(I)}\varphi_k(x)\langle f,2^{-\frac{k}{2}}\textrm{Tr}_{y+c(I)}\varphi_k\rangle\,dy
\end{align}
$$
In the third identity we used $\overline{\tilde{\varphi}} = \varphi$, where $\tilde{\varphi}(x) = \varphi(-x)$, because $\hat{\varphi}$ is real. In the last term, let's define $\textrm{Tr}_{c(I)}\textrm{Dil}_{2^{-k}}^2\varphi = \varphi_s$, where $s$ denotes the tile $(c(I)+[-2^{-k-1},2^{-k-1}])\times [-2^k,0]$. We rewrite then the last integral as the average
$$
\varphi_k*\varphi_k*f(x) = \frac{1}{2Y}\int_{-Y}^{Y}\sum_{\vert I\vert= 2^{-k}}\textrm{Tr}_y\varphi_s(x)\langle f,\textrm{Tr}_y\varphi_s\rangle\,dy,
$$
where $Y = 2^{-k-1}$; however, you can modify the argument above to see that actually you can take the limit $Y\to\infty$. Summing up we have
$$
P_-f(x) = \lim_{Y\to\infty}\frac{1}{2Y}\int_{-Y}^Y\sum_s \langle f,\textrm{Tr}_y\varphi_s\rangle\textrm{Tr}_y\varphi_s(x)\,dy,
$$
where the tiles $s$ are those here constructed.
We left open the assumption $\varphi_k*\varphi_{k'}=0$ unless $k=k'$. In the paper they took $\hat{\varphi}$ supported in an interval of length $\frac{1}{4}$, but it is then impossible, I think, to get $\sum_k\hat{\varphi}_k = 1_{(-\infty,0]}$. I suspect the the average in dilation helps to solve the problem here; in fact, I would try to find a function $\varphi$ such that $\sum_k \int_2^4\hat{\varphi}_k^2(t\xi)\frac{dt}{t} = 1_{(-\infty,0]}$, and such that the supports of $\hat{\varphi}_k$ and $\hat{\varphi}_{k'}$ are disjoint, but not sure. The square $\hat{\varphi}_k^2$ is to use the same trick $P^2_-$.
In any case, we are reduced to the operator
$$
Tf := \sum_s\langle f,\varphi_s\rangle\varphi_s.
$$
Now if we try to use it in the Carleson operator, we have to deal with
$$
\vert T(e^{iN\cdot}f)\vert = \vert\sum_s\langle f,\textrm{Mod}_{-N}\varphi_s\rangle\textrm{Mod}_{-N}\varphi_s\vert.
$$
Now there is another technicality, for the frequency support of $\textrm{Mod}_{-N}\varphi_s$ doesn't belong to the mesh $(j2^k,(j+1)2^k)$. The way out from this nuisance is to average over translations of the mesh by using $\textrm{Mod}_\xi$, much like we did above for the intervals $I$. 
Excuse me if I do not complete all the computations, but I think the idea it is more or less clear. The point is that the operator $Q_\xi$ allows to get rid from many technicalities, but mainly with problems with the relative position of the mesh $\mathcal{D}$.
