Normal modes themselves are only the patterns in which the system naturally oscillates. You don't use initial conditions to determine the normal modes themselves; instead, initial conditions would tell you how the system's energy is distributed among the normal modes.
As far as how to continue...
Yep, that's the idea.
A good way to think about it is that the components of \mathbf{p} are multiplicative operators. Just like the differential operator \frac{\partial}{\partial x} takes a function f(x) and turns it into f'(x), a multiplicative operator x takes a function f(x) and turns it...
This isn't really a dot product, though, it's a dot "application." That is, it still means that \mathbf{a}\cdot\mathbf{b} = a_x b_x + a_y b_y + a_z b_z, but you can't assume that e.g. a_x b_x means "a_x multiplied by b_x," because you are dealing with things that don't get multiplied. Think...
Like Mute said, binomial expansion is a special case of Taylor expansion, so you could say they're the same thing in that sense.
It's standard practice to say "binomial expansion" to mean any expansion of this sort, even if it doesn't use the binomial theorem for nonnegative integer exponents.
If I remember this stuff correctly, then yes. You're just using a unit vector in the direction of the electric field rather than in the z direction. Alternatively, you could rotate your coordinate system so that the electric field points in the z direction, solve the problem, and then rotate...
Think about this: equal probability of emerging in any direction means the probability distribution is uniform over a two-dimensional sphere. You could write this as ΔP=f(θ,φ)ΔθΔφ. What would f(θ,φ) be in that case? Remember to think the properties of spherical coordinates.
You wound up with the right answer here, but your work is a mess because you keep reusing the same variables for different indices. It makes it almost impossible to follow. Now, if you can do this for yourself without getting confused, then it's fine, but for anything you're going to be sharing...
Oh wait, I think I may have slightly misled you - I forgot about the fact that \epsilon^{ab} is antisymmetric. That automatically takes care of the inversion of the transformation, so you don't need to make one of them \phi^a \to \phi^a - \epsilon^{ab}\phi^b; you can just use \phi^a \to \phi^a +...
That's right, \lambda(\phi^a\phi^a)^2 = \lambda\times(\phi^a\phi^a)\times(\phi^a\phi^a) \sim \lambda \Phi^4. Pretty much every Lagrangian you encounter in basic QFT will be a perturbative expansion of the form \sum_n(\text{const.})(\text{fields})^n, but you don't get arbitrary functions.
Now...
\lambda is a constant, if that helps, so it doesn't really act in any nontrivial manner... other than that, I'll come back and take a closer look at this shortly.