How to properly make sense of the $\mathcal{S}$-matrix as a correlator on a sphere? In the book "Lectures on the Infrared Structure of Gravity and Gauge Theories" by Andrew Strominger, the author discusses in Chapter 3 the idea of "The $\mathcal{S}$-matrix as a Celestial Correlator".
He says:

The alternate description of scattering is as a type of correlation function on a sphere, depicted in figure 16. In the massless case, "in" and "out" massless particles are labelled by operators denoted $$\mathcal{O}_k(z,\bar{z}).\tag{3.0.1}$$ Here, $$z = \frac{x^1+ i x^2}{r+x^3}\tag{3.0.2}$$ denotes the point on the sphere at past or future null infinity where the particle of type $k$ enters or exits spacetime.

After a brief discussion about the massive case, which doesn't affect my doubt, the author continues:

It is then natural to express the $n$-particle scattering amplitudes in the form of a celestial correlator on $\mathcal{CS}^2$: $$\langle \text{out}|\mathcal{S}|\text{in}\rangle \to \langle \mathcal{O}_1(z_1,\bar{z}_1),\dots, \mathcal{O}_n(z_n,\bar{z}_n)\tag{3.0.3}\rangle$$ Note that no assumptions are being made here, we are simply rewriting the $\mathcal{S}$-matrix in a different notation.

Now my question is on making sense of Eq. (3.0.3) precisely.
What really are the $\mathcal{O}_k$ operators? My intuition tells me that they are just the creation/annihilation operators $a(\mathbf{p}),a^\dagger(\mathbf{p})$ expressed in terms of the asymptotic fields at $\mathcal{I}^\pm$ and with $\mathbf{p}$ written in terms of $(z,\bar{z})$. Is that the case?
Also what is the state with respect to which the mean value in Eq. (3.0.3) is taken? Is it the same vacuum on the asymptotic Fock space?
In summary, what is the precise way to translate the LHS of Eq. (3.0.3) to the RHS?
 A: The following answer is based on my interpretation of the article you cited, and on some of its references.
"Motivation" for the LHS $\rightarrow$ RHS step
What's the point of re-writing $\mathcal{S}$ into $\mathcal{O}_k$? It's probably just to break down a potentially complex $\mathcal{S}$ into "manageable" single-particle/excitation operators $\mathcal{O}_k$. Rendering the calculation easier. And more importantly disentangling an $n-$body problem into $n$ $1-$body problems.
This is similar to working out the mathematical expression of perturbative amplitudes in terms of simple(r) step-by-step Feynman rules for each line & vertex.
What are the $\mathcal{O}_k$ operators?
For this paper (one of the references in yours actually) I found the following half-decent definition:

The $n$-particle scattering amplitudes $\mathcal{A}_n$ of any
  four-dimensional quantum field thoery (QFT$_4$) can be described as a
  collection of $n$-point correlation functions on the two-sphere
  ($S^2$) with coordinates $(z,\bar{z})$: 
  $$ \mathcal{A}_n = \langle O_1(E_1, z_1, \bar{z}_1) \dots O_n(E_n, z_n, \bar{z}_n) \rangle,$$
  where $O_k$ creates (if $E_k<0$) or annihilates (if $E_k>0$) an
  asymptotic particle with energy $|E_k|$ at the point $(z_k,
> \bar{z}_k)$ where the particle crosses the asymptotic $S^2$ at null
  infinity ($\mathcal{I}$).

Based on the use of the word "asymptotic", it seems to me that your suggestion might be right. You are ignoring the specifics of the interactions, and just looking at $t\rightarrow -\infty$ when the particle was created ($\mathcal{O}_1(z_1, \bar{z}_1)$) and at $t\rightarrow \infty$ when the particle was destroyed ($\mathcal{O}_1(z_2, \bar{z}_2)$). I used the same subscript $1$ for the same praticle, though the creation and destruction occur at distinct angles $(z, \bar{z})$.
What is the state wrt which the mean is taken?
Because they are talking about "asymptotic" particles, and because they want to make $n$ $1-$body problems out of an $n-$body one, the state is the vacuum.
I.e., you start with nothing and create your single particle $1$ from it. Then you destroy so it returns to nothing.
You can use the non-interacting vacuum $|0\rangle$, but since it talks about non-Abelian stuff immediately afterwards you might want to use the interacting vacuum $|\Omega\rangle$.
A: Understanding the ${\cal S}$-matrix as a celestial correlator can be done in two different levels. In the first level we actually just introduce some suggestive notation. The so-defined correlator, however, still won't transform as a correlation function of quasi-primary operators (after all, where are the scaling dimensions?). This leads to the second level, where we actually consider a new basis for the asymptotic states in which indeed the ${\cal S}$-matrix behaves as a conformal correlator. It is then that one boldly postulates that there exists a two-dimensional conformal field theory whose correlation functions reproduce the ${\cal S}$-matrix in such a basis. This is where we start to get into the subject known as celestial holography. But let us explain a little more the details.
Let us recall that the starting point of the LSZ procedure to evaluate the ${\cal S}$-matrix is to write it in terms of Heisenberg picture creation/annihilation operators
\begin{eqnarray}
 \langle\text{out}|\text{in}\rangle &=& \langle 0;{\text out}|a_{\rm out}(p_1')\cdots a_{\rm out}(p_n')a^\dagger_{\rm in}(p_1)\cdots a^\dagger_{\rm in}(p_m)|0;\text{in}\rangle,
\end{eqnarray}
where, in distinction to standard QFT textbooks, we have allowed for the possibility of vacuum degeneracy, since this is in fact something clearly established at this point in the lectures under consideration. The operators $a_{\rm out}(p_i')$ and $a_{\rm in}^\dagger(p_j)$ can all be written in terms of in/out fields using the Klein-Gordon inner product in the asymptotic regions of spacetime. Equivalently one may introduce the in/out fields which have the standard plane wave expansion like free fields, and are such that the full interacting fields asymptote to them in the in/out regions of spacetime. This is done, e.g. for the photon, in equation (2.8.16) of the lectures:
\begin{eqnarray}
    A_z^{(0)}(u,z,\bar z) &=& -\dfrac{i}{8\pi^2}\dfrac{\sqrt{2}e}{1+z\bar z}\int_0^\infty d\omega\left[a_{\rm out}^+(\omega\hat{x}(z,\bar z))e^{-i\omega u}-a_{\rm out}^{-}(\omega\hat{x}(z,\bar z))^\dagger e^{i\omega u}\right] \tag{2.8.16},
\end{eqnarray}
where the upperscripts $\pm$ indicate the photon helicities. In that case, either by using the KG inner product, or by constructing an asymptotic expansion like the above and inverting it to extract the creation and annihilation operators, we observe that they will be extracted from bulk fields with the direction of the momentum vectors corresponding to points on the celestial sphere. In that case, it is natural to view $a_{\rm out}(p_i')$ and $a_{\rm in}^\dagger(p_j)$ as operators ${\cal O}^\pm_\omega(z,\bar z)$ defined on the celestial sphere, where $\pm$ picks out or in. We are just introducing a new suggestive notation for $a_{\rm out}(p_i')$ and $a_{\rm in}^\dagger(p_j)$ at this level: they are the same operators. In that sense, the ${\cal S}$-matrix is indeed
\begin{eqnarray}
 \langle\text{out}|\text{in}\rangle &=& \langle 0;{\text {out}}|{\cal O}_{\omega_1'}^+(z_1',\bar z_1')\cdots {\cal O}_{\omega_n'}^+(z_n',\bar z_n'){\cal O}_{\omega_1}^-(z_1,\bar z_1)\cdots {\cal O}_{\omega_m}^-(z_m,\bar z_m)|0;\text{in}\rangle.
\end{eqnarray}
Now, the motivation behind all of this is clearly a holographic perspective to gravity in asymptotically flat spacetimes. As shown in the lectures before getting to this point, quantum field theory soft theorems are really Ward identities of asymptotic symmetries, and when written in this language, they pretty much resemble Ward identities for symmetries in a two-dimensional conformal field theory living on the celestial sphere. In particular, quite remarkably, the subleading soft graviton theorem gives rise to an identity obeyed by the ${\cal S}$-matrix which looks like the Ward identity of a stress-tensor in a two-dimensional CFT (see section 2.6 of https://arxiv.org/abs/2107.02075 for a nice review). The only puzzling thing is that the conformal weights $(h,\bar h)$ are operators. Therefore this change of notation is not the full story yet behind the ${\cal S}$-matrix as a celestial correlator. The idea now is that if one changes the asymptotic basis so as to diagonalize these $(h,\bar h)$ operators one would have the standard Ward identity of a stress tensor.
A companion perspective is that Lorentz ${\rm SL}(2,\mathbb{C})$ transformations act on the celestial sphere as global conformal transformations, and one hope would be that the ${\cal S}$-matrix would transform as a conformal correlator under its action. But still, if the ${\cal S}$-matrix elements are to be conformal correlators, something is still missing. The candidate operators ${\cal O}_{\omega}^\pm(z,\bar z)$ may be local operators on the celestial sphere, but they depend on this parameter $\omega$ which is acted upon by Lorentz. Moreover, where are the scaling dimensions?
The answer to all of these questions is the introduction of the so-called conformal primary bases. The full details are nicely reviewed in https://arxiv.org/abs/2107.02075, but focusing on massless particles, one defines the states in the one particle Hilbert space
\begin{eqnarray}
 |\Delta,z,\bar z\rangle &=& \int_0^\infty d\omega \omega^{\Delta-1} |\omega,z,\bar z\rangle,
\end{eqnarray}
where $|\omega,z,\bar z\rangle$ are the standard momentum eigenstates with momentum $p = \omega \hat{q}(z,\bar z)$ where $\hat{q}(z,\bar z)$ is the map from points on the celestial sphere to null directions $$\hat{q}(z,\bar z)=(1+z\bar z,z+\bar z,-i(z-\bar z),1-z\bar z).$$ The ${\cal S}$-matrix in this basis can again be computed as a correlator of celestial operators ${\cal O}_\Delta^\pm(z,\bar z)$. But now, it turns out that by definition (see the aforementioned reference to understand the details of why this is true), it transforms indeed under Lorentz as a two-dimensional conformal correlator of quasi-primary operators. The candidate stress tensor Ward identity now takes the form of an actual Ward identity with numeric $(h,\bar h)$ and moreover one may show that these candidate operators are actually primaries with respect to said stress tensor, see https://arxiv.org/abs/1906.10149.
The motivation behind all of this is certainly holography. Now one boldly proposes that there exists a two-dimensional celestial conformal field theory whose correlation functions give the ${\cal S}$-matrix in this basis. One would then hope that if one finds such a theory one may be able to non-perturbatively define quantum gravity in asymptotically flat spacetimes. This is an active field of research known as celestial holography today and up to this point many properties such theories should obey have been found by translating properties known from quantum field theories. But it has also been able to give insights on the bulk already! For example the recently discovered $w_{1+\infty}$ symmetry of gravity (https://arxiv.org/abs/2105.14346) and its identification in a study of Einstein's equations (https://arxiv.org/abs/2112.15573).
So answering the question: equation (3.0.3) is really just the change of notation mentioned here. The operators are just the creation/annihilation operators written in a new notation and everything is happening in the Fock space from QFT. But this is not the end of the story: the really make sense of the ${\cal S}$-matrix as a celestial correlator one must change the basis to a conformal primary basis. Then one might imagine that there exists such a celestial CFT and that one is computing correlators in it.
