What happens when a single photon travels through an dielectric medium? Let's consider a single photon traveling along the $z$ direction, we can describe such a photon with the quantum state
$$
\vert f\rangle
=
\int\frac{dp_z}{(2\pi)\sqrt{2p_z}}
f(p_z)
e^{ip_z(t-z)}
\hat{a}^\dagger(0,0,p_z)
\vert0\rangle
\tag{1}
$$
in natural units and $f(p_z)$ is the spectrum/wave function of the photon.
For instance, $f(p_z)\propto e^{-(p_z-k_0)^2/4\sigma^2}$ would correspond to a Gaussian wave packet with momentum centered at $k_0$ with spread $\sigma$ in momentum space.
What happens now, when we place a dielectric medium of length $L$ with dielectric constant $\varepsilon$ in the path of the photon?
We know that the energy density of the electric field inside a medium with dielectric constant $\varepsilon$ is $\varepsilon E^2$, therefore, we can write the Hamiltonian density
$$
\mathcal{H}
=
\frac{1}{2}
\left(
  \varepsilon(z)E(t,z)^2
  +
  B(t,z)^2
\right)
\tag{2}
$$
wherein $E,B$ are the electric and magnetic fields and
$$
\varepsilon(z)
=
\begin{cases}
\varepsilon & 0\leq z\leq L
\\
0 & \text{otherwise}
\tag{3}
\end{cases}.
$$
Inserting the mode decomposition of the electric and magnetic field
$$
\begin{aligned}
E(t,z)
&=
\int\frac{dp_z}{(2\pi)\sqrt{2p_z}}\left(-ip_z\right)
\left\{
\hat{a}(0,0,p_z)e^{-ip_z(t-z)}
-
\text{h.c.}
\right\}
\\
B(t,z)
&=
\int\frac{dp_z}{(2\pi)\sqrt{2p_z}}\left(+ip_z\right)
\left\{
\hat{a}(0,0,p_z)e^{-ip_z(t-z)}
-
\text{h.c.}
\right\}
\end{aligned}
\tag{4}
$$
I believe (?) we find for the Hamiltonian
$$
H
=
\int dz\ \mathcal{H}
=
-
\frac{1}{2}
\int\frac{dp_z}{2\pi}
\omega(p_z)
\hat{a}^\dagger(0,0,p_z)
\hat{a}(0,0,p_z)
\tag{5}
$$
wherein
$$
\omega(p_z)
=
\begin{cases}
p_z & \text{outside the dielectric}
\\
(1+\varepsilon)p_z & \text{inside the dielectric}
\end{cases}
\tag{6}.
$$
So the presence of the dielectric rescales the eigenfrequencies of the field (I think this statement is true although I am not sure if the explicit form eq. (6) of my calculation is correct).
But what does that mean to the coherent state? How does it change while propagating through the dielectric?
My idea is to switch to the Schroedinger picture where the state is time-independent and evolve the state with the time-evolution operator
$$
U
=
\exp\left\{
-i\int_0^{L/v}dt^\prime
H(t^\prime)
\right\}
=
\exp\left\{
-i\omega(p_z)L/v
\right\}
\tag{7}
$$
where we used the $\omega(p_z)$ for inside the cavity. The coherent state leaving the dielectric should then be given by $U\vert f\rangle$.
What do you think? Is this reasonable?
Update 1:
Building onto the answer given by Tales Rick Perche lets take the interaction Hamiltonian to be
$$
H_\text{int}(t)
=
\int_{-L/2}^{+L/2}dz\colon E(t,z)^2\colon
=
\frac{1}{4\pi^2}
\int_{-L/2}^{+L/2}dz
\int dp
\int dq
\left\{
e^{-i(p+q)(t-z)}
\hat{a}_p\hat{a}_q
+
e^{-i(p-q)(t-z)}
\hat{a}_p\hat{a}_q^\dagger
+
e^{+i(p-q)(t-z)}
\hat{a}_p^\dagger\hat{a}_q
+
e^{-i(p+q)(t-z)}
\hat{a}_p^\dagger\hat{a}_q^\dagger
\right\}
.\tag{8}
$$
In the optical domain, we commonly apply the rotating wave approximation which lets us drop the nonlinear terms
$$
H_\text{int}(t)
=
\frac{1}{4\pi^2}
\int_{-L/2}^{+L/2}dz
\int dp
\int dq
\left\{
e^{-i(p-q)(t-z)}
\hat{a}_p\hat{a}_q^\dagger
+
e^{+i(p-q)(t-z)}
\hat{a}_p^\dagger\hat{a}_q
\right\}
.\tag{9}
$$
Now using
$$
\int_{-L/2}^{+L/2}dz\ e^{iaz}
=
\frac{2}{a}
\sin\left(\frac{aL}{2}\right)
\tag{10}
$$
we find
$$
H_\text{int}(t)
=
\frac{1}{2\pi^2}
\int dp
\int dq
\frac{\sin\left(\frac{p-q}{2}L\right)}{p-q}
\left\{
e^{-i(p-q)t}
\hat{a}_p\hat{a}_q^\dagger
-
e^{+i(p-q)t}
\hat{a}_p^\dagger\hat{a}_q
\right\}
\tag{11}
$$
Using the first term of the Magnus expansion leads us to an approximation of the interaction-picture time-evolution operator
$$
\begin{aligned}
U_\text{int}(t_1,t_0)
&=
\exp\left\{
-i\int_{t_0}^{t_1}dt^\prime
H_\text{int}(t^\prime)
\right\}
\\
&=
\exp\left\{
\frac{1}{2\pi^2}
\int dp
\int dq
\frac{\sin\left(\frac{p-q}{2}L\right)}{p-q}
\left\{
\frac{e^{+i(p-q)t_0}-e^{+i(p-q)t_1}}{p-q}
\hat{a}_p\hat{a}_q^\dagger
+
\frac{e^{-i(p-q)t_0}-e^{-i(p-q)t_1}}{p-q}
\hat{a}_p^\dagger\hat{a}_q
\right\}
\right\}
\end{aligned}
\tag{12}
$$
where $t_1=t_0+L/v_{gr}$ where $v_{gr}$ is the group velocity inside the medium. Defining $t_0=0$ to be the time the wave packet hits the dielectricum at $z=-L/2$, eq. (12) simplifies to
$$
\begin{aligned}
U_\text{int}(T)
&=
\exp\left\{
\frac{1}{2\pi^2}
\int dp
\int dq
\frac{\sin\left(\frac{p-q}{2}L\right)}{p-q}
\left\{
\frac{1-e^{+i(p-q)T}}{p-q}
\hat{a}_p\hat{a}_q^\dagger
+
\frac{1-e^{-i(p-q)T}}{p-q}
\hat{a}_p^\dagger\hat{a}_q
\right\}
\right\}
\end{aligned}
\tag{13}
$$
with $T=L/v_{gr}$ being the interaction (or transmit) time.
The argument of eq. (13) is strongly enhanced at the resonance at $p=q$ where the momenta match. At the same time eq. (13) has some similarity with the generator of beam splitter transformations.
Update 2:
I found the answer to the questions in Macroscopic quantum electrodynamics — concepts and applications by Scheel and Buhmann for the lossless beam splitter (p. 11) and the absorptive beam splitter (compatible with statistical physics) (p. 32).
Additionally, QED in Dispersing and Absorbing Media by Knoell, Scheel and Welsch show a general method using the method of Green tensors (p. 20) which is valid in classical electrodynamics but applies also to the complex Fourier amplitudes.
Update 3:
The sinc is our interaction Hamiltonian is also known as phase-matching function, see Very nonlinear quantum optics by Mejia. p. 33
 A: You are using a very simplified treatment for the problem that ignores perhaps the most important part of it: the boundary conditions for the electromagnetic field in contact with the dielectric.
However, considering this treatment, your expression (5) for the Hamiltonian is incorrect. The Hamiltonian is not position dependent, and in your expression you expand it in momentum basis, so $\omega(p)$ cannot have different values depending on the position: $\omega$ does not depend on the position, otherwise the Hamiltonian would.
In this simplified formulation of the problem, the correct way of approaching this problem is by considering the Hamiltonian, given by
$$H = \int d z \left(\epsilon(z) :E(t,z)^2: + :B(t,z)^2:\right).$$
The normal ordering is necessary to prevent divergences, and is present even in the standard case, with no dielectric medium. Notice that $\epsilon(z)$ should be defned by
\begin{equation}
   \epsilon(z) = \begin{cases}\epsilon, 0\leq z\leq L\\
1, \text{otherwise}\end{cases}.
\end{equation}
After plugging in the values for $E$ and $B$ given by your expression (4), you will obtain a result that will look like
$$H = \frac{1}{4\pi}\int dp \:p\: \hat{a}^\dagger_p \hat{a}_p +(\epsilon-1)\int_0^L d z :E(t,z)^2:,$$
so that the modification to the dynamics of the system will come from the second term on the right hand side of the equation above. However, you would not be able to use expressions for the Dirac delta or other functions in that expression due to the finite interval, and it would end up looking quite bad:
$$\begin{aligned}
\int_0^L dz :E(t,z)^2: = \frac{1}{8\pi^2} \int_0^L dz \int dp \int dp'\Big(e^{-i(p+p')(t-z)}\hat{a}_p\hat{a}_{p'}+e^{-i(p-p')(t-z)}\hat{a}_{p'}^\dagger\hat{a}_{p}\\+e^{i(p-p')(t-z)}\hat{a}_p^\dagger\hat{a}_{p'}+e^{i(p+p')(t-z)}\hat{a}_p^\dagger\hat{a}_{p'}^\dagger\Big).
\end{aligned}$$
Using
$$\int_0^L dz e^{az} = \frac{1}{a}(e^{aL}-1)$$
in each one of the above integrals you would get your result. (It looks a bit nicer making the cavity from $-L$ to $L$, so that everything becomes $\text{sinc}$ functions). In essence, there is no way of handling the quadratic terms on the creation and annihilation operators, which will make it impossible to obtain a simple solution for your state. Perturbation theory could be an option for finite times...
In any case, the appropriate way of handling this problem is to impose boundary conditions for the electric field inside the cavity. Under these conditions the field modes inside it would be discrete, while the modes outside would be continuous. However, it is very important to notice that the notion of photon wavepacket you have defined there depends explicitly on whether your field is within a cavity, or in infinite space. This is a general feature of particle excitations in quantum field theory: they are not local, and depend on the mode decomposition chosen for the field, which depends on the boundary conditions. In fact, defining particle states for the case of fields with general boundary conditions is not an easy problem to tackle.
