Quoting the definition on Wikipedia,
Let $(\Omega_1, \Sigma_1, \mu_1)$ and $(\Omega_2, \Sigma_2, \mu_2)$ be $\sigma$-finite measure spaces. Suppose $1 \leq p_0 \leq p_1 \leq \infty$, $1 \leq q_0 \leq q_1 \leq \infty$, and let $T : L^{p_0}(\mu_1) + L^{p_1}(\mu_2) \to L^{q_0}(\mu_2) + L^{q_1}(\mu_2)$ be a linear operator that maps $L^{p_0}(\mu_1)$ (resp. $L^{p_1}(\mu_1)$) boundedly into $L^{q_0}(\mu_2)$ (resp. $L^{q_1}(\mu_2)$). For $0 < \theta < 1$, let $p_\theta$, $q_\theta$ be defined as above. Then $T$ maps $L^{p_\theta}(\mu_1)$ boundedly into $L^{q_\theta}(\mu_2)$ and satisifies the operator norm estimate
$$ \|T\|_{L^{p_\theta} \to L^{q_\theta}} \leq \|T\|^{1-\theta}_{L^{p_0} \to L^{q_0}} \|T\|^\theta_{L^{p_1} \to L^{q_1}}.$$
I am having some difficulty visualizing the statement of this theorem, and therefore do not understand why this is an "interpolation" theorem.
What is the intuition behind this theorem? What is being interpolated here?
Riesz-Thorin interpolation (RTI) bounds the norms of linear maps acting between $L^p$ spaces.
Unlike Marcinkiewicz interpolation, RTI only works for strong type operators.
Overview of interpolation: to understand how to obtain control on the expression $||Tf||_{L^q}$ for operator $T$ & function $f,$ one would divide $f$ into two (or more) components (eg, into one with large $f$ & small $f,$ or where $f$ is oscillating with high frequency or only varying with low frequency). Each component would be estimated using a chosen combination of the extreme estimates available. Finally optimising over these choices & summing up we hope to get a good estimate on the original expression.
The general intuition behind interpolation is as follows: Let’s say we have two upper bound estimates of the form
\begin{equation*} X_0\leq Y_0,~X_1\leq Y_1. \end{equation*}
We are looking to establish boundedness estiamtes such as
\begin{equation*} ||Tf||_{L^q}\leq C||f||_{L^p},~\frac{1}{p}+\frac{1}{q}=1. \end{equation*}
We want to reduce the task of proving such estimates to that of proving various endpoint versions of these estimates.
So we conclude a family of intermediate estimates $X_{\theta}\leq Y_{\theta},~0<\theta<1.$
To get the formula in your theorem, let $p_0,p_1\in\mathbb{R}$ such that $0<p_0<p_1\leq \infty.$ Then for $0<\theta<1$ define $p_{\theta}$ by
\begin{equation*} \frac{1}{p_{\theta}}=1-\frac{\theta}{p_0}+\frac{\theta}{p_1} \end{equation*}
Now split the function $f$ in the $L^{p_{\theta}}$ space as the product
\begin{equation*} |f|=|f|^{1-\theta}|f|^{\theta}. \end{equation*}
Now we increase its integrability via Holder’s inequality (to the $p_{\theta}$ power) to get
\begin{equation*} ||f||_{p_{\theta}}\leq ||f||^{1-\theta}_{p_0}||f||^{\theta}_{p_1},~f\in L^{p_0}\cap L^{p_1}. \end{equation*}
I’m happy to expand on any points, or go into interpolation in more detail if you like. Hope this helps!