Drive frequency for second order quantum transitions Summary
$
\newcommand{\ket}[1]{\left \lvert #1 \right \rangle}
\newcommand{\bra}[1]{\left \langle #1 \right \rvert}
\newcommand{\braket}[2]{\left \langle #1 | #2 \right \rangle}
\newcommand{\bbraket}[3]{\left \langle #1 | #2 | #3 \right \rangle}
%
$How do we show, from first principles, that the drive frequency $\omega_d$ needed to excite a two-photon transition from state $\ket{i}$ to state $\ket{f}$ is $\omega_d = (\omega_f - \omega_i) / 2$?
Setup
Consider a Hamiltonian
$$H(t) = H_0 + V(t)$$
where the time dependent $H_0$ is considered "unperturbed" and $V(t)$ is considered a perturbation.
Let us further suppose that
$$V(t) = V_0 \, v(t) \, \mathcal{O} \, $$
where $V_0$ captures the magnitude of the perturbation, $v(t)$ captures the time dependence, and $\mathcal{O}$ is a time-indepdendent operator.
It's often useful to work in the interaction picture.
In that picture, the Schrodinger equation becomes
$$ i \hbar \partial_t \ket{\Psi'(t)} = V'(t) \ket{\Psi'(t)}$$
where $V'(t) \equiv U^{-1}(t) V(t) U(t)$, $U(t)$ is the propagator associated with $H_0$, and primes on the states indicates that the interaction picture states are related to the Schrodinger picture ones via $\ket{n'} = U^{-1}(t)\ket{n(t)}$.
Dyson series
We can solve the Schrodinger equation formally as
$$\ket{\Psi'(t)} = \ket{\Psi'(0)} -\frac{i}{\hbar}\int_0^t dt' V'(t') \ket{\Psi'(t')} \, .$$
Plugging this equation into itself gives
$$
\ket{\Psi'(t)} =
\ket{\Psi'(0)}
-\frac{i}{\hbar} \int_0^t dt' \, V'(t')
\left[
  \ket{\Psi'(0)} - \frac{i}{\hbar} \int_0^{t'} dt'' \, V'(t'') \ket{\Psi'(t'')}
\right] \, .
$$
So far, this equation is exact.
Now suppose we iterate again, but this time make an approximation by keeping only the term where $\ket{\Psi'(t'')} = \ket{\Psi'(0)}$.
Let's also hit the whole thing with $\bra{f}$ to compute the transition amplitude to $\ket{f}$.
The result is
$$
\braket{f}{\Psi'(t)} \approx
\underbrace{\langle f \ket{\Psi'(0)}}_{0^\text{th}\text{ order}}
- \underbrace{
  \frac{i}{\hbar} \int_0^t dt' \, \langle f | V'(t') \ket{\Psi'(0)}
}_{1^\text{st} \text{ order}}
- \underbrace{\frac{1}{\hbar^2} \int_0^t \int_0^{t'} dt' \, dt'' \,
  \langle f | V'(t') V'(t'') \ket{\Psi'(0)}
}_{2^\text{nd} \text{ order}} \, .
$$
Using the following facts


*

*$V'(t) = U^{-1}(t)V(t)U(t)$

*For an eigenstate $\ket{n}$ of $H_0$ we have $U(t) \ket{n} = \exp(-i \omega_n t) \ket{n}$

*$\sum_n \ket{n} \bra{n} = \text{identity}$
we can reduce the second order term to
\begin{align}
\braket{f}{\Psi'(t)} = &
  - \left( \frac{V_0}{\hbar} \right)^2 \sum_n \int_0^t \int_0^{t'} dt' \, dt'' \\
  & v(t') v(t'') \bbraket{f}{\mathcal{O}}{n} \bbraket{n}{\mathcal{O}}{i} \\
  & \exp \left[
    i \left(
      \omega_f t' - \omega_n t' + \omega_n t'' - \omega_i t''
    \right)
  \right]
\end{align}
where we've used a shortcut notation $\ket{i} \equiv \ket{\Psi'(0)}$.
Question
I've heard that if $v(t)$ is a sinusoidal drive, e.g. $v(t) = \cos(\omega_d t)$, then the second order term we've computed has large magnitude when $\omega_d = (\omega_f - \omega_i)/2$.
How can we show this?
Is there some way we can do either the sum or the integrals in some meaningful limits, e.g. when $t \rightarrow \infty$?
 A: $
\newcommand{\ket}[1]{\left \lvert #1 \right \rangle}
\newcommand{\bra}[1]{\left \langle #1 \right \rvert}
\newcommand{\braket}[2]{\left \langle #1 | #2 \right \rangle}
\newcommand{\bbraket}[3]{\left \langle #1 | #2 | #3 \right \rangle}
%
$Your final expression factorizes quite nicely as
\begin{align}
\braket{f}{\Psi'(t)} = &
  - \left( \frac{V_0}{\hbar} \right)^2 
\sum_n \bbraket{f}{\mathcal{O}}{n} \bbraket{n}{\mathcal{O}}{i}
\\ & \qquad \times
\int_0^t  \mathrm dt'
v(t') \exp\left(i(\omega_f-\omega_n)t'\right)
\\ & \qquad \quad \times
\int_0^{t'} \mathrm dt''  
v(t'')  \exp\left(i(\omega_n-\omega_i)t''\right)
,
\end{align}
and here (i) the first integral is essentially trivial, yielding some exponentials times a sinc function, while (ii) the second one is slightly more involved but still yields immediately to symbolic integration in Mathematica.
I don't see much point in repeating the full result here, but the crucial aspect is that if you fix the system frequencies ($\omega_i$, $\omega_f$ and $\omega_n$) and you scan over the pump frequency $\omega_p$, then you will get some form of bump around $(\omega_i+\omega_f)/2$ and this bump will become sharper and sharper in the $t\gg 2\pi/\omega$ limit:

(In this sample, $\omega_i=0$ and $\omega_f=1$, without loss of generality, and then $\omega_p=0.55+0.001i$, with the imaginary part added just to regularize the calculation. The sidebands are due to the precise location of $\mathrm{Re}(\omega_p)$, but the bright center at $\omega_p=0.5$ is the two-photon resonance process.)
If you want a more enlightening analytical result, it's better to drop the counter-rotating terms and just use an exponential form for the potential, $v(t) = e^{-i\omega_p t}$, under which the integrals will simplify significantly to
$$
\frac{1}{\omega_i-\omega_n+\omega_p}
\left[
\frac{e^{i(\omega_f-\omega_i-2\omega_p)t}-1}{\omega_f-\omega_i-2\omega_p}
-
\frac{e^{i(\omega_f-\omega_n-\omega_p)t}-1}{\omega_f-\omega_n-\omega_p}
\right]
,
$$
which contains an explicit term in $\operatorname{sinc}((\omega_f-\omega_i-2\omega_p)t/2)$; this will converge at $t\to\infty$ to a Dirac delta on precisely the energy-conservation condition $2\omega_p=\omega_f-\omega_i$.

Finally, if you have multiple intermediate levels $|n⟩$ contributing to the process, then the total amplitude will be the interference of all the individual amplitudes, but the energy-conserving term in $\operatorname{sinc}((\omega_f-\omega_i-2\omega_p)t/2)$ will factor out, and the process you're interested in will have an amplitude that reads
\begin{align}
\braket{f}{\Psi'(t)} \approx &
  - \left( \frac{V_0}{\hbar} \right)^2 
\frac{e^{i(\omega_f-\omega_i-2\omega_p)t}-1}{\omega_f-\omega_i-2\omega_p}
\sum_n \frac{\bbraket{f}{\mathcal{O}}{n} \bbraket{n}{\mathcal{O}}{i}}{\omega_i-\omega_n+\omega_p}
,
\end{align}
in a form that should be pretty recognizable from multiple other places in perturbation theory. Here you get a sum of two-photon matrix elements over various intermediate states, divided by their detuning, with a global sinc-like energy-conservation factor multiplying the whole sum.
A: Emilio has already put a nice answer (it makes me feel that in my comment I made the problem sound much more complicated then it needs to be!). When I read the question, I understood it as asking, when you are driving a two-photon transition, if it is better to use $\omega_d \approx \frac{\omega_f - \omega_i}{2}$ than use other combinations such as, e.g.  $\omega_{d_1} \approx \frac{\omega_f - \omega_i}{8}$ and $\omega_{d_2} \approx \frac{3(\omega_f-\omega_i)}{8}$, and I don't think Emilio addressed this point explicitly so I would like to add my explanation. I don't think much of what I write below (borrowing heavily from the book "Atom-Photon Interactions") will be different from what Emilio has already said, but maybe the presentation looks more formal. 

(Starting from page 28, equation (21) from API), we start with the expression for the second-order contribution to the transition amplitude $\mathcal{P}_{fi}^{(2)}$:

$\mathcal{P}_{fi}^{(2)} = \left(\frac{1}{i\hbar}\right)^2 \underset{-T/2 \leq \tau_1 \leq \tau_2 \leq T/2}{\int d\tau_1 d\tau_2} \sum_{k} V_{fk}V_{ki}\,\, e^{i(E_f - E_k) \tau_2 
 / \hbar} \, e^{i(E_k - E_i)\tau_1 / \hbar}$

The sum represents sum over intermediate states. To eliminate the ordering $\tau_1 \leq \tau_2$, we use the integral representation of Heaviside theta function:

$e^{-i E_k (\tau_2 - \tau_1) /\hbar} \theta (\tau_2 - \tau_1) = \underset{\eta \to 0_{+}}{\lim} -\frac{1}{2\pi i} \int_{-\infty}^{\infty} \frac{e^{-i E (\tau_2 - \tau_1) / \hbar} }{E + i\eta - E_k} dE$

Using this expression, we can re-write $\mathcal{P}_{fi}^{(2)}$ as

$\mathcal{P}_{fi}^{(2)} = \left(\frac{1}{i\hbar}\right)^2 \left(\frac{-1}{2\pi i}\right) \int_{-T/2}^{T/2} d\tau_2 \int_{-T/2}^{T/2} d\tau_1 \int_{-\infty}^{\infty} dE \,\, e^{i(E_f - E)\tau_2 / \hbar} e^{-i(E_i - E)\tau_1 / \hbar} \underset{\eta \to 0_{+}}{\lim} \sum_{k} \frac{V_{fk}V_{ki}}{E+i\eta - E_k}$

Now we need to generalize this to the situation when the interaction $\hat{V}$ has some time-dependence. The matrix element $V_{fk}$ will be written as $V_{fk}(\tau_2)$ and $V_{ki}$ as $V_{ki}(\tau_1)$. In the problem we consider, the explicit time dependence comes from the envelop $v(\tau)$. So $V_{fk}(\tau_1) = V_{fk}v(\tau_1)$ and so on. 
For a sinusoidal drive $v(\tau) = A_1\sin (\omega_{d_1} \tau) + A_2 \sin (\omega_{d_2} \tau)$ we obtain sinc functions after we carry out the time integrals. As Emilio said, we can ignore the counter-rotating terms ($e^{+i\omega\tau}$). And I will also focus on contribution from both $\omega_{d_1}$ and $\omega_{d_2}$, assuming that both of them and their harmonics are detuned from resonance. With this approximation we get

$\left(\frac{2 \pi \hbar}{i\hbar}\right)^2 \left(\frac{-1}{2\pi i} \right) \int dE \,\delta^{(T)} (E-E_i - \hbar\omega_{d_1}) \delta^{(T)}(E+\hbar\omega_{d_2}-E_f)\underset{\eta \to 0_{+}}{\lim} \sum_{k} A_1 A_2\frac{V_{fk}V_{ki}}{E+i\eta - E_k} $

where the function $\delta^{(T)} (E_a - E_b)$ being

$\delta^{(T)} (E_a - E_b) = \frac{1}{2\pi \hbar} \int_{-T/2}^{T/2} d\tau \, e^{i(E_a - E_b) \tau / \hbar} = \frac{1}{\pi} \frac{\sin(E_f - E_i) T/ 2 \hbar}{E_f - E_i}$

obviously $\delta^{(T)}$ approximates a Dirac delta function at $T \to \infty$ and hence the notation.
Here we take a break from re-writing expressions and make some physical interpretations...The strength of the transition will be determined by the magnitude of the matrix elements and the energy widths of $W_{fi}(E) = \underset{\eta \to 0_{+}}{\lim} \sum_{k} A_1 A_2\frac{V_{fk}V_{ki}}{E+i\eta - E_k}$ and $\delta^{(T)}(E)$'s. Let's assume that there is no "special" intermediate state with large matrix element, so that we focus on the effects of the energy widths of those aforementioned functions. The $\delta^{(T)}(E)$'s are practically zero when $E > \hbar / T$. So we see that we get major contribution  to the energy integral when $E$ lies in the range such that $E_i + \hbar\omega_{d_1} - E $ and $E_f-\hbar\omega_{d_2}-E$ are within $\hbar / T$ and $W_{fi} (E)$ varies slowly in the same range. We see that the more dissimilar $\omega_{d_1}$ and $\omega_{d_2}$ are, the harder it is to satisfy the overlap problem. 
To complete the analysis, once you assume $W_{fi}(E)$ is slowly varying, you approximate it as $W_{fi}(E_m)$, where $E_m$ gives you the strongest overlap. If $\omega_{d_1} = \omega_{d_2}$, then $E_m$ is simply $\hbar\omega_d$. The integral of the product of delta function gives you another delta function that enforces the conservation of energy for two  -photon transition (basically reduces to what Emilio wrote at the end). 
