What's the "effective potential" for photons in $X$-ray diffraction? The slickest way to introduce $X$-ray diffraction is to invoke scattering theory in quantum mechanics. One treats the incoming photon as just another particle in a scattering problem; by Fermi's golden rule and the Born approximation, the scattering rate from $\mathbf{k}$ to $\mathbf{k}'$ is
$$\Gamma(\mathbf{k}', \mathbf{k}) \propto |\langle \mathbf{k}'| V | \mathbf{k} \rangle|^2$$
where $V$ is the potential experienced by the photon. Since $V$ has the periodicity of the lattice, it follows that the scattering rate vanishes unless $\mathbf{k} - \mathbf{k}'$ is a reciprocal lattice vector.
However, most sources do not say much about this "photon potential" or justify its form. For example, David Tong's lecture notes simply avoid the issue:

Firing a beam of particles — whether neutrons, electrons or photons in the X-ray spectrum — at the solid reveals a characteristic diffraction pattern. [...] Our starting point is the standard asymptotic expression describing a wave scattering off a central potential.

Tong declines to comment on the potential at all. Steve Simon's solids textbook says a bit more:

If we think of the incoming wave as being a particle, then we should think of the sample as being some potential $V(r)$ that the particle experiences. [...] X-rays scatter from the electrons in a system. As a result, the scattering potential is proportional to the electron density.

That is, the potential is the electron density, and nuclei don't contribute because they are heavier.
This sounds plausible, but I have no idea how to derive this from first principles. From the QFT side, I imagine such a calculation would start from the QED interaction
$$\mathcal{L}_{\text{int}} = \bar{\psi} \gamma^\mu A_\mu \psi.$$
Treating the field $\psi$ as a classical, static background field we have
$$\mathcal{L}_{\text{int}} = j^\mu A_\mu = \rho A_0$$
since the components $j^i$ are zero. But I'm not sure how this is supposed to be a "potential for the photon"; it seems to only affect one component of $A_\mu$. And it's unclear where the mass of the fermion is going to come in, to make the electrons count more than the nuclei. 
On the other end, I suppose one could start from classical electromagnetism. In that case we're talking about Thomson scattering, and heavier particles indeed contribute less. The challenge is then exporting quantum mechanical ideas, such as partial waves and the Born approximation, to this classical context. Maybe this is manageable, but I've never seen this done anywhere either.
The derivation of X-ray diffraction peaks from quantum scattering theory is very slick, but how is it justified in detail?
 A: The comments by KF Gauss are on the right track: non-relativistic quantum electrodynamics (NRQED) is a good foundation. The question is then how to recover the "effective potential" description of X-ray diffraction from NRQED.
This isn't my specialty, but I reviewed the NRQED derivation of elastic X-ray scattering from a crystal, using the references listed below. At lowest order in the small-coupling expansion, for elastic scattering (ignoring resonances), the scattering cross-section can be written essentially as shown in the OP. An explicit expression for the effective potential $V$ is shown at the end of the "Result" section, below.
The rest of this answer outlines the derivation, starting from NRQED. This outline is based on the following sources:


*

*Sections 2, 3, and 5 in Rohringer, "Introduction to the theory of x-ray matter interaction", https://arxiv.org/abs/1811.12052

*Chapter 2 in Slowik's 2015 dissertation "Quantum effects in nonresonant x-ray scattering", https://pure.mpg.de/rest/items/item_2283270/component/file_2283267/content
For additional context, see 


*

*Bennett et al (2018), "Monitoring molecular nonadiabatic dynamics with femtosecond X-ray diffraction", PNAS vol 115 pages 6538–6547, https://www.pnas.org/content/115/26/6538

*Bennett's 2017 dissertation "Quantum-Field Effects in Ultrafast Nonlinear X-Ray and Optical Spectroscopy of Molecules", https://escholarship.org/content/qt98n675c7/qt98n675c7.pdf
 The model 
I'll start by specifying the model (NRQED). The ingredients are:


*

*The EM field, expressed in terms of a gauge field $\mathbf{A}$ in the Coulomb gauge.

*A non-relativistic fermion field $\psi$ for the electrons.

*An external potential to account for the lattice of nuclei, whose dynamics can be neglected.
I'll write $H_\chi$ for the the electron-electron and electron-nucleus Coulomb interaction terms in the Hamiltonian. These terms define the structure of the crystal (hence the subscript $\chi$), but we won't need to use them explicitly because we're trying to understand what the crystal does to the X-rays, not trying to understand the structure of the crystal itself. The total Hamiltonian is
$$
 H = H_{A} + H_{A\psi} + H_{\chi}
$$
with
\begin{align}
H_A &= \int d^3x \frac{\mathbf{E}^2+\mathbf{B}^2}{2} + \text{const}
\\
H_{A\psi} &= \int d^3x\ \left(\frac{(\mathbf{D}\psi)^\dagger
  \cdot (\mathbf{D}\psi)}{2m}
 +\frac{q}{2m}\mathbf{B}\cdot\psi^\dagger\sigma\psi\right)
\end{align}
and
\begin{align*}
 \mathbf{B} &= \nabla\times \mathbf{A} \\
 \mathbf{D} &= \nabla+iq\mathbf{A},
\end{align*}
with commutation relations
\begin{gather*}
 \big\{\psi_a(x),\psi_b^\dagger(y)\big\}
  =\delta_{ab}\delta^3(x-y)
\\
 \big[A_j(y),E_k(x)\big]=i\big(\delta_{jk}\delta^3(x-y)\big)^\perp
\end{gather*}
where $\perp$ indicates the transverse part (which would be $\delta_{jk}-\mathbf{k}_j\mathbf{k}_k/\mathbf{k}^2$ in the Fourier-transformed domain). In addition to the equations of motion implied by the Hamiltonian, we also have the Gauss-law constraint
$$
 \nabla\cdot\mathbf{E}\propto \psi^\dagger\psi
$$
that defines the longitudinal electric field operator in terms of the charge density operator.
 The quantity to be calculated 
The initial state consists of a crystal and an incoming photon. We can construct such a state by applying an appropriate combination of electron creation operators $\psi^\dagger$ to the vacuum state, and by also applying a photon creation operator with the desired wavenumber $\mathbf{k}$. I'll write this state as 
$$
|I\rangle = a^\dagger(\mathbf{k})|\chi\rangle,
$$
 where $|\chi\rangle$ denotes the state of the electrons in the crystal. To construct the final state, we can use 
$$
|F\rangle = a^\dagger(\mathbf{k}')|\chi\rangle,
$$
 with a different wavenumber $\mathbf{k}'$ for the photon. By using the same $|\chi\rangle$ for the initial and final state of the crystal, we're assuming elastic scattering. Momentum is not conserved in this model because the nuclei are represented by an external potential.
The goal is to calculate the transition amplitude $\langle F|I\rangle$. We can do this by working to second order in the electron charge $q$, except for the $H_\chi$ terms that will implicitly be included to all orders.
 Which terms are important? 
The electron-photon interaction terms are contained in $H_{A\psi}$, which can be written
$$
 H_{A\psi} = H_\psi + H_1 + H_2
$$
where $H_1$ is linear in $q\mathbf{A}$ and $H_2$ is quadratic in $q\mathbf{A}$. We'll treat $H_1+H_2$ as the interaction Hamiltonian $H_\text{int}$ in perturbation theory. Schematically, the transition amplitude to second order in $H_\text{int}$ is
\begin{align*}
 \langle F|I\rangle
 &\sim
 \delta(E_F-E_I)
 \left(
 \langle F|H_\text{int}|I\rangle
 \phantom{\frac{1}{2}}\right.
\\
 &+
 \left.
 \sum_M\frac{\langle F|H_\text{int}|M\rangle \langle M|H_\text{int}|I\rangle}{
  E_I-E_M}
 \right).
\end{align*}
This is an abbreviation of equation (25) in Rohringer's paper. We immediately see that the term of first order in $q$ does not contribute (because it would change the number of photons, but our initial and final states have the same number of photons), so the lowest non-zero order is $q^2$. After discarding higher-order terms, we get
\begin{align*}
 \langle F|I\rangle
 &\sim
 \delta(E_F-E_I)
 \left(
 \langle F|H_2|I\rangle
 \phantom{\frac{1}{2}}\right.
\\
 &+
 \left.\sum_M\frac{\langle F|H_1|M\rangle \langle M|H_1|I\rangle}{
  E_I-E_M}
 \right).
\end{align*}
The second term, which involves $H_1$, is expected to be important near a resonance, or at any photon-energy for which the denominator is close to zero. For elastic scattering with photon energies far away from any of the energies $E_M$, we can neglect the second term. (This seems plausible, but I'm trusting the literature here, because I haven't actually checked the relative magnitudes of these terms under realistic conditions.) 
 The result 
With the approximations outlined above, we are left with
\begin{align*}
 \langle F|I\rangle
 &\sim
 \delta(E_F-E_I)
 \langle F|H_2|I\rangle
\\
 &\sim
 \delta(E_F-E_I)
 \int d^3x\ 
 \langle \chi|\psi^\dagger(x)\psi(x)|\chi\rangle
\\
 &\phantom{\sim\delta(E_F-E_I)\int d^3x\ }
 \times
 \frac{1}{m}\langle 0|a(\mathbf{k}')\mathbf{A}^2(x)a^\dagger(\mathbf{k})|0\rangle
\\
 &\sim
 \delta(E_F-E_I)
 \frac{1}{m}\int d^3x\ 
 \rho(\mathbf{x})
 \,
 e^{i(\mathbf{k}'-\mathbf{k})\cdot\mathbf{x}}
 \times \text{polarization factor}
\end{align*}
where $\rho(\mathbf{x})$ is the crystal's electron density
$$
\rho(\mathbf{x})
 := \langle \chi|\psi^\dagger(\mathbf{x})\psi(\mathbf{x})|\chi\rangle.
$$
To relate this to the form shown in the OP, write it as
$$
 \langle F|I\rangle
 \sim
 \delta(E_F-E_I)
 \,
 \langle \mathbf{k}'|V |\mathbf{k}\rangle
$$
with
$$
 |\mathbf{k}\rangle := a^\dagger(\mathbf{k})|0\rangle
\hskip2cm
 V := \frac{1}{m}\int d^3x\ \rho(\mathbf{x})\mathbf{A}^2(x).
$$
This is the "effective potential" for photons in elastic diffraction. In words, it is the square of the gauge field operator (in the Coulomb gauge, so only the physical transverse components are involved) smeared against the electron density of the crystal.
 X-ray scattering from nuclei 
The model considered above treated the nuclei as non-dynamic sources of an external potential for the electrons. That model doesn't conserve momentum, and it doesn't allow for X-ray scattering from the nuclei. To fix those defects, we can add another term
$$
H_{A\varphi} = \int d^3x\ \frac{(\mathbf{D}\varphi)^\dagger
  \cdot (\mathbf{D}\varphi)}{2M}
$$
(plus a magnetic-moment term, not shown) where $\varphi^\dagger(x)$ and $\varphi(x)$ are the creation/annihilation operators for a given species of nucleus with mass $M$, and we can replace the external-potental term implicit in $H_\chi$ with an electron-nucleus Coulomb interaction term. The latter replacement restores momentum conservation, and adding $H_{A\varphi}$ allows for X-ray scattering from nuclei.
These changes affect the result in two ways:


*

*We can no longer use strictly the same crystal-state $|\chi\rangle$ in both the initial and final states, because now the crystal recoils when it scatters the photon. However, because the recoil is so slight, the initial and final crystal-states $|\chi_I\rangle$ and $|\chi_F\rangle$ will still be approximately proportional to each other (as vectors in the Hilbert space), so 
they might as well be equal.

*The "effective potential" $V$ now has two terms, namely
$$
 V := \int d^3x\ 
 \left(\frac{\rho(\mathbf{x})}{m}
 +\frac{\rho'(\mathbf{x})}{M}
 \right)\mathbf{A}^2(x)
$$
where $\rho'$ is the density of nuclei in the crystal, analogous to the electron density $\rho$ but involving the nucleus field operator $\varphi$ in place of $\psi$. Because of the masses in the denominoators, the nucleus term is suppressed by a factor $m/M< 1/2000$ compared to the electron term, so for most purposes we can neglect the contribution from the nuclei.
 Deriving NRQED from QED 
This is the really difficult part. If the EM field were just a background classical field, then it would be easier: Starting with the Dirac equation, use the first couple of orders of the Foldy-Wouthuysen transformation to get the non-relativistic Pauli equation. However, if we treated EM field as just a background classical field, then X-rays wouldn't scatter.
As soon as we promote the EM field to a quantum field, deriving the non-relativistic approximation becomes vastly more difficult. This is part of a whole circle of difficult problems in QED that stems from the fact that we don't know what combinations of field operators annihilate the (true) vacuum state. As a result, we don't know what operators create (true) single-particle states, and so on. 
(In relativistic QFT, perturbation theory is great for correlation functions, but it's useless for state-vectors, for reasons related to Haag's theorem. Discretizing everything evades Haag's theorem in the strict sense, invalidating superficial claims that QED "doesn't exist," but Haag's theorem still causes trouble: turning on the slightest interaction makes the new and old vacuum states practically orthogonal to each other, making perturbation theory useless for state-vectors.)
We can get around that problem by cheating a little. We can make an educated guess about what degrees of freedom are important in the NR approximation, define a new model with a quantum field for each of those degrees of freedom, and make another educated guess about how these degrees of freedom are related to the original fields in the original model, at least in the context of correlation functions where perturbation theory works. Once we've made these guesses, we can calculate a bunch of corresponding things in both theories, choosing the coefficients in the NR theory so that its predictions match those of the original theory, order-by-order in the appropriate expansion parameter(s). This is the NR Effective Field Theory approach. I listed a few references about this in my answer to 
Recovering nonrelativistic quantum mechanics from quantum field theory
But because of the initial guesswork, this isn't really a derivation of NRQED from QED. I wish I had a better answer.
By the way, if nuclear magnetic moments were important (but I don't think they are for X-ray diffraction), then relativistic QED wouldn't be the right foundation. QED + QCD would be better. This isn't an issue if we just start with NRQED, because then the magnetic moment term has its own independent coefficient that we can tune as needed.
