Limitations of Ehrenfest's theorem and non-stationary states I am trying to look at the cases where Ehrenfest's theorem would fail.
So for instance, the formal derivation makes use of the hermiticity of the Hamiltonian operator and thus I think the theorem would fail for non-hermitian hamiltonians (having noise terms included; think an open system, etc)
Also, I have not been able to find a good source that talks about the limitations of this theorem in detail.
Is this theorem fairly general or does one have to be careful while applying this theorem?
What about the case of discontinuous potentials, tunneling phenomena and time-dependent potentials?
Can the theorem be used in any way in such cases?
EDIT : Consider the following gif$-$

Are there any parts of this gif that cannot be explained by Ehrenfest's Theorem?
I am thinking, in particular, of the part where the probability density sneaks past the steep climbing potential (red curve) and shows the weird behavior (probably due to interference of the wavefunction ?).
 A: Ehrenfest's theorem is derived from the Schrödinger equation, but the derivation does not introduce any new approximation. Therefore, it must hold at least in any system correctly described by the Schrödinger equation (with a hermitian Hamiltonian). In particular, it should hold for discontinuous and/or time-dependent potentials
For systems described by other time-evolution equations, it is usually possible to generalize Ehrenfest's theorem. For example, for a non-hermitian Hamiltonian, we can write :
$$i\hbar \frac{\text d}{\text dt}|\psi\rangle = (H+i\Gamma)|\psi\rangle$$
with $H$ and $\Gamma$ hermitian operators. Then, the usual derivation of Ehrenfest's theorem yields, for any operator $A$ :
$$\frac{\text d}{\text dt}\langle A \rangle = \left\langle \frac{\partial A}{\partial t}\right\rangle + \frac{1}{i\hbar}\langle [A,H]\rangle + \frac1\hbar \langle\{A,\Gamma\}\rangle$$
As another system, open systems (with noise and dissipation) can sometimes be described by the Lindblad equation :
$$\dot\rho = -\frac i\hbar[H,\rho] + \sum_{j} \big(2 L_i\rho L_i^\dagger - L_i^\dagger L_i \rho - \rho L_i^\dagger L_i\big)$$
Then, for an operator $A$, its expectation value $\langle A \rangle = \operatorname{Tr}(A\rho)$ evolves according to :
\begin{align}
\frac{\text d}{\text dt}\langle A \rangle &= \operatorname{Tr}\left(\frac{\partial A}{\partial t}\rho  + A\dot \rho\right)\\
&=\left\langle \frac{\partial A}{\partial t}\right\rangle + \operatorname{Tr}\left(A \left(-\frac i\hbar[H,\rho] + \sum_{j} \big(2 L_i\rho L_i^\dagger - L_i^\dagger L_i \rho - \rho L_i^\dagger L_i\big)\right)\right) \\
&= \left\langle \frac{\partial A}{\partial t}\right\rangle - \frac{i}{\hbar}\operatorname{Tr}\left(AH\rho - A\rho H\right) + \operatorname{Tr}\left(\sum_{j} \big(2 AL_i\rho L_i^\dagger - AL_i^\dagger L_i \rho - A\rho L_i^\dagger L_i\big)\right) \\
&= \left\langle \frac{\partial A}{\partial t}\right\rangle - \frac{i}{\hbar}\big\langle [A,H]\big\rangle + \sum_{j}\left\langle 2L_i^\dagger AL_i- AL_i^\dagger L_i  -L_i^\dagger L_i A \right\rangle
\end{align}
Edits

*

*For a time-dependent Hamiltonian $H(t)$, the Ehrenfest theorem is unchanged :
$$\frac{\text d}{\text dt} \langle A\rangle = \left\langle\frac{\partial A}{\partial t}\right\rangle + \frac{1}{i\hbar} \langle [A,H(t)]\rangle$$


*The Ehrenfest theorem is equivalent to the Schrödinger equation, but in most cases the Schrödinger equation is itself an approximation, namely considering the system as isolated (ie neglecting noise, dissipation and entanglement with the environment). Moreover, this equivalence requires considering all possible initial states (to remove the expectation values). This derivation does not show how to compute (or "explain") the time evolution of one particular wave-function.


*The Ehrenfest theorem is (as far as I know) most useful (as in "most exploitable") when dealing with a well localized wave-packet, in which we can approximate $\langle \nabla V(\hat x)\rangle \simeq \nabla V(\langle \hat x\rangle)$. For a potential with terms of order $\geq 3$, knowing $\langle x\rangle, \langle p\rangle$ at $t=0$ is not enough to determine their further evolution. In your gif, it is easy to "predict/explain" the first 2 seconds using the formulae from Ehrenfest's theorem (for $x$ and $p$), while the further evolution (where the wave packets spreads a lot and interferes with itself) cannot.


*For a discontinuous potential the easiest way is to do the calculations like for distributions. Eg, for $H = \frac{p^2}{2m} + V\theta(x)$ with $\theta(x)$ the Heaviside step function, we find :
\begin{align}
\frac{\text d}{\text dt} \langle x \rangle &= \frac{p}{m} \\
\frac{\text d}{\text dt} \langle p\rangle &= \frac{1}{i\hbar}\langle [p,H]\rangle = \frac{V}{i\hbar} \langle [p,\theta(x)]\rangle
\end{align}
We can compute the commutator : for any states $|\psi\rangle,|\phi\rangle$, we have :
\begin{align}
\left\langle \psi ~\middle|~ \frac{1}{i\hbar}p\theta(x) ~\middle|~\phi\right\rangle  &= \int \psi^*(x) \frac{\text d}{\text dx} ( \theta(x)\phi(x))\text dx \\
&= \int \psi^*(x) ( \delta(x) \phi(x) + \theta(x) \phi'(x))\text dx \\
&= \psi^*(0)\phi(0) + \left\langle \psi ~\middle|~ \frac{1}{i\hbar}\theta(x)p ~\middle|~ \phi\right\rangle
\end{align}
So the Ehrenfest formulae read :
\begin{align}
\frac{\text d}{\text dt} \langle x \rangle &= \frac{p}{m} \\
\frac{\text d}{\text dt} \langle p\rangle &= V\langle \delta(x)\rangle = V|\psi(0)|^2
\end{align}
Edit 2
The Ehrenfest theorem is the statement :

For any initial state and any observable $A$, the expectation value satisfies :
$$\frac{\text d}{\text dt}\langle A \rangle = \left\langle \frac{\partial A}{\partial t}\right\rangle + \frac{1}{i\hbar}\langle [A,H]\rangle$$

This is equivalent (assuming that time-evolution is linear and unitary) to the Schrödinger equation :

For any initial state $|\psi(t=0)\rangle$, we have :
$$i\hbar \frac{\text d}{\text dt} |\psi\rangle = H|\psi\rangle$$

This equivalence is very general : with arbitrary potential, in finite dimensional systems, many-particles systems with spins, etc.
In the case of a higher order potential, the Ehrenfest equations are true, but hard to use to compute anything useful, as soon as the wave-packet is no well localized.
A: 
On the same lines, for cases of discontinuous potentials doesn't the right-hand side (expectation of derivative of potential) become undefined in the equation given by Ehrenfest theorem?

If $V$ is discontinuous, then its derivative can be understood in the distributional sense.  The idea is that even if $f$ does not have a well-defined derivative in the ordinary sense of a difference quotient, there may exist a (possibly generalized) function $g$ such that, for all smooth and rapidly decreasing functions $\psi$, we have that
$$\int f(x) \psi'(x) dx = -\int g(x) \psi(x) dx$$
We then call $g$ the distributional derivative of $f$; if $g$ is a genuine function (not e.g. a delta function or some other generalized function), it is called a weak derivative of $f$.  Clearly if $f$ is differentiable then $g=f'$ via integration by parts; however, if e.g. $f$ is the Heaviside function
$$f(x) = \begin{cases} 0 & x<0 \\ 1 & x\geq 0 \end{cases}
$$
then one finds that $g(x) = \delta(x)$.  In that sense, the delta function is the distributional derivative of the Heaviside function.
In precisely this way, a discontinuous function may have a weak or distributional derivative, and since we integrate when taking expected values, we find
$$\langle V'(\hat X)\rangle_\psi := \int \psi^*(x) V'(x) \psi(x)  dx = - \int V(x) \frac{d}{dx}\big[\psi^*(x)\psi(x)\big] dx$$

For instance, in the gif given in the question it is straight forward to see (for most part) how the mean value follows the Ehrenfest theorem (loosely) gaining speed where the curve steeps down etc. But at the left end where the wave reaches a sharply increasing potential it no longer is clear to me how I can apply Ehrenfest there.

Ehrenfest's theorem is a statement about the expectation values of position and momentum for a particular state.  If the wavefunction is e.g. a Gaussian wavepacket, then these expected values can be understood easily and intuitively as the position and momentum of a semi-classical point particle and in this case Ehrenfest's theorem looks like the corresponding classical equations of motion.
On the other hand, if the wavefunction is not a simple wavepacket, then these interpretations are no longer straightforward and obvious.  Any arbitrary square integrable function can be written as a superposition of momentum eigenstates via Fourier transform, but if it's some rapidly oscillating mess then you're going to need to sit down and actually calculate precisely what those expected values are - and when you do, you'll see (barring technicalities like domain issues, the existence of all the relevant expectation values, etc) that the results will obey Ehrenfest's theorem (though that fact may yield relatively little physical intuition).

Also, this question was motivated from an assignment which had given the above gif and asked to tell the part where the Ehrenfest theorem 'fails'.

My response would be to ask for further clarification from my instructor, essentially just as I asked for clarification from you in my comment. Perhaps they mean that a naive application of the theorem wherein one simply imagines wavepackets always evolve into different wavepackets fails; perhaps they mean for you to understand that $\langle V(\hat X)\rangle \neq V(\langle \hat X\rangle )$ in general, and that the naive expectation that these quantities are equal fails dramatically at some point in the gif.


Is it almost always possible to treat the derivative of V in the distributional sense?

It depends on what you mean by almost always.  Bounded functions which have a countable number of finite jumps have well-defined $n^{th}$ derivatives for all $n$, for example.  Functions such as
$$f(x) = \begin{cases} 0 & x\leq 0 \\ \frac{1}{x} & x > 0\end{cases}$$
do not, because there are smooth, rapidly-decreasing functions $\psi$ such that $\int f(x) \psi'(x) dx$ does not exist; more generally, it is sufficient that $f$ be locally integrable, meaning that for all intervals $[a,b]$, the integral $\int_a^b f(x) dx$ is finite.

I am interested in finding out about the technical issues you mention "(barring technicalities like domain issues, the existence of all the relevant expectation values, etc)"

A few potential technicalities are:

*

*If your Hilbert space is $L^2(\mathbb R)$, unbounded operators like $\hat X$ and $\hat P$ are not defined on the entire Hilbert space.  There exist square-integrable functions $f$ such that $\hat X f$ is not square-integrable; such functions are not in the domain of $\hat X$.  An example is the square-integrable function
$$f(x) = \begin{cases}0 & x <1 \\ \frac{1}{x} & x\geq 1\end{cases}$$

*Even if $f\in \mathrm{dom}(\hat X)$, its possible that $\langle \hat X\rangle_f$ may not be defined.  An example is the function
$$f(x) = \begin{cases}0 & x \leq 0 \\ \frac{1}{x} & x >0 \wedge x < 1 \\ 0 & x > 1\end{cases}$$
Therefore, when you see an expression like
$$m\frac{d}{dt}\langle \hat X\rangle_f = \langle \hat P\rangle_f$$
then you should implicitly recognize that $f$ has to be in the domain of both $\hat X$ and $\hat P$, and must be such that the relevant expectation values exist.
