Do Quantum Averages Obey Newton’s Laws? A Visual Guide to the Ehrenfest Theorem

Last Updated on July 27, 2025 by Sushanta Barman

In quantum mechanics, particles are described by wave functions, not classical paths. However, the Ehrenfest theorem provides a clear connection between quantum and classical physics.

Simply put, it tells us that the average behavior of quantum particles can closely follow Newton’s classical laws under certain conditions.

In this article, we will learn what the Ehrenfest theorem is, why it matters, and how to verify it numerically using split-step Fourier method (SSFM) simulation.

Ehrenfest Theorem

The Ehrenfest theorem states that, in quantum mechanics, the time evolution of the expectation values (averages) of position and momentum for a particle follows equations that closely resemble Newton’s laws of motion for classical particles.

Mathematically, it is expressed as:

\[\frac{d}{dt}\langle x (t) \rangle = \frac{\langle p (t) \rangle}{m}, \quad \frac{d}{dt}\langle p (t) \rangle = -\left\langle \frac{\partial V}{\partial x} \right\rangle\]

where \(\langle x (t) \rangle\) is the average (expectation value) of the particle’s position, \(\langle p (t) \rangle\) is the average momentum, \(m\) is the mass, and \(V(x)\) is the potential energy as a function of position. The notation \(\langle \rangle\) means the quantum expectation value, calculated using the wave function of the particle.

This theorem is important because it bridges quantum mechanics with classical physics. It shows that, under suitable conditions, the average motion of a quantum particle behaves just like a classical particle following Newton’s laws.

Example: An Electron in a Constant Electric Field

Let’s take an example to see how the Ehrenfest theorem works in practice.

We consider an electron moving in one dimension under the influence of a constant electric field \(\vec{E}\), as shown in Fig. 1. The electron has charge \(q=−e\) and mass \(m\).

Potential energy: The potential energy of the electron in this field is given by:

\[V(x) = -qEx.\]

Classical equations of motion: According to Newton’s laws, a classical particle in this potential would follow:

\[x_{\text{cl}}(t) = x_0 + \frac{p_0}{m}t + \frac{qE}{2m}t^2, \quad p_{\text{cl}}(t) = p_0 + qEt,\]

where \(x_0\) is the initial position and \(p_0\) is the initial momentum.

Quantum mechanical expressions: For the quantum version, we focus on the expectation values of position and momentum:

• The average position: \(\langle x(t) \rangle = \int x\,|\psi(x, t)|^2 dx\), where \(\psi(x,t)\) is the wave function of the particle.
• The average momentum: \(\langle p (t)\rangle = \int \psi^*(x, t)\left(-i\hbar \frac{\partial}{\partial x}\right)\psi(x, t)\, dx\)

According to the Ehrenfest theorem, these quantum averages should satisfy:

\[\frac{d}{dt}\langle x (t)\rangle = \frac{\langle p (t) \rangle}{m}, \quad \frac{d}{dt}\langle p (t)\rangle = qE.\]

How we compare in simulation: In our simulation, we shall calculate both the classical solutions and the quantum expectation values as functions of time. We then compare:

• The quantum average position \(\langle x (t) \rangle\) with the classical trajectory \(x_{\text{cl}}(t)\)
• The quantum average momentum \(\langle p (t) \rangle\) with the classical momentum \(p_{\text{cl}}(t)\)

If the Ehrenfest theorem holds, these quantum and classical results should closely match—especially in the case of a constant electric field, where the potential is linear.

Simulating Quantum Motion

In this section, we explain how to solve the time-dependent Schrödinger equation for our system using the split-step Fourier method (SSFM). We also discuss how to calculate the quantum averages of position, momentum, and force from the simulated wave function.

The Time-Dependent Schrödinger Equation (TDSE) for the System

For an electron in a constant electric field \(\vec{E}\) along \(x-\)axis, the potential energy is linear in position: \(V(x) = -qEx\). So the TDSE for our system becomes:

\[i\hbar \frac{\partial \psi(x, t)}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2 \psi(x, t)}{\partial x^2} – qEx\,\psi(x, t)\]

where, \(\psi(x, t)\) is the wave function, describing the quantum state of the electron, \(m\) is the mass of the electron, \(q\) is the charge of the electron (with \(q = -e\)), and \(\hbar\) is the reduced Planck constant.

Our goal is to numerically solve this equation to observe how the quantum state of the electron evolves over time in the presence of a constant electric field.

Split-Step Fourier Method

The split-step Fourier method is a powerful and efficient technique to solve the TDSE numerically. Here’s how it works:

1. Discretization:

• Space: Divide the spatial domain into a grid with \(N\) points. Each point is spaced by \(\Delta x\), covering the range \(x \in [0, L_x)\).
• Time: Choose a small time step \(\Delta t\) and advance the solution in discrete time steps.

2. Initial Wave Function: Start with an initial wave packet, typically a Gaussian centered at position \(x_0\) with momentum \(p_0\):
\[\psi(x, 0) = A \exp\left(-\frac{(x – x_0)^2}{2\sigma_x^2}\right) \exp\left(\frac{i p_0 (x – x_0)}{\hbar}\right),\]
where, \(x_0\) is the initial center position of the wave packet, \(\sigma_x\) is the spatial width, and \(p_0\) is the average momentum. The factor \(A\) is a normalization constant chosen so that the total probability \(\int |\psi(x, 0)|^2 dx\) is equal to 1. The first exponential term ensures the wave packet is localized around \(x_0\), while the second term adds a phase corresponding to an average momentum \(p_0\).

This Gaussian wave packet is a standard choice for simulating quantum evolution because it is both localized and has a defined mean momentum.

3. Evolution Steps: At each time step, the wave function is evolved by alternately applying the potential and kinetic operators. This is done in three main sub-steps for each small time step \(\Delta t\):

a) Potential Step (First Half): Multiply the wave function by an exponential factor involving the potential for half a time step:

\[\psi(x) \leftarrow \psi(x) \cdot \exp\left(-\frac{i}{2\hbar} V(x) \Delta t \right)\]

b) Kinetic Step (Full): Transform the wave function to momentum (Fourier) space using a Fast Fourier Transform (FFT). In this space, apply the kinetic energy operator:

\[\tilde{\psi}(k) \leftarrow \tilde{\psi}(k) \cdot \exp\left(-\frac{i\hbar k^2}{2m} \Delta t\right)\]

Then, transform back to real space using the inverse FFT.

c) Potential Step (Second Half): Again, multiply by the potential exponential for the remaining half-time step:

\[\psi(x) \leftarrow \psi(x) \cdot \exp\left(-\frac{i}{2\hbar} V(x) \Delta t \right)\]

• Repeat these steps (a-c) for each time step to evolve the wave function forward in time.
• This splitting keeps the numerical method both accurate and efficient.

Advantages:

• The method is fast due to the efficiency of FFTs.
• It preserves the normalization of the wave function (that is, the total probability \(\int |\psi(x, t)|^2 dx = 1\) should remain constant during the simulation).
• It is stable for long-time simulations when using small enough time steps.

Quantum Averages of Position, Momentum, and Force

Once the wave function is obtained at each time step, we can compute key quantum observables:

1. Average Position (\(\langle x (t)\rangle\)):

\[\langle x (t) \rangle = \int x\, |\psi(x, t)|^2 dx\]

In the code, this is computed as a sum over all grid points:

\[\langle x (t) \rangle \approx \sum_{j} x_j\, |\psi(x_j, t)|^2\, \Delta x\]

2. Average Momentum (\(\langle p (t) \rangle\)):

\[\langle p (t) \rangle = \hbar \int \mathrm{Im}\left(\psi^*(x, t)\frac{\partial \psi(x, t)}{\partial x}\right) dx\]

Numerically, the derivative is calculated using spectral methods (Fourier techniques), making it both accurate and efficient.

3. Average Force (\(\langle F \rangle\)):
The force on the particle is given by the negative gradient of the potential:

\[F(x) = -\frac{dV}{dx}\]

For a constant electric field, \(F(x)\) is simply the constant \(qE\).
For more general potentials, you would compute the quantum average force as:

\[\langle F \rangle = – \int |\psi(x, t)|^2\, \frac{dV}{dx}\, dx\]

These quantum averages let us compare the quantum results directly with the classical predictions and check the validity of the Ehrenfest theorem.

Table 1: Values of the simulation parameters.

Python Code to visualize the Ehrenfest Theorem

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.ticker import ScalarFormatter
from dataclasses import dataclass

# -------------------------------
# Global style
# -------------------------------
plt.rcParams.update({
    "font.size": 14,
    "axes.titlesize": 18,
    "axes.labelsize": 16,
    "legend.fontsize": 13,
    "xtick.labelsize": 13,
    "ytick.labelsize": 13
})

BLUE = "#1f77b4"
RED  = "#d62728"

# -------------------------------
# Physical constants (SI)
# -------------------------------
HBAR = 1.054571817e-34     # J*s
M_E  = 9.1093837015e-31    # kg
E_CH = 1.602176634e-19     # C

@dataclass
class SimParams:
    Nx: int = 1024
    Lx: float = 400e-9         # 400 nm
    dt: float = 2e-16          # 0.2 fs
    steps: int = 2000          # 400 fs total
    frame_skip: int = 5        # ~400 frames
    mass: float = M_E
    charge: float = -E_CH
    Efield: float = 1e7        # V/m
    x0: float = 80e-9
    lambda_dB: float = 1.5e-9
    sigma_x: float = 10e-9

    @property
    def p0(self):
        return 2*np.pi*HBAR/self.lambda_dB  # h/lambda

def make_grid(params):
    dx = params.Lx / params.Nx
    x = np.linspace(0.0, params.Lx, params.Nx, endpoint=False)
    k = 2*np.pi * np.fft.fftfreq(params.Nx, d=dx)
    return x, k, dx

def initial_gaussian(x, x0, p0, sigma_x):
    norm = (1.0/(np.pi*sigma_x**2))**0.25
    psi = norm * np.exp(-(x-x0)**2/(2*sigma_x**2)) * np.exp(1j*p0*(x-x0)/HBAR)
    dx = x[1]-x[0]
    psi /= np.sqrt(np.sum(np.abs(psi)**2)*dx)
    return psi

def linear_potential(x, F0):
    return -F0*x

def expectation_x(x, psi, dx):
    return np.real(np.sum(x*np.abs(psi)**2)*dx)

def expectation_p_spectral(psi, k, dx):
    psi_k = np.fft.fft(psi)
    dpsi_dx = np.fft.ifft(1j*k*psi_k)
    integrand = np.imag(np.conj(psi)*dpsi_dx)
    return HBAR*np.sum(integrand)*dx

def finite_diff(arr, dt):
    deriv = np.zeros_like(arr)
    deriv[1:-1] = (arr[2:]-arr[:-2])/(2*dt)
    deriv[0] = (arr[1]-arr[0])/dt
    deriv[-1] = (arr[-1]-arr[-2])/dt
    return deriv

def simulate(params):
    x, k, dx = make_grid(params)
    dt = params.dt
    F0 = params.charge*params.Efield
    Vx = linear_potential(x, F0)
    expV = np.exp(-1j*Vx*dt/(2*HBAR))
    expK = np.exp(-1j*(HBAR*k**2)*dt/(2*params.mass))

    psi = initial_gaussian(x, params.x0, params.p0, params.sigma_x)

    nframes = params.steps//params.frame_skip + 1
    psi_list = np.zeros((nframes, params.Nx))
    times = np.zeros(nframes)
    x_exp = np.zeros(nframes)
    p_exp = np.zeros(nframes)

    # Classical COM
    v0 = params.p0/params.mass
    full_times = np.arange(params.steps+1)*dt
    x_class = params.x0 + v0*full_times + 0.5*(F0/params.mass)*full_times**2
    x_class_f = x_class[::params.frame_skip]

    # record initial
    f = 0
    psi_list[f] = np.abs(psi)**2
    times[f] = 0.0
    x_exp[f] = expectation_x(x, psi, dx)
    p_exp[f] = expectation_p_spectral(psi, k, dx)

    for n in range(1, params.steps+1):
        psi *= expV
        psi_k = np.fft.fft(psi)
        psi_k *= expK
        psi = np.fft.ifft(psi_k)
        psi *= expV
        if n % params.frame_skip == 0:
            f += 1
            psi_list[f] = np.abs(psi)**2
            times[f] = n*dt
            x_exp[f] = expectation_x(x, psi, dx)
            p_exp[f] = expectation_p_spectral(psi, k, dx)

    # Normalize panel 1
    psi_list /= psi_list.max()

    dt_frame = params.frame_skip*dt
    dxdt = finite_diff(x_exp, dt_frame)
    dpdt = finite_diff(p_exp, dt_frame)

    # Clean any NaNs/Infs
    dxdt = np.nan_to_num(dxdt)
    dpdt = np.nan_to_num(dpdt)

    F_array = np.full_like(times, F0)

    return x, times, psi_list, x_exp, p_exp, x_class_f, dxdt, dpdt, F_array, params

def animate():
    x, times, psi_list, x_exp, p_exp, x_class, dxdt, dpdt, F_array, params = simulate(SimParams())

    fig, axs = plt.subplots(2, 2, figsize=(8, 7))  # 16:9 aspect
    ax1, ax2, ax3, ax4 = axs.flatten()

    # Panel 1: |ψ|^2 (blue)
    line1, = ax1.plot([], [], color=BLUE, lw=2)
    ax1.set_xlim(0.0, params.Lx*1e9)
    ax1.set_ylim(0, 1.1)
    ax1.set_xlabel("x (nm)")
    ax1.set_ylabel("|ψ|²")
    ax1.set_title("Wavepacket evolution")

    # Panel 2: ⟨x⟩ vs classical x
    line2a, = ax2.plot([], [], color=BLUE, lw=2, label="⟨x⟩(t)")
    line2b, = ax2.plot([], [], color=RED,  lw=1.8, linestyle="--", label="Classical x(t)")
    mark2a, = ax2.plot([], [], 'o', color=BLUE, markersize=8)
    mark2b, = ax2.plot([], [], 'o', color=RED,  markersize=4)
    ax2.set_xlim(times[0]*1e15, times[-1]*1e15)
    ymin = min(x_exp.min(), x_class.min())*1e9
    ymax = max(x_exp.max(), x_class.max())*1e9
    margin = 0.05*(ymax-ymin if ymax>ymin else 1.0)
    ax2.set_ylim(ymin-margin, ymax+margin)
    ax2.set_xlabel("Time (fs)")
    ax2.set_ylabel("x (nm)")
    ax2.set_title("⟨x⟩ vs classical")
    ax2.legend(loc="best")

    # Panel 3: d⟨x⟩/dt vs ⟨p⟩/m
    line3a, = ax3.plot([], [], color=BLUE, lw=2, label="d⟨x⟩/dt")
    line3b, = ax3.plot([], [], color=RED,  lw=1.8, linestyle="--", label="⟨p⟩/m")
    mark3a, = ax3.plot([], [], 'o', color=BLUE, markersize=8)
    mark3b, = ax3.plot([], [], 'o', color=RED,  markersize=4)
    ax3.set_xlim(times[0]*1e15, times[-1]*1e15)
    vel = np.concatenate([dxdt, p_exp/params.mass])
    vmin, vmax = vel.min(), vel.max()
    vmargin = 0.05*(vmax - vmin if vmax>vmin else 1.0)
    ax3.set_ylim(vmin - vmargin, vmax + vmargin)
    ax3.set_xlabel("Time (fs)")
    ax3.set_ylabel("Velocity (m/s)")
    ax3.set_title("Velocity check")
    fmt3 = ScalarFormatter(useMathText=True)
    fmt3.set_powerlimits((-2, 2))
    ax3.yaxis.set_major_formatter(fmt3)
    ax3.ticklabel_format(style='sci', axis='y', scilimits=(-2, 2))
    ax3.legend(loc="best")

    # Panel 4: d⟨p⟩/dt vs -∂V/∂x
    line4a, = ax4.plot([], [], color=BLUE, lw=2, label="d⟨p⟩/dt")
    line4b, = ax4.plot([], [], color=RED,  lw=1.8, linestyle="--", label="-∂V/∂x = -eE")
    mark4a, = ax4.plot([], [], 'o', color=BLUE, markersize=6)
    mark4b, = ax4.plot([], [], 'o', color=RED,  markersize=4)
    ax4.set_xlim(times[0]*1e15, times[-1]*1e15)
    F0 = F_array[0]
    span_min = max(0.1*abs(F0), 1e-14)
    resid = np.max(np.abs(dpdt - F0))
    span = max(span_min, 5*resid)
    ax4.set_ylim(F0 - span, F0 + span)
    ax4.set_xlabel("Time (fs)")
    ax4.set_ylabel("Force (N)")
    ax4.set_title("Force check")
    fmt4 = ScalarFormatter(useMathText=True)
    fmt4.set_powerlimits((-2, 2))
    ax4.yaxis.set_major_formatter(fmt4)
    ax4.ticklabel_format(style='sci', axis='y', scilimits=(-2, 2))
    ax4.legend(loc="best")

    def init():
        artists = [line1, line2a, line2b, mark2a, mark2b, line3a, line3b, mark3a, mark3b, line4a, line4b, mark4a, mark4b]
        for ln in artists:
            ln.set_data([], [])
        return artists

    def update(frame):
        tfs = times[:frame+1]*1e15
        t_now = times[frame]*1e15

        line1.set_data(x*1e9, psi_list[frame])
        line2a.set_data(tfs, x_exp[:frame+1]*1e9)
        line2b.set_data(tfs, x_class[:frame+1]*1e9)
        mark2a.set_data([t_now], [x_exp[frame]*1e9])
        mark2b.set_data([t_now], [x_class[frame]*1e9])

        line3a.set_data(tfs, dxdt[:frame+1])
        line3b.set_data(tfs, (p_exp[:frame+1]/params.mass))
        mark3a.set_data([t_now], [dxdt[frame]])
        mark3b.set_data([t_now], [(p_exp[frame]/params.mass)])

        line4a.set_data(tfs, dpdt[:frame+1])
        line4b.set_data(tfs, F_array[:frame+1])
        mark4a.set_data([t_now], [dpdt[frame]])
        mark4b.set_data([t_now], [F_array[frame]])

        return [line1, line2a, line2b, mark2a, mark2b, line3a, line3b, mark3a, mark3b, line4a, line4b, mark4a, mark4b]

    ani = FuncAnimation(fig, update, frames=len(times), init_func=init, blit=False, interval=50)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    animate()
            

Below is a brief explanation of the above Pyhton script:

1. Imports and Styling
   It uses NumPy for arrays, Matplotlib for plotting/animation, and a simple @dataclass to hold simulation parameters. Global plot styles ensure clear, consistent figures.
2. Simulation Parameters
   The SimParams class defines the grid size (\(N_x\), \(L_x\)), time step (\(dt\) and total steps), physical constants (electron mass and charge) field strength, and initial wave‑packet properties (\(x_0\), \(\sigma_x\), and de Broglie wavelength).
3. Helper Functions
   • make_grid builds the spatial (\(x\)) and spectral (\(k\)) grids.
   • initial_gaussian creates and normalizes a Gaussian wave packet with mean momentum \(p_0\).
   • linear_potential returns \(V(x)=-qEx\).
   • expectation_x and expectation_p_spectral compute \(\langle x\rangle\) and \(\langle p\rangle\) from the wave function.
   • finite_diff estimates time derivatives via central differences.
4. Simulate Function
   • Precomputes split‑step operators for half‑step potential and full‑step kinetic evolution.
   • Iteratively applies these operators to advance the wave function, recording \(|\psi|^2\), \(\langle x\rangle\), and \(\langle p\rangle\) every few steps.
   • Also computes the classical trajectory \(x_{\rm cl}(t)\) for comparison.
   • Finally, it normalizes densities for plotting and computes \(d\langle x\rangle/dt\) and \(d\langle p\rangle/dt\).
5. animate Function
   • Sets up a \(2\times2\) figure: Panel 1 shows \(|\psi|^2\); Panel 2 plots \(\langle x\rangle\) vs. classical \(x_{\rm cl}\); Panel 3 compares \(d\langle x\rangle/dt\) to \(\langle p\rangle/m\); Panel 4 compares \(d\langle p\rangle/dt\) to the constant force \(qE\).
   • Uses Matplotlib’s FuncAnimation to update these plots over time with moving markers highlighting the current frame.

Running the script displays an animation that demonstrates how quantum averages follow Newton’s laws under a linear potential.

Figure 2 Final frame from a 1D split‑step Fourier simulation of an electron wave packet in a constant electric field, illustrating the Ehrenfest theorem. (a) Normalized probability density \(|\psi(x,t)|^{2}\) at \(t=400\,\mathrm{fs}\) on a domain \(L_{x}=400\,\mathrm{nm}\) with \(N_{x}=1024\). (b) Position: quantum expectation \(\langle x\rangle(t)\) (blue, solid) versus classical trajectory \(x_{\mathrm{cl}}(t)=x_{0}+ (p_{0}/m)t+(qE/2m)t^{2}\) (red, dashed). (c) Velocity check: \(d\langle x\rangle/dt\) (blue, solid) versus \(\langle p\rangle/m\) (red, dashed), where \(\langle p\rangle\) is computed spectrally from \(\psi\). (d) Force check: \(d\langle p\rangle/dt\) (blue, solid) versus the constant force \(qE\) (red, dashed). Parameters: electron mass \(m\), charge \(q=-e\), \(E=10^{7}\,\mathrm{V/m}\), \(\Delta t=0.2\,\mathrm{fs}\), \(x_{0}=80\,\mathrm{nm}\), \(\sigma_{x}=10\,\mathrm{nm}\), de Broglie wavelength \(\lambda_{\mathrm{dB}}=1.5\,\mathrm{nm}\) (\(p_{0}=2\pi\hbar/\lambda_{\mathrm{dB}}\)). Agreement between blue and red curves in (b)–(d) confirms \(d\langle x\rangle/dt=\langle p\rangle/m\) and \(d\langle p\rangle/dt=qE\), demonstrating the Ehrenfest relations under a linear potential.

Discussion of Simulation Results

The results from the simulation are shown in Fig. 2(a)–(d) and in the animation below (Video 1). Each panel of the figure is explained below.

• Density evolution (Panel (a)): The probability density \(|\psi(x,t)|^2\) first moves in the \(+x\) direction (because the initial momentum \(p_0>0\)), then decelerates since the constant force \(qE\) points toward \(-x\) for an electron (\(q=-e\), \(E>0\)). It subsequently accelerates in the \(-x\) direction with minimal distortion, as expected under a uniform force. The total probability is conserved throughout the run, i.e., \(\int |\psi(x,t)|^2\,dx \approx 1\) (numerical unitarity).
• Position vs classical trajectory (Panel (b)): The quantum average position \(\langle x\rangle(t)\) (blue) overlaps the classical path \(x_{\rm cl}(t)=x_0+(p_0/m)t+(qE/2m)t^2\) (red, dashed) within plotting resolution, including the expected quadratic curvature set by the constant acceleration \(a=qE/m\). This directly confirms the position‑level correspondence.
• Velocity identity (Panel (c)): The numerically differentiated \(\frac{d}{dt}\langle x\rangle\) tracks \(\langle p\rangle/m\) at all times, verifying the first Ehrenfest relation. Momentum is computed via a spectral derivative, ensuring high spatial accuracy.
• Force identity (Panel (d)). The quantity \(\frac{d}{dt}\langle p\rangle\) is flat and equal to the constant force \(qE\) (red, dashed), validating the second Ehrenfest relation for a linear potential \(V(x)=-qEx\). Small residuals arise from discrete time stepping and finite‑difference differentiation, not from physics.

Therefore, panels (b)–(d) of Fig. 2 demonstrate that the two Ehrenfest identities are satisfied to within numerical tolerance.

Animation of the Simulation Results

Conclusion

This article demonstrated, both analytically and numerically, how the Ehrenfest theorem links quantum dynamics to Newtonian motion. Using a 1D split–step Fourier (SSFM) simulation of an electron in a constant electric field \(E\), we tracked the quantum expectation values \(\langle x\rangle(t)\) and \(\langle p\rangle(t)\) and compared them with their classical counterparts. The results confirm the theorem’s two identities to numerical tolerance: \(d\langle x\rangle/dt \approx \langle p\rangle/m\) and \(d\langle p\rangle/dt \approx qE\), while \(\langle x\rangle(t)\) overlaps the classical trajectory \(x_{\rm cl}(t)\). Small deviations are attributable to discretization and finite‑difference derivatives, not a failure of the theorem. The linear potential \(V(x)=-qEx\) provides an ideal test case; for nonlinear or strongly varying forces, \(\langle F(x)\rangle \neq F(\langle x\rangle)\) and departures are expected. The framework here extends naturally to harmonic, anharmonic, and time‑dependent fields, enabling systematic studies of quantum–classical correspondence beyond the constant‑force regime.

References

1. G. Friesecke and M. Koppen, “On the Ehrenfest theorem of quantum mechanics,” J. Math. Phys. 50, 082102 (2009). doi:10.1063/1.3191679.
2. L. E. Ballentine, Y. Yang, and J. P. Zibin, “Inadequacy of Ehrenfest’s theorem to characterize the classical regime,” Phys. Rev. A 50, 2854 (1994). doi: 10.1103/PhysRevA.50.2854
3. M. D. Feit, J. A. Fleck, Jr., and A. Steiger, “Solution of the Schrödinger equation by a spectral method,” J. Comput. Phys. 47, 412–433 (1982). doi:10.1016/0021-9991(82)90091-2
4. J. A. C. Weideman and B. M. Herbst, “Split-step methods for the solution of the nonlinear Schrödinger equation,” SIAM J. Numer. Anal. 23, 485–507 (1986). doi:10.1137/0723033
5. M. R. Hermann and J. A. Fleck, Jr., “Split-operator spectral method for solving the time-dependent Schrödinger equation in spherical coordinates,” Phys. Rev. A 38, 6000 (1988). doi:10.1103/PhysRevA.38.6000

Sushanta Barman

I am a senior research scholar in the Department of Physics at IIT Kanpur. My work focuses on ion-beam optics and matter-wave phenomena. I am also interested in emerging matter-wave technologies.

sites.google.com/view/sushanta-barman/home?authuser=0