The proof is relatively long, so I focus on displaying the reasoning and major steps.

We work on a Black-Scholes model. Without loss of generality, we focus on an option with strike $P$ to buy at $t_e$ a European call option expiring at $T$, written on a stock $S$. Expectations are always taken with respect to the risk-neutral measure $Q$ unless otherwise stated, and we write $E_t(\cdot):=E(\cdot|\mathscr{F}_t)$.  

The value $C_t$ of a European call option is equal to:
$$C_t=E_t\left(e^{-r(T-t)}\max\{S_T-K\}\right)$$

The value $O_t$ of an option over a European call option is given by:
$$\begin{align}
O_t
&=E_t\left(e^{-r(t_e-t)}\max\{C_{t_e}-P,0\}\right)
\\&=e^{-r(t_e-t)}E_t\left(
\max\left\{E_{t_e}\left(e^{-r(T-t_e)}\max\{S_T-K,0\}\right)-P,0\right\}\right)
\\&=e^{-r(T-t)}E_t
\left(1_{S_{t_e}\geq S^\star}\left(E_{t_e}\left(1_{S_T\geq K}\left(S_T-K\right)\right)-\tilde{P}\right)\right)
\\&=e^{-r(T-t)}E_t\left(E_{t_e}
\left(1_{S_{t_e}\geq S^\star}1_{S_T\geq K}\left(S_T-K\right)-1_{S_{t_e}\geq S^\star}\tilde{P}\right)\right)
\\\tag{1}&=e^{-r(T-t)}\left(
E_t\left(1_{S_{t_e}\geq S^\star, S_T\geq K}S_T\right)
-E_t\left(1_{S_{t_e}\geq S^\star, S_T\geq K}\right)K
-E_t\left(1_{S_{t_e}\geq S^\star}\right)\tilde{P}
\right)
\end{align}$$
where the last inequality stems from the law of iterated expectations, $\tilde{P}:=e^{r(T-t_e)}P$ is the compounded strike, and $S^\star$ is the value of $s$ that solves the following equation:
$$\tag{2}c(s,T-t_e)-\tilde{P}=0$$
where $c$ is the undiscounted Black-Scholes price for a European call option:
$$c(s,\tau):=se^{r\tau}\Phi\left(\frac{\ln\frac{s}{K}+\left(r+\frac{\sigma^2}{2}\right)\tau}{\sigma\sqrt{\tau}}\right)-K\Phi\left(\frac{\ln\frac{s}{K}+\left(r-\frac{\sigma^2}{2}\right)\tau}{\sigma\sqrt{\tau}}\right)$$

In practice, the quantity $S^\star$ can be computed by solving numerically Equation $(2)$, for example by Newton-Raphson. Now, the third term in Equation $(1)$ is merely the probability of a log-normal variable being above $S^\star$. By analogy with classical Black-Scholes, this is equal to:
$$\tag{3.a}E_t\left(1_{S_{t_e}\geq S^\star}\right)=\Phi\left(\frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}}\right)$$

To compute the second term in $(1)$, we define $Z$ and $Y$ as two independent normal random variables with zero mean and unit variance in order to represent the Brownian increments from $t$ to $t_e$ and from there to $T$. Note that:
$$\begin{align}
1_{\{S_{t_e}\geq S^\star,\ S_T\geq K\}}
&=1_{\left\{S_te^{\left(r-\frac{\sigma^2}{2}\right)(t_e-t)+\sigma \sqrt{t_e-t}Z}\geq S^\star,\
S_te^{\left(r-\frac{\sigma^2}{2}\right)(T-t)+\sigma\left(\sqrt{t_e-t}Z+\sqrt{T-t_e}Y\right)}\geq K\right\}}
\\\tag{4}
&=1_{\left\{Z\leq \frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\
X\leq \frac{\ln\frac{S_t}{K}+\left(r-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}}
\end{align}$$
where $X$ is a third standard normal variable with the following correlation with $Z$:
$$\rho:=\frac{\text{Cov}(\sqrt{t_e-t}Z+\sqrt{T-t_e}Y,Z)}{\sqrt{V(\sqrt{T-t_e}Z+\sqrt{t_e-t}Y)V(Z)}}=\sqrt{\frac{t_e-t}{T-t}}$$ 
Hence the second term in Equation $(1)$ is the cumulative bi-variate normal probability over $Z$ and $X$ parameterised by their correlation $\rho$:
$$\begin{align}
&E_t\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}\right)
\\[6pt]
&\quad=\Phi_\rho\left(\frac{\ln\frac{S_t}{S^\star}+\left(r-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\frac{\ln\frac{S_t}{K}+\left(r-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right)
\end{align}$$

For the first term in Equation $(1)$, we change the measure of the expectation to the stock measure where the numéraire is the stock $S$, which we will write $\mathcal{S}$:
$$\begin{align}
E_t^\mathcal{Q}\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right)
&=E_t^\mathcal{S}\left(e^{r(T-t)}\frac{S_t}{S_T}1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right)
\\
&=e^{r(T-t)}S_tE_t^\mathcal{S}\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}\right)
\end{align}$$

The Radon-Nikodym process implied by this change of measure is:
$$\begin{align}
\left.\frac{d\mathcal{Q}}{d\mathcal{S}}\right|_{\mathscr{F}_t}
&=e^{r(T-t)}\frac{S_t}{S_T}
\\
&=e^{\frac{\sigma^2}{2}(T-t)-\sigma W^\mathcal{Q}_t}
\end{align}$$
The "Novikov process" for this change of measure is therefore $\theta_t:=-\sigma t$. Thus the following process is the Brownian Motion under the stock measure:
$$W^\mathcal{S}_t=W^\mathcal{Q}_t-\sigma t$$

This implies that the drift of the stock under this new measure increases by $\sigma^2$. We can leverage Equation $(4)$ but this time on standard normal variables $Z'$, $Y'$ and $X'$ under the stock measure $\mathcal{S}$:
$$\begin{align}
1_{\{S_{t_e}\geq S^\star,\ S_T\geq K\}}
&=1_{\left\{Z'\leq \frac{\ln\frac{S_t}{S^\star}+\left((r+\sigma^2)-\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\
X'\leq \frac{\ln\frac{S_t}{K}+\left((r+\sigma^2)-\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}}
\\
&=1_{\left\{Z'\leq \frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\
X'\leq \frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right\}}
\end{align}$$
That is:
$$\begin{align}
&E_t\left(1_{S_{t_e}\geq S^\star,\ S_T\geq K}S_T\right) =
\\[6pt]
\tag{3.c}
& \quad e^{r(T-t)}S_t \Phi_\rho\left(\frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}},\frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}\right)
\end{align}$$
Defining:
$$\begin{align}
d_1 & := \frac{\ln\frac{S_t}{K}+\left(r+\frac{\sigma^2}{2}\right)(T-t)}{\sigma\sqrt{T-t}}
\\
d_1^\star & := \frac{\ln\frac{S_t}{S^\star}+\left(r+\frac{\sigma^2}{2}\right)(t_e-t)}{\sigma\sqrt{t_e-t}}
\\[12pt]
d_2 & := d_1-\sigma\sqrt{T-t}
\\[18pt]
d_2^\star & := d_1^\star-\sigma\sqrt{T-t}
\end{align}$$
Then combining Equations $(1)$, $(3.a)$, $(3.b)$ and $(3.c)$, we obtain the desired result:
$$O_t=S_t\Phi_\rho\left(d_1^\star,d_1\right)-e^{-r(T-t)} K\Phi_\rho\left(d_2^\star,d_2\right)-e^{-r(t_e-t)}P\Phi(d_2^\star) \quad \square$$

As you can see, it is very similar to the Black-Scholes Equation for a call option:

- The third term of the valuation formula for $O_t$ is equivalent to the second term in the Black-Scholes formula, namely the discounted strike $P$ times the probability that the option on the option will be exercised;
- The first two terms combined look very close to the Black-Scholes value of a call option: this is to be expected, because the option is written on a call option. However, the probabilities $\Phi_\rho(d_1^\star,d_1)$ and $\Phi_\rho(d_2^\star,d_2)$ account for the value of the stock price at $t_e$ and $T$. This is because when you enter into a vanilla call option, you know the stock price at trade date $t$ but not at expiry $T$. In this case, you will enter the option at a future time $t_e>t$, hence you do not know the value of the underlying neither at inception $t_e$ nor at expiry $T$: this heightened uncertainty is captured by the bi-variate normal distribution.