The problem of using the free fields inside a dielectric has been discussed elsewhere as well, for example in 'Why you should not use the electric field to quantize in nonlinear optics' by Quesada and Sipe. They present a simple and general argument for why Faraday's law does not hold when quantizing the dielectric Hamiltonian using the free field.
In this case, the nonlinear Hamiltonian of order $N$ will be some $N+1$ order polynomial of the creation and annihilation operators:
$$
H = \mathrm{poly}_{N+1}(\hat{a},\hat{a}^{\dagger})
$$
while the electric field is (by construction) a first order polynomial of these operations, and in a non-magnetic material the magnetic field is as well. Faraday's law is
$$
\tag{1}
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}.
$$
Using Heisenberg's equation of motion the r.h.s. can be written as
$$
\frac{\partial \mathbf{B}}{\partial t} = -\frac{i}{\hbar}[\mathbf{B},H] =
[\mathrm{poly}_{1}(\hat{a},\hat{a}^{\dagger}),\mathrm{poly}_{N+1}(\hat{a},\hat{a}^{\dagger})]
$$
Using the relations $[\hat{a},f(\hat{a},\hat{a}^{\dagger})] = \frac{\partial f}{\partial \hat{a}}$ and $[\hat{a}^{\dagger},f(\hat{a},\hat{a}^{\dagger})] = -\frac{\partial f}{\partial \hat{a}^{\dagger}}$ we can infer that
$$
\frac{\partial \mathbf{B}}{\partial t} = \mathrm{poly}_{N}(\hat{a},\hat{a}^{\dagger})
$$
but this is in contradiction with
$$
\nabla \times \mathbf{E} = \mathrm{poly}_{1}(\hat{a},\hat{a}^{\dagger})
$$
and thus (1) does not hold.
In practice this leads to several problems. A straightforward one, pointed out by Hillery, is that if you go through the trouble of deriving the wave equation for the free fields inside the dielectric you get
$$
\nabla^2 \mathbf{E} = \frac{1}{n^2c^2}\frac{\partial^2 \mathbf{E}}{\partial t}
$$
and the solutions to this equation are waves that propagate at $nc > c$. As discussed by Quesada et al. there are also problems in nonlinear interactions, for example in nonlinear self-interactions.
When is the approach using the free field justified? One could say "never", since it is simply wrong, but people still use it to get correct results. One reason for this is that parametric down-conversion is often described in the interaction picture, and the time-dependence of the field operators is usually introduced explicitly. For example, Loudon (2000) introduces the electric-field operator
$$
\hat{E}^+ = i\int_0^{\infty} \left(\frac{\hbar\omega}{4\pi\epsilon_0 c A \eta(\omega)} \right)^{1/2} \hat{a}\ \mathrm{exp}[-i\omega t + ik(\omega)z].
$$
This operator corresponds to waves propagating at the correct speed by construction, and the problem of the Hamiltonian actually generating a different time-evolution is swept under the rug. As for the problems with nonlinear interactions, these don't really appear in the low-gain regime of PDC, and it's worth keeping in mind that most perturbative approaches to PDC completely ignore time ordering as well (which does not work when self-interaction becomes significant).
Finally, regarding your question about the QFT formulation, it is not actually necessary to start from a Lagrangian density and quantize the field "properly". An alternative approach is presented in the resources you linked. The correct approach is to quantize the $D$-field instead of the $E$-field. The reason for this is that the $E$-field is not fundamental inside a dielectric medium, and in a macroscopic approach to quantization the fundamental excitations inside the medium are in fact combinations of light and matter excitations.
Hillery goes through the steps of how to re-write the nonlinear dielectric Hamiltonian in terms of the $D$-field and the inverse susceptibilities $\beta$, though this process is a bit tedious. After introducing the appropriate creation operator, you end up with the field operators
$$
\mathbf{D}(\mathbf{r},t) =
\sum_{\mathbf{k},\alpha}
i\sqrt{\frac{\hbar \epsilon_0 \omega_k }{2 V}}
\hat{e}_{\mathbf{k},\alpha} \Big(
\hat{a}_{\mathbf{k},\alpha}(t)
e^{i\mathbf{k}\cdot \mathbf{r}}
- \hat{a}^{\dagger}_{\mathbf{k},\alpha}(t)
e^{-i\mathbf{k}\cdot \mathbf{r}}\Big)
$$
$$
\mathbf{B}(\mathbf{r},t) = \sum_{\mathbf{k},\alpha}
i\sqrt{\frac{\hbar}{2 \epsilon_0 \omega_k V}}
(\mathbf{k}\times\hat{e}_{\mathbf{k},\alpha}) \Big(
\hat{a}_{\mathbf{k},\alpha}(t)
e^{i\mathbf{k}\cdot \mathbf{r}}
- \hat{a}^{\dagger}_{\mathbf{k},\alpha}(t)
e^{-i\mathbf{k}\cdot \mathbf{r}}\Big) \\
$$
If you take the corresponding free Hamiltonian $H = \int \left(\frac{1}{2\mu_0}\mathbf{B}^2 + \frac{1}{2}\beta^{(1)}\mathbf{D}^2 \right)dr^3$
and apply Heiseberg's equation of motion, you will see that the $D$-field obeys
$$
\nabla^2 \mathbf{D} =
\frac{n^2}{c^2}\frac{\partial^2 \mathbf{D}}{\partial t^2}
$$
and the dynamics are therefore correct. In conclusion, you can use the "standard" description by simply choosing to work with the displacement field instead of the electric field. You do also see this approach taken in some of the literature. However, when working in the interaction picture with the field operators introduced by hand, the answer for most experimentally relevant quantities are the same.
