Diffraction in Electromagnetic Waves
You seem to be asking two questions: (1) what is the medium for the EM field is and (2) what is the mechanism for diffraction.
There is no material medium - this is what the Michelson-Morley experiment was all about. You could say that the separate harmonic oscillators in the second quantized EM field and the photon field they represent have replaced the "aether" sought in the MM experiment. The difference is that the second quantized field can be made Lorentz invariant (respecting special relativity), well almost (see footnote). The Wikipedia entry on this topic uses a gauge that is not Lorentz invariant, but this problem can be overcome and the Wikipedia page gives you a good idea of what's (messily!) involved.
The simplest answer to what is the mechanism of diffraction I can think of is this: consider a field on a plane, say $z = 0$ and split it up using Fourier decomposition into constituent plane waves. Each constituent plane wave has a different direction defined by the wavevector $\left(k_x, k_y, k_z\right)$ with $k^2 = k_x^2 + k_y^2 + k_z^2$, i.e. all the wavevectors have the same magnitude but different directions. So, when we ask what the field looks like at a different value of $z$, we build the field up from our plane wave constituents at this point (use an inverse Fourier transform). However, now, because the wavevectors are all in different directions, the plane waves have all undergone different phase delays in reaching the new value of $z$ (even though their phase advances by $k$ radians per unit length in the direction of the respective wave vector). Therefore, the field's configuration gets scrambled by all these different phase delays.
You can try this thinking with your slit - set the analysis up in two dimensions so that there are only 2D wave vectors for simplicity. The screen with the slit is in the $z = 0$ plane and the one orthogonal direction is the $x$ axis. All the Cartesian components of the fields fulfil the same (Helmholtz) equation, so we can discuss the principles by just looking at one scalar field $\psi$ (say, the electric field's $x$-component). Each plane wave has the form $\psi(k_x) = \exp\left(i (k_x x + k_z z)\right)$ The Fourier transform of the field output from the slit is then $\frac{1}{k_x}\sin\left(\frac{w k_x}{2}\right)$ where $w$ is the slit width, and unless the slit is very wide, the Fourier transform has a wide spread of frequencies. This means that for $z > 0$ ("downstream of the the slit's output) the field is the superposition $\int\limits_{-\infty}^\infty \frac{1}{k_x}\sin\left(\frac{w k_x}{2}\right) \exp\left(i (k_x x + k_z z)\right) d k_x$. When we plug $z = 0$ in, the integral is simply the inverse FT and we get our original slit field. But now put some nonzero value of $z$ in: because $k_x^2 + k_z^2 = k^2$, we have $k_z = \sqrt{k^2 - k_x^2}$ (assuming the field is running in the $+z$ direction), we get $\int\limits_{-\infty}^\infty \frac{1}{k_x}\sin\left(\frac{w k_x}{2}\right) \exp\left(i (k_x x + \sqrt{k^2 - k_x^2} z)\right) d k_x$. You can see the "scrambling", $k_x$-dependent phase factor $\exp\left(i \sqrt{k^2 - k_x^2} z\right) = \exp\left(i k \cos\theta_x\right)$ (where $\theta_x$ is the angle that the plane wave with wavevector $(k_x, k_z)$ makes with the $z$-axis) will yield the complicated scrambling you see with the slit. Various approximations, notably Fraunhofer and Fresnel, are applied to this integral. Interestingly, the $thinner$ your slit, the wider the direction spread and the more noticeable your diffraction. Wide slits with coherent light produce nondiverging beams (like a laser beam). This is related to the uncertainty principle: the space of momentum eigenstates is found by Fourier transforming the space of position eigenstates, so that the more spread out in one domain a wave is, the more tightly packed in the other it is.
Footnote: A second quantized EM field comprises an infinite number of oscillators - countably so in a box, uncountably so in $\mathbb{R}^3$. The ground state therefore has an energy $\frac{h \nu}{2}$ for each oscillator, and so holds infinite energy. Depending on the approach, this may or may not be a problem - there are procedures for subtracting this infinity out of calculations. However, one may want to "regularize" calculations by assuming that the oscillators don't go up to infinite frequency, but instead are "cut off" at some high frequency (presumably well above the frequencies of cosmic ray $\gamma$ rays). So, although the second quantized EM field can be made Lorentz invariant, I don't believe there is any way one can make this cutoff process truly Lorentz invariant (the cutoff frequency, unless infinite, will depend on observers).